basis/math/intervals/intervals.factor

   1 ! Copyright (C) 2007, 2008 Slava Pestov, Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 ! Based on Slate's src/unfinished/interval.slate by Brian Rice.
   4 USING: accessors kernel sequences arrays math math.order
   5 combinators generic layouts ;
   6 IN: math.intervals
   7
   8 SYMBOL: empty-interval
   9
  10 SYMBOL: full-interval
  11
  12 TUPLE: interval { from read-only } { to read-only } ;
  13
  14 : <interval> ( from to -- interval )
  15     2dup [ first ] bi@ {
  16         { [ 2dup > ] [ 2drop 2drop empty-interval ] }
  17         { [ 2dup = ] [
  18             2drop 2dup [ second ] both?
  19             [ interval boa ] [ 2drop empty-interval ] if
  20         ] }
  21         [ 2drop interval boa ]
  22     } cond ;
  23
  24 : open-point ( n -- endpoint ) f 2array ;
  25
  26 : closed-point ( n -- endpoint ) t 2array ;
  27
  28 : [a,b] ( a b -- interval )
  29     [ closed-point ] dip closed-point <interval> ; foldable
  30
  31 : (a,b) ( a b -- interval )
  32     [ open-point ] dip open-point <interval> ; foldable
  33
  34 : [a,b) ( a b -- interval )
  35     [ closed-point ] dip open-point <interval> ; foldable
  36
  37 : (a,b] ( a b -- interval )
  38     [ open-point ] dip closed-point <interval> ; foldable
  39
  40 : [a,a] ( a -- interval )
  41     closed-point dup <interval> ; foldable
  42
  43 : [-inf,a] ( a -- interval ) -1./0. swap [a,b] ; inline
  44
  45 : [-inf,a) ( a -- interval ) -1./0. swap [a,b) ; inline
  46
  47 : [a,inf] ( a -- interval ) 1./0. [a,b] ; inline
  48
  49 : (a,inf] ( a -- interval ) 1./0. (a,b] ; inline
  50
  51 : [-inf,inf] ( -- interval ) full-interval ; inline
  52
  53 : compare-endpoints ( p1 p2 quot -- ? )
  54     [ 2dup [ first ] bi@ ] dip call [
  55         2drop t
  56     ] [
  57         2dup [ first ] bi@ = [
  58             [ second ] bi@ not or
  59         ] [
  60             2drop f
  61         ] if
  62     ] if ; inline
  63
  64 : endpoint< ( p1 p2 -- ? ) [ < ] compare-endpoints ;
  65
  66 : endpoint<= ( p1 p2 -- ? ) [ endpoint< ] 2keep = or ;
  67
  68 : endpoint> ( p1 p2 -- ? ) [ > ] compare-endpoints ;
  69
  70 : endpoint>= ( p1 p2 -- ? ) [ endpoint> ] 2keep = or ;
  71
  72 : endpoint-min ( p1 p2 -- p3 ) [ endpoint< ] most ;
  73
  74 : endpoint-max ( p1 p2 -- p3 ) [ endpoint> ] most ;
  75
  76 : interval>points ( int -- from to )
  77     [ from>> ] [ to>> ] bi ;
  78
  79 : points>interval ( seq -- interval )
  80     dup [ first fp-nan? ] any?
  81     [ drop [-inf,inf] ] [
  82         dup first
  83         [ [ endpoint-min ] reduce ]
  84         [ [ endpoint-max ] reduce ]
  85         2bi <interval>
  86     ] if ;
  87
  88 : (interval-op) ( p1 p2 quot -- p3 )
  89     [ [ first ] [ first ] [ call ] tri* ]
  90     [ drop [ second ] both? ]
  91     3bi 2array ; inline
  92
  93 : interval-op ( i1 i2 quot -- i3 )
  94     {
  95         [ [ from>> ] [ from>> ] [ ] tri* (interval-op) ]
  96         [ [ to>>   ] [ from>> ] [ ] tri* (interval-op) ]
  97         [ [ to>>   ] [ to>>   ] [ ] tri* (interval-op) ]
  98         [ [ from>> ] [ to>>   ] [ ] tri* (interval-op) ]
  99     } 3cleave 4array points>interval ; inline
 100
 101 : do-empty-interval ( i1 i2 quot -- i3 )
 102     {
 103         { [ pick empty-interval eq? ] [ 2drop ] }
 104         { [ over empty-interval eq? ] [ drop nip ] }
 105         { [ pick full-interval eq? ] [ 2drop ] }
 106         { [ over full-interval eq? ] [ drop nip ] }
 107         [ call ]
 108     } cond ; inline
 109
 110 : interval+ ( i1 i2 -- i3 )
 111     [ [ + ] interval-op ] do-empty-interval ;
 112
 113 : interval- ( i1 i2 -- i3 )
 114     [ [ - ] interval-op ] do-empty-interval ;
 115
 116 : interval-intersect ( i1 i2 -- i3 )
 117     {
 118         { [ over empty-interval eq? ] [ drop ] }
 119         { [ dup empty-interval eq? ] [ nip ] }
 120         { [ over full-interval eq? ] [ nip ] }
 121         { [ dup full-interval eq? ] [ drop ] }
 122         [
 123             [ interval>points ] bi@ swapd
 124             [ [ swap endpoint< ] most ]
 125             [ [ swap endpoint> ] most ] 2bi*
 126             <interval>
 127         ]
 128     } cond ;
 129
 130 : intervals-intersect? ( i1 i2 -- ? )
 131     interval-intersect empty-interval eq? not ;
 132
 133 : interval-union ( i1 i2 -- i3 )
 134     {
 135         { [ over empty-interval eq? ] [ nip ] }
 136         { [ dup empty-interval eq? ] [ drop ] }
 137         { [ over full-interval eq? ] [ drop ] }
 138         { [ dup full-interval eq? ] [ nip ] }
 139         [ [ interval>points 2array ] bi@ append points>interval ]
 140     } cond ;
 141
 142 : interval-subset? ( i1 i2 -- ? )
 143     dupd interval-intersect = ;
 144
 145 : interval-contains? ( x int -- ? )
 146     dup empty-interval eq? [ 2drop f ] [
 147         dup full-interval eq? [ 2drop t ] [
 148             [ from>> first2 [ >= ] [ > ] if ]
 149             [ to>>   first2 [ <= ] [ < ] if ]
 150             2bi and
 151         ] if
 152     ] if ;
 153
 154 : interval-zero? ( int -- ? )
 155     0 swap interval-contains? ;
 156
 157 : interval* ( i1 i2 -- i3 )
 158     [ [ [ * ] interval-op ] do-empty-interval ]
 159     [ [ interval-zero? ] either? ]
 160     2bi [ 0 [a,a] interval-union ] when ;
 161
 162 : interval-1+ ( i1 -- i2 ) 1 [a,a] interval+ ;
 163
 164 : interval-1- ( i1 -- i2 ) -1 [a,a] interval+ ;
 165
 166 : interval-neg ( i1 -- i2 ) -1 [a,a] interval* ;
 167
 168 : interval-bitnot ( i1 -- i2 ) interval-neg interval-1- ;
 169
 170 : interval-sq ( i1 -- i2 ) dup interval* ;
 171
 172 : special-interval? ( interval -- ? )
 173     { empty-interval full-interval } memq? ;
 174
 175 : interval-singleton? ( int -- ? )
 176     dup special-interval? [
 177         drop f
 178     ] [
 179         interval>points
 180         2dup [ second ] both?
 181         [ [ first ] bi@ = ]
 182         [ 2drop f ] if
 183     ] if ;
 184
 185 : interval-length ( int -- n )
 186     {
 187         { [ dup empty-interval eq? ] [ drop 0 ] }
 188         { [ dup full-interval eq? ] [ drop 1/0. ] }
 189         [ interval>points [ first ] bi@ swap - ]
 190     } cond ;
 191
 192 : interval-closure ( i1 -- i2 )
 193     dup [ interval>points [ first ] bi@ [a,b] ] when ;
 194
 195 : interval-integer-op ( i1 i2 quot -- i3 )
 196     [
 197         2dup [ interval>points [ first integer? ] both? ] both?
 198     ] dip [ 2drop [-inf,inf] ] if ; inline
 199
 200 : interval-shift ( i1 i2 -- i3 )
 201     #! Inaccurate; could be tighter
 202     [
 203         [
 204             [ interval-closure ] bi@
 205             [ shift ] interval-op
 206         ] interval-integer-op
 207     ] do-empty-interval ;
 208
 209 : interval-shift-safe ( i1 i2 -- i3 )
 210     [
 211         dup to>> first 100 > [
 212             2drop [-inf,inf]
 213         ] [
 214             interval-shift
 215         ] if
 216     ] do-empty-interval ;
 217
 218 : interval-max ( i1 i2 -- i3 )
 219     #! Inaccurate; could be tighter
 220     [ [ interval-closure ] bi@ [ max ] interval-op ] do-empty-interval ;
 221
 222 : interval-min ( i1 i2 -- i3 )
 223     #! Inaccurate; could be tighter
 224     [ [ interval-closure ] bi@ [ min ] interval-op ] do-empty-interval ;
 225
 226 : interval-interior ( i1 -- i2 )
 227     dup special-interval? [
 228         interval>points [ first ] bi@ (a,b)
 229     ] unless ;
 230
 231 : interval-division-op ( i1 i2 quot -- i3 )
 232     {
 233         { [ 0 pick interval-closure interval-contains? ] [ 3drop [-inf,inf] ] }
 234         { [ pick interval-zero? ] [ call 0 [a,a] interval-union ] }
 235         [ call ]
 236     } cond ; inline
 237
 238 : interval/ ( i1 i2 -- i3 )
 239     [ [ [ / ] interval-op ] interval-division-op ] do-empty-interval ;
 240
 241 : interval/-safe ( i1 i2 -- i3 )
 242     #! Just a hack to make the compiler work if bootstrap.math
 243     #! is not loaded.
 244     \ integer \ / method [ interval/ ] [ 2drop f ] if ;
 245
 246 : interval/i ( i1 i2 -- i3 )
 247     [
 248         [
 249             [
 250                 [ interval-closure ] bi@
 251                 [ /i ] interval-op
 252             ] interval-integer-op
 253         ] interval-division-op
 254     ] do-empty-interval ;
 255
 256 : interval/f ( i1 i2 -- i3 )
 257     [ [ [ /f ] interval-op ] interval-division-op ] do-empty-interval ;
 258
 259 : (interval-abs) ( i1 -- i2 )
 260     interval>points [ first2 [ abs ] dip 2array ] bi@ 2array ;
 261
 262 : interval-abs ( i1 -- i2 )
 263     {
 264         { [ dup empty-interval eq? ] [ ] }
 265         { [ dup full-interval eq? ] [ drop 0 [a,inf] ] }
 266         { [ 0 over interval-contains? ] [ (interval-abs) { 0 t } suffix points>interval ] }
 267         [ (interval-abs) points>interval ]
 268     } cond ;
 269
 270 : interval-mod ( i1 i2 -- i3 )
 271     #! Inaccurate.
 272     [
 273         [
 274             nip interval-abs to>> first [ neg ] keep (a,b)
 275         ] interval-division-op
 276     ] do-empty-interval ;
 277
 278 : interval-rem ( i1 i2 -- i3 )
 279     #! Inaccurate.
 280     [
 281         [
 282             nip interval-abs to>> first 0 swap [a,b)
 283         ] interval-division-op
 284     ] do-empty-interval ;
 285
 286 : interval-recip ( i1 -- i2 ) 1 [a,a] swap interval/ ;
 287
 288 : interval-2/ ( i1 -- i2 ) -1 [a,a] interval-shift ;
 289
 290 SYMBOL: incomparable
 291
 292 : left-endpoint-< ( i1 i2 -- ? )
 293     [ swap interval-subset? ]
 294     [ nip interval-singleton? ]
 295     [ [ from>> ] bi@ = ]
 296     2tri and and ;
 297
 298 : right-endpoint-< ( i1 i2 -- ? )
 299     [ interval-subset? ]
 300     [ drop interval-singleton? ]
 301     [ [ to>> ] bi@ = ]
 302     2tri and and ;
 303
 304 : (interval<) ( i1 i2 -- i1 i2 ? )
 305     2dup [ from>> ] bi@ endpoint< ;
 306
 307 : interval< ( i1 i2 -- ? )
 308     {
 309         { [ 2dup [ special-interval? ] either? ] [ incomparable ] }
 310         { [ 2dup interval-intersect empty-interval eq? ] [ (interval<) ] }
 311         { [ 2dup left-endpoint-< ] [ f ] }
 312         { [ 2dup right-endpoint-< ] [ f ] }
 313         [ incomparable ]
 314     } cond 2nip ;
 315
 316 : left-endpoint-<= ( i1 i2 -- ? )
 317     [ from>> ] dip to>> = ;
 318
 319 : right-endpoint-<= ( i1 i2 -- ? )
 320     [ to>> ] dip from>> = ;
 321
 322 : interval<= ( i1 i2 -- ? )
 323     {
 324         { [ 2dup [ special-interval? ] either? ] [ incomparable ] }
 325         { [ 2dup interval-intersect empty-interval eq? ] [ (interval<) ] }
 326         { [ 2dup right-endpoint-<= ] [ t ] }
 327         [ incomparable ]
 328     } cond 2nip ;
 329
 330 : interval> ( i1 i2 -- ? )
 331     swap interval< ;
 332
 333 : interval>= ( i1 i2 -- ? )
 334     swap interval<= ;
 335
 336 : interval-bitand-pos ( i1 i2 -- ? )
 337     [ to>> first ] bi@ min 0 swap [a,b] ;
 338
 339 : interval-bitand-neg ( i1 i2 -- ? )
 340     dup from>> first 0 < [ drop ] [ nip ] if
 341     0 swap to>> first [a,b] ;
 342
 343 : interval-nonnegative? ( i -- ? )
 344     from>> first 0 >= ;
 345
 346 : interval-bitand ( i1 i2 -- i3 )
 347     #! Inaccurate.
 348     [
 349         {
 350             {
 351                 [ 2dup [ interval-nonnegative? ] both? ]
 352                 [ interval-bitand-pos ]
 353             }
 354             {
 355                 [ 2dup [ interval-nonnegative? ] either? ]
 356                 [ interval-bitand-neg ]
 357             }
 358             [ 2drop [-inf,inf] ]
 359         } cond
 360     ] do-empty-interval ;
 361
 362 : interval-bitor ( i1 i2 -- i3 )
 363     #! Inaccurate.
 364     [
 365         2dup [ interval-nonnegative? ] both?
 366         [
 367             [ interval>points [ first ] bi@ ] bi@
 368             4array supremum 0 swap >integer next-power-of-2 [a,b]
 369         ] [ 2drop [-inf,inf] ] if
 370     ] do-empty-interval ;
 371
 372 : interval-bitxor ( i1 i2 -- i3 )
 373     #! Inaccurate.
 374     interval-bitor ;
 375
 376 : interval-log2 ( i1 -- i2 )
 377     {
 378         { empty-interval [ empty-interval ] }
 379         { full-interval [ 0 [a,inf] ] }
 380         [
 381             to>> first 1 max dup most-positive-fixnum >
 382             [ drop full-interval interval-log2 ]
 383             [ 1+ >integer log2 0 swap [a,b] ]
 384             if
 385         ]
 386     } case ;
 387
 388 : assume< ( i1 i2 -- i3 )
 389     dup special-interval? [ drop ] [
 390         to>> first [-inf,a) interval-intersect
 391     ] if ;
 392
 393 : assume<= ( i1 i2 -- i3 )
 394     dup special-interval? [ drop ] [
 395         to>> first [-inf,a] interval-intersect
 396     ] if ;
 397
 398 : assume> ( i1 i2 -- i3 )
 399     dup special-interval? [ drop ] [
 400         from>> first (a,inf] interval-intersect
 401     ] if ;
 402
 403 : assume>= ( i1 i2 -- i3 )
 404     dup special-interval? [ drop ] [
 405         from>> first [a,inf] interval-intersect
 406     ] if ;
 407
 408 : integral-closure ( i1 -- i2 )
 409     dup special-interval? [
 410         [ from>> first2 [ 1+ ] unless ]
 411         [ to>> first2 [ 1- ] unless ]
 412         bi [a,b]
 413     ] unless ;