! Copyright (C) 2024 Giftpflanze.
! See https://factorcode.org/license.txt for BSD license.
USING: combinators dragonbox.lookup-table kernel math
math.bitwise math.floating-point math.functions math.order ;
IN: dragonbox

! https://github.com/jk-jeon/dragonbox/blob/master/other_files/Dragonbox.pdf

CONSTANT: q 64
CONSTANT: p 52
CONSTANT: E_min -1022
CONSTANT: bias 1023
CONSTANT: κ 2
CONSTANT: Q 128

: ⌊nlog10_2⌋ ( n -- m ) 315653 * -20 shift ; inline
: ⌊nlog2_10⌋ ( n -- m ) 1741647 * -19 shift ; inline
: ⌊nlog10_2-log10_4/3⌋ ( n -- m )
    631305 * 261663 - -21 shift ; inline

: k ( E -- k ) ⌊nlog10_2⌋ neg κ + ; inline
: β ( E k -- β ) ⌊nlog2_10⌋ + ; inline

MACRO: (i) ( quot quot -- quot: ( F φ β -- (i)' int? ) )
    '[
        [ 2 * @ ] [ * ] [ Q - ] tri*
        [ shift @ ] [ neg bits zero? ] 2bi
    ] ;

: x(i) ( F φ β -- x-odd? x-int? ) [ 1 - ] [ odd? ] (i) ; inline
: y(i) ( F φ β -- y-odd? y-int? ) [     ] [ odd? ] (i) ; inline
: z(i) ( F φ β -- z(i)   z-int? ) [ 1 + ] [      ] (i) ; inline

: s ( z(i) -- s ) 2361183241434822607 * -71 shift ; inline
: r ( z(i) s -- r ) 1000 * - ; inline
: s/r ( z(i) -- s r ) [ s ] keep over r ; inline

: δ(i) ( φ β -- δ(i) ) Q - 1 + shift ; inline

: strip-zeroes ( s -- s' d )
    0 [
        over 10 /mod zero?
        [ [ nipd swap 1 + ] [ drop ] if ] keep
    ] loop ; inline

:: normal-interval ( F E -- f e )
    E k :> k
    k φ :> φ
    E k β :> β
    F φ β z(i) :> ( z(i) z-int? )
    z(i) s/r :> ( s r )
    s r :> ( s̃! r̃! )
    φ β δ(i) :> δ(i)
    F even? :> w∈I?
    r δ(i) <=> {
        { +gt+ [ f ] }
        { +lt+ [
            w∈I? r zero? z-int? and and
            [ [ 1000 r̃! s 1 - s̃! ] when ] keep not
        ] }
        { +eq+ [
            F φ β x(i) :> ( x-odd? x-int? )
            x-odd? [ w∈I? not x-int? or ] unless*
        ] }
    } case [
        s strip-zeroes k - κ + 1 +
    ] [
        r̃ 50 + δ(i) 2/ - :> D
        D 100 /i :> t
        s̃ 10 * t + {
            { [ t 100 * D number= ] [ ] }
            { [
                F φ β y(i) :> ( y-odd? y-int? )
                D 50 - even? y-odd? eq? w∈I? not y-int? and or
            ] [ 1 - ] }
            [ ]
        } cond κ k -
    ] if ; inline

: k0 ( E -- k0 ) ⌊nlog10_2-log10_4/3⌋ neg ; inline

:: x(i)' ( φ β -- x(i) )
    φ q Q - shift
    φ q Q - p - 2 - shift -
    p q - β + 1 + shift ; inline

:: z(i)' ( φ β -- z(i) )
    φ q Q - shift
    φ q Q - p - 1 - shift +
    p q - β + 1 + shift ; inline

! x̃(i) = x(i)   if x integer and x ∈ 10^k0 I
!        x(i)+1 otherwise
! x∈I             ⇔ F even
! x(i) is integer ⇔ 2 ≤ e ≤ 2 + ⌊log2 10^(d1+1)/3⌋ = 3
: x̃(i) ( x(i) F E -- x̃(i) )
    [ even? ] [ 2 3 between? ] bi* and [ 1 + ] unless ;

! z̃(i) = z(i)   if z not integer or z ∈ 10^k0 I
!        z(i)-1 otherwise
! z∈I             ⇔ F even
! z(i) is integer ⇔ 0 ≤ e ≤ 2 + ⌊log2 10^(d2+1)/3⌋ = 3
: z̃(i) ( z(i) F E -- z̃(i) )
    [ even? ] [ 0 3 between? not ] bi* or [ 1 - ] unless ;

: y(ru) ( φ β -- y(ru) ) p + 2 + Q - shift 1 + 2/ ;

:: shorter-interval ( F E -- f e )
    E k0 :> k0
    k0 φ :> φ
    E k0 β :> β
    ! x(i), z(i)  (5.2.1)
    ! x̃(i), z̃(i)  (5.2.3)
    φ β x(i)' F E x̃(i) :> x(i)
    φ β z(i)' F E z̃(i) :> z(i)
    ! ⌊z̃(i)/10⌋
    z(i) 10 /i :> z*
    ! x̃(i) ≤ ⌊z̃(i)/10⌋⋅10: ⌊z̃(i)/10⌋⋅10^(-k0+1) strip-zeroes
    x(i) z* 10 * <= [ z* strip-zeroes k0 - 1 + ] [
        ! otw. y(ru)  (5.2.2)
        φ β y(ru) :> y(ru)
        y(ru) y(ru) -77 number= [
            ! if tie (5.2.4): y(r)⋅10^k0
            y(ru) odd? [ 1 - ] when
        ] [
            ! otw. y(ru) ≥ x̃(i) : y(ru)⋅10^k0 ? (y(ru)+1)⋅10^k0
            y(ru) x(i) >= [ 1 + ] unless
        ] if k0
    ] if ; inline

: double>s/f/e ( d -- s f e )
    >double< swap [
        [ E_min ] [ [ p set-bit ] [ bias - ] bi* ] if-zero p -
    ] 2keep [ zero? ] [ 1 > ] bi* and
    [ shorter-interval ] [ normal-interval ] if ;