math.intervals: tighter interval arithmetic for intervals with infinities

basis/math/intervals/intervals-tests.factor

@@ -348,6 +348,10 @@ comparison-ops [
 
 [ t ] [ -10 10 [a,b] interval-abs 0 10 [a,b] = ] unit-test
 
+[ t ] [ full-interval interval-abs [0,inf] = ] unit-test
+
+[ t ] [ [0,inf] interval-sq [0,inf] = ] unit-test
+
 ! Test that commutative interval ops really are
 : random-interval-or-empty ( -- obj )
     10 random 0 = [ empty-interval ] [ random-interval ] if ;

basis/math/intervals/intervals.factor

@@ -94,21 +94,25 @@ MEMO: array-capacity-interval ( -- interval )
 : interval>points ( int -- from to )
     [ from>> ] [ to>> ] bi ;
 
-: points>interval ( seq -- interval )
+: points>interval ( seq -- interval nan? )
-    dup [ first fp-nan? ] any?
+    [ first fp-nan? not ] partition
-    [ drop [-inf,inf] ] [
+    [
-        dup first
+        [ [ ] [ endpoint-min ] map-reduce ]
-        [ [ endpoint-min ] reduce ]
+        [ [ ] [ endpoint-max ] map-reduce ] bi
-        [ [ endpoint-max ] reduce ]
+        <interval>
-        2bi <interval>
+    ]
-    ] if ;
+    [ empty? not ]
+    bi* ;
+
+: nan-ok ( interval nan? -- interval ) drop ; inline
+: nan-not-ok ( interval nan? -- interval ) [ drop full-interval ] when ; inline
 
 : (interval-op) ( p1 p2 quot -- p3 )
     [ [ first ] [ first ] [ call ] tri* ]
     [ drop [ second ] both? ]
     3bi 2array ; inline
 
-: interval-op ( i1 i2 quot -- i3 )
+: interval-op ( i1 i2 quot -- i3 nan? )
     {
         [ [ from>> ] [ from>> ] [ ] tri* (interval-op) ]
         [ [ to>>   ] [ from>> ] [ ] tri* (interval-op) ]
@@ -126,10 +130,10 @@ MEMO: array-capacity-interval ( -- interval )
     } cond ; inline
 
 : interval+ ( i1 i2 -- i3 )
-    [ [ + ] interval-op ] do-empty-interval ;
+    [ [ + ] interval-op nan-ok ] do-empty-interval ;
 
 : interval- ( i1 i2 -- i3 )
-    [ [ - ] interval-op ] do-empty-interval ;
+    [ [ - ] interval-op nan-ok ] do-empty-interval ;
 
 : interval-intersect ( i1 i2 -- i3 )
     {
@@ -154,7 +158,7 @@ MEMO: array-capacity-interval ( -- interval )
         { [ dup empty-interval eq? ] [ drop ] }
         { [ over full-interval eq? ] [ drop ] }
         { [ dup full-interval eq? ] [ nip ] }
-        [ [ interval>points 2array ] bi@ append points>interval ]
+        [ [ interval>points 2array ] bi@ append points>interval nan-not-ok ]
     } cond ;
 
 : interval-subset? ( i1 i2 -- ? )
@@ -173,7 +177,7 @@ MEMO: array-capacity-interval ( -- interval )
     0 swap interval-contains? ;
 
 : interval* ( i1 i2 -- i3 )
-    [ [ [ * ] interval-op ] do-empty-interval ]
+    [ [ [ * ] interval-op nan-ok ] do-empty-interval ]
     [ [ interval-zero? ] either? ]
     2bi [ 0 [a,a] interval-union ] when ;
 
@@ -220,7 +224,7 @@ MEMO: array-capacity-interval ( -- interval )
     [
         [
             [ interval-closure ] bi@
-            [ shift ] interval-op
+            [ shift ] interval-op nan-not-ok
         ] interval-integer-op
     ] do-empty-interval ;
 
@@ -235,11 +239,11 @@ MEMO: array-capacity-interval ( -- interval )
 
 : interval-max ( i1 i2 -- i3 )
     #! Inaccurate; could be tighter
-    [ [ interval-closure ] bi@ [ max ] interval-op ] do-empty-interval ;
+    [ [ interval-closure ] bi@ [ max ] interval-op nan-not-ok ] do-empty-interval ;
 
 : interval-min ( i1 i2 -- i3 )
     #! Inaccurate; could be tighter
-    [ [ interval-closure ] bi@ [ min ] interval-op ] do-empty-interval ;
+    [ [ interval-closure ] bi@ [ min ] interval-op nan-not-ok ] do-empty-interval ;
 
 : interval-interior ( i1 -- i2 )
     dup special-interval? [
@@ -254,7 +258,7 @@ MEMO: array-capacity-interval ( -- interval )
     } cond ; inline
 
 : interval/ ( i1 i2 -- i3 )
-    [ [ [ / ] interval-op ] interval-division-op ] do-empty-interval ;
+    [ [ [ / ] interval-op nan-not-ok ] interval-division-op ] do-empty-interval ;
 
 : interval/-safe ( i1 i2 -- i3 )
     #! Just a hack to make the compiler work if bootstrap.math
@@ -266,13 +270,13 @@ MEMO: array-capacity-interval ( -- interval )
         [
             [
                 [ interval-closure ] bi@
-                [ /i ] interval-op
+                [ /i ] interval-op nan-not-ok
             ] interval-integer-op
         ] interval-division-op
     ] do-empty-interval ;
 
 : interval/f ( i1 i2 -- i3 )
-    [ [ [ /f ] interval-op ] interval-division-op ] do-empty-interval ;
+    [ [ [ /f ] interval-op nan-not-ok ] interval-division-op ] do-empty-interval ;
 
 : (interval-abs) ( i1 -- i2 )
     interval>points [ first2 [ abs ] dip 2array ] bi@ 2array ;
@@ -281,10 +285,13 @@ MEMO: array-capacity-interval ( -- interval )
     {
         { [ dup empty-interval eq? ] [ ] }
         { [ dup full-interval eq? ] [ drop [0,inf] ] }
-        { [ 0 over interval-contains? ] [ (interval-abs) { 0 t } suffix points>interval ] }
+        { [ 0 over interval-contains? ] [ (interval-abs) { 0 t } suffix points>interval nan-not-ok ] }
-        [ (interval-abs) points>interval ]
+        [ (interval-abs) points>interval nan-not-ok ]
     } cond ;
 
+: interval-absq ( i1 -- i2 )
+    interval-abs interval-sq ;
+
 : interval-recip ( i1 -- i2 ) 1 [a,a] swap interval/ ;
 
 : interval-2/ ( i1 -- i2 ) -1 [a,a] interval-shift ;