basis/math/intervals/intervals.factor

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