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