basis/math/intervals/intervals.factor

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