core/math/integers/integers.factor

   1 ! Copyright (C) 2004, 2010 Slava Pestov.
   2 ! Copyright (C) 2008, Doug Coleman.
   3 ! See https://factorcode.org/license.txt for BSD license.
   4 USING: combinators kernel kernel.private math math.order
   5 math.private ;
   6 IN: math.integers
   7
   8 <PRIVATE
   9
  10 : fixnum-min ( x y -- z ) [ fixnum< ] most ; foldable
  11 : fixnum-max ( x y -- z ) [ fixnum> ] most ; foldable
  12
  13 M: integer numerator ; inline
  14 M: integer denominator drop 1 ; inline
  15 M: integer >fraction 1 ; inline
  16
  17 M: fixnum >fixnum ; inline
  18 M: fixnum >bignum fixnum>bignum ; inline
  19 M: fixnum >integer ; inline
  20 M: fixnum >float fixnum>float ; inline
  21 M: fixnum integer>fixnum ; inline
  22 M: fixnum integer>fixnum-strict ; inline
  23
  24 M: fixnum hashcode* nip ; inline
  25 M: fixnum equal? over bignum? [ >bignum bignum= ] [ 2drop f ] if ; inline
  26 M: fixnum number= eq? ; inline
  27
  28 M: fixnum < fixnum< ; inline
  29 M: fixnum <= fixnum<= ; inline
  30 M: fixnum > fixnum> ; inline
  31 M: fixnum >= fixnum>= ; inline
  32
  33 M: fixnum u< fixnum< ; inline
  34 M: fixnum u<= fixnum<= ; inline
  35 M: fixnum u> fixnum> ; inline
  36 M: fixnum u>= fixnum>= ; inline
  37
  38 M: fixnum min over fixnum? [ fixnum-min ] [ call-next-method ] if ; inline
  39 M: fixnum max over fixnum? [ fixnum-max ] [ call-next-method ] if ; inline
  40
  41 M: fixnum + fixnum+ ; inline
  42 M: fixnum - fixnum- ; inline
  43 M: fixnum * fixnum* ; inline
  44 M: fixnum /i fixnum/i ; inline
  45
  46 M: fixnum mod fixnum-mod ; inline
  47
  48 M: fixnum /mod fixnum/mod ; inline
  49
  50 M: fixnum bitand fixnum-bitand ; inline
  51 M: fixnum bitor fixnum-bitor ; inline
  52 M: fixnum bitxor fixnum-bitxor ; inline
  53 M: fixnum shift integer>fixnum fixnum-shift ; inline
  54
  55 M: fixnum bitnot fixnum-bitnot ; inline
  56
  57 : fixnum-bit? ( x n -- ? )
  58     { fixnum fixnum } declare
  59     dup 0 >= [ neg shift even? not ] [ 2drop f ] if ;
  60
  61 M: fixnum bit? integer>fixnum-strict fixnum-bit? ; inline
  62
  63 : fixnum-log2 ( x -- n )
  64     { fixnum } declare
  65     0 swap [ dup 1 eq? ] [
  66         [ 1 fixnum+fast ] [ 2/ ] bi*
  67     ] until drop ;
  68
  69 M: fixnum (log2) fixnum-log2 { fixnum } declare ; inline
  70
  71 M: bignum >fixnum bignum>fixnum ; inline
  72 M: bignum >bignum ; inline
  73 M: bignum integer>fixnum bignum>fixnum ; inline
  74 M: bignum integer>fixnum-strict bignum>fixnum-strict ; inline
  75
  76 M: bignum hashcode* nip bignum>fixnum ;
  77
  78 M: bignum equal?
  79     over bignum? [ bignum= ] [
  80         swap dup fixnum? [ >bignum bignum= ] [ 2drop f ] if
  81     ] if ; inline
  82
  83 M: bignum number= bignum= ; inline
  84
  85 M: bignum < bignum< ; inline
  86 M: bignum <= bignum<= ; inline
  87 M: bignum > bignum> ; inline
  88 M: bignum >= bignum>= ; inline
  89
  90 M: bignum u< bignum< ; inline
  91 M: bignum u<= bignum<= ; inline
  92 M: bignum u> bignum> ; inline
  93 M: bignum u>= bignum>= ; inline
  94
  95 M: bignum + bignum+ ; inline
  96 M: bignum - bignum- ; inline
  97 M: bignum * bignum* ; inline
  98 M: bignum /i bignum/i ; inline
  99 M: bignum mod bignum-mod ; inline
 100
 101 M: bignum /mod bignum/mod ; inline
 102
 103 M: bignum bitand bignum-bitand ; inline
 104 M: bignum bitor bignum-bitor ; inline
 105 M: bignum bitxor bignum-bitxor ; inline
 106 M: bignum shift integer>fixnum bignum-shift ; inline
 107
 108 M: bignum bitnot bignum-bitnot ; inline
 109 M: bignum bit? bignum-bit? ; inline
 110 M: bignum (log2) bignum-log2 ; inline
 111
 112 ! Converting ratios to floats. Based on FLOAT-RATIO from
 113 ! sbcl/src/code/float.lisp, which has the following license:
 114
 115 ! "The software is in the public domain and is
 116 ! provided with absolutely no warranty."
 117
 118 ! First step: pre-scaling
 119 ! As an optimization to minimize the size of the operands of the bignum
 120 ! divisions we will do, we start by stripping any trailing zeros from
 121 ! the denominator and move it into the scale factor.
 122 ! We want a 54 bit result, starting with leading 1, followed by
 123 ! the 52 bit mantissa and then a guard bit: 1mmmmmmmmmm...mmmmmmmmmmmg
 124 ! So we shift the numerator to get the result of the integer division
 125 ! "num/den" in the range ]2^54; 2^53]; Our shift is only a guess
 126 ! based on the magnitude of the inputs, so it
 127 ! will actually give results in the range ]2^55; 2^53].
 128 ! Note: epsilon is used for rounding in step 3.
 129 : twos ( x -- y ) dup 1 - bitxor log2 ; inline
 130
 131 : scale-denonimator ( den -- scaled-den scale' )
 132     dup twos neg [ shift ] keep ; inline
 133
 134 : (epsilon?) ( num shift -- ? )
 135     dup neg? [ neg 2^ 1 - bitand zero? not ] [ 2drop f ] if ; inline
 136
 137 : scale-numerator ( num den -- epsilon? num' scale )
 138     over [ log2 ] bi@ - [
 139         54 + [ (epsilon?) ] [ shift ] 2bi
 140     ] keep ; inline
 141
 142 : pre-scale ( num den -- epsilon? num' den' scale )
 143     scale-denonimator [
 144         [ scale-numerator ] keep swap
 145     ] dip swap - ; inline
 146
 147 ! Second step: compute mantissa
 148 ! "num/den" would be in the range ]2^55; 2^53]. After this step
 149 ! it will be in the range ]2^54; 2^53]. Compute "num/den" and the
 150 ! reminder used for rounding
 151 ! For subnormals, after we know the final value of the exponent,
 152 ! we shift the numerator again to get the correct precision.
 153 ! We do it before rounding so that subnormals are correctly rounded.
 154 : (2/-with-epsilon) ( epsilon? num -- epsilon?' num' )
 155     [ 1 bitand zero? not or ] [ 2/ ] bi ; inline
 156
 157 : (shift-with-epsilon) ( epsilon? num den scale -- epsilon?' num' den scale )
 158     [
 159         nip 1021 +
 160         [ neg 2^ 1 - bitand zero? not or ] [ shift ] 2bi
 161     ] 2keep ; inline
 162
 163 : mantissa-and-guard ( epsilon? num den scale -- epsilon?' mantissa-and-guard rem scale' )
 164     2over /i log2 53 >
 165     [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] when
 166     ! At this point, the scale value is the exponent minus 1.
 167     dup -1021 < [ (shift-with-epsilon) ] when
 168     [ /mod ] dip ; inline
 169
 170 ! Third step: rounding
 171 !
 172 ! if the guard bit is 0, round down
 173 ! else if the guard bit is 1 and (rem != 0 or epsilon is true), round up
 174 ! else break the tie by alternating rounding down or up to avoid accumulating errors
 175 !
 176 ! The epsilon trick works because epsilon is true if numerator bits were discarded.
 177 ! Mathematically, (num+epsilon)/denom = (num/denum) + (epsilon/denom)
 178 ! We have actually computed the "num/denum" part and use the "epsilon/denom"
 179 ! to choose the correct rounding.
 180 !
 181 ! Note that rounding down means doing nothing because we will
 182 ! discard the guard bit after this
 183 : round-to-nearest ( epsilon? mantissa-and-guard rem -- mantissa-and-guard' )
 184     over odd?
 185     [
 186         zero? [
 187             dup 2 bitand zero? not rot or [ 1 + ] when
 188         ] [ nip 1 + ] if
 189     ] [ drop nip ] if ; inline
 190
 191 ! Fourth step: post-scaling
 192 ! Because of rounding, our mantissa with guard bit may have overflowed
 193 ! the 54 bit precision to 2^54 so we have to handle it specially.
 194 ! For subnormals, the rounding may also have overflowed the precision,
 195 ! but the overflowed value is actually the correct value by chance
 196 ! (even in the case when the biggest subnormal is rounded up to
 197 ! the smallest normal float) because we interpret it directly
 198 ! as the bits of the resulting double.
 199 : scale-float ( mantissa scale -- float' )
 200     ! the scale value is the exponent minus 1.
 201     {
 202         { [ dup 1024 > ] [ 2drop 1/0. ] }
 203         { [ dup -1021 < ] [ drop bits>double ] } ! subnormals and underflow
 204         [ [ 52 2^ 1 - bitand ] dip 1022 + 52 shift bitor bits>double ]
 205     } cond ; inline
 206
 207 : post-scale ( mantissa scale -- n )
 208     [ 2/ ] dip ! drop guard bit
 209     over 53 2^ = [ [ 2/ ] [ 1 + ] bi* ] when
 210     scale-float ; inline
 211
 212 ! Main word
 213 : /f-abs ( m n -- f )
 214     over zero? [ nip zero? 0/0. 0.0 ? ] [
 215         [ drop 1/0. ] [
 216             pre-scale
 217             mantissa-and-guard
 218             [ round-to-nearest ] dip
 219             post-scale
 220         ] if-zero
 221     ] if ; inline
 222
 223 : bignum/f ( m n -- f )
 224     [ [ abs ] bi@ /f-abs ] [ [ 0 < ] bi@ xor ] 2bi [ neg ] when ; inline
 225
 226 M: bignum /f { bignum bignum } declare bignum/f ;
 227
 228 CONSTANT: bignum/f-threshold 0x20,0000,0000,0000
 229
 230 : fixnum/f ( m n -- m/n )
 231     [ >float ] bi@ float/f ; inline
 232
 233 M: fixnum /f
 234     { fixnum fixnum } declare
 235     2dup [ abs bignum/f-threshold >= ] either?
 236     [ bignum/f ] [ fixnum/f ] if ; inline
 237
 238 : bignum>float ( bignum -- float )
 239     { bignum } declare 1 >bignum bignum/f ;
 240
 241 M: bignum >float bignum>float ; inline
 242
 243 PRIVATE>