source: project/release/4/numbers/trunk/numbers.scm @ 32321

Last change on this file since 32321 was 32321, checked in by sjamaan, 6 years ago

numbers: Fix calculation of "l" in Burnikel-Ziegler, so it doesn't split the numbers in the wrong way (this is pretty unclear from the paper)

File size: 77.1 KB
Line 
1;;;; numbers.scm
2;
3; Copyright (c) 2008-2015 The CHICKEN Team
4; Copyright (c) 2000-2007, Felix L. Winkelmann
5; All rights reserved.
6;
7; Redistribution and use in source and binary forms, with or without
8; modification, are permitted provided that the following conditions
9; are met:
10; 1. Redistributions of source code must retain the above copyright
11;    notice, this list of conditions and the following disclaimer.
12; 2. Redistributions in binary form must reproduce the above copyright
13;    notice, this list of conditions and the following disclaimer in the
14;    documentation and/or other materials provided with the distribution.
15; 3. The name of the authors may not be used to endorse or promote products
16;    derived from this software without specific prior written permission.
17;
18; THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
19; IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21; IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT,
22; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23; NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24; DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25; THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27; THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28;
29; Abbreviations of paper and book titles used in comments are:
30; [Knuth] Donald E. Knuth, "The Art of Computer Programming", Volume 2
31; [MpNT]  Tiplea at al., "MpNT: A Multi-Precision Number Theory Package"
32; [MCA]   Richard P. Brent & Paul Zimmermann, "Modern Computer Arithmetic"
33
34(declare 
35  (no-bound-checks)
36  (no-procedure-checks))
37
38(module numbers 
39    (+ - * / = > < >= <= eqv?
40       add1 sub1 signum number->string string->number integer-length bit-set?
41       bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift
42       equal? ; From scheme. Structural & bytevector comparisons Just Work
43       exp log sin cos tan atan acos asin conj
44       expt sqrt exact-integer-sqrt exact-integer-nth-root
45       quotient modulo remainder quotient&modulo quotient&remainder
46       floor/ floor-quotient floor-remainder
47       truncate/ truncate-quotient truncate-remainder
48       numerator denominator
49       abs max min fxgcd fpgcd gcd lcm
50       positive? negative? odd? even? zero?
51       exact? inexact?
52       rationalize
53       random randomize
54       floor ceiling truncate round
55       inexact->exact inexact exact->inexact exact
56       number? complex? real? rational? integer? exact-integer?
57       make-rectangular make-polar real-part imag-part magnitude angle
58       bignum? ratnum? cflonum? rectnum? compnum? cintnum? cplxnum?
59       nan? finite? infinite?
60
61       ;; Specialization hooks, to be removed when integrating into core
62       @fixnum-2-plus @integer-2-plus @basic-2-plus
63       @fixnum-negate @integer-negate @basic-negate
64       @fixnum-abs @integer-abs @basic-abs
65       @fixnum-signum @flonum-signum @integer-signum @basic-signum
66       @fixnum-2-minus @integer-2-minus @basic-2-minus
67       @fixnum-2-times @integer-2-times @basic-2-times
68       @fixnum-quotient @flonum-quotient @integer-quotient @basic-quotient
69       @fixnum-remainder @flonum-remainder @integer-remainder @basic-remainder
70       @fixnum-quotient&remainder @flonum-quotient&remainder
71       @integer-quotient&remainder @basic-quotient&remainder
72       @fixnum-2-gcd @flonum-2-gcd
73       @basic-positive? @integer-positive? @fixnum-positive?
74       @basic-negative? @integer-negative? @fixnum-negative?
75       @basic-2-= @integer-2-= @basic-2-< @integer-2-<
76       @basic-2-<= @integer-2-<= @basic-2-> @integer-2->
77       @basic-2->= @integer-2->=
78       @basic-even? @basic-odd? @integer-even? @integer-odd?
79       @basic-nan? @flonum-nan? @basic-zero?
80       @basic-finite? @flonum-finite? @basic-infinite? @flonum-infinite?
81       @basic-number->string @fixnum->string @flonum->string @integer->string
82
83       @integer-2-bitwise-and @integer-2-bitwise-ior @integer-2-bitwise-xor
84       @integer-arithmetic-shift @fixnum-length @integer-length
85       @integer-bit-set?
86
87       ;; These must stay exported because they're called as hooks from C.
88       ;; We might later use them in the types db, but it won't help much?
89       @extended-2-plus @extended-2-minus @extended-2-times
90       @bignum-2-times-karatsuba @bignum-2-divrem-burnikel-ziegler
91       @extended-abs @extended-signum @extended-negate
92       @integer->string/recursive @extended-number->string
93       @error-hook)
94
95  (import (except scheme
96                  + - * / = > < >= <= eqv?
97                  number->string string->number 
98                  exp log sin cos tan atan acos asin expt sqrt
99                  quotient modulo remainder
100                  numerator denominator
101                  abs max min gcd lcm
102                  positive? negative? odd? even? zero?
103                  exact? inexact?
104                  rationalize
105                  floor ceiling truncate round
106                  inexact->exact exact->inexact
107                  number? complex? real? rational? integer?
108                  make-rectangular make-polar real-part imag-part magnitude angle)
109          (except chicken add1 sub1 random randomize conj signum finite?
110                  bitwise-and bitwise-ior bitwise-xor bitwise-not
111                  arithmetic-shift bit-set?)
112          foreign)
113
114(foreign-declare "#include \"numbers-c.c\"")
115
116;;; THESE SERVE AS SPECIALIZATION ENTRIES.
117;;; Once this is integrated into core, they can be killed and the
118;;; C functions inlined directly.  Remember to fix the allocations!
119(define @basic-2-plus (##core#primitive "C_2_basic_plus"))
120(define (@fixnum-2-plus a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_plus" 6) a b))
121(define @integer-2-plus (##core#primitive "C_u_2_integer_plus"))
122
123(define @basic-abs (##core#primitive "C_basic_abs"))
124(define (@fixnum-abs x) (##core#inline_allocate ("C_a_u_i_fixnum_abs" 6) x))
125(define @integer-abs (##core#primitive "C_u_integer_abs"))
126
127(define @basic-signum (##core#primitive "C_basic_signum"))
128(define (@fixnum-signum x) (##core#inline "C_u_i_fixnum_signum" x))
129(define (@flonum-signum x) (##core#inline_allocate ("C_a_u_i_flonum_signum" 4) x))
130(define (@integer-signum x) (##core#inline "C_u_i_integer_signum" x))
131
132(define @basic-negate (##core#primitive "C_basic_negate"))
133(define (@fixnum-negate x) (##core#inline_allocate ("C_a_u_i_fixnum_negate" 6) x))
134(define @integer-negate (##core#primitive "C_u_integer_negate"))
135
136(define @basic-2-minus (##core#primitive "C_2_basic_minus"))
137(define (@fixnum-2-minus a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_minus" 6) a b))
138(define @integer-2-minus (##core#primitive "C_u_2_integer_minus"))
139
140(define @basic-2-times (##core#primitive "C_2_basic_times"))
141(define (@fixnum-2-times a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_times" 7) a b))
142(define @integer-2-times (##core#primitive "C_u_2_integer_times"))
143
144(define @basic-quotient (##core#primitive "C_basic_quotient"))
145(define (@fixnum-quotient a b) (##core#inline_allocate ("C_a_u_i_fixnum_quotient_checked" 6) a b))
146(define (@flonum-quotient a b) (##core#inline_allocate ("C_a_i_flonum_actual_quotient_checked" 4) a b))
147(define @integer-quotient (##core#primitive "C_u_integer_quotient"))
148
149(define @basic-remainder (##core#primitive "C_basic_remainder"))
150(define (@fixnum-remainder a b) (##core#inline "C_u_i_fixnum_remainder_checked" a b))
151(define (@flonum-remainder a b) (##core#inline_allocate ("C_a_i_flonum_remainder_checked" 4) a b))
152(define @integer-remainder (##core#primitive "C_u_integer_remainder"))
153
154(define @basic-quotient&remainder (##core#primitive "C_basic_divrem"))
155(define (@fixnum-quotient&remainder a b)
156  (values (##core#inline_allocate ("C_a_u_i_fixnum_quotient_checked" 6) a b)
157          (##core#inline "C_u_i_fixnum_remainder_checked" a b)))
158(define (@flonum-quotient&remainder a b)
159  (values (##core#inline_allocate ("C_a_i_flonum_actual_quotient_checked" 4) a b)
160          (##core#inline_allocate ("C_a_i_flonum_remainder_checked" 4) a b)))
161(define @integer-quotient&remainder (##core#primitive "C_u_integer_divrem"))
162
163;; There is no specialization for flonums because of the integer check.
164;; Perhaps that should change?
165(define (@fixnum-2-gcd a b) (##core#inline "C_u_i_2_fixnum_gcd" a b))
166(define (@flonum-2-gcd a b) (##core#inline_allocate ("C_a_u_i_2_flonum_gcd" 4) a b))
167
168(define (@basic-2-= a b) (##core#inline "C_i_2_basic_equalp" a b))
169(define (@integer-2-= a b) (##core#inline "C_u_i_2_integer_equalp" a b))
170
171(define (@basic-2-< a b) (##core#inline "C_i_2_basic_lessp" a b))
172(define (@integer-2-< a b) (##core#inline "C_u_i_2_integer_lessp" a b))
173
174(define (@basic-2-<= a b) (##core#inline "C_i_2_basic_less_or_equalp" a b))
175(define (@integer-2-<= a b) (##core#inline "C_u_i_2_integer_less_or_equalp" a b))
176
177(define (@basic-2-> a b) (##core#inline "C_i_2_basic_greaterp" a b))
178(define (@integer-2-> a b) (##core#inline "C_u_i_2_integer_greaterp" a b))
179
180(define (@basic-2->= a b) (##core#inline "C_i_2_basic_greater_or_equalp" a b))
181(define (@integer-2->= a b) (##core#inline "C_u_i_2_integer_greater_or_equalp" a b))
182
183(define (@basic-even? x) (##core#inline "C_i_basic_evenp" x))
184(define (@basic-odd? x) (##core#inline "C_i_basic_oddp" x))
185(define (@integer-even? x) (##core#inline "C_u_i_integer_evenp" x))
186(define (@integer-odd? x) (##core#inline "C_u_i_integer_oddp" x))
187(define (@basic-positive? x) (##core#inline "C_i_basic_positivep" x))
188(define (@basic-negative? x) (##core#inline "C_i_basic_negativep" x))
189(define (@integer-positive? x) (##core#inline "C_u_i_integer_positivep" x))
190(define (@integer-negative? x) (##core#inline "C_u_i_integer_negativep" x))
191(define (@fixnum-positive? x) (##core#inline "C_u_i_fixnum_positivep" x))
192(define (@fixnum-negative? x) (##core#inline "C_u_i_fixnum_negativep" x))
193
194(define @integer-2-bitwise-and (##core#primitive "C_u_2_integer_bitwise_and"))
195(define @integer-2-bitwise-ior (##core#primitive "C_u_2_integer_bitwise_ior"))
196(define @integer-2-bitwise-xor (##core#primitive "C_u_2_integer_bitwise_xor"))
197(define @integer-arithmetic-shift (##core#primitive "C_u_integer_shift"))
198
199(define (@fixnum-length x) (##core#inline "C_u_i_fixnum_length" x))
200(define (@integer-length x) (##core#inline "C_u_i_integer_length" x))
201(define (@integer-bit-set? n i) (##core#inline "C_u_i_integer_bit_setp" n i))
202
203;; Listing these predicates here may seem silly, but in core's
204;; types.db it can improve things considerably
205(define (@flonum-nan? x) (##core#inline "C_u_i_flonum_nanp" x))
206(define (@basic-nan? x) (##core#inline "C_i_nanp" x))
207(define (@flonum-finite? x) (##core#inline "C_u_i_flonum_finitep" x))
208(define (@basic-finite? x) (##core#inline "C_i_numbers_finitep" x))
209(define (@flonum-infinite? x) (##core#inline "C_u_i_flonum_infinitep" x))
210(define (@basic-infinite? x) (##core#inline "C_i_numbers_infinitep" x))
211(define (@basic-zero? x) (##core#inline "C_i_numbers_zerop" x))
212
213(define @basic-number->string (##core#primitive "C_basic_number_to_string"))
214(define @fixnum->string (##core#primitive "C_u_fixnum_to_string"))
215(define @flonum->string (##core#primitive "C_u_flonum_to_string"))
216(define @integer->string (##core#primitive "C_u_integer_to_string"))
217
218(define-compiler-syntax exact-integer?
219  (syntax-rules ()
220    ((_ x) (or (##core#inline "C_fixnump" x) (##core#inline "C_i_bignump" x)))))
221
222;;; Error handling
223
224(define (bad-number loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a number" x))
225(define (bad-real loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a real number" x))
226(define (bad-ratnum loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a rational number" x))
227(define (bad-integer loc x) (##sys#signal-hook #:type-error loc "bad argument type - not an integer" x))
228(define (bad-natural loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an nonnegative integer" x))
229(define (bad-positive loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be a positive (non-zero) integer" x))
230(define (bad-complex/o loc x) (##sys#signal-hook #:type-error loc "bad argument type - complex number has no ordering" x))
231(define (bad-base loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a valid base" x))
232(define (bad-inexact loc x) (##sys#signal-hook #:type-error loc "bad argument type - inexact number has no exact representation" x))
233(define (bad-exact loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an exact number" x))
234(define (log0 loc x) (##sys#signal-hook #:arithmetic-error loc "log of exact 0 is undefined" x))
235(define (expt0 loc x y) (##sys#signal-hook #:arithmetic-error loc "exponent of exact 0 with complex argument is undefined" x y))
236(define (div/0 loc x y) (##sys#signal-hook #:arithmetic-error loc "division by zero" x y))
237
238;; Ugly override because we add our own codes
239(define (@error-hook code loc . args)
240  (case code
241    ((99948) (bad-real loc (car args)))
242    ((99949) (bad-complex/o loc (car args)))
243    (else (apply ##sys#error-hook code loc args))))
244
245(define-inline (%init-tags tagvec) (##core#inline "init_tags" tagvec))
246
247;;; Primitives
248
249(define-inline (fp/ x y) (##core#inline_allocate ("C_a_i_flonum_quotient" 4) x y))
250(define-inline (fp* x y) (##core#inline_allocate ("C_a_i_flonum_times" 4) x y))
251(define-inline (fp= x y) (##core#inline "C_flonum_equalp" x y))
252
253(define-inline (compnum-real c) (##sys#slot c 1))
254(define-inline (compnum-imag c) (##sys#slot c 2))
255(define-inline (%make-complex r i) (##sys#make-structure 'compnum r i))
256
257(define-inline (ratnum-numerator c) (##sys#slot c 1))
258(define-inline (ratnum-denominator c) (##sys#slot c 2))
259(define-inline (%make-ratnum r i) (##sys#make-structure 'ratnum r i))
260
261;;; Setup
262
263(%init-tags
264 (vector 'bignum                        ; BIG_TAG
265         'ratnum                        ; RAT_TAG
266         'compnum))                     ; COMP_TAG
267
268(##sys#gc #f)                           ; move tag-vector into 2nd generation
269
270
271;;; Comparisons:
272
273(define-inline (%= a b) (##core#inline "C_i_2_basic_equalp" a b))
274(define-inline (%> a b) (##core#inline "C_i_2_basic_greaterp" a b))
275(define-inline (%< a b) (##core#inline "C_i_2_basic_lessp" a b))
276(define-inline (%>= a b) (##core#inline "C_i_2_basic_greater_or_equalp" a b))
277(define-inline (%<= a b) (##core#inline "C_i_2_basic_less_or_equalp" a b))
278
279(define (eqv? a b) (##core#inline "C_i_numbers_eqvp" a b))
280(define = (##core#primitive "C_numbers_nequalp"))
281(define > (##core#primitive "C_numbers_greaterp"))
282(define >= (##core#primitive "C_numbers_greater_or_equal_p"))
283(define < (##core#primitive "C_numbers_lessp"))
284(define <= (##core#primitive "C_numbers_less_or_equal_p"))
285
286;;; Basic arithmetic:
287
288(define (+ . args)
289  (if (null? args) 
290      0
291      (let ((x (##sys#slot args 0))
292            (args (##sys#slot args 1)))
293        (if (null? args)
294            (if (number? x) x (bad-number '+ x))
295            (let loop ((args (##sys#slot args 1))
296                       (x (%+ x (##sys#slot args 0))))
297              (if (null? args)
298                  x
299                  (loop (##sys#slot args 1) (%+ x (##sys#slot args 0))) ) )  ) ) ) )
300
301(define %+ (##core#primitive "C_2_basic_plus"))
302
303(define (@extended-2-plus x y)
304  (cond ((or (cplxnum? x) (cplxnum? y))
305         ;; Just add real and imag parts together
306         (let ((r (%+ (real-part x) (real-part y)))
307               (i (%+ (imag-part x) (imag-part y))) )
308           (make-complex r i) ))
309        ((ratnum? x)
310         (if (ratnum? y)
311             (rat+/- '+ %+ x y)
312             ;; a/b + c/d = (a*d + b*c)/(b*d)  [with d = 1]
313             (let* ((b (ratnum-denominator x))
314                    (numerator (%+ (ratnum-numerator x) (%* b y))))
315               (if (##core#inline "C_i_flonump" numerator)
316                   (%/ numerator b)
317                   (%make-ratnum numerator b)))))
318        ((ratnum? y)
319         ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
320         (let* ((d (ratnum-denominator y))
321                (numerator (%+ (%* x d) (ratnum-numerator y))))
322           (if (##core#inline "C_i_flonump" numerator)
323               (%/ numerator d)
324               (%make-ratnum numerator d))))
325        (else (bad-number '+ y)) ) )
326
327(define (- arg1 . args)
328  (if (null? args)
329      ((##core#primitive "C_basic_negate") arg1)
330      (let loop ((args (##sys#slot args 1)) (x (%- arg1 (##sys#slot args 0))))
331        (if (null? args)
332            x
333            (loop (##sys#slot args 1) (%- x (##sys#slot args 0))) ) ) ) )
334
335(define (@extended-negate x)
336  (cond ((ratnum? x)
337         (%make-ratnum ((##core#primitive "C_u_integer_negate")
338                        (ratnum-numerator x))
339                       (ratnum-denominator x)))
340        ((cplxnum? x)
341         (%make-complex ((##core#primitive "C_basic_negate") (compnum-real x))
342                        ((##core#primitive "C_basic_negate") (compnum-imag x))))
343        (else (bad-number '- x)) ) ) ; loc?
344
345(define %- (##core#primitive "C_2_basic_minus"))
346
347(define (@extended-2-minus x y)
348  (cond ((or (cplxnum? x) (cplxnum? y))
349         ;; Just subtract real and imag parts from eachother
350         (let ((r (%- (real-part x) (real-part y)))
351               (i (%- (imag-part x) (imag-part y))))
352           (make-complex r i) ))
353        ((ratnum? x)
354         (if (ratnum? y)
355             (rat+/- '- %- x y)
356             ;; a/b - c/d = (a*d - b*c)/(b*d)  [with d = 1]
357             (let* ((b (ratnum-denominator x))
358                    (numerator (%- (ratnum-numerator x) (%* b y))))
359               (if (##core#inline "C_i_flonump" numerator)
360                   (%/ numerator b)
361                   (%make-ratnum numerator b)))))
362        ((ratnum? y)
363         ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
364         (let* ((d (ratnum-denominator y))
365                (numerator (%- (%* x d) (ratnum-numerator y))))
366           (if (##core#inline "C_i_flonump" numerator)
367               (%/ numerator d)
368               (%make-ratnum numerator d))))
369        (else (bad-number '- y)) ) )
370
371(define (* . args)
372  (if (null? args) 
373      1
374      (let ((x (##sys#slot args 0))
375            (args (##sys#slot args 1)))
376        (if (null? args)
377            (if (number? x) x (bad-number '+ x))
378            (let loop ((args (##sys#slot args 1))
379                       (x (%* x (##sys#slot args 0))))
380              (if (null? args)
381                  x
382                  (loop (##sys#slot args 1) (%* x (##sys#slot args 0))) ) )  ) ) ) )
383
384(define %* (##core#primitive "C_2_basic_times"))
385
386(define (@extended-2-times x y)
387  (define (nonrat*rat x y)
388    ;; a/b * c/d = a*c / b*d  [with b = 1]
389    ;;  =  ((a / g) * c) / (d / g)
390    ;; With   g = gcd(a, d)   and  a = x   [Knuth, 4.5.1]
391    (let* ((d (ratnum-denominator y))
392           (g (%gcd '* x d)))
393      (%ratnum (%* (quotient x g) (ratnum-numerator y))
394               (quotient d g))))
395
396  (cond ((or (cplxnum? x) (cplxnum? y))
397         (let* ((a (real-part x)) (b (imag-part x))
398                (c (real-part y)) (d (imag-part y))
399                (r (%- (%* a c) (%* b d)))
400                (i (%+ (%* a d) (%* b c))) )
401           (make-complex r i) ) )
402        ((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
403         ;; This may be incorrect when one is a ratnum consisting of bignums
404         (fp* (exact->inexact y) (exact->inexact x))) ; loc?
405        ((ratnum? x)
406         (if (ratnum? y)
407             ;; a/b * c/d = a*c / b*d  [generic]
408             ;;   = ((a / g1) * (c / g2)) / ((b / g2) * (d / g1))
409             ;; With  g1 = gcd(a, d)  and   g2 = gcd(b, c) [Knuth, 4.5.1]
410             (let* ((a (ratnum-numerator x)) (b (ratnum-denominator x))
411                    (c (ratnum-numerator y)) (d (ratnum-denominator y))
412                    (g1 (%integer-gcd a d))
413                    (g2 (%integer-gcd b c)))
414               (%ratnum (%* (quotient a g1) (quotient c g2))
415                        (%* (quotient b g2) (quotient d g1))))
416             (nonrat*rat y x)))
417        ((ratnum? y) (nonrat*rat x y))
418        (else (bad-number '* x))))
419
420(define-inline (bignum-digit-count b) (##core#inline "C_u_i_bignum_size" b))
421(define %extract-digits (##core#primitive "C_u_bignum_extract_digits"))
422
423;; Karatsuba multiplication: invoked from C when the two numbers are
424;; large enough to make it worthwhile.  Complexity is O(n^log2(3)),
425;; where n is max(len(x), len(y)).  The description in [Knuth, 4.3.3]
426;; leaves a lot to be desired.  [MCA, 1.3.2] and [MpNT, 3.2] are a bit
427;; easier to understand.  We assume that length(x) <= length(y).
428(define (@bignum-2-times-karatsuba x y)
429  (let* ((same? (eqv? x y))             ; Check before calling (abs)
430         (rs (fx* (##core#inline "C_u_i_integer_signum" x)
431                  (##core#inline "C_u_i_integer_signum" y)))
432         (x ((##core#primitive "C_u_integer_abs") x))
433         (n (bignum-digit-count y))
434         (n/2 (fxshr n 1))
435         (bits (fx* n/2 (foreign-value "C_BIGNUM_DIGIT_LENGTH" int)))
436         (x-hi (%extract-digits x n/2 #f))
437         (x-lo (%extract-digits x 0 n/2)))
438    (if same?              ; This looks pointless, but reduces garbage
439        (let* ((a  (%* x-hi x-hi))
440               (b  (%* x-lo x-lo))
441               (ab (%- x-hi x-lo))
442               (c  (%* ab ab)))
443          (%+ (arithmetic-shift a (fxshl bits 1))
444              (%+ (arithmetic-shift (%+ b (%- a c)) bits)
445                  b)))
446        (let* ((y ((##core#primitive "C_u_integer_abs") y))
447               (y-hi (%extract-digits y n/2 #f))
448               (y-lo (%extract-digits y 0 n/2))
449               (a  (%* x-hi y-hi))
450               (b  (%* x-lo y-lo))
451               (c  (%* (%- x-hi x-lo)
452                       (%- y-hi y-lo))))
453          (%* rs (%+ (arithmetic-shift a (fxshl bits 1))
454                     (%+ (arithmetic-shift (%+ b (%- a c)) bits)
455                         b)))))))
456
457(define (/ arg1 . args)
458  (if (null? args) 
459      (%/ 1 arg1)
460      (let loop ((args (##sys#slot args 1)) (x (%/ arg1 (##sys#slot args 0))))
461        (if (null? args)
462            x
463            (loop (##sys#slot args 1) (%/ x (##sys#slot args 0))) ) ) ) )
464
465(define (%/ x y)
466  (when (eq? y 0) (div/0 '/ x y))
467  (cond ((and (exact-integer? x) (exact-integer? y))
468         (let ((g (%integer-gcd x y)))
469           (%ratnum ((##core#primitive "C_u_integer_quotient") x g)
470                    ((##core#primitive "C_u_integer_quotient") y g))))
471        ;; Compnum *must* be checked first
472        ((or (cplxnum? x) (cplxnum? y))
473         (let* ((a (real-part x)) (b (imag-part x))
474                (c (real-part y)) (d (imag-part y))
475                (r (%+ (%* c c) (%* d d)))
476                (x (%/ (%+ (%* a c) (%* b d)) r))
477                (y (%/ (%- (%* b c) (%* a d)) r)) )
478           (make-complex x y) ))
479        ((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
480         ;; This may be incorrect when one is a ratnum consisting of bignums
481         (fp/ (exact->inexact x) (exact->inexact y)))
482        ((ratnum? x)
483         (if (ratnum? y)
484             ;; a/b / c/d = a*d / b*c  [generic]
485             ;;   = ((a / g1) * (d / g2) * sign(a)) / abs((b / g2) * (c / g1))
486             ;; With   g1 = gcd(a, c)   and    g2 = gcd(b, d) [Knuth, 4.5.1 ex. 4]
487             (let* ((a (ratnum-numerator x)) (b (ratnum-denominator x))
488                    (c (ratnum-numerator y)) (d (ratnum-denominator y))
489                    (g1 (%integer-gcd a c))
490                    (g2 (%integer-gcd b d)))
491               (%ratnum (%* (quotient a g1) (quotient d g2))
492                        (%* (quotient b g2) (quotient c g1))))
493             ;; a/b / c/d = a*d / b*c  [with d = 1]
494             ;;   = ((a / g) * sign(a)) / abs(b * (c / g))
495             ;; With   g = gcd(a, c)   and  c = y  [Knuth, 4.5.1 ex. 4]
496             (let* ((a (ratnum-numerator x))
497                    (g (%gcd '/ a y))
498                    (num (quotient a g))
499                    (denom (%* (ratnum-denominator x) (quotient y g))))
500               (if (##core#inline "C_i_flonump" denom)
501                   (%/ num denom)
502                   (%ratnum num denom)))))
503        ((ratnum? y)
504         ;; a/b / c/d = a*d / b*c  [with b = 1]
505         ;;   = ((a / g1) * d * sign(a)) / abs(c / g1)
506         ;; With   g1 = gcd(a, c)   and   a = x  [Knuth, 4.5.1 ex. 4]
507         (let* ((c (ratnum-numerator y))
508                (g (%gcd '/ x c))
509                (num (%* (quotient x g) (ratnum-denominator y)))
510                (denom (quotient c g)))
511           (if (##core#inline "C_i_flonump" denom)
512               (%/ num denom)
513               (%ratnum num denom))))
514        ((not (number? x)) (bad-number '/ x))
515        (else (bad-number '/ y))) )
516
517;; Burnikel-Ziegler recursive division: Split high number (x) in three
518;; or four parts and divide by the lowest number (y), split in two
519;; parts.  There are descriptions in [MpNT, 4.2], [MCA, 1.4.3] and the
520;; paper "Fast Recursive Division" by Christoph Burnikel & Joachim
521;; Ziegler is freely available.  There is also a description in Karl
522;; Hasselstrom's thesis "Fast Division of Integers".
523;;
524;; The complexity of this is supposedly O(r*s^{log(3)-1} + r*log(s)),
525;; where s is the length of a, and r is the length of b (in digits).
526;;
527;; TODO: See if it's worthwhile to implement "division without remainder"
528;; from the Burnikel-Ziegler paper.
529(define (@bignum-2-divrem-burnikel-ziegler x y return-quot? return-rem?)
530  (define-inline (digit-bits n)
531    (fx* (foreign-value "C_BIGNUM_DIGIT_LENGTH" int) n))
532  (define DIV-LIMIT (foreign-value "C_BURNIKEL_ZIEGLER_THRESHOLD" int))
533
534  ;; Here and in 2n/1n we pass both b and [b1, b2] to avoid splitting
535  ;; up the number more than once.
536  (define (burnikel-ziegler-3n/2n a12 a3 b b1 b2 n)
537    (receive (q^ r1)
538        (if (%< (arithmetic-shift a12 (fxneg (digit-bits n))) b1)
539            (let* ((n/2 (fxshr n 1))
540                   (b11 (%extract-digits b1 n/2 #f))
541                   (b12 (%extract-digits b1 0 n/2)))
542              (burnikel-ziegler-2n/1n a12 b1 b11 b12 n))
543            (let ((base*n (digit-bits n)))
544              (values (%- (arithmetic-shift 1 base*n) 1) ; B^n-1
545                      (%+ (%- a12 (arithmetic-shift b1 base*n)) b1))))
546      (let ((r1a3 (%+ (arithmetic-shift r1 (digit-bits n)) a3)))
547        (let lp ((r^ (%- r1a3 (%* q^ b2)))
548                 (q^ q^))
549          (if (negative? r^)
550              (lp (%+ r^ b) (%- q^ 1))
551              (values q^ r^))))))
552
553  (define (burnikel-ziegler-2n/1n a b b1 b2 n)
554    (if (or (fxodd? n) (fx< n DIV-LIMIT))
555        ((##core#primitive "C_u_integer_divrem") a b)
556        (let* ((n/2 (fxshr n 1))
557               ;; Split a and b into n-sized parts a1,..,a4 and b1,b2
558               (a12 (%extract-digits a n #f))
559               (a3  (%extract-digits a n/2 n))
560               (a4  (%extract-digits a 0 n/2)))
561          (receive (q1 r1) (burnikel-ziegler-3n/2n a12 a3 b b1 b2 n/2)
562            (receive (q2 r) (burnikel-ziegler-3n/2n r1 a4 b b1 b2 n/2)
563              (values (%+ (arithmetic-shift q1 (digit-bits n/2)) q2) r))))))
564
565  ;; The caller will ensure that abs(x) > abs(y)
566  (let* ((q-neg? (if (negative? y) (not (negative? x)) (negative? x)))
567         (r-neg? (negative? x))
568         (x (abs x))
569         (y (abs y))
570         (s (bignum-digit-count y))
571         ;; Define m as min{2^k|(2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
572         ;; This ensures we shift as little as possible (less pressure
573         ;; on the GC) while maintaining a power of two until we drop
574         ;; below the threshold, so we can always split N in half.
575         (m (fxshl 1 (integer-length (fx/ s DIV-LIMIT))))
576         (j (fx/ (fx+ s (fx- m 1)) m))  ; j = s/m, rounded up
577         (n (fx* j m))
578         (norm-shift (fx- (digit-bits n) (integer-length y)))
579         (x (arithmetic-shift x norm-shift))
580         (y (arithmetic-shift y norm-shift))
581         ;; l needs to be the smallest value so that a < base^{l*n}/2
582         (l (fx/ (fx+ (bignum-digit-count x) n) n))
583         (l (if (fx= (digit-bits l) (integer-length x)) (fx+ l 1) l))
584         (t (fxmax l 2))
585         (y-hi (%extract-digits y (fxshr n 1) #f))
586         (y-lo (%extract-digits y 0 (fxshr n 1))))
587    (let lp ((zi (arithmetic-shift x (fxneg (digit-bits (fx* (fx- t 2) n)))))
588             (i (fx- t 2))
589             (quot 0))
590      (receive (qi ri) (burnikel-ziegler-2n/1n zi y y-hi y-lo n)
591        (let ((quot (%+ (arithmetic-shift quot (digit-bits n)) qi)))
592          (if (fx> i 0)
593              (let ((zi-1 (let* ((base*n*i-1 (fx* n (fx- i 1)))
594                                 (base*n*i   (fx* n i))
595                                 (xi-1 (%extract-digits x base*n*i-1 base*n*i)))
596                            (%+ (arithmetic-shift ri (digit-bits n)) xi-1))))
597                (lp zi-1 (fx- i 1) quot))
598              (let ((rem (if (or (not return-rem?) (eq? 0 norm-shift))
599                             ri
600                             (arithmetic-shift ri (fxneg norm-shift)))))
601                ;; Return requested values (quot, rem or both) with correct sign:
602                (cond ((and return-quot? return-rem?)
603                       (values (if q-neg? (- quot) quot)
604                               (if r-neg? (- rem) rem)))
605                      (return-quot? (if q-neg? (- quot) quot))
606                      (else (if r-neg? (- rem) rem))))))))))
607
608
609;;; Complex numbers
610
611(define (make-complex r i)
612  (if (or (eq? i 0) (and (##core#inline "C_i_flonump" i) (fp= i 0.0)))
613      r
614      (%make-complex (if (inexact? i) (exact->inexact r) r)
615                     (if (inexact? r) (exact->inexact i) i)) ) )
616
617(define (make-rectangular r i)
618  (unless (real? r) (bad-real 'make-rectangular r))
619  (unless (real? i) (bad-real 'make-rectangular i))
620  (make-complex r i) )
621
622(define (make-polar r phi)
623  (unless (real? r) (bad-real 'make-rectangular r))
624  (unless (real? phi) (bad-real 'make-rectangular phi))
625  (let ((fphi (exact->inexact phi)))
626    (make-complex (%* r (##core#inline_allocate ("C_a_i_cos" 4) fphi))
627                  (%* r (##core#inline_allocate ("C_a_i_sin" 4) fphi)))))
628
629(define (real-part x)
630  (cond ((cplxnum? x) (compnum-real x))
631        ((number? x) x)
632        (else (bad-number 'real-part x))))
633
634(define (imag-part x)
635  (cond ((cplxnum? x) (compnum-imag x))
636        ((##core#inline "C_i_flonump" x) 0.0)
637        ((number? x) 0)
638        (else (bad-number 'imag-part x))))
639
640(define (magnitude x)
641  (cond ((cplxnum? x)
642         (let ((r (compnum-real x))
643               (i (compnum-imag x)) )
644           (%sqrt 'magnitude (%+ (%* r r) (%* i i))) ))
645        ;; Avoid bad error message (location is hardcoded in abs)
646        ((number? x) ((##core#primitive "C_basic_abs") x))
647        (else (bad-number 'magnitude x))))
648
649(define (angle x)
650  (if (number? x)
651      (##core#inline_allocate
652       ("C_a_i_atan2" 4)
653       (exact->inexact (imag-part x))
654       (exact->inexact (real-part x)))
655      (else (bad-number 'angle x))))
656
657;;; Rationals
658
659(define (%ratnum m n)
660  (cond
661   ((eq? n 1) m)
662   ((eq? n -1) ((##core#primitive "C_u_integer_negate") m))
663   ((negative? n)
664    (%make-ratnum ((##core#primitive "C_u_integer_negate") m)
665                  ((##core#primitive "C_u_integer_negate") n)))
666   (else (%make-ratnum m n))))
667
668;; Knuth, 4.5.1
669(define (rat+/- loc op x y)
670  (let ((a (ratnum-numerator x)) (b (ratnum-denominator x))
671        (c (ratnum-numerator y)) (d (ratnum-denominator y)))
672    (let ((g1 (%integer-gcd b d)))
673      (cond
674       ((eq? g1 1) (%make-ratnum (op (%* a d) (%* b c)) (%* b d)))
675       ;; Save a quotient and multiplication if the gcd is equal
676       ;; to one of the denominators since quotient of b or d and g1 = 1
677       ((%= g1 b) (let* ((t (op (%* a (quotient d g1)) c))
678                         (g2 (%integer-gcd t g1)))
679                    (%ratnum (quotient t g2) (quotient d g2))))
680       ((%= g1 d) (let* ((t (op a (%* c (quotient b g1))))
681                         (g2 (%integer-gcd t g1)))
682                    (%ratnum (quotient t g2) (quotient b g2))))
683       (else (let* ((b/g1 (quotient b g1))
684                    (t (op (%* a (quotient d g1))
685                           (%* c b/g1)))
686                    (g2 (%integer-gcd t g1)))
687               (%make-ratnum (quotient t g2)
688                             (%* b/g1 (quotient d g2)))))))))
689
690(define (numerator x)
691  (cond ((exact-integer? x) x)
692        ((##core#inline "C_i_flonump" x)
693         (cond ((not (finite? x)) (bad-inexact 'numerator x))
694               ((##core#inline "C_u_i_fpintegerp_fixed" x) x)
695               (else (exact->inexact (numerator (inexact->exact x))))))
696        ((ratnum? x) (ratnum-numerator x))
697        (else (bad-ratnum 'numerator x))))
698
699(define (denominator x)
700  (cond ((exact-integer? x) 1)
701        ((##core#inline "C_i_flonump" x)
702         (cond ((not (finite? x)) (bad-inexact 'denominator x))
703               ((##core#inline "C_u_i_fpintegerp_fixed" x) 1.0)
704               (else (exact->inexact (denominator (inexact->exact x))))))
705        ((ratnum? x) (ratnum-denominator x))
706        (else (bad-ratnum 'denominator x))))
707
708;;; Enhanced versions of other standard procedures
709
710(define (@extended-abs x)
711  (cond ((ratnum? x)
712         (%make-ratnum
713          ((##core#primitive "C_u_integer_abs") (ratnum-numerator x))
714          (ratnum-denominator x)))
715        ((cplxnum? x)
716         (##sys#signal-hook
717          #:type-error 'abs
718          "can not compute absolute value of complex number" x))
719        (else (bad-number 'abs x))))
720
721(define abs (##core#primitive "C_basic_abs"))
722
723(define (number? x) (##core#inline "C_i_numbers_numberp" x))
724(define (integer? x) (##core#inline "C_i_numbers_integerp" x))
725
726(set! ##sys#integer? integer?)
727
728(set! ##sys#number? number?)
729(define complex? number?)
730
731;; All numbers are real, except for compnums
732(define (real? x)
733  (and (##core#inline "C_i_numbers_numberp" x)
734       (not (##sys#structure? x 'compnum)) ) )
735
736(define (rational? x) (and (real? x) (finite? x)))
737
738(define (exact-integer? x)
739  (or (##core#inline "C_fixnump" x) (##core#inline "C_i_bignump" x)) )
740
741
742(define (exact? x)
743  (cond ((##core#inline "C_fixnump" x))
744        ((##core#inline "C_i_bignump" x))
745        ((##core#inline "C_i_flonump" x) #f)
746        ((ratnum? x))
747        ((cplxnum? x)
748         (and (exact? (compnum-real x)) (exact? (compnum-imag x))))
749        (else (bad-number 'exact? x))))
750
751(define ##sys#exact? exact?)
752
753(define (inexact? x)
754  (cond ((##core#inline "C_fixnump" x) #f)
755        ((##core#inline "C_i_bignump" x) #f)
756        ((##core#inline "C_i_flonump" x))
757        ((ratnum? x) #f)
758        ((cplxnum? x)
759         (or (inexact? (compnum-real x)) (inexact? (compnum-imag x))))
760        (else (bad-number 'inexact? x))))
761
762(define ##sys#inexact? inexact?)
763
764(define (positive? x) (##core#inline "C_i_basic_positivep" x))
765(define (negative? x) (##core#inline "C_i_basic_negativep" x))
766
767(define (zero? x) (##core#inline "C_i_numbers_zerop" x))
768
769(define (odd? x) (##core#inline "C_i_basic_oddp" x))
770(define (even? x) (##core#inline "C_i_basic_evenp" x))
771
772(define (max x1 . xs)
773  (let loop ((i (##core#inline "C_i_flonump" x1)) (m x1) (xs xs))
774    (if (null? xs)
775        (if i (exact->inexact m) m)
776        (let ((h (##sys#slot xs 0)))
777          (loop (or i (##core#inline "C_i_flonump" h))
778                (if (%> h m) h m)
779                (##sys#slot xs 1)) ) ) )  )
780
781(define (min x1 . xs)
782  (let loop ((i (##core#inline "C_i_flonump" x1)) (m x1) (xs xs))
783    (if (null? xs)
784        (if i (exact->inexact m) m)
785        (let ((h (##sys#slot xs 0)))
786          (loop (or i (##core#inline "C_i_flonump" h))
787                (if (%< h m) h m)
788                (##sys#slot xs 1)) ) ) )  )
789
790(define quotient (##core#primitive "C_basic_quotient"))
791(define truncate-quotient quotient)
792
793(define remainder (##core#primitive "C_basic_remainder"))
794(define truncate-remainder remainder)
795
796;; Modulo's sign follows y  (whereas remainder's sign follows x)
797(define (modulo x y)
798  (let ((r (remainder x y)))
799    (if (positive? y)
800        (if (negative? r) (%+ r y) r)
801        (if (positive? r) (%+ r y) r))))
802
803(define floor-remainder modulo)
804
805(define quotient&remainder (##core#primitive "C_basic_divrem"))
806(define truncate/ quotient&remainder)
807
808;; Modulo's sign follows y  (whereas remainder's sign follows x)
809(define (quotient&modulo x y)
810  (receive (div rem) (quotient&remainder x y)
811    (if (positive? y)
812        (if (negative? rem)
813            (values div (%+ rem y))
814            (values div rem))
815        (if (positive? rem)
816            (values div (%+ rem y))
817            (values div rem)))))
818
819;; Same as above, but quotient gets adjusted along with the remainder
820(define (floor/ x y)
821  (receive (div rem) (quotient&remainder x y)
822    (if (positive? y)
823        (if (negative? rem)
824            (values (%- div 1) (%+ rem y))
825            (values div rem))
826        (if (positive? rem)
827            (values (%- div 1) (%+ rem y))
828            (values div rem)))))
829
830;; XXX This is a bad bad bad definition; very inefficient But to
831;; improve it we would need to provide another implementation of the
832;; quotient procedure which floors instead of truncates.
833(define (floor-quotient x y)
834  (receive (div rem) (floor/ x y) div))
835
836
837(define (%flo->rat x)
838  ;; Try to multiply by two until we reach an integer
839  (define (float-fraction-length x)
840    (do ((x x (fp* x 2.0))
841         (i 0 (fx+ i 1)))
842        ((##core#inline "C_u_i_fpintegerp_fixed" x) i)))
843
844  (define (deliver y d)
845    (let* ((q (%integer-power 2 (float-fraction-length y)))
846           (scaled-y (%* y (exact->inexact q))))
847      (if (finite? scaled-y)          ; Shouldn't this always be true?
848          (%/ (%/ ((##core#primitive "C_u_flo_to_int") scaled-y) q) d)
849          (bad-inexact 'inexact->exact x))))
850
851  (if (and (%< x 1.0)         ; Watch out for denormalized numbers
852           (%> x -1.0))       ; XXX: Needs a test, it seems pointless
853      (deliver (%* x (expt 2.0 flonum-precision))
854               ;; Can be bignum (is on 32-bit), so must wait until after init.
855               ;; We shouldn't need to calculate this every single time, tho..
856               (%integer-power 2 flonum-precision))
857      (deliver x 1)))
858
859(define (inexact->exact x)
860  (cond ((exact? x) x)
861        ((##core#inline "C_i_flonump" x)
862         (cond ((##core#inline "C_u_i_fpintegerp_fixed" x)
863                ((##core#primitive "C_u_flo_to_int") x))
864               ((##core#inline "C_u_i_flonum_finitep" x)
865                (%flo->rat x))
866               (else (bad-inexact 'inexact->exact x))))
867        ((cplxnum? x)
868         (make-complex (inexact->exact (compnum-real x))
869                       (inexact->exact (compnum-imag x))))
870        (else (bad-number 'inexact->exact x))))
871
872(define ##sys#inexact->exact inexact->exact)
873(define exact inexact->exact)
874
875;; Exponent of the lowest allowed flonum; if we get any lower we get zero.
876;; In other words, this is the first (smallest) flonum after 0.
877;; Equal to (expt 2.0 (- flonum-minimum-exponent flonum-precision))
878(define minimum-denorm-flonum-expt (fx- flonum-minimum-exponent flonum-precision))
879
880(define (exact->inexact x)
881  (cond ((##core#inline "C_fixnump" x)
882         (##core#inline_allocate ("C_a_i_fix_to_flo" 4) x))
883        ((##core#inline "C_i_flonump" x) x)
884        ((##core#inline "C_i_bignump" x)
885         (##core#inline_allocate ("C_a_u_i_big_to_flo" 4) x))
886        ((ratnum? x)
887         ;; This tries to keep the numbers within representable ranges
888         ;; and tries to drop as few significant digits as possible by
889         ;; bringing the two numbers to within the same powers of two.
890         ;; See algorithms M & N in Knuth, 4.2.1
891         (let* ((n1 (ratnum-numerator x))
892                (an ((##core#primitive "C_u_integer_abs") n1))
893                (d1 (ratnum-denominator x))
894                ;; Approximate distance between the numbers in powers
895                ;; of 2 ie, 2^e-1 < n/d < 2^e+1 (e is the *un*biased
896                ;; value of e_w in M2)
897                ;; XXX: What if b != 2 (ie, flonum-radix is not 2)?
898                (e (fx- (integer-length an) (integer-length d1)))
899                (rnd (lambda (n d e) ; Here, 1 <= n/d < 2  (normalized) [N5]
900                       ;; Cannot shift above the available precision,
901                       ;; and can't have an exponent that's below the
902                       ;; minimum flonum exponent.
903                       (let* ((s (min (fx- flonum-precision 1)
904                                      (fx- e minimum-denorm-flonum-expt)))
905                              (normalized (%/ (arithmetic-shift n s) d))
906                              (r (round normalized))
907                              (fraction (exact->inexact r))
908                              (exp (fx- e s)))
909                         (let ((res (fp* fraction (expt 2.0 exp))))
910                           (if (negative? n1) (%- 0 res) res)))))
911                (scale (lambda (n d)  ; Here, 1/2 <= n/d < 2   [N3]
912                         (if (%< n d) ; n/d < 1?
913                             ;; Scale left [N3]; only needed once (see note in M3)
914                             (rnd (arithmetic-shift n 1) d (fx- e 1))
915                             ;; Already normalized
916                             (rnd n d e)))))
917           ;; After this step, which shifts the smaller number to
918           ;; align with the larger, "f" in algorithm N is represented
919           ;; in the procedures above by n/d.
920           (if (negative? e)
921               (scale (arithmetic-shift an (%- 0 e)) d1)
922               (scale an (arithmetic-shift d1 e)))))
923        ((cplxnum? x)
924         (make-complex (exact->inexact (compnum-real x))
925                       (exact->inexact (compnum-imag x))))
926        (else (bad-number 'exact->inexact x))))
927
928(define inexact exact->inexact)
929(define ##sys#exact->inexact exact->inexact)
930
931(define (fxgcd a b)
932  (##core#inline "C_u_i_2_fixnum_gcd" a b))
933
934(define (fpgcd a b)
935  (##core#inline_allocate ("C_a_u_i_2_flonum_gcd" 4) a b))
936
937(define (%integer-gcd a b)
938  (define k fixnum-precision) ; Can be anything between 2 and min(F, B).
939  (define k/2 (fx/ k 2))      ; F is fixnum precision and B bits in a big digit
940  (define-inline (len x) (##core#inline "C_u_i_integer_length" x))
941
942  ;; This is Lehmer's GCD algorithm with Jebelean's quotient test, as
943  ;; it is presented in the paper "An Analysis of Lehmer’s Euclidean
944  ;; GCD Algorithm", by J. Sorenson.  Fuck the ACM and their goddamn
945  ;; paywall; you can currently find the paper here:
946  ;; http://www.csie.nuk.edu.tw/~cychen/gcd/An%20analysis%20of%20Lehmer%27s%20Euclidean%20GCD%20algorithm.pdf
947  ;; If that URI fails, it's also explained in [MpNT, 5.2]
948  ;;
949  ;; The basic idea is to avoid divisions which yield only small
950  ;; quotients, in which the remainder won't reduce the numbers by
951  ;; much.  This can be detected by dividing only the leading k bits.
952  (define (lehmer-gcd u v)
953    (let ((-h (fxneg (fx- (len u) k))))
954      (let lp ((i-even? #t)
955               (u^ (arithmetic-shift u -h))
956               (v^ (arithmetic-shift v -h))
957               (x-prev 1) (y-prev 0)
958               (x-curr 0) (y-curr 1))
959        (let* ((q^ (fx/ u^ v^))     ; Estimated quotient for this step
960               (x-next (fx- x-prev (fx* q^ x-curr)))
961               (y-next (fx- y-prev (fx* q^ y-curr))))
962          ;; Euclidian GCD swap on u^ and v^
963          (let ((u^ v^)
964                (v^ (fx- u^ (fx* q^ v^))))
965            (let ((done? (if i-even?
966                             (or (fx< v^ (fxneg y-next))
967                                 (fx< (fx- u^ v^) (fx- x-next x-curr)))
968                             (or (fx< v^ (fxneg x-next))
969                                 (fx< (fx- u^ v^) (fx- y-next y-curr))))))
970              (if done?
971                  (values (+ (* x-prev u) (* y-prev v))
972                          (+ (* x-curr u) (* y-curr v)))
973                  (lp (not i-even?) u^ v^ x-curr y-curr x-next y-next))))))))
974
975  ;; This implements the basic Euclidian GCD algorithm, with a
976  ;; conditional call to Lehmer's GCD algorithm when the length
977  ;; difference between a and b is at most one halfdigit.
978  ;; The complexity of the whole thing is supposedly O(n^2/log n)
979  ;; where n is the number of bits in a and b.
980  (let* ((a (abs a)) (b (abs b))   ; Enforce loop invariant on input:
981         (swap? (%< a b)))         ; both must be positive, and a >= b
982    (let lp ((a (if swap? b a))
983             (b (if swap? a b)))
984      (cond ((eq? b 0) a)
985            ((fx< (fx- (len a) (len b)) k/2)
986             (receive (a b) (lehmer-gcd a b)
987               (if (eq? b 0)
988                   a
989                   (lp b ((##core#primitive "C_u_integer_remainder") a b)))))
990            ((fixnum? a)      ; b MUST be fixnum due to loop invariant
991             (##core#inline "C_u_i_2_fixnum_gcd" a b))
992            (else
993             (lp b ((##core#primitive "C_u_integer_remainder") a b)))))))
994
995(define (%gcd loc a b)
996  (cond ((exact-integer? a)
997         (cond ((exact-integer? b)
998                (%integer-gcd a b))
999               ((and (##core#inline "C_i_flonump" b)
1000                     (##core#inline "C_u_i_fpintegerp_fixed" b))
1001                (exact->inexact (%integer-gcd a (inexact->exact b))))
1002               (else (bad-integer loc b))))
1003        ((and (##core#inline "C_i_flonump" a)
1004              (##core#inline "C_u_i_fpintegerp_fixed" a))
1005         (cond ((##core#inline "C_i_flonump" b)
1006                (##core#inline_allocate ("C_a_u_i_2_flonum_gcd" 4) a b))
1007               ((exact-integer? b)
1008                (exact->inexact (%integer-gcd (inexact->exact a) b)))
1009               (else (bad-integer loc b))))
1010        (else (bad-integer loc a))))
1011
1012(define (gcd . ns)
1013  (if (eq? ns '())
1014      0
1015      (let loop ((head (##sys#slot ns 0))
1016                 (next (##sys#slot ns 1)))
1017        (if (null? next)
1018            (if (integer? head)
1019                ((##core#primitive "C_basic_abs") head)
1020                (bad-integer 'gcd head))
1021            (let ((n2 (##sys#slot next 0)))
1022              (loop (%gcd 'gcd head n2)
1023                    (##sys#slot next 1)) ) )  ) ) )
1024
1025(define (lcm . ns)
1026  (if (null? ns)
1027      1
1028      (let loop ((head (##sys#slot ns 0))
1029                 (next (##sys#slot ns 1)))
1030        (if (null? next)
1031            (if (integer? head)
1032                ((##core#primitive "C_basic_abs") head)
1033                (bad-integer 'gcd head))
1034            (let ((n2 (##sys#slot next 0)))
1035              (loop (quotient (%* head n2)
1036                     (%gcd 'lcm head n2))
1037                    (##sys#slot next 1)) ) )  ) ) )
1038
1039(define (floor x)
1040  (cond ((exact-integer? x) x)
1041        ((##core#inline "C_i_flonump" x) (fpfloor x))
1042        ;; (floor x) = greatest integer <= x
1043        ((ratnum? x) (let* ((n (ratnum-numerator x))
1044                            (q (quotient n (ratnum-denominator x))))
1045                       (if (%>= n 0) q (%- q 1))))
1046        (else (bad-real 'floor x))))
1047
1048(define (ceiling x)
1049  (cond ((exact-integer? x) x)
1050        ((##core#inline "C_i_flonump" x) (fpceiling x))
1051        ;; (ceiling x) = smallest integer >= x
1052        ((ratnum? x) (let* ((n (ratnum-numerator x))
1053                            (q (quotient n (ratnum-denominator x))))
1054                       (if (%>= n 0) (%+ q 1) q)))
1055        (else (bad-real 'ceiling x))) )
1056
1057(define (truncate x)
1058  (cond ((exact-integer? x) x)
1059        ((##core#inline "C_i_flonump" x) (fptruncate x))
1060        ;; (rational-truncate x) = integer of largest magnitude <= (abs x)
1061        ((ratnum? x) (quotient (ratnum-numerator x)
1062                                (ratnum-denominator x)))
1063        (else (bad-real 'truncate x))))
1064
1065(define (round x)
1066  (cond ((exact-integer? x) x)
1067        ((##core#inline "C_i_flonump" x)
1068         (##core#inline_allocate ("C_a_i_flonum_round_proper" 4) x))
1069        ((ratnum? x) (let* ((x+1/2 (%+ x (%make-ratnum 1 2)))
1070                            (r (floor x+1/2)))
1071                       (if (and (%= r x+1/2) (odd? r)) (%- r 1) r)))
1072        (else (bad-real 'round x))))
1073
1074(define (find-ratio-between x y)
1075  (define (sr x y)
1076    (let ((fx (inexact->exact (floor x))) 
1077          (fy (inexact->exact (floor y))))
1078      (cond ((not (%< fx x)) (list fx 1))
1079            ((%= fx fy) 
1080             (let ((rat (sr (%/ 1 (%- y fy)) (%/ 1 (%- x fx)))))
1081               (list (%+ (cadr rat) (%* fx (car rat))) (car rat))))
1082            (else (list (%+ 1 fx) 1)))))
1083  (cond ((%< y x) (find-ratio-between y x))
1084        ((not (%< x y)) (list x 1))
1085        ((positive? x) (sr x y))
1086        ((negative? y) (let ((rat (sr (%- 0 y) (%- 0 x))))
1087                         (list (%- 0 (car rat)) (cadr rat))))
1088        (else '(0 1))))
1089
1090(define (find-ratio x e) (find-ratio-between (%- x e) (%+ x e)))
1091
1092(define (rationalize x e)
1093  (let ((result (apply %/ (find-ratio x e))))
1094    (if (or (inexact? x) (inexact? e))
1095        (exact->inexact result)
1096        result)))
1097
1098(define (exp n)
1099  (cond ((not (number? n)) (bad-number 'exp n))
1100        ((cplxnum? n)
1101         (%* (##core#inline_allocate ("C_a_i_exp" 4) (compnum-real n))
1102             (let ((p (compnum-imag n)))
1103               (make-complex
1104                (##core#inline_allocate ("C_a_i_cos" 4) p)
1105                (##core#inline_allocate ("C_a_i_sin" 4) p) ) ) ))
1106        (else (##core#inline_allocate ("C_a_i_exp" 4) (exact->inexact n)))))
1107
1108(define (%log x)
1109  (cond
1110   ((eq? x 0)                        ; Exact zero?  That's undefined
1111    (log0 'log x))
1112   ;; avoid calling inexact->exact on X here (to avoid overflow?)
1113   ((or (cplxnum? x) (negative? x)) ; General case
1114    (%+ (%log (magnitude x)) (* (make-complex 0 1) (angle x))))
1115   (else ; Real number case (< already ensured the argument type is a number)
1116    (##core#inline_allocate ("C_a_i_log" 4) (exact->inexact x)))))
1117
1118(define (log a #!optional b)
1119  (if b (%/ (%log a) (%log b)) (%log a)))
1120
1121;; These can be removed after integration into core and bootstrapping,
1122;; when the compiler can write these objects natively.
1123(define %i (%make-complex 0 1))
1124(define %-i (%make-complex 0 -1))
1125(define %i2 (%make-complex 0 2))
1126
1127(define (sin n)
1128  (cond ((not (number? n)) (bad-number 'sin n))
1129        ((cplxnum? n) (let ((in (%* %i n)))
1130                        (%/ (%- (exp in) (exp (%- 0 in))) %i2)))
1131        (else (##core#inline_allocate ("C_a_i_sin" 4) (exact->inexact n)))))
1132
1133(define (cos n)
1134  (cond ((not (number? n)) (bad-number 'cos n))
1135        ((cplxnum? n) (let ((in (%* %i n)))
1136                        (%/ (%+ (exp in) (exp (%- 0 in))) 2) ))
1137        (else (##core#inline_allocate ("C_a_i_cos" 4) (exact->inexact n)))))
1138
1139(define (tan n)
1140  (cond ((not (number? n)) (bad-number 'tan n))
1141        ((cplxnum? n) (%/ (sin n) (cos n)))
1142        (else (##core#inline_allocate ("C_a_i_tan" 4) (exact->inexact n)))))
1143
1144;; General case: sin^{-1}(z) = -i\ln(iz + \sqrt{1-z^2})
1145(define (asin n)
1146  (cond ((not (number? n)) (bad-number 'asin n))
1147        ((and (##core#inline "C_i_flonump" n) (fp>= n -1.0) (fp<= n 1.0))
1148         (##core#inline_allocate ("C_a_i_asin" 4) n))
1149        ((and (##core#inline "C_fixnump" n) (fx>= n -1) (fx<= n 1))
1150         (##core#inline_allocate ("C_a_i_asin" 4)
1151                                 (##core#inline_allocate
1152                                  ("C_a_i_fix_to_flo" 4) n)))
1153        ;; General definition can return compnums
1154        (else (%* %-i (%log (%+ (%* %i n) (%sqrt 'asin (%- 1 (%* n n)))))))))
1155
1156;; General case:
1157;; cos^{-1}(z) = 1/2\pi + i\ln(iz + \sqrt{1-z^2}) = 1/2\pi - sin^{-1}(z) = sin(1) - sin(z)
1158(define acos
1159  (let ((asin1 (##core#inline_allocate ("C_a_i_asin" 4) 1)))
1160    (lambda (n)
1161      (cond ((not (number? n)) (bad-number 'acos n))
1162            ((and (##core#inline "C_i_flonump" n) (fp>= n -1.0) (fp<= n 1.0))
1163             (##core#inline_allocate ("C_a_i_acos" 4) n))
1164            ((and (##core#inline "C_fixnump" n) (fx>= n -1) (fx<= n 1))
1165             (##core#inline_allocate ("C_a_i_acos" 4)
1166                                     (##core#inline_allocate
1167                                      ("C_a_i_fix_to_flo" 4) n)))
1168            ;; General definition can return compnums
1169            (else (%- asin1 (asin n)))))))
1170
1171(define (atan n #!optional b)
1172  (cond ((not (number? n)) (bad-number 'atan n))
1173        ((cplxnum? n) (if b
1174                          (bad-real 'atan n)
1175                          (let ((in (%* %i n)))
1176                            (%/ (%- (%log (%+ 1 in)) (%log (%- 1 in))) %i2))))
1177        (else (if b
1178                  (##core#inline_allocate
1179                   ("C_a_i_atan2" 4) (exact->inexact n) (exact->inexact b))
1180                  (##core#inline_allocate
1181                   ("C_a_i_atan" 4) (exact->inexact n))))))
1182
1183;; This is "Karatsuba Square Root" as described by Paul Zimmermann,
1184;; which is 3/2K(n) + O(n log n) for an input of 2n words, where K(n)
1185;; is the number of operations performed by Karatsuba multiplication.
1186(define (%exact-integer-sqrt a)
1187  ;; Because we assume a3b+a2 >= b^2/4, we must check a few edge cases:
1188  (if (and (fixnum? a) (fx<= a 4))
1189      (case a
1190        ((0 1) (values a 0))
1191        ((2)   (values 1 1))
1192        ((3)   (values 1 2))
1193        ((4)   (values 2 0))
1194        (else  (error "this should never happen")))
1195      (let*-values
1196          (((len/4) (fxshr (fx+ (integer-length a) 1) 2))
1197           ((len/2) (fxshl len/4 1))
1198           ((s^ r^) (%exact-integer-sqrt (arithmetic-shift a (fxneg len/2))))
1199           ((mask)  (%- (arithmetic-shift 1 len/4) 1))
1200           ((a0)    (bitwise-and a mask))
1201           ((a1)    (bitwise-and (arithmetic-shift a (fxneg len/4)) mask))
1202           ((q u)   (quotient&remainder (%+ (arithmetic-shift r^ len/4) a1)
1203                                        (arithmetic-shift s^ 1)))
1204           ((s)     (%+ (arithmetic-shift s^ len/4) q))
1205           ((r)     (%+ (arithmetic-shift u len/4) (%- a0 (%* q q)))))
1206        (if (negative? r)
1207            (values (%- s 1) (%- (%+ r (arithmetic-shift s 1)) 1))
1208            (values s r)))))
1209
1210(define (exact-integer-sqrt x)
1211  (if (or (not (exact-integer? x)) (negative? x))
1212      (bad-natural 'exact-integer-sqrt x)
1213      (%exact-integer-sqrt x)))
1214
1215;; This procedure is so large because it tries very hard to compute
1216;; exact results if at all possible.
1217(define (%sqrt loc n)
1218  (cond ((cplxnum? n)     ; Must be checked before we call "negative?"
1219         (let ((p (%/ (angle n) 2))
1220               (m (##core#inline_allocate ("C_a_i_sqrt" 4) (magnitude n))) )
1221           (make-complex (%* m (cos p)) (%* m (sin p)) ) ))
1222        ((negative? n)
1223         (make-complex
1224          0.0
1225          (##core#inline_allocate
1226           ("C_a_i_sqrt" 4)
1227           (exact->inexact ((##core#primitive "C_basic_negate") n)))))
1228        ((exact-integer? n)
1229         (receive (s^2 r)
1230             (%exact-integer-sqrt n)
1231           (if (eq? 0 r)
1232               s^2
1233               (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))
1234        ((ratnum? n)                    ; Try to compute exact sqrt
1235         (receive (ns^2 nr)             ; We already know n is positive!
1236             (%exact-integer-sqrt (ratnum-numerator n))
1237           (if (eq? nr 0)
1238               (receive (ds^2 dr) (%exact-integer-sqrt (ratnum-denominator n))
1239                 (if (eq? dr 0)
1240                     (%/ ns^2 ds^2)
1241                     (%sqrt loc (exact->inexact n))))
1242               (%sqrt loc (exact->inexact n)))))
1243        (else (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))
1244
1245(define (sqrt x) (%sqrt 'sqrt x))
1246
1247;; Not really used
1248(define (square x) (%* x x))
1249
1250;; Generalized Newton's algorithm for positive integers, with a little help
1251;; from Wikipedia ;)  https://en.wikipedia.org/wiki/Nth_root_algorithm
1252(define (%exact-integer-nth-root loc k n)
1253  (cond
1254   ((or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exact-integer-sqrt on n=2?
1255    (values k 0))
1256   ((not (positive? n))
1257    (bad-positive loc n))
1258   (else
1259    (let ((len (integer-length k)))
1260      (if (%< len n)          ; Idea from Gambit: 2^{len-1} <= k < 2^{len}
1261          (values 1 (%- k 1)) ; Since we know x >= 2, we know x^{n} can't exist
1262          ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n)
1263          (let* ((shift-amount (inexact->exact (ceiling (/ (fx+ len 1) n))))
1264                 (g0 (arithmetic-shift 1 shift-amount))
1265                 (n-1 (%- n 1)))
1266            (let lp ((g0 g0)
1267                     (g1 (quotient (%+ (%* n-1 g0) (quotient k (%integer-power g0 n-1))) n)))
1268              (if (%< g1 g0)
1269                  (lp g1 (quotient (%+ (%* n-1 g1) (quotient k (%integer-power g1 n-1))) n))
1270                  (values g0 (%- k (%integer-power g0 n)))))))))))
1271
1272(define (exact-integer-nth-root k n)
1273  (cond ((or (not (exact-integer? k))
1274             (##core#inline "C_u_i_2_integer_lessp" k 0))
1275         (bad-natural 'exact-integer-nth-root k))
1276        ((not (exact-integer? n))
1277         (bad-integer 'exact-integer-nth-root n))
1278        (else (%exact-integer-nth-root 'exact-integer-nth-root k n))))
1279
1280;; This requires e to be an integer, but base does not need to be
1281(define (%integer-power base e)
1282  (if (negative? e)
1283      (%/ 1 (%integer-power base ((##core#primitive "C_u_integer_negate") e)))
1284      (let lp ((res 1) (e2 e))
1285        (cond
1286         ((eq? e2 0) res)
1287         ((even? e2) ; recursion is faster than iteration here
1288          (%* res (square (lp 1 (arithmetic-shift e2 -1)))))
1289         (else
1290          (lp (%* res base) (%- e2 1)))))))
1291
1292(define (expt a b)
1293  (define (log-expt a b)
1294    (exp (%* b (%log a))))
1295  (define (slow-expt a b)
1296    (if (eq? 0 a)
1297        (expt0 'expt a b)
1298        (exp (%* b (%log a)))))
1299  (cond ((not (number? a)) (bad-number 'expt a))
1300        ((not (number? b)) (bad-number 'expt b))
1301        ((and (ratnum? a) (not (inexact? b)))
1302         ;; (n*d)^b = n^b * d^b = n^b * x^{-b}  | x = 1/b
1303         ;; Hopefully faster than integer-power
1304         (%* (expt (ratnum-numerator a) b)
1305             (expt (ratnum-denominator a)
1306                   ((##core#primitive "C_basic_negate") b))))
1307        ((ratnum? b)
1308         ;; x^{a/b} = (x^{1/b})^a
1309         (cond
1310          ((exact-integer? a)
1311           (if (negative? a)
1312               (log-expt (exact->inexact a) (exact->inexact b))
1313               (receive (ds^n r)
1314                   (%exact-integer-nth-root 'expt a (ratnum-denominator b))
1315                 (if (eq? r 0)
1316                     (%integer-power ds^n (ratnum-numerator b))
1317                     ((##core#primitive "C_expt")
1318                      (exact->inexact a) (exact->inexact b))))))
1319          ((##core#inline "C_i_flonump" a)
1320           (log-expt a (exact->inexact b)))
1321          (else (slow-expt a b))))
1322        ((or (cplxnum? b) (and (cplxnum? a) (not (integer? b))))
1323         (slow-expt a b))
1324        ((and (##core#inline "C_i_flonump" b)
1325              (not (##core#inline "C_u_i_fpintegerp_fixed" b)))
1326         (if (negative? a)
1327             (log-expt (exact->inexact a) (exact->inexact b))
1328             ((##core#primitive "C_expt") (exact->inexact a) (exact->inexact b))))
1329        ((##core#inline "C_i_flonump" a)
1330         ((##core#primitive "C_expt") (exact->inexact a) (exact->inexact b)))
1331        ;; this doesn't work that well, yet...
1332        ;; (XXX: What does this mean? why not? I do know this is ugly... :P)
1333        (else (if (or (inexact? a) (inexact? b))
1334                  (exact->inexact (%integer-power a (inexact->exact b)))
1335                  (%integer-power a b)))) )
1336
1337(define (conj x)                        ; Nonstandard, not part of CHICKEN
1338  (cond ((not (number? x)) (bad-number 'conj x))
1339        ((cplxnum? x)
1340         (make-complex (compnum-real x)
1341                       ((##core#primitive "C_basic_negate") (compnum-imag x))))
1342        (else x)) )
1343
1344(define (add1 n) (%+ n 1))
1345(define (sub1 n) (%- n 1))
1346
1347(define (@extended-signum x)
1348  (cond ((ratnum? x) (##core#inline "C_u_i_integer_signum" (ratnum-numerator x)))
1349        ((cplxnum? x) (make-polar 1 (angle x)))
1350        (else (bad-number 'signum x))))
1351
1352(define signum (##core#primitive "C_basic_signum"))
1353
1354;; This is incompatible with SRFI-33 or SRFI-60! (arg order swapped)
1355(define (bit-set? n i)
1356  (unless (exact-integer? n) (bad-exact 'bit-set? n))
1357  (unless (exact-integer? i) (bad-exact 'bit-set? i))
1358  (##core#inline "C_u_i_integer_bit_setp" n i))
1359
1360;; From SRFI-33 (oddly not in Chicken core).  We could move the whole
1361;; thing into C, but that doesn't buy us much since we want to keep an
1362;; inlineable variant anyway.
1363(define (integer-length x)
1364  (unless (exact-integer? x) (bad-exact 'integer-length x))
1365  (##core#inline "C_u_i_integer_length" x))
1366
1367(define (bitwise-and . xs)
1368  (if (null? xs)
1369      -1
1370      (let ((x1 (##sys#slot xs 0)))
1371        (unless (exact-integer? x1) (bad-exact 'bitwise-and x1))
1372        (let loop ((x x1) (xs (##sys#slot xs 1)))
1373          (if (null? xs)
1374              x
1375              (let ((xi (##sys#slot xs 0)))
1376                (unless (exact-integer? xi) (bad-exact 'bitwise-and xi))
1377                (loop
1378                 ((##core#primitive "C_u_2_integer_bitwise_and") x xi)
1379                 (##sys#slot xs 1) ) ) ) ))) )
1380
1381(define (bitwise-ior . xs)
1382  (if (null? xs)
1383      0
1384      (let ((x1 (##sys#slot xs 0)))
1385        (unless (exact-integer? x1) (bad-exact 'bitwise-ior x1))
1386        (let loop ((x x1) (xs (##sys#slot xs 1)))
1387          (if (null? xs)
1388              x
1389              (let ((xi (##sys#slot xs 0)))
1390                (unless (exact-integer? xi) (bad-exact 'bitwise-ior xi))
1391                (loop
1392                 ((##core#primitive "C_u_2_integer_bitwise_ior") x xi)
1393                 (##sys#slot xs 1) ) ) ) ))) )
1394
1395(define (bitwise-xor . xs)
1396  (if (null? xs)
1397      0
1398      (let ((x1 (##sys#slot xs 0)))
1399        (unless (exact-integer? x1) (bad-exact 'bitwise-ior x1))
1400        (let loop ((x x1) (xs (##sys#slot xs 1)))
1401          (if (null? xs)
1402              x
1403              (let ((xi (##sys#slot xs 0)))
1404                (unless (exact-integer? xi) (bad-exact 'bitwise-xor xi))
1405                (loop
1406                 ((##core#primitive "C_u_2_integer_bitwise_xor") x xi)
1407                 (##sys#slot xs 1) ) ) ) ))) )
1408
1409(define (bitwise-not n)
1410  (unless (exact-integer? n) (bad-exact 'bitwise-not n))
1411  ((##core#primitive "C_u_2_integer_minus") -1 n) )
1412
1413(define (arithmetic-shift n m)
1414  (unless (exact-integer? n) (bad-exact 'arithmetic-shift n))
1415  ;; Strictly speaking, shifting *right* is okay for any number
1416  ;; (ie, shifting by a negative bignum would just result in 0 or -1)...
1417  (unless (##core#inline "C_fixnump" m)
1418    (##sys#signal-hook #:type-error 'arithmetic-shift
1419                       "can only shift by fixnum amounts" n m))
1420  ((##core#primitive "C_u_integer_shift") n m))
1421
1422;; This simple enough idea is from
1423;; http://www.numberworld.org/y-cruncher/algorithms/radix-conversion.html
1424(define (@integer->string/recursive n base expected-string-size)
1425  (let*-values (((halfsize) (fxshr (fx+ expected-string-size 1) 1))
1426                ((b^M/2) (%integer-power base halfsize))
1427                ((hi lo) (quotient&remainder n b^M/2))
1428                ((strhi) (number->string hi base))
1429                ((strlo) (number->string (abs lo) base)))
1430    (string-append strhi
1431                   ;; Fix up any leading zeroes that were stripped from strlo
1432                   (make-string (fx- halfsize (string-length strlo)) #\0)
1433                   strlo)))
1434
1435(define @extended-number->string
1436  (let ((string-append string-append))
1437    (lambda (n base)
1438      (cond
1439       ((ratnum? n) (string-append (number->string (ratnum-numerator n) base)
1440                                   "/"
1441                                   (number->string (ratnum-denominator n) base)))
1442       ;; What about bases that include an "i"?  That could lead to
1443       ;; ambiguous results.
1444       ((cplxnum? n) (let ((r (compnum-real n))
1445                           (i (compnum-imag n)) )
1446                       (string-append
1447                        (number->string r base)
1448                        ;; The infinities and NaN always print their sign
1449                        (if (and (finite? i) (positive? i)) "+" "")
1450                        (number->string i base) "i") ))
1451       (else (bad-number 'number->string n)))  ) ) )
1452
1453(define number->string (##core#primitive "C_basic_number_to_string"))
1454(define ##sys#number->string number->string) ; for printer
1455
1456;; We try to prevent memory exhaustion attacks by limiting the
1457;; maximum exponent value.  Perhaps this should be a parameter?
1458(define-constant +maximum-allowed-exponent+ 10000)
1459
1460;; From "Easy Accurate Reading and Writing of Floating-Point Numbers"
1461;; by Aubrey Jaffer.
1462(define (mantexp->dbl mant point)
1463  (if (not (negative? point))
1464      (exact->inexact (%* mant (%integer-power 10 point)))
1465      (let* ((scl (%integer-power 10 (abs point)))
1466             (bex (fx- (fx- (integer-length mant) (integer-length scl))
1467                       flonum-precision)))
1468        (if (fx< bex 0)
1469            (let* ((num (arithmetic-shift mant (fxneg bex)))
1470                   (quo (round-quotient num scl)))
1471              (cond ((%> (integer-length quo) flonum-precision)
1472                     ;; Too many bits of quotient; readjust
1473                     (set! bex (fx+ 1 bex))
1474                     (set! quo (round-quotient num (%* scl 2)))))
1475              (ldexp (exact->inexact quo) bex))
1476            ;; Fall back to exact calculation in extreme cases
1477            (%* mant (%integer-power 10 point))))))
1478
1479(define ldexp (foreign-lambda double "ldexp" double int))
1480
1481;; Should we export this?
1482(define (round-quotient n d)
1483  (let ((q ((##core#primitive "C_u_integer_quotient") n d)))
1484    (if ((if (even? q) > >=) (%* (abs (remainder n d)) 2) (abs d))
1485        (%+ q (if (eqv? (negative? n) (negative? d)) 1 -1))
1486        q)))
1487
1488;; Shorthand
1489(define-inline (%subchar s i) (##core#inline "C_subchar" s i))
1490(define (%string->compnum radix str offset exactness)
1491  (define (go-inexact!)
1492    ;; Go inexact unless exact was requested (with #e prefix)
1493    (unless (eq? exactness 'e) (set! exactness 'i)))
1494  (define (safe-exponent value e)
1495    (and e (cond
1496            ((not value) 0)
1497            ((%> e +maximum-allowed-exponent+)
1498             (and (eq? exactness 'i)
1499                  (cond ((zero? value) 0.0)
1500                        ((%> value 0.0) +inf.0)
1501                        (else -inf.0))))
1502            ((%< e (fxneg +maximum-allowed-exponent+))
1503             (and (eq? exactness 'i) +0.0))
1504            ((eq? exactness 'i) (mantexp->dbl value e))
1505            (else (%* value (%integer-power 10 e))))))
1506  (let* ((len (##sys#size str))
1507         (0..r (integer->char (fx+ (char->integer #\0) (fx- radix 1))))
1508         (a..r (integer->char (fx+ (char->integer #\a) (fx- radix 11))))
1509         (A..r (integer->char (fx+ (char->integer #\A) (fx- radix 11))))
1510         ;; Ugly flag which we need (note that "exactness" is mutated too!)
1511         ;; Since there is (almost) no backtracking we can do this.
1512         (seen-hashes? #f)
1513         ;; All these procedures return #f or an object consed onto an end
1514         ;; position.  If the cdr is false, that's the end of the string.
1515         ;; If just #f is returned, the string contains invalid number syntax.
1516         (scan-digits
1517          (lambda (start)
1518            (let lp ((i start))
1519              (if (fx= i len)
1520                  (and (fx> i start) (cons i #f))
1521                  (let ((c (%subchar str i)))
1522                    (if (fx<= radix 10)
1523                        (if (and (char>=? c #\0) (char<=? c 0..r))
1524                            (lp (fx+ i 1))
1525                            (and (fx> i start) (cons i i)))
1526                        (if (or (and (char>=? c #\0) (char<=? c #\9))
1527                                (and (char>=? c #\a) (char<=? c a..r))
1528                                (and (char>=? c #\A) (char<=? c A..r)))
1529                            (lp (fx+ i 1))
1530                            (and (fx> i start) (cons i i)))))))))
1531         (scan-hashes
1532          (lambda (start)
1533            (let lp ((i start))
1534              (if (fx= i len)
1535                  (and (fx> i start) (cons i #f))
1536                  (let ((c (%subchar str i)))
1537                    (if (eq? c #\#)
1538                        (lp (fx+ i 1))
1539                        (and (fx> i start) (cons i i))))))))
1540         (scan-digits+hashes
1541          (lambda (start neg? all-hashes-ok?)
1542            (let* ((digits (and (not seen-hashes?) (scan-digits start)))
1543                   (hashes (if digits
1544                               (and (cdr digits) (scan-hashes (cdr digits)))
1545                               (and all-hashes-ok? (scan-hashes start))))
1546                   (end (or hashes digits)))
1547              (and-let* ((end)
1548                         (num ((##core#primitive "C_digits_to_integer")
1549                               str start (car end) radix neg?)))
1550                (when hashes            ; Eeewww. Feeling dirty yet?
1551                  (set! seen-hashes? #t)
1552                  (go-inexact!))
1553                (cons num (cdr end))))))
1554         (scan-exponent
1555          (lambda (start)
1556            (and (fx< start len)
1557                 (let ((sign (case (%subchar str start)
1558                               ((#\+) 'pos) ((#\-) 'neg) (else #f))))
1559                   (and-let* ((start (if sign (fx+ start 1) start))
1560                              (end (scan-digits start)))
1561                     (go-inexact!)
1562                     (cons ((##core#primitive "C_digits_to_integer")
1563                            str start (car end) radix (eq? sign 'neg))
1564                           (cdr end)))))))
1565         (scan-decimal-tail             ; The part after the decimal dot
1566          (lambda (start neg? decimal-head)
1567            (and (fx< start len)
1568                 (let* ((tail (scan-digits+hashes start neg? decimal-head))
1569                        (next (if tail (cdr tail) start)))
1570                   (and (or decimal-head (not next)
1571                            (fx> next start)) ; Don't allow empty "."
1572                        (case (and next (%subchar str next))
1573                          ((#\e #\s #\f #\d #\l
1574                            #\E #\S #\F #\D #\L)
1575                           (and-let* (((fx> len next))
1576                                      (ee (scan-exponent (fx+ next 1)))
1577                                      (e (car ee))
1578                                      (h (safe-exponent decimal-head e)))
1579                             (let* ((te (and tail (fx- e (fx- (cdr tail) start))))
1580                                    (num (and tail (car tail)))
1581                                    (t (safe-exponent num te)))
1582                               (cons (if t (%+ h t) h) (cdr ee)))))
1583                          (else (let* ((last (or next len))
1584                                       (te (and tail (fx- start last)))
1585                                       (num (and tail (car tail)))
1586                                       (t (safe-exponent num te))
1587                                       (h (or decimal-head 0)))
1588                                  (cons (if t (%+ h t) h) next)))))))))
1589         (scan-ureal
1590          (lambda (start neg?)
1591            (if (and (fx> len (fx+ start 1)) (eq? radix 10)
1592                     (eq? (%subchar str start) #\.))
1593                (begin
1594                  (go-inexact!)
1595                  (scan-decimal-tail (fx+ start 1) neg? #f))
1596                (and-let* ((end (scan-digits+hashes start neg? #f)))
1597                  (case (and (cdr end) (%subchar str (cdr end)))
1598                    ((#\.)
1599                     (go-inexact!)
1600                     (and (eq? radix 10)
1601                          (if (fx> len (fx+ (cdr end) 1))
1602                              (scan-decimal-tail (fx+ (cdr end) 1) neg? (car end))
1603                              (cons (car end) #f))))
1604                    ((#\e #\s #\f #\d #\l
1605                      #\E #\S #\F #\D #\L)
1606                     (and-let* (((eq? radix 10))
1607                                ((fx> len (cdr end)))
1608                                (ee (scan-exponent (fx+ (cdr end) 1)))
1609                                (num (car end))
1610                                (val (safe-exponent num (car ee))))
1611                       (cons val (cdr ee))))
1612                    ((#\/)
1613                     (set! seen-hashes? #f) ; Reset flag for denominator
1614                     (and-let* (((fx> len (cdr end)))
1615                                (d (scan-digits+hashes (fx+ (cdr end) 1) #f #f))
1616                                (num (car end))
1617                                (denom (car d)))
1618                       (if (not (eq? denom 0))
1619                           (cons (%/ num denom) (cdr d))
1620                           ;; Hacky: keep around an inexact until we decide we
1621                           ;; *really* need exact values, then fail at the end.
1622                           (and (not (eq? exactness 'e))
1623                                (case (signum num)
1624                                  ((-1) (cons -inf.0 (cdr d)))
1625                                  ((0)  (cons +nan.0 (cdr d)))
1626                                  ((+1) (cons +inf.0 (cdr d))))))))
1627                    (else end))))))
1628         (scan-real
1629          (lambda (start)
1630            (and (fx< start len)
1631                 (let* ((sign (case (%subchar str start)
1632                                ((#\+) 'pos) ((#\-) 'neg) (else #f)))
1633                        (next (if sign (fx+ start 1) start)))
1634                   (and (fx< next len)
1635                        (case (%subchar str next)
1636                          ((#\i #\I)
1637                           (or (and sign
1638                                    (cond
1639                                     ((fx= (fx+ next 1) len) ; [+-]i
1640                                      (cons (if (eq? sign 'neg) -1 1) next))
1641                                     ((and (fx<= (fx+ next 5) len)
1642                                           (string-ci=? (substring str next (fx+ next 5)) "inf.0"))
1643                                      (go-inexact!)
1644                                      ;; XXX TODO: Make this +inf/-inf
1645                                      (cons (fp/ (if (eq? sign 'neg) -1.0 1.0) 0.0)
1646                                            (and (fx< (fx+ next 5) len)
1647                                                 (fx+ next 5))))
1648                                     (else #f)))
1649                               (scan-ureal next (eq? sign 'neg))))
1650                          ((#\n #\N)
1651                           (or (and sign
1652                                    (fx<= (fx+ next 5) len)
1653                                    (string-ci=? (substring str next (fx+ next 5)) "nan.0")
1654                                    (begin (go-inexact!)
1655                                           ;; XXX TODO: Make this +nan.0
1656                                           ;; The atof() hack is used to ensure
1657                                           ;; the nan is positive (otherwise we get
1658                                           ;; into trouble with numbers-syntax on old
1659                                           ;; CHICKENs)
1660                                           (cons (foreign-value "atof(\"+nan\")" float)
1661                                                 (and (fx< (fx+ next 5) len)
1662                                                      (fx+ next 5)))))
1663                               (scan-ureal next (eq? sign 'neg))))
1664                          (else (scan-ureal next (eq? sign 'neg)))))))))
1665         (number (and-let* ((r1 (scan-real offset)))
1666                   (case (and (cdr r1) (%subchar str (cdr r1)))
1667                     ((#f) (car r1))
1668                     ((#\i #\I) (and (fx= len (fx+ (cdr r1) 1))
1669                                     (or (eq? (%subchar str offset) #\+) ; ugh
1670                                         (eq? (%subchar str offset) #\-))
1671                                     (make-rectangular 0 (car r1))))
1672                     ((#\+ #\-)
1673                      (set! seen-hashes? #f) ; Reset flag for imaginary part
1674                      (and-let* ((r2 (scan-real (cdr r1)))
1675                                 ((cdr r2))
1676                                 ((fx= len (fx+ (cdr r2) 1)))
1677                                 ((or (eq? (%subchar str (cdr r2)) #\i)
1678                                      (eq? (%subchar str (cdr r2)) #\I))))
1679                        (make-rectangular (car r1) (car r2))))
1680                     ((#\@)
1681                      (set! seen-hashes? #f) ; Reset flag for angle
1682                      (and-let* ((r2 (scan-real (fx+ (cdr r1) 1)))
1683                                 ((not (cdr r2))))
1684                        (make-polar (car r1) (car r2))))
1685                     (else #f)))))
1686    (and number (if (eq? exactness 'i)
1687                    (exact->inexact number)
1688                    ;; Ensure we didn't encounter +inf.0 or +nan.0 with #e
1689                    (and (finite? number) number)))))
1690
1691(define (string->number str #!optional (base 10))
1692  (##sys#check-string str 'string->number)
1693  (unless (and (##core#inline "C_fixnump" base)
1694               (fx< 1 base) (fx< base 37)) ; We only have 0-9 and the alphabet!
1695    (bad-base 'string->number base))
1696  (let scan-prefix ((i 0)
1697                    (exness #f)
1698                    (radix #f)
1699                    (len (##sys#size str)))
1700    (if (and (fx< (fx+ i 2) len) (eq? (%subchar str i) #\#))
1701        (case (%subchar str (fx+ i 1))
1702          ((#\i #\I) (and (not exness) (scan-prefix (fx+ i 2) 'i radix len)))
1703          ((#\e #\E) (and (not exness) (scan-prefix (fx+ i 2) 'e radix len)))
1704          ((#\b #\B) (and (not radix) (scan-prefix (fx+ i 2) exness 2 len)))
1705          ((#\o #\O) (and (not radix) (scan-prefix (fx+ i 2) exness 8 len)))
1706          ((#\d #\D) (and (not radix) (scan-prefix (fx+ i 2) exness 10 len)))
1707          ((#\x #\X) (and (not radix) (scan-prefix (fx+ i 2) exness 16 len)))
1708          (else #f))
1709        (%string->compnum (or radix base) str i exness))))
1710
1711;;; Reader hook
1712(define (##sys#string->number str #!optional (radix 10) exactness)
1713  (%string->compnum radix str 0 exactness))
1714
1715(define (randomize . n)
1716  (let ((nn (if (null? n)
1717                (##sys#flo2fix (fp/ (current-seconds) 1000.0)) ; wall clock time
1718                (car n))))
1719    (unless (exact-integer? nn) (bad-exact 'randomize nn))
1720    (##core#inline "C_u_i_integer_randomize" nn) ) )
1721
1722(define (random n)
1723  (unless (exact-integer? n) (bad-exact 'random n))
1724  (if (eq? n 0)
1725      0
1726      ((##core#primitive "C_u_integer_random") n) ) )
1727
1728
1729;;; Non-standard type procedures
1730
1731(define (bignum? x) (##core#inline "C_i_bignump" x))
1732
1733(define (nan? x) (##core#inline "C_i_nanp" x))
1734(define (infinite? x) (##core#inline "C_i_numbers_infinitep" x))
1735(define (finite? x) (##core#inline "C_i_numbers_finitep" x))
1736
1737(define (ratnum? x) (##sys#structure? x 'ratnum)) ; rational number
1738
1739;; XXX THE ONES BELOW ARE EXTREMELY CONFUSING
1740;; Especially considering the type tag in a complex number is "compnum"!
1741;; Best to rename cplxnum? to compnum? and ditch the rest.
1742;; A user can easily define them themselves
1743(define (cplxnum? x) (##sys#structure? x 'compnum)) ; complex number
1744
1745(define (rectnum? x)    ; "exact" complex number
1746  (and (cplxnum? x)
1747       (integer? (compnum-real x))
1748       (integer? (compnum-imag x))))
1749
1750(define (compnum? x)    ; inexact complex number
1751  (and (cplxnum? x)
1752       (inexact? (compnum-real x))
1753       (inexact? (compnum-imag x))))
1754
1755(define (cintnum? x)    ; integer number
1756  (or (integer? x) (rectnum? x)) )
1757
1758(define (cflonum? x)    ; floatingpoint number
1759  (or (##core#inline "C_i_flonump" x) (compnum? x)) )
1760
1761;;; What we provide
1762
1763(register-feature! #:full-numeric-tower)
1764
1765)
Note: See TracBrowser for help on using the repository browser.