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

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

numbers: Fix error reporting when supplied non-numeric arguments

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