factor/extra/ryu/ryu.factor

! Copyright (C) 2018 Alexander Ilin.
! See http://factorcode.org/license.txt for BSD license.
USING: formatting kernel locals math math.bitwise math.functions
math.order ryu.data sequences shuffle strings vectors ;

IN: ryu

<PRIVATE

: mul-shift ( x mul shift -- y )
    [ first2 rot [ * ] keep swapd * -64 shift + ] [ 64 - neg ] bi* shift ;

: mul-shift-all ( mmShift m mul shift -- vm vp vr )
    [ 4 * ] 2dip
    [ [ 1 - swap - ] 2dip mul-shift ]
    [ [ 2 +        ] 2dip mul-shift ]
    [                     mul-shift ] 3tri ;

:: pow-5-factor ( x -- y )
    x :> value!
    f 0 [ 2dup x <= swap not and ] [
        value 5 /mod zero? [ value! 1 + ] [ nipd swap ] if
    ] while nip ; inline

: multiple-of-power-of-5 ( p value -- ? )
    pow-5-factor <= ;

: double-pow-5-bits ( n -- m )
    [ 1 ] [
        DOUBLE_LOG2_5_NUMERATOR * DOUBLE_LOG2_5_DENOMINATOR + 1 -
        DOUBLE_LOG2_5_DENOMINATOR /i
    ] if-zero ; inline

: decimal-length ( m -- n )
    {
        10
        100
        1000
        10000
        100000
        1000000
        10000000
        100000000
        1000000000
        10000000000
        100000000000
        1000000000000
        10000000000000
        100000000000000
        1000000000000000
        10000000000000000
        100000000000000000
        1000000000000000000
    } [ dupd >= ] find-last [ 2 + ] [ drop 1 ] if nip ; inline

CONSTANT: mantissaBits 52
CONSTANT: exponentBits 11
CONSTANT: offset 1023 ! (1 << (exponentBits - 1)) - 1

:: unpack-bits ( value -- e2 m2 acceptBounds ieeeExponent<=1? neg? string/f )
    value double>bits
    dup mantissaBits exponentBits + bit? :> sign
    dup mantissaBits bits :> ieeeMantissa
    mantissaBits neg shift exponentBits bits :> ieeeExponent
    0 :> m2!
    0 :> e2!
    exponentBits on-bits ieeeExponent = [
        ieeeMantissa zero? [ sign "-Inf" "Inf" ? ] [ "NaN" ] if
    ] [
        ieeeExponent [
            ieeeMantissa [ sign "-0e0" "0e0" ? ] [
                m2!
                -1 offset - mantissaBits - e2!
                f
            ] if-zero
        ] [
            offset - mantissaBits - 2 - e2!
            ieeeMantissa mantissaBits set-bit m2!
            f
        ] if-zero
    ] if [ e2 m2 dup even? ieeeExponent 1 <= sign ] dip ; inline

