extra/ryu/ryu.factor

   1 ! Copyright (C) 2018 Alexander Ilin.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: formatting kernel locals math math.bitwise math.functions
   4 math.order ryu.data sequences shuffle strings vectors ;
   5
   6 IN: ryu
   7
   8 <PRIVATE
   9
  10 : mul-shift ( x mul shift -- y )
  11     [ first2 rot [ * ] keep swapd * -64 shift + ] [ 64 - neg ] bi* shift ;
  12
  13 : mul-shift-all ( mmShift m mul shift -- vm vp vr )
  14     [ 4 * ] 2dip
  15     [ [ 1 - swap - ] 2dip mul-shift ]
  16     [ [ 2 +        ] 2dip mul-shift ]
  17     [                     mul-shift ] 3tri ;
  18
  19 :: pow-5-factor ( x -- y )
  20     x :> value!
  21     f 0 [ 2dup x <= swap not and ] [
  22         value 5 /mod zero? [ value! 1 + ] [ nipd swap ] if
  23     ] while nip ; 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! vplength acceptBounds vmIsTrailingZeros! vrIsTrailingZeros! vr! vm! -- vplength' output )
  84     ! vr is converted into the output
  85     0 vplength
  86     ! the if has this stack-effect: ( lastRemovedDigit vplength -- lastRemovedDigit' vplength' 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 [ over zero? vrIsTrailingZeros! ] when
  93             vr 10 /mod -roll vr! nip ! lastRemovedDigit!
  94             1 - ! vplength!
  95         ] while 2drop
  96         vmIsTrailingZeros [
  97             [ vm dup 10 /i dup 10 * swapd = ] [
  98                 vm!
  99                 vrIsTrailingZeros [ over zero? vrIsTrailingZeros! ] when
 100                 vr 10 /mod -roll vr! nip ! lastRemovedDigit!
 101                 vp 10 /i vp!
 102                 1 - ! vplength!
 103             ] while drop ! Drop (vm 10 /i) result from the while condition.
 104         ] when
 105         vrIsTrailingZeros [
 106             over 5 = [
 107                 vr even? [ 4 -rot nip ] when ! 4 lastRemovedDigit!
 108             ] when
 109         ] when
 110         vr pick 5 >= [ 1 + ] [
 111             dup vm = [
 112                 acceptBounds vmIsTrailingZeros and not [ 1 + ] when
 113             ] when
 114         ] if
 115     ] [
 116         ! common
 117         [ vp 10 /i vm 10 /i 2dup > ] [
 118             vm! vp!
 119             vr 10 /mod -roll vr! nip ! lastRemovedDigit!
 120             1 - ! vplength!
 121         ] while 2drop
 122         vr dup vm = [ 1 + ] [
 123             pick 5 >= [ 1 + ] when
 124         ] if
 125     ] if nipd ; inline
 126
 127 : write-char ( index seq char -- index+1 seq' )
 128     -rot [ tuck ] dip [ set-nth 1 + ] keep ; inline
 129
 130 : write-exp ( exp index result -- result' )
 131     CHAR: e write-char
 132     pick neg? [
 133         CHAR: - write-char [ neg ] 2dip
 134     ] when
 135     pick dup 100 >= [
 136         100 /i CHAR: 0 + write-char
 137         [ 100 mod 2 * ] 2dip
 138         pick DIGIT_TABLE nth write-char
 139         [ 1 + DIGIT_TABLE nth ] 2dip [ set-nth ] keep
 140     ] [
 141         10 >= [
 142             [ 2 * ] 2dip
 143             pick DIGIT_TABLE nth write-char
 144             [ 1 + DIGIT_TABLE nth ] 2dip [ set-nth ] keep
 145         ] [
 146             [ CHAR: 0 + ] 2dip [ set-nth ] keep
 147         ] if
 148     ] if ; inline
 149
 150 :: produce-output ( exp sign olength output2! -- string )
 151     25 <vector> 0 :> ( result i! )
 152     0 sign [ CHAR: - swap result set-nth 1 ] when :> index!
 153     [ output2 10000 >= ] [
 154         output2 dup 10000 /i dup output2! 10000 * - :> c
 155         index olength + i - 1 - :> res-index
 156         c 100 mod 2 *
 157         dup DIGIT_TABLE nth res-index result set-nth
 158         1 + DIGIT_TABLE nth res-index 1 + result set-nth
 159         c 100 /i 2 *
 160         dup DIGIT_TABLE nth res-index 2 - result set-nth
 161         1 + DIGIT_TABLE nth res-index 1 - result set-nth
 162         i 4 + i!
 163     ] while
 164     output2 100 >= [
 165         output2 dup 100 /i dup output2! 100 * - 2 * :> c
 166         index olength + i - :> res-index
 167         c DIGIT_TABLE nth res-index 1 - result set-nth
 168         c 1 + DIGIT_TABLE nth res-index result set-nth
 169         i 2 + i!
 170     ] when
 171     output2 10 >= [
 172         output2 2 * :> c
 173         index olength + i - :> res-index
 174         c 1 + DIGIT_TABLE nth res-index result set-nth
 175         c DIGIT_TABLE nth index result set-nth
 176     ] [ CHAR: 0 output2 + index result set-nth ] if
 177     index 1 + index!
 178     olength 1 > [
 179         CHAR: . index result set-nth
 180         index olength + index!
 181     ] when exp index result write-exp >string ; inline
 182
 183 PRIVATE>
 184
 185 :: print-float ( value -- string )
 186     value >float unpack-bits [
 187         [ 5drop ] dip
 188     ] [| e2 m2 acceptBounds ieeeExponent<=1 sign |
 189         m2 4 * :> mv
 190         mantissaBits 2^ m2 = not ieeeExponent<=1 or 1 0 ? :> mmShift
 191         f f 0 0 0 :> ( vmIsTrailingZeros! vrIsTrailingZeros! e10! vr! vm! )
 192         ! After the following loop vp is left on stack.
 193         e2 0 >= [
 194             e2 DOUBLE_LOG10_2_NUMERATOR * DOUBLE_LOG10_2_DENOMINATOR /i 0 max :> q
 195             q e10!
 196             q double-pow-5-bits DOUBLE_POW5_INV_BITCOUNT + 1 - :> k
 197             q k + e2 - :> i
 198             mmShift m2 q DOUBLE_POW5_INV_SPLIT nth i mul-shift-all vr! swap vm! ! vp on stack
 199             q 21 <= [
 200                 mv 5 divisor? [
 201                     q mv multiple-of-power-of-5 vrIsTrailingZeros!
 202                 ] [
 203                     acceptBounds [
 204                         q mv mmShift - 1 - multiple-of-power-of-5 vmIsTrailingZeros!
 205                     ] [
 206                         q mv 2 + multiple-of-power-of-5 1 0 ? - ! vp!
 207                     ] if
 208                 ] if
 209             ] when
 210         ] [ ! e2 < 0
 211             e2 neg DOUBLE_LOG10_5_NUMERATOR * DOUBLE_LOG10_5_DENOMINATOR /i 1 [-] :> q
 212             q e2 + e10!
 213             e2 neg q - :> i
 214             i double-pow-5-bits DOUBLE_POW5_BITCOUNT - :> k
 215             q k - :> j
 216             mmShift m2 i DOUBLE_POW5_SPLIT nth j mul-shift-all vr! swap vm! ! vp on stack
 217             q 1 <= [
 218                 mv 1 bitand bitnot q >= vrIsTrailingZeros!
 219                 acceptBounds [
 220                     mv 1 - mmShift - bitnot 1 bitand q >= vmIsTrailingZeros!
 221                 ] [ 1 - ] if ! vp!
 222             ] [
 223                 q 63 < [
 224                     q 1 - on-bits mv bitand zero? vrIsTrailingZeros!
 225                 ] when
 226             ] if
 227         ] if
 228         dup decimal-length ! vp vplength
 229         dup e10 + 1 - sign 2swap ! exp and sign for produce-output
 230         acceptBounds vmIsTrailingZeros vrIsTrailingZeros vr vm
 231         prepare-output produce-output
 232     ] if* ;
 233
 234 ALIAS: d2s print-float