diff --git a/interval.tal b/interval.tal
new file mode 100644
index 0000000..aefd10c
--- /dev/null
+++ b/interval.tal
@@ -0,0 +1,78 @@
+( 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
+
+@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
+
+( [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 ORA ?&positive
+    ( negative ) x16-negate SWP2 !x16-negate
+    &positive JMP2r
+    &cross SWP2 x16-negate x16-max #0000 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 GTH2 ?&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** )
+
+~fix16.tal