extra/ryu/ryu.factor

   1 ! Copyright (C) 2018 Alexander Ilin.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: combinators.smart kernel math math.bitwise
   4 math.functions math.order math.parser ryu.data sequences
   5 sequences.private ;
   6
   7 IN: ryu
   8
   9 <PRIVATE
  10
  11 : mul-shift ( x mul shift -- y )
  12     [ first2 rot [ * ] keep swapd * -64 shift + ] [ 64 - neg ] bi* shift ;
  13
  14 : mul-shift-all ( mmShift m mul shift -- vm vp vr )
  15     [ 4 * ] 2dip
  16     [ [ 1 - swap - ] 2dip mul-shift ]
  17     [ [ 2 +        ] 2dip mul-shift ]
  18     [                     mul-shift ] 3tri ;
  19
  20 :: pow-5-factor ( x -- y )
  21     x f 0 [ 2dup x > or ] [
  22         [ 5 /mod ] 2dip rot zero? [ 1 + ] [ nip dupd ] if
  23     ] until 2nip ; inline
  24
  25 : multiple-of-power-of-5 ( p value -- ? )
  26     pow-5-factor <= ;
  27
  28 : double-pow-5-bits ( n -- m )
  29     [ 1 ] [
  30         DOUBLE_LOG2_5_NUMERATOR * DOUBLE_LOG2_5_DENOMINATOR + 1 -
  31         DOUBLE_LOG2_5_DENOMINATOR /i
  32     ] if-zero ; inline
  33
  34 : decimal-length ( m -- n )
  35     {
  36         10
  37         100
  38         1000
  39         10000
  40         100000
  41         1000000
  42         10000000
  43         100000000
  44         1000000000
  45         10000000000
  46         100000000000
  47         1000000000000
  48         10000000000000
  49         100000000000000
  50         1000000000000000
  51         10000000000000000
  52         100000000000000000
  53         1000000000000000000
  54     } [ dupd >= ] find-last [ 2 + ] [ drop 1 ] if nip ; inline
  55
  56 CONSTANT: mantissaBits 52
  57 CONSTANT: exponentBits 11
  58 CONSTANT: offset 1023 ! (1 << (exponentBits - 1)) - 1
  59
  60 :: unpack-bits ( value -- e2 m2 acceptBounds ieeeExponent<=1? neg? string/f )
  61     value double>bits
  62     dup mantissaBits exponentBits + bit? :> sign
  63     dup mantissaBits bits :> ieeeMantissa
  64     mantissaBits neg shift exponentBits bits :> ieeeExponent
  65     0 :> m2!
  66     0 :> e2!
  67     exponentBits on-bits ieeeExponent = [
  68         ieeeMantissa zero? [ sign "-Inf" "Inf" ? ] [ "NaN" ] if
  69     ] [
  70         ieeeExponent [
  71             ieeeMantissa [ sign "-0e0" "0e0" ? ] [
  72                 m2!
  73                 -1 offset - mantissaBits - e2!
  74                 f
  75             ] if-zero
  76         ] [
  77             offset - mantissaBits - 2 - e2!
  78             ieeeMantissa mantissaBits set-bit m2!
  79             f
  80         ] if-zero
  81     ] if [ e2 m2 dup even? ieeeExponent 1 <= sign ] dip ; inline
  82
  83 :: prepare-output ( vp! acceptBounds vmIsTrailingZeros! vrIsTrailingZeros! vr! vm! -- output )
  84     ! vr is converted into the output
  85     0
  86     ! the if has this stack-effect: ( lastRemovedDigit -- lastRemovedDigit' output )
  87     vmIsTrailingZeros vrIsTrailingZeros or [
  88         ! rare
  89         [ vp 10 /i vm 10 /i 2dup > ] [
  90             vm! vp!
  91             vmIsTrailingZeros [ vm 10 divisor? vmIsTrailingZeros! ] when
  92             vrIsTrailingZeros [ dup zero? vrIsTrailingZeros! ] when
  93             vr 10 /mod swap vr! nip ! lastRemovedDigit!
  94         ] while 2drop
  95         vmIsTrailingZeros [
  96             [ vm dup 10 /i dup 10 * swapd = ] [
  97                 vm!
  98                 vrIsTrailingZeros [ dup zero? vrIsTrailingZeros! ] when
  99                 vr 10 /mod swap vr! nip ! lastRemovedDigit!
 100                 vp 10 /i vp!
 101             ] while drop ! Drop (vm 10 /i) result from the while condition.
 102         ] when
 103         vrIsTrailingZeros [
 104             dup 5 = [
 105                 vr even? [ drop 4 ] when ! 4 lastRemovedDigit!
 106             ] when
 107         ] when
 108         vr over 5 >= [ 1 + ] [
 109             dup vm = [
 110                 acceptBounds vmIsTrailingZeros and not [ 1 + ] when
 111             ] when
 112         ] if
 113     ] [
 114         ! common
 115         [ vp 10 /i vm 10 /i 2dup > ] [
 116             vm! vp!
 117             vr 10 /mod swap vr! nip ! lastRemovedDigit!
 118         ] while 2drop
 119         vr dup vm = [ 1 + ] [
 120             over 5 >= [ 1 + ] when
 121         ] if
 122     ] if nip ; inline
 123
 124 :: produce-output ( exp sign output -- string )
 125     [
 126         sign "-" f ?
 127         output number>string 1 cut-slice dup empty? f "." ? swap
 128         "e"
 129         exp number>string
 130     ] "" append-outputs-as ; inline
 131
 132 PRIVATE>
 133
 134 :: print-float ( value -- string )
 135     value >float unpack-bits
 136     :> ( e2 m2 acceptBounds ieeeExponent<=1 sign output )
 137     output [
 138         m2 4 * :> mv
 139         mantissaBits 2^ m2 = not ieeeExponent<=1 or 1 0 ? :> mmShift
 140         f f 0 0 0 :> ( vmIsTrailingZeros! vrIsTrailingZeros! e10! vr! vm! )
 141         ! After the following loop vp is left on stack.
 142         e2 0 >= [
 143             e2 DOUBLE_LOG10_2_NUMERATOR * DOUBLE_LOG10_2_DENOMINATOR /i 0 max :> q
 144             q e10!
 145             q double-pow-5-bits DOUBLE_POW5_INV_BITCOUNT + 1 - :> k
 146             q k + e2 - :> i
 147             mmShift m2 q DOUBLE_POW5_INV_SPLIT nth-unsafe i mul-shift-all vr! swap vm! ! vp on stack
 148             q 21 <= [
 149                 mv 5 divisor? [
 150                     q mv multiple-of-power-of-5 vrIsTrailingZeros!
 151                 ] [
 152                     acceptBounds [
 153                         q mv mmShift - 1 - multiple-of-power-of-5 vmIsTrailingZeros!
 154                     ] [
 155                         q mv 2 + multiple-of-power-of-5 1 0 ? - ! vp!
 156                     ] if
 157                 ] if
 158             ] when
 159         ] [ ! e2 < 0
 160             e2 neg DOUBLE_LOG10_5_NUMERATOR * DOUBLE_LOG10_5_DENOMINATOR /i 1 [-] :> q
 161             q e2 + e10!
 162             e2 neg q - :> i
 163             i double-pow-5-bits DOUBLE_POW5_BITCOUNT - :> k
 164             q k - :> j
 165             mmShift m2 i DOUBLE_POW5_SPLIT nth-unsafe j mul-shift-all vr! swap vm! ! vp on stack
 166             q 1 <= [
 167                 mv 1 bitand bitnot q >= vrIsTrailingZeros!
 168                 acceptBounds [
 169                     mv 1 - mmShift - bitnot 1 bitand q >= vmIsTrailingZeros!
 170                 ] [ 1 - ] if ! vp!
 171             ] [
 172                 q 63 < [
 173                     q 1 - on-bits mv bitand zero? vrIsTrailingZeros!
 174                 ] when
 175             ] if
 176         ] if
 177         [ decimal-length e10 + 1 - sign ] keep ! exp sign vp
 178         acceptBounds vmIsTrailingZeros vrIsTrailingZeros vr vm
 179         prepare-output produce-output
 180     ] unless* ;
 181
 182 ALIAS: d2s print-float