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

Last change on this file since 33133 was 33133, checked in by sjamaan, 3 years ago

numbers: Fix reference URL which moved, update copyrights

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