basis/math/intervals/intervals.factor

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