From 4b53916a12cb02595340f881db3d6cb369ad87c5 Mon Sep 17 00:00:00 2001 From: Slava Pestov Date: Wed, 19 Aug 2009 02:32:18 -0500 Subject: [PATCH] math.intervals: tighter interval arithmetic for intervals with infinities --- basis/math/intervals/intervals-tests.factor | 4 ++ basis/math/intervals/intervals.factor | 49 ++++++++++++--------- 2 files changed, 32 insertions(+), 21 deletions(-) diff --git a/basis/math/intervals/intervals-tests.factor b/basis/math/intervals/intervals-tests.factor index 07c3d8fae7..a2bdf6d98f 100644 --- a/basis/math/intervals/intervals-tests.factor +++ b/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 ; diff --git a/basis/math/intervals/intervals.factor b/basis/math/intervals/intervals.factor index 8ea28b2235..99997ab8cb 100755 --- a/basis/math/intervals/intervals.factor +++ b/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 ) - dup [ first fp-nan? ] any? - [ drop [-inf,inf] ] [ - dup first - [ [ endpoint-min ] reduce ] - [ [ endpoint-max ] reduce ] - 2bi - ] if ; +: points>interval ( seq -- interval nan? ) + [ first fp-nan? not ] partition + [ + [ [ ] [ endpoint-min ] map-reduce ] + [ [ ] [ endpoint-max ] map-reduce ] bi + + ] + [ 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 ] } - [ (interval-abs) points>interval ] + { [ 0 over interval-contains? ] [ (interval-abs) { 0 t } suffix points>interval nan-not-ok ] } + [ (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 ;