## Copyright 2020 Alexander Bolz
##
## Distributed under the Boost Software License, Version 1.0.
## (See accompanying file LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
## --------------------------------------------------------------------------------------------------
## This file contains an implementation of the Schubfach algorithm as described in
##
## [1] Raffaello Giulietti, "The Schubfach way to render doubles",
## https://drive.google.com/open?id=1luHhyQF9zKlM8yJ1nebU0OgVYhfC6CBN
## --------------------------------------------------------------------------------------------------
import std/private/digitsutils
when defined(nimPreviewSlimSystem):
import std/assertions
template sf_Assert(x: untyped): untyped =
assert(x)
## ==================================================================================================
##
## ==================================================================================================
type
ValueType = float32
BitsType = uint32
Single {.bycopy.} = object
bits: BitsType
const
significandSize: int32 = 24
MaxExponent = 128
exponentBias: int32 = MaxExponent - 1 + (significandSize - 1)
maxIeeeExponent: BitsType = BitsType(2 * MaxExponent - 1)
hiddenBit: BitsType = BitsType(1) shl (significandSize - 1)
significandMask: BitsType = hiddenBit - 1
exponentMask: BitsType = maxIeeeExponent shl (significandSize - 1)
signMask: BitsType = not (not BitsType(0) shr 1)
proc constructSingle(bits: BitsType): Single {.constructor.} =
result.bits = bits
proc constructSingle(value: ValueType): Single {.constructor.} =
result.bits = cast[typeof(result.bits)](value)
proc physicalSignificand(this: Single): BitsType {.noSideEffect.} =
return this.bits and significandMask
proc physicalExponent(this: Single): BitsType {.noSideEffect.} =
return (this.bits and exponentMask) shr (significandSize - 1)
proc isFinite(this: Single): bool {.noSideEffect.} =
return (this.bits and exponentMask) != exponentMask
proc isInf(this: Single): bool {.noSideEffect.} =
return (this.bits and exponentMask) == exponentMask and
(this.bits and significandMask) == 0
proc isNaN(this: Single): bool {.noSideEffect.} =
return (this.bits and exponentMask) == exponentMask and
(this.bits and significandMask) != 0
proc isZero(this: Single): bool {.noSideEffect.} =
return (this.bits and not signMask) == 0
proc signBit(this: Single): int {.noSideEffect.} =
return int((this.bits and signMask) != 0)
## ==================================================================================================
## Returns floor(x / 2^n).
##
## Technically, right-shift of negative integers is implementation defined...
## Should easily be optimized into SAR (or equivalent) instruction.
proc floorDivPow2(x: int32; n: int32): int32 {.inline.} =
return x shr n
## Returns floor(log_10(2^e))
## static inline int32_t FloorLog10Pow2(int32_t e)
## {
## SF_ASSERT(e >= -1500);
## SF_ASSERT(e <= 1500);
## return FloorDivPow2(e * 1262611, 22);
## }
## Returns floor(log_10(3/4 2^e))
## static inline int32_t FloorLog10ThreeQuartersPow2(int32_t e)
## {
## SF_ASSERT(e >= -1500);
## SF_ASSERT(e <= 1500);
## return FloorDivPow2(e * 1262611 - 524031, 22);
## }
## Returns floor(log_2(10^e))
proc floorLog2Pow10(e: int32): int32 {.inline.} =
sf_Assert(e >= -1233)
sf_Assert(e <= 1233)
return floorDivPow2(e * 1741647, 19)
const
kMin: int32 = -31
kMax: int32 = 45
g: array[kMax - kMin + 1, uint64] = [0x81CEB32C4B43FCF5'u64, 0xA2425FF75E14FC32'u64,
0xCAD2F7F5359A3B3F'u64, 0xFD87B5F28300CA0E'u64, 0x9E74D1B791E07E49'u64,
0xC612062576589DDB'u64, 0xF79687AED3EEC552'u64, 0x9ABE14CD44753B53'u64,
0xC16D9A0095928A28'u64, 0xF1C90080BAF72CB2'u64, 0x971DA05074DA7BEF'u64,
0xBCE5086492111AEB'u64, 0xEC1E4A7DB69561A6'u64, 0x9392EE8E921D5D08'u64,
0xB877AA3236A4B44A'u64, 0xE69594BEC44DE15C'u64, 0x901D7CF73AB0ACDA'u64,
0xB424DC35095CD810'u64, 0xE12E13424BB40E14'u64, 0x8CBCCC096F5088CC'u64,
0xAFEBFF0BCB24AAFF'u64, 0xDBE6FECEBDEDD5BF'u64, 0x89705F4136B4A598'u64,
0xABCC77118461CEFD'u64, 0xD6BF94D5E57A42BD'u64, 0x8637BD05AF6C69B6'u64,
0xA7C5AC471B478424'u64, 0xD1B71758E219652C'u64, 0x83126E978D4FDF3C'u64,
0xA3D70A3D70A3D70B'u64, 0xCCCCCCCCCCCCCCCD'u64, 0x8000000000000000'u64,
0xA000000000000000'u64, 0xC800000000000000'u64, 0xFA00000000000000'u64,
0x9C40000000000000'u64, 0xC350000000000000'u64, 0xF424000000000000'u64,
0x9896800000000000'u64, 0xBEBC200000000000'u64, 0xEE6B280000000000'u64,
0x9502F90000000000'u64, 0xBA43B74000000000'u64, 0xE8D4A51000000000'u64,
0x9184E72A00000000'u64, 0xB5E620F480000000'u64, 0xE35FA931A0000000'u64,
0x8E1BC9BF04000000'u64, 0xB1A2BC2EC5000000'u64, 0xDE0B6B3A76400000'u64,
0x8AC7230489E80000'u64, 0xAD78EBC5AC620000'u64, 0xD8D726B7177A8000'u64,
0x878678326EAC9000'u64, 0xA968163F0A57B400'u64, 0xD3C21BCECCEDA100'u64,
0x84595161401484A0'u64, 0xA56FA5B99019A5C8'u64, 0xCECB8F27F4200F3A'u64,
0x813F3978F8940985'u64, 0xA18F07D736B90BE6'u64, 0xC9F2C9CD04674EDF'u64,
0xFC6F7C4045812297'u64, 0x9DC5ADA82B70B59E'u64, 0xC5371912364CE306'u64,
0xF684DF56C3E01BC7'u64, 0x9A130B963A6C115D'u64, 0xC097CE7BC90715B4'u64,
0xF0BDC21ABB48DB21'u64, 0x96769950B50D88F5'u64, 0xBC143FA4E250EB32'u64,
0xEB194F8E1AE525FE'u64, 0x92EFD1B8D0CF37BF'u64, 0xB7ABC627050305AE'u64,
0xE596B7B0C643C71A'u64, 0x8F7E32CE7BEA5C70'u64, 0xB35DBF821AE4F38C'u64]
proc computePow10Single(k: int32): uint64 {.inline.} =
## There are unique beta and r such that 10^k = beta 2^r and
## 2^63 <= beta < 2^64, namely r = floor(log_2 10^k) - 63 and
## beta = 2^-r 10^k.
## Let g = ceil(beta), so (g-1) 2^r < 10^k <= g 2^r, with the latter
## value being a pretty good overestimate for 10^k.
## NB: Since for all the required exponents k, we have g < 2^64,
## all constants can be stored in 128-bit integers.
sf_Assert(k >= kMin)
sf_Assert(k <= kMax)
return g[k - kMin]
proc lo32(x: uint64): uint32 {.inline.} =
return cast[uint32](x)
proc hi32(x: uint64): uint32 {.inline.} =
return cast[uint32](x shr 32)
when defined(sizeof_Int128):
proc roundToOdd(g: uint64; cp: uint32): uint32 {.inline.} =
let p: uint128 = uint128(g) * cp
let y1: uint32 = lo32(cast[uint64](p shr 64))
let y0: uint32 = hi32(cast[uint64](p))
return y1 or uint32(y0 > 1)
elif defined(vcc) and defined(cpu64):
proc umul128(x, y: uint64, z: ptr uint64): uint64 {.importc: "_umul128", header: "<intrin.h>".}
proc roundToOdd(g: uint64; cpHi: uint32): uint32 {.inline.} =
var p1: uint64 = 0
var p0: uint64 = umul128(g, cpHi, addr(p1))
let y1: uint32 = lo32(p1)
let y0: uint32 = hi32(p0)
return y1 or uint32(y0 > 1)
else:
proc roundToOdd(g: uint64; cp: uint32): uint32 {.inline.} =
let b01: uint64 = uint64(lo32(g)) * cp
let b11: uint64 = uint64(hi32(g)) * cp
let hi: uint64 = b11 + hi32(b01)
let y1: uint32 = hi32(hi)
let y0: uint32 = lo32(hi)
return y1 or uint32(y0 > 1)
## Returns whether value is divisible by 2^e2
proc multipleOfPow2(value: uint32; e2: int32): bool {.inline.} =
sf_Assert(e2 >= 0)
sf_Assert(e2 <= 31)
return (value and ((uint32(1) shl e2) - 1)) == 0
type
FloatingDecimal32 {.bycopy.} = object
digits: uint32 ## num_digits <= 9
exponent: int32
proc toDecimal32(ieeeSignificand: uint32; ieeeExponent: uint32): FloatingDecimal32 {.
inline.} =
var c: uint32
var q: int32
if ieeeExponent != 0:
c = hiddenBit or ieeeSignificand
q = cast[int32](ieeeExponent) - exponentBias
if 0 <= -q and -q < significandSize and multipleOfPow2(c, -q):
return FloatingDecimal32(digits: c shr -q, exponent: 0'i32)
else:
c = ieeeSignificand
q = 1 - exponentBias
let isEven: bool = (c mod 2 == 0)
let lowerBoundaryIsCloser: bool = (ieeeSignificand == 0 and ieeeExponent > 1)
## const int32_t qb = q - 2;
let cbl: uint32 = 4 * c - 2 + uint32(lowerBoundaryIsCloser)
let cb: uint32 = 4 * c
let cbr: uint32 = 4 * c + 2
## (q * 1262611 ) >> 22 == floor(log_10( 2^q))
## (q * 1262611 - 524031) >> 22 == floor(log_10(3/4 2^q))
sf_Assert(q >= -1500)
sf_Assert(q <= 1500)
let k: int32 = floorDivPow2(q * 1262611 - (if lowerBoundaryIsCloser: 524031 else: 0), 22)
let h: int32 = q + floorLog2Pow10(-k) + 1
sf_Assert(h >= 1)
sf_Assert(h <= 4)
let pow10: uint64 = computePow10Single(-k)
let vbl: uint32 = roundToOdd(pow10, cbl shl h)
let vb: uint32 = roundToOdd(pow10, cb shl h)
let vbr: uint32 = roundToOdd(pow10, cbr shl h)
let lower: uint32 = vbl + uint32(not isEven)
let upper: uint32 = vbr - uint32(not isEven)
## See Figure 4 in [1].
## And the modifications in Figure 6.
let s: uint32 = vb div 4
## NB: 4 * s == vb & ~3 == vb & -4
if s >= 10:
let sp: uint32 = s div 10
## = vb / 40
let upInside: bool = lower <= 40 * sp
let wpInside: bool = 40 * sp + 40 <= upper
## if (up_inside || wp_inside) // NB: At most one of u' and w' is in R_v.
if upInside != wpInside:
return FloatingDecimal32(digits: sp + uint32(wpInside), exponent: k + 1)
let uInside: bool = lower <= 4 * s
let wInside: bool = 4 * s + 4 <= upper
if uInside != wInside:
return FloatingDecimal32(digits: s + uint32(wInside), exponent: k)
let mid: uint32 = 4 * s + 2
## = 2(s + t)
let roundUp: bool = vb > mid or (vb == mid and (s and 1) != 0)
return FloatingDecimal32(digits: s + uint32(roundUp), exponent: k)
## ==================================================================================================
## ToChars
## ==================================================================================================
proc printDecimalDigitsBackwards(buf: var openArray[char]; pos: int; output: uint32): int32 {.inline.} =
var output = output
var pos = pos
var tz: int32 = 0
## number of trailing zeros removed.
var nd: int32 = 0
## number of decimal digits processed.
## At most 9 digits remaining
if output >= 10000:
let q: uint32 = output div 10000
let r: uint32 = output mod 10000
output = q
dec(pos, 4)
if r != 0:
let rH: uint32 = r div 100
let rL: uint32 = r mod 100
utoa2Digits(buf, pos, rH)
utoa2Digits(buf, pos + 2, rL)
tz = trailingZeros2Digits(if rL == 0: rH else: rL) + (if rL == 0: 2 else: 0)
else:
tz = 4
nd = 4
if output >= 100:
let q: uint32 = output div 100
let r: uint32 = output mod 100
output = q
dec(pos, 2)
utoa2Digits(buf, pos, r)
if tz == nd:
inc(tz, trailingZeros2Digits(r))
inc(nd, 2)
if output >= 100:
let q2: uint32 = output div 100
let r2: uint32 = output mod 100
output = q2
dec(pos, 2)
utoa2Digits(buf, pos, r2)
if tz == nd:
inc(tz, trailingZeros2Digits(r2))
inc(nd, 2)
sf_Assert(output >= 1)
sf_Assert(output <= 99)
if output >= 10:
let q: uint32 = output
dec(pos, 2)
utoa2Digits(buf, pos, q)
if tz == nd:
inc(tz, trailingZeros2Digits(q))
else:
let q: uint32 = output
sf_Assert(q >= 1)
sf_Assert(q <= 9)
dec(pos)
buf[pos] = chr(uint32('0') + q)
return tz
proc decimalLength(v: uint32): int32 {.inline.} =
sf_Assert(v >= 1)
sf_Assert(v <= 999999999'u)
if v >= 100000000'u:
return 9
if v >= 10000000'u:
return 8
if v >= 1000000'u:
return 7
if v >= 100000'u:
return 6
if v >= 10000'u:
return 5
if v >= 1000'u:
return 4
if v >= 100'u:
return 3
if v >= 10'u:
return 2
return 1
proc formatDigits(buffer: var openArray[char]; pos: int; digits: uint32; decimalExponent: int32;
forceTrailingDotZero: bool = false): int {.inline.} =
const
minFixedDecimalPoint: int32 = -4
maxFixedDecimalPoint: int32 = 9
var pos = pos
assert(minFixedDecimalPoint <= -1, "internal error")
assert(maxFixedDecimalPoint >= 1, "internal error")
sf_Assert(digits >= 1)
sf_Assert(digits <= 999999999'u)
sf_Assert(decimalExponent >= -99)
sf_Assert(decimalExponent <= 99)
var numDigits: int32 = decimalLength(digits)
let decimalPoint: int32 = numDigits + decimalExponent
let useFixed: bool = minFixedDecimalPoint <= decimalPoint and
decimalPoint <= maxFixedDecimalPoint
## Prepare the buffer.
## Avoid calling memset/memcpy with variable arguments below...
for i in 0..<32: buffer[pos+i] = '0'
assert(minFixedDecimalPoint >= -30, "internal error")
assert(maxFixedDecimalPoint <= 32, "internal error")
var decimalDigitsPosition: int32
if useFixed:
if decimalPoint <= 0:
## 0.[000]digits
decimalDigitsPosition = 2 - decimalPoint
else:
## dig.its
## digits[000]
decimalDigitsPosition = 0
else:
## dE+123 or d.igitsE+123
decimalDigitsPosition = 1
var digitsEnd = pos + decimalDigitsPosition + numDigits
let tz: int32 = printDecimalDigitsBackwards(buffer, digitsEnd, digits)
dec(digitsEnd, tz)
dec(numDigits, tz)
## decimal_exponent += tz; // => decimal_point unchanged.
if useFixed:
if decimalPoint <= 0:
## 0.[000]digits
buffer[pos+1] = '.'
pos = digitsEnd
elif decimalPoint < numDigits:
## dig.its
for i in countdown(7, 0):
buffer[i + decimalPoint + 1] = buffer[i + decimalPoint]
buffer[pos+decimalPoint] = '.'
pos = digitsEnd + 1
else:
## digits[000]
inc(pos, decimalPoint)
if forceTrailingDotZero:
buffer[pos] = '.'
buffer[pos+1] = '0'
inc(pos, 2)
else:
buffer[pos] = buffer[pos+1]
if numDigits == 1:
## dE+123
inc(pos)
else:
## d.igitsE+123
buffer[pos+1] = '.'
pos = digitsEnd
let scientificExponent: int32 = decimalPoint - 1
## SF_ASSERT(scientific_exponent != 0);
buffer[pos] = 'e'
buffer[pos+1] = if scientificExponent < 0: '-' else: '+'
inc(pos, 2)
let k: uint32 = cast[uint32](if scientificExponent < 0: -scientificExponent else: scientificExponent)
if k < 10:
buffer[pos] = chr(uint32('0') + k)
inc pos
else:
utoa2Digits(buffer, pos, k)
inc(pos, 2)
return pos
proc float32ToChars*(buffer: var openArray[char]; v: float32; forceTrailingDotZero = false): int {.
inline.} =
let significand: uint32 = physicalSignificand(constructSingle(v))
let exponent: uint32 = physicalExponent(constructSingle(v))
var pos = 0
if exponent != maxIeeeExponent:
## Finite
buffer[pos] = '-'
inc(pos, signBit(constructSingle(v)))
if exponent != 0 or significand != 0:
## != 0
let dec: auto = toDecimal32(significand, exponent)
return formatDigits(buffer, pos, dec.digits, dec.exponent, forceTrailingDotZero)
else:
buffer[pos] = '0'
buffer[pos+1] = '.'
buffer[pos+2] = '0'
buffer[pos+3] = ' '
inc(pos, if forceTrailingDotZero: 3 else: 1)
return pos
if significand == 0:
buffer[pos] = '-'
inc(pos, signBit(constructSingle(v)))
buffer[pos] = 'i'
buffer[pos+1] = 'n'
buffer[pos+2] = 'f'
buffer[pos+3] = ' '
return pos + 3
else:
buffer[pos] = 'n'
buffer[pos+1] = 'a'
buffer[pos+2] = 'n'
buffer[pos+3] = ' '
return pos + 3