1 ! Copyright (C) 2004, 2010 Slava Pestov.
2 ! Copyright (C) 2008, Doug Coleman.
3 ! See http://factorcode.org/license.txt for BSD license.
4 USING: combinators kernel kernel.private math math.order
10 : fixnum-min ( x y -- z ) [ fixnum< ] most ; foldable
11 : fixnum-max ( x y -- z ) [ fixnum> ] most ; foldable
13 M: integer numerator ; inline
14 M: integer denominator drop 1 ; inline
15 M: integer >fraction 1 ; inline
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
24 M: fixnum hashcode* nip ; inline
25 M: fixnum equal? over bignum? [ >bignum bignum= ] [ 2drop f ] if ; inline
26 M: fixnum number= eq? ; inline
28 M: fixnum < fixnum< ; inline
29 M: fixnum <= fixnum<= ; inline
30 M: fixnum > fixnum> ; inline
31 M: fixnum >= fixnum>= ; inline
33 M: fixnum u< fixnum< ; inline
34 M: fixnum u<= fixnum<= ; inline
35 M: fixnum u> fixnum> ; inline
36 M: fixnum u>= fixnum>= ; inline
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
41 M: fixnum + fixnum+ ; inline
42 M: fixnum - fixnum- ; inline
43 M: fixnum * fixnum* ; inline
44 M: fixnum /i fixnum/i ; inline
46 M: fixnum mod fixnum-mod ; inline
48 M: fixnum /mod fixnum/mod ; inline
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
55 M: fixnum bitnot fixnum-bitnot ; inline
57 : fixnum-bit? ( x n -- ? )
58 { fixnum fixnum } declare
59 dup 0 >= [ neg shift even? not ] [ 2drop f ] if ;
61 M: fixnum bit? integer>fixnum-strict fixnum-bit? ; inline
63 : fixnum-log2 ( x -- n )
65 0 swap [ dup 1 eq? ] [
66 [ 1 fixnum+fast ] [ 2/ ] bi*
69 M: fixnum (log2) fixnum-log2 { fixnum } declare ; inline
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
76 M: bignum hashcode* nip bignum>fixnum ;
79 over bignum? [ bignum= ] [
80 swap dup fixnum? [ >bignum bignum= ] [ 2drop f ] if
83 M: bignum number= bignum= ; inline
85 M: bignum < bignum< ; inline
86 M: bignum <= bignum<= ; inline
87 M: bignum > bignum> ; inline
88 M: bignum >= bignum>= ; inline
90 M: bignum u< bignum< ; inline
91 M: bignum u<= bignum<= ; inline
92 M: bignum u> bignum> ; inline
93 M: bignum u>= bignum>= ; inline
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
101 M: bignum /mod bignum/mod ; inline
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
108 M: bignum bitnot bignum-bitnot ; inline
109 M: bignum bit? bignum-bit? ; inline
110 M: bignum (log2) bignum-log2 ; inline
112 ! Converting ratios to floats. Based on FLOAT-RATIO from
113 ! sbcl/src/code/float.lisp, which has the following license:
115 ! "The software is in the public domain and is
116 ! provided with absolutely no warranty."
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 result in ]2^54;2^53] to find the mantissa
123 ! in ]2^53,2^52] with 1 extra "guard" bit for rounding.
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
131 : scale-denonimator ( den -- scaled-den scale' )
132 dup twos neg [ shift ] keep ; inline
134 : (epsilon?) ( num shift -- ? )
135 dup neg? [ neg 2^ 1 - bitand zero? not ] [ 2drop f ] if ; inline
137 : scale-numerator ( num den -- epsilon? num' scale )
138 over [ log2 ] bi@ - [
139 54 + [ (epsilon?) ] [ shift ] 2bi
142 : pre-scale ( num den -- epsilon? num' den' scale )
144 [ scale-numerator ] keep swap
145 ] dip swap - ; inline
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 : (2/-with-epsilon) ( epsilon? num -- epsilon?' num' )
152 [ 1 bitand zero? not or ] [ 2/ ] bi ; inline
154 : mantissa-and-guard ( epsilon? num den scale -- epsilon?' mantissa-and-guard rem scale' )
156 [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] when
157 [ /mod ] dip ; inline
159 ! Third step: rounding
161 ! if the guard bit is 0, round down
162 ! else if the guard bit is 1 and (rem != 0 or epsilon is true), round up
163 ! else break the tie by alternating rounding down or up to avoid accumulating errors
165 ! The epsilon trick works because epsilon is true if numerator bits were discarded.
166 ! Mathematically, (num+epsilon)/denom = (num/denum) + (epsilon/denom)
167 ! We have actually computed the "num/denum" part and use the "epsilon/denom"
168 ! to choose the correct rounding.
170 ! Note that rounding down means doing nothing because we will
171 ! discard the guard bit after this
172 : round-to-nearest ( epsilon? mantissa-and-guard rem -- mantissa-and-guard' )
176 dup 2 bitand zero? not rot or [ 1 + ] when
178 ] [ drop nip ] if ; inline
180 ! Fourth step: post-scaling
181 ! Because of rounding, our mantissa with guard bit is now in the
182 ! range [2^54;2^53], so we have to handle 2^54 specially.
183 : scale-float ( mantissa scale -- float' )
184 ! At this point, the scale value is the exponent minus 1.
186 { [ dup 1024 > ] [ 2drop 1/0. ] }
187 { [ dup -1021 < ] [ 1021 + shift bits>double ] } ! subnormals and underflow
188 [ [ 52 2^ 1 - bitand ] dip 1022 + 52 shift bitor bits>double ]
191 : post-scale ( mantissa scale -- n )
192 [ 2/ ] dip ! drop guard bit
193 over 53 2^ = [ [ 2/ ] [ 1 + ] bi* ] when
197 : /f-abs ( m n -- f )
198 over zero? [ nip zero? 0/0. 0.0 ? ] [
202 [ round-to-nearest ] dip
207 : bignum/f ( m n -- f )
208 [ [ abs ] bi@ /f-abs ] [ [ 0 < ] bi@ xor ] 2bi [ neg ] when ; inline
210 M: bignum /f ( m n -- f ) { bignum bignum } declare bignum/f ;
212 CONSTANT: bignum/f-threshold 0x20,0000,0000,0000
214 : fixnum/f ( m n -- m/n )
215 [ >float ] bi@ float/f ; inline
218 { fixnum fixnum } declare
219 2dup [ abs bignum/f-threshold >= ] either?
220 [ bignum/f ] [ fixnum/f ] if ; inline
222 : bignum>float ( bignum -- float )
223 { bignum } declare 1 >bignum bignum/f ;
225 M: bignum >float bignum>float ; inline