( interval.tal )
( )
( interval arithmetic over 8.8 fixed point numbers. )
( )
( intervals are closed, i.e. [lower, upper]. points )
( are represented as [x, x], and the empty interval )
( is not representable. )

@iv-point ( x* -> x** ) DUP2 JMP2r

@iv-is-negative ( x** -> x<0^ ) NIP2 #0000 !x16-lt
@iv-is-positive ( x** -> x<0^ ) POP2 #0000 !x16-gt
@iv-not-positive ( x** -> x<0^ ) NIP2 #0000 !x16-lteq
@iv-not-negative ( x** -> x<0^ ) POP2 #0000 !x16-gteq

( +x / [-a, +b] -> [ -x/a, +x/b] )
( -x / [-a, +b] -> [ -x/b, +x/a] )
( +x / [-a,  0] -> [-256x, -x/a] )
( -x / [-a,  0] -> [  x/a, 256x] )
( +x / [ 0, +b] -> [  x/b, 256x] )
( -x / [ 0, +b] -> [-256x, -x/b] )
( +x / [+a, +b] -> [  x/b,  x/a] )
( -x / [+a, +b] -> [ -x/a, -x/b] )
( +x / [-a, -b] -> [ -x/a, -x/b] )
( -x / [-a, -b] -> [  x/b,  x/a] )
(  x / [ 0,  0] -> error         )

(  0 / y        -> [0, 0] )
( +x / [-a, -b] -> [-x/b, -x/a] )
( +x / [+a, +b] -> [+x/b, +x/a] )
( -x / [-a, -b] -> [+x/a, +x/b] )
( -x / [+a, +b] -> [-x/b, -x/a] )


@iv-crosses-zero ( x** -> bool^ )
    #0000 SWP2 OVR2 x16-lteq    ( x0* 0* x1<0^ )
    ?&no !x16-lt                ( x0<0^ )
    &no POP2 POP2 #00 JMP2r     ( 0^ )

@iv-is-point ( x** -> bool^ )
    EQU2 JMP2r

( )
@iv-emit ( x** -> )
    [ LIT "[ #18 DEO ]
    SWP2 x16-emit
    [ LIT ", #18 DEO #20 #18 DEO ]
    x16-emit
    [ LIT "] #18 DEO ]
    JMP2r

( [x0, x1] | [y0, y1] -> [min[x0,y0], max[x1,y1]] )
@iv-union ( x** y** -> x|y** )
    ROT2 x16-max STH2
    x16-min STH2r JMP2r

( [x0, x1] intersects [y0, y1] if x1 >= y0 and x0 <= y1. )
@iv-intersects ( x** y** -> bool^ )
    iv-disjoint #00 EQU JMP2r

( [x0, x1] is disjoint with [y0, y1] if x1 < y0 or x0 > y1. )
@iv-disjoint ( x** y** -> bool^ )
    STH2 ROT2                   ( x1* y0* x0* [y1*] )
    STH2r x16-gt STH            ( x1* y0* [x0>y1^] )
    x16-lt STHr ORA JMP2r       ( x1<y0|x0>y1^ )

( [x0, x1] -> [-x1, -x0] )
@iv-negate ( x** -> -x** )
    x16-negate SWP2 !x16-negate

@iv-abs ( x** -> |x|** )
    OVR2 OVR2 iv-crosses-zero ?&cross
    DUP2 #8000 LTH2 ?&positive
    ( negative ) x16-negate SWP2 !x16-negate
    &positive JMP2r
    &cross SWP2 x16-negate x16-max #0000 SWP2 JMP2r

@iv-scalar-add ( x** y* -> x+y** )
    STH2k x16-add SWP2
    STH2r x16-add SWP2 JMP2r

( [x0, x1] + [y0, y1] -> [x0+y0, x1+y1] )
@iv-add ( x** y** -> x+y** )
    ROT2 x16-add STH2           ( x0* y0* [x1+y1*] )
    x16-add STH2r JMP2r         ( x0+y0* x1+y1* )

( [x0, x1] - [y0, y1] -> [x0-y1, x1-y0] )
@iv-sub ( x** y** -> x-y** )
    STH2 x16-sub SWP2           ( x1-y0* x0* [y1*] )
    STH2r x16-sub SWP2 JMP2r    ( x0-y1* x1-y0* )

( [x0, x1] * (+y) -> [x0y, x1y] )
( [x0, x1] * (-y) -> [x1y, x0y] )
( [x0, x1] * 0    -> [0, 0]     )
@iv-scalar-mul ( x** y* -> xy** )
    STH2k #0000 x16-gt ?&positive ( x0* x1* [y*] )
    SWP2 &positive                ( a* b* [y*] )
    STH2kr x16-mul SWP2           ( by* a* [y*] )
    STH2r  x16-mul SWP2 JMP2r     ( ay* by* )

( [x0, x1] * [y0, y1] -> [min[x0y0, x0y1], max[x1y0, x1y1]] )
@iv-mul ( x** y** -> xy** )
    ROT2k                       ( x0* x1* y0* y1* y0* y1* x1* )
    iv-scalar-mul               ( x0* x1* y0* y1* a0* a1* )
    STH2 STH2                   ( x0* x1* y0* y1* [a1* a0*] )
    ROT2 POP2 ROT2              ( y0* y1* x0* [a1* a0*] )
    iv-scalar-mul STH2r STH2r   ( b0* b1* a0* a1* )
    !iv-union                   ( a|b** )

@iv-square ( x** -> xx** )
    iv-abs OVR2 OVR2 !iv-mul

@iv-pow ( x** k^ -> x^k** )
    DUP #00 EQU ?&one
    DUP #01 EQU ?&id
    DUP #02 EQU ?&square
    DUP #01 AND ?&odd
    #01 SFT iv-pow OVR2 OVR2 !iv-mul
    &one POP POP2 POP2 #0100 DUP2 JMP2r
    &id POP JMP2r
    &square POP !iv-square
    &odd #01 SUB STH OVR2 OVR2 STHr iv-pow !iv-mul

( FIXME: not reliable )
@iv-sin ( x** -> sin[x]** )
    x16-sin SWP2 x16-sin
    OVR2 OVR2 x16-lt JMP SWP2 JMP2r

( FIXME: not reliable )
@iv-cos ( x** -> sin[x]** )
    x16-cos SWP2 x16-cos
    OVR2 OVR2 x16-lt JMP SWP2 JMP2r

~fix16.tal