:: prepare-output ( vp! vplength acceptBounds vmIsTrailingZeros! vrIsTrailingZeros! vr! vm! -- vplength' output )
    ! vr is converted into the output
    0 vplength
    ! the if has this stack-effect: ( lastRemovedDigit vplength -- lastRemovedDigit' vplength' output )
    vmIsTrailingZeros vrIsTrailingZeros or [
        ! rare
        [ vp 10 /i vm 10 /i 2dup > ] [
            vm! vp!
            vmIsTrailingZeros [ vm 10 divisor? vmIsTrailingZeros! ] when
            vrIsTrailingZeros [ over zero? vrIsTrailingZeros! ] when
            vr 10 /mod -roll vr! nip ! lastRemovedDigit!
            1 - ! vplength!
        ] while 2drop
        vmIsTrailingZeros [
            [ vm dup 10 /i dup 10 * swapd = ] [
                vm!
                vrIsTrailingZeros [ over zero? vrIsTrailingZeros! ] when
                vr 10 /mod -roll vr! nip ! lastRemovedDigit!
                vp 10 /i vp!
                1 - ! vplength!
            ] while drop ! Drop (vm 10 /i) result from the while condition.
        ] when
        vrIsTrailingZeros [
            over 5 = [
                vr even? [ 4 -rot nip ] when ! 4 lastRemovedDigit!
            ] when
        ] when
        vr pick 5 >= [ 1 + ] [
            dup vm = [
                acceptBounds vmIsTrailingZeros and not [ 1 + ] when
            ] when
        ] if
    ] [
        ! common
        [ vp 10 /i vm 10 /i 2dup > ] [
            vm! vp!
            vr 10 /mod -roll vr! nip ! lastRemovedDigit!
            1 - ! vplength!
        ] while 2drop
        vr dup vm = [ 1 + ] [
            pick 5 >= [ 1 + ] when
        ] if
    ] if nipd ; inline

: write-char ( index seq char -- index+1 seq' )
    -rot [ tuck ] dip [ set-nth 1 + ] keep ; inline

: write-exp ( exp index result -- result' )
    CHAR: e write-char
    pick neg? [
        CHAR: - write-char [ neg ] 2dip
    ] when
    pick dup 100 >= [
        100 /i CHAR: 0 + write-char
        [ 100 mod 2 * ] 2dip
        pick DIGIT_TABLE nth write-char
        [ 1 + DIGIT_TABLE nth ] 2dip [ set-nth ] keep
    ] [
        10 >= [
            [ 2 * ] 2dip
            pick DIGIT_TABLE nth write-char
            [ 1 + DIGIT_TABLE nth ] 2dip [ set-nth ] keep
        ] [
            [ CHAR: 0 + ] 2dip [ set-nth ] keep
        ] if
    ] if ; inline

:: produce-output ( exp sign olength output2! -- string )
    25 <vector> 0 :> ( result i! )
    0 sign [ CHAR: - swap result set-nth 1 ] when :> index!
    [ output2 10000 >= ] [
        output2 dup 10000 /i dup output2! 10000 * - :> c
        index olength + i - 1 - :> res-index
        c 100 mod 2 *
        dup DIGIT_TABLE nth res-index result set-nth
        1 + DIGIT_TABLE nth res-index 1 + result set-nth
        c 100 /i 2 *
        dup DIGIT_TABLE nth res-index 2 - result set-nth
        1 + DIGIT_TABLE nth res-index 1 - result set-nth
        i 4 + i!
    ] while
    output2 100 >= [
        output2 dup 100 /i dup output2! 100 * - 2 * :> c
        index olength + i - :> res-index
        c DIGIT_TABLE nth res-index 1 - result set-nth
        c 1 + DIGIT_TABLE nth res-index result set-nth
        i 2 + i!
    ] when
    output2 10 >= [
        output2 2 * :> c
        index olength + i - :> res-index
        c 1 + DIGIT_TABLE nth res-index result set-nth
        c DIGIT_TABLE nth index result set-nth
    ] [ CHAR: 0 output2 + index result set-nth ] if
    index 1 + index!
    olength 1 > [
        CHAR: . index result set-nth
        index olength + index!
    ] when exp index result write-exp >string ; inline

PRIVATE>

:: print-float ( value -- string )
    value >float unpack-bits [
        [ 5drop ] dip
    ] [| e2 m2 acceptBounds ieeeExponent<=1 sign |
        m2 4 * :> mv
        mantissaBits 2^ m2 = not ieeeExponent<=1 or 1 0 ? :> mmShift
        f f 0 0 0 :> ( vmIsTrailingZeros! vrIsTrailingZeros! e10! vr! vm! )
        ! After the following loop vp is left on stack.
        e2 0 >= [
            e2 DOUBLE_LOG10_2_NUMERATOR * DOUBLE_LOG10_2_DENOMINATOR /i 0 max :> q
            q e10!
            q double-pow-5-bits DOUBLE_POW5_INV_BITCOUNT + 1 - :> k
            q k + e2 - :> i
            mmShift m2 q DOUBLE_POW5_INV_SPLIT nth i mul-shift-all vr! swap vm! ! vp on stack
            q 21 <= [
                mv 5 divisor? [
                    q mv multiple-of-power-of-5 vrIsTrailingZeros!
                ] [
                    acceptBounds [
                        q mv mmShift - 1 - multiple-of-power-of-5 vmIsTrailingZeros!
                    ] [
                        q mv 2 + multiple-of-power-of-5 1 0 ? - ! vp!
                    ] if
                ] if
            ] when
        ] [ ! e2 < 0
            e2 neg DOUBLE_LOG10_5_NUMERATOR * DOUBLE_LOG10_5_DENOMINATOR /i 1 - 0 max :> q
            q e2 + e10!
            e2 neg q - :> i
            i double-pow-5-bits DOUBLE_POW5_BITCOUNT - :> k
            q k - :> j
            mmShift m2 i DOUBLE_POW5_SPLIT nth j mul-shift-all vr! swap vm! ! vp on stack
            q 1 <= [
                mv 1 bitand bitnot q >= vrIsTrailingZeros!
                acceptBounds [
                    mv 1 - mmShift - bitnot 1 bitand q >= vmIsTrailingZeros!
                ] [ 1 - ] if ! vp!
            ] [
                q 63 < [
                    q 1 - 2^ 1 - mv bitand zero? vrIsTrailingZeros!
                ] when
            ] if
        ] if
        dup decimal-length ! vp vplength
        dup e10 + 1 - sign 2swap ! exp and sign for produce-output
        acceptBounds vmIsTrailingZeros vrIsTrailingZeros vr vm
        prepare-output produce-output
    ] if* ;

ALIAS: d2s print-float