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

Last change on this file since 31412 was 31412, checked in by sjamaan, 7 years ago

numbers: Convert remainder procedures to core naming conventions

File size: 77.2 KB
Line 
1;;;; numbers.scm
2;
3; Copyright (c) 2008-2014 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
30(declare 
31  (no-bound-checks)
32  (no-procedure-checks))
33
34(module numbers 
35    (+ - * / = > < >= <= eqv?
36       add1 sub1 signum number->string string->number integer-length
37       bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift
38       equal? ; From scheme. Structural & bytevector comparisons Just Work
39       exp log sin cos tan atan acos asin conj
40       expt sqrt exact-integer-sqrt exact-integer-nth-root
41       quotient modulo remainder quotient&modulo quotient&remainder
42       floor/ floor-quotient floor-remainder
43       truncate/ truncate-quotient truncate-remainder
44       numerator denominator
45       abs max min gcd lcm
46       positive? negative? odd? even? zero?
47       exact? inexact?
48       rationalize
49       random randomize
50       floor ceiling truncate round
51       inexact->exact inexact exact->inexact exact
52       number? complex? real? rational? integer? exact-integer?
53       make-rectangular make-polar real-part imag-part magnitude angle
54       bignum? ratnum? cflonum? rectnum? compnum? cintnum? cplxnum?
55       nan? finite? infinite?)
56
57  (import (except scheme
58                  + - * / = > < >= <= eqv?
59                  number->string string->number 
60                  exp log sin cos tan atan acos asin expt sqrt
61                  quotient modulo remainder
62                  numerator denominator
63                  abs max min gcd lcm
64                  positive? negative? odd? even? zero?
65                  exact? inexact?
66                  rationalize
67                  floor ceiling truncate round
68                  inexact->exact exact->inexact
69                  number? complex? real? rational? integer?
70                  make-rectangular make-polar real-part imag-part magnitude angle)
71          (except chicken add1 sub1 random randomize conj signum finite?
72                  bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift)
73          foreign)
74
75(foreign-declare "#include \"numbers-c.c\"")
76
77(define-foreign-variable FIX integer)
78(define-foreign-variable FLO integer)
79(define-foreign-variable BIG integer)
80(define-foreign-variable RAT integer)
81(define-foreign-variable COMP integer)
82(define-foreign-variable NONE integer)
83
84;;; Error handling
85
86(define (bad-number loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a number" x))
87(define (bad-real loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a real number" x))
88(define (bad-ratnum loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a rational number" x))
89(define (bad-integer loc x) (##sys#signal-hook #:type-error loc "bad argument type - not an integer" x))
90(define (bad-natural loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an nonnegative integer" x))
91(define (bad-positive loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be a positive (non-zero) integer" x))
92(define (bad-complex/o loc x) (##sys#signal-hook #:type-error loc "bad argument type - complex number has no ordering" x))
93(define (bad-base loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a valid base" x))
94(define (bad-inexact loc x) (##sys#signal-hook #:type-error loc "bad argument type - inexact number has no exact representation" x))
95(define (bad-exact loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an exact number" x))
96(define (log0 loc x) (##sys#signal-hook #:arithmetic-error loc "log of exact 0 is undefined" x))
97(define (expt0 loc x y) (##sys#signal-hook #:arithmetic-error loc "exponent of exact 0 with complex argument is undefined" x y))
98(define (div/0 loc x y) (##sys#signal-hook #:arithmetic-error loc "division by zero" x y))
99
100(define-inline (%init-tags tagvec) (##core#inline "init_tags" tagvec))
101(define-inline (%check-number x) (##core#inline "check_number" x))
102
103(define-inline (assert-number x loc)
104  (when (eq? NONE (%check-number x))
105    (bad-number loc x) ) )
106
107(define-inline (fix-div/0 x y loc)
108  (if (eq? y 0)
109      (div/0 loc x y)
110      y) )
111
112;;; Primitives
113
114(define-inline (fp/ x y) (##core#inline_allocate ("C_a_i_flonum_quotient" 4) x y))
115(define-inline (fp+ x y) (##core#inline_allocate ("C_a_i_flonum_plus" 4) x y))
116(define-inline (fp- x y) (##core#inline_allocate ("C_a_i_flonum_difference" 4) x y))
117(define-inline (fp* x y) (##core#inline_allocate ("C_a_i_flonum_times" 4) x y))
118(define-inline (fp= x y) (##core#inline "C_flonum_equalp" x y))
119(define-inline (fp> x y) (##core#inline "C_flonum_greaterp" x y))
120(define-inline (fp< x y) (##core#inline "C_flonum_lessp" x y))
121
122(define-inline (%flonum? x) (##core#inline "C_i_flonump" x))
123(define-inline (%flo-integer? x) (##core#inline "C_u_i_fpintegerp" x))
124
125(define-inline (complex-real c) (##sys#slot c 1))
126(define-inline (complex-imag c) (##sys#slot c 2))
127(define-inline (%make-complex r i) (##sys#make-structure 'compnum r i))
128
129(define-inline (rat-numerator c) (##sys#slot c 1))
130(define-inline (rat-denominator c) (##sys#slot c 2))
131(define-inline (%make-rat r i) (##sys#make-structure 'ratnum r i))
132
133(define-inline (%fix->flo n) (##core#inline_allocate ("C_a_i_fix_to_flo" 4) n))
134(define-inline (%big->flo n) (##core#inline_allocate ("C_a_u_i_big_to_flo" 4) n))
135
136(define %fix+fix (##core#primitive "C_u_2_fixnum_plus"))
137(define %fix+big (##core#primitive "C_u_fixnum_plus_bignum"))
138(define %big+big (##core#primitive "C_u_2_bignum_plus"))
139
140(define %big-neg (##core#primitive "C_u_bignum_negate"))
141;; Can't use fxneg because that breaks in the edge case of negating
142;; the most negative fixnum.  Yes, 2's complement is fun!
143(define %fix-neg (##core#primitive "C_u_fixnum_neg"))
144
145(define %fix-big (##core#primitive "C_u_fixnum_minus_bignum"))
146(define %big-fix (##core#primitive "C_u_bignum_minus_fixnum"))
147(define %big-big (##core#primitive "C_u_2_bignum_minus"))
148
149(define %fix*fix (##core#primitive "C_u_2_fixnum_times"))
150(define %fix*big (##core#primitive "C_u_fixnum_times_bignum"))
151(define %big*big (##core#primitive "C_u_2_bignum_times"))
152
153(define %big-quotient-fix (##core#primitive "C_u_bignum_quotient_fixnum"))
154(define %big-quotient-big (##core#primitive "C_u_bignum_quotient_bignum"))
155
156(define %big-remainder-fix (##core#primitive "C_u_bignum_remainder_fixnum"))
157(define %big-remainder-big (##core#primitive "C_u_bignum_remainder_bignum"))
158
159(define %big-divrem-fix (##core#primitive "big_divrem_fix"))
160(define %big-divrem-big (##core#primitive "big_divrem_big"))
161
162;; This one should really be part of Chicken, hence the name
163(define (fxgcd x y) (##core#inline "C_u_i_2_fixnum_gcd" x y))
164(define biggcd (##core#primitive "big_gcd"))
165(define (fpgcd x y) (##core#inline_allocate ("C_a_u_i_2_flonum_gcd" 4) x y))
166
167(define-inline (%big-cmp x y) (##core#inline "C_u_i_bignum_cmp" x y))
168
169(define %big-abs (##core#primitive "C_u_bignum_abs"))
170
171(define-inline (%big-odd? x) (##core#inline "C_u_i_bignum_oddp" x))
172(define-inline (%big-negative? x) (##core#inline "C_u_i_bignum_negativep" x))
173
174(define %%expt-0 (##core#primitive "C_expt"))
175
176(define (%expt-0 a b)
177  (if (and (negative? a) (not (##sys#integer? b)))
178      (%exp (%* b (%log a)))
179      (%%expt-0 a b)))
180
181(define %quotient-0 (##core#primitive "C_quotient"))
182
183(define %flo->integer (##core#primitive "C_u_flo_to_int"))
184
185(define %int-bitwise-int (##core#primitive "C_u_2_integer_bitwise_op"))
186(define %int-shift-fix (##core#primitive "C_u_int_shift_fix"))
187(define-inline (%int-length i) (##core#inline "C_u_i_int_length" i))
188
189(define number->string-0 (##core#primitive "C_number_to_string"))
190
191(define %big->string (##core#primitive "C_u_bignum_to_digits"))
192(define %digits->integer (##core#primitive "C_digits_to_integer"))
193
194(define-inline (%subchar s i) (##core#inline "C_subchar" s i))
195
196(define-inline (%fix-randomize n) (##core#inline "C_randomize" n))
197(define-inline (%big-randomize n) (##core#inline "big_randomize" n))
198
199(define-inline (%fix-random n) (##core#inline "C_random_fixnum" n))
200(define %big-random (##core#primitive "C_u_bignum_random"))
201
202
203;;; Support macros
204
205(define-syntax switchq
206  (syntax-rules (else)
207    ((_ "aux" _) (##core#undefined))
208    ((_ "aux" _ (else body ...))
209     (begin body ...))
210    ((_ "aux" tmp (val body ...) more ...)
211     (if (eq? tmp val)
212         (begin body ...)
213         (switchq "aux" tmp more ...)))
214    ((_ exp body ...)
215     (let ((tmp exp))
216       (switchq "aux" tmp body ...)))))
217
218
219;;; Setup
220
221(%init-tags
222 (vector 'bignum                        ; BIG_TAG
223         'ratnum                        ; RAT_TAG
224         'compnum))                     ; COMP_TAG
225
226(##sys#gc #f)                           ; move tag-vector into 2nd generation
227
228
229;;; Basic arithmetic:
230
231(define (+ . args)
232  (if (null? args) 
233      0
234      (let ((x (##sys#slot args 0))
235            (rest (##sys#slot args 1)))
236        (cond ((null? rest)
237               (assert-number x '+) 
238               x)
239              (else
240               (let loop ((args rest) (x x))
241                 (if (null? args)
242                     x
243                     (loop (##sys#slot args 1) (%+ x (##sys#slot args 0))) ) ) ) ) ) ) )
244
245(define (%+ x y)
246  (switchq (%check-number x)
247    [FIX 
248     (switchq (%check-number y)
249       [FIX (%fix+fix x y)]
250       [FLO (fp+ (%fix->flo x) y)]
251       [BIG (%fix+big x y)]
252       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
253       [RAT (let ((d (rat-denominator y)))
254              (%/ (%+ (%* x d) (rat-numerator y)) d))]
255       [COMP (%comp+comp (%make-complex x 0) y)]
256       [else (bad-number '+ y)] ) ]
257    [FLO 
258     (switchq (%check-number y)
259       [FIX (fp+ x (%fix->flo y))]
260       [FLO (fp+ x y)]
261       [BIG (fp+ x (%big->flo y))]
262       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
263       [RAT (let ((d (rat-denominator y)))
264              (%/ (%+ (%* x d) (rat-numerator y)) d))]
265       [COMP (%comp+comp (%make-complex x 0) y)]
266       [else (bad-number '+ y)] ) ]
267    [BIG
268     (switchq (%check-number y)
269       [FIX (%fix+big y x)]
270       [FLO (fp+ (%big->flo x) y)]
271       [BIG (%big+big x y)]
272       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
273       [RAT (let ((d (rat-denominator y)))
274              (%/ (%+ (%* x d) (rat-numerator y)) d))]
275       [COMP (%comp+comp (%make-complex x 0) y)]
276       [else (bad-number '+ y)] ) ]
277    [RAT
278     (switchq (%check-number y)
279       [RAT (rat+/- '+ %+ x y)]
280       [COMP (%comp+comp (%make-complex x 0) y)]
281       [NONE (bad-number '+ y)]
282       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with d = 1]
283       [else (let ((b (rat-denominator x)))
284              (%/ (%+ (rat-numerator x) (%* b y)) b))] ) ]
285    [COMP
286     (switchq (%check-number y)
287       [COMP (%comp+comp x y)]
288       [NONE (bad-number '+ y)]
289       [else (%comp+comp x (%make-complex y 0))] ) ]
290    [else (bad-number '+ x)] ) )
291
292(define (%comp+comp x y)
293  (let ([r (%+ (complex-real x) (complex-real y))]
294        [i (%+ (complex-imag x) (complex-imag y))] )
295    (make-complex r i) ) )
296
297(define (- arg1 . args)
298  (if (null? args) 
299      (switchq (%check-number arg1)
300        [FIX (%fix-neg arg1)]
301        [FLO (fpneg arg1)]
302        [BIG (%big-neg arg1)]
303        [RAT (%make-rat (- (rat-numerator arg1)) (rat-denominator arg1))]
304        [COMP (%make-complex (%- 0 (complex-real arg1)) (%- 0 (complex-imag arg1)))]
305        [else (bad-number '- arg1)] )
306      (let loop ([args (##sys#slot args 1)] [x (%- arg1 (##sys#slot args 0))])
307        (if (null? args)
308            x
309            (loop (##sys#slot args 1) (%- x (##sys#slot args 0))) ) ) ) )
310
311(define (%- x y)
312  (switchq (%check-number x)
313    [FIX 
314     (switchq (%check-number y)
315       [FIX (let ((n (%fix-neg y)))
316              (if (= (%check-number n) BIG) ;; fix-neg(most negative) => bignum
317                  (%fix+big x n)
318                  (%fix+fix x n)))]
319       [FLO (fp- (%fix->flo x) y)]
320       [BIG (%fix-big x y)]
321       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
322       [RAT (let ((d (rat-denominator y)))
323              (%/ (%- (%* x d) (rat-numerator y)) d))]
324       [COMP (%comp-comp (%make-complex x 0) y)]
325       [else (bad-number '- y)] ) ]
326    [FLO 
327     (switchq (%check-number y)
328       [FIX (fp- x (%fix->flo y))]
329       [FLO (fp- x y)]
330       [BIG (fp- x (%big->flo y))]
331       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
332       [RAT (let ((d (rat-denominator y)))
333              (%/ (%- (%* x d) (rat-numerator y)) d))]
334       [COMP (%comp-comp (%make-complex x 0) y)]
335       [else (bad-number '- y)] ) ]
336    [BIG
337     (switchq (%check-number y)
338       [FIX (%big-fix x y)]
339       [FLO (fp- (%big->flo x) y)]
340       [BIG (%big-big x y)]             
341       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
342       [RAT (let ((d (rat-denominator y)))
343              (%/ (%- (%* x d) (rat-numerator y)) d))]
344       [COMP (%comp-comp (%make-complex x 0) y)]
345       [else (bad-number '- y)] ) ]
346    [RAT
347     (switchq (%check-number y)
348       [RAT (rat+/- '- %- x y)]
349       [COMP (%comp-comp (%make-complex x 0) y)]
350       [NONE (bad-number '- y)]
351       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with d = 1]
352       [else (let ((b (rat-denominator x)))
353               (%/ (%- (rat-numerator x) (%* b y)) b))] ) ]
354    [COMP
355     (switchq (%check-number y)
356       [COMP (%comp-comp x y)]
357       [NONE (bad-number '- y)] 
358       [else (%comp-comp x (%make-complex y 0))] ) ]
359    [else (bad-number '- x)] ) )
360
361(define (%comp-comp x y)
362  (let ([r (%- (complex-real x) (complex-real y))]
363        [i (%- (complex-imag x) (complex-imag y))] )
364    (make-complex r i) ) )
365
366(define (* . args)
367  (if (null? args) 
368      1
369      (let ((x (##sys#slot args 0))
370            (rest (##sys#slot args 1)))
371        (cond ((null? rest)
372               (assert-number x '*) 
373               x)
374              (else
375               (let loop ((args rest) (x x))
376                 (if (null? args)
377                     x
378                     (loop (##sys#slot args 1) (%* x (##sys#slot args 0))) ) ) ) ) ) ) )
379
380(define-inline (%nonrat*rat x y)
381  ;; a/b * c/d = a*c / b*d  [with b = 1]
382  ;;  =  ((a / g) * c) / (d / g)
383  ;; With   g = gcd(a, d)   and  a = x   [Knuth, 4.5.1]
384  (let* ((d (rat-denominator y))
385         (g (%gcd-0 '* x d)))
386    (ratnum (%* (%quotient '* x g) (rat-numerator y))
387            (%quotient '* d g))))
388
389(define (%* x y)
390  (switchq (%check-number x)
391    [FIX
392     (switchq (%check-number y)
393       [FIX (%fix*fix x y)]
394       [FLO (fp* (%fix->flo x) y)]
395       [BIG (%fix*big x y)]
396       [RAT (%nonrat*rat x y)]
397       [COMP (%comp*comp (%make-complex x 0) y)]
398       [else (bad-number '* y)] ) ]
399    [FLO
400     (switchq (%check-number y)
401       [FIX (fp* x (%fix->flo y))]
402       [FLO (fp* x y)]
403       [BIG (fp* x (%big->flo y))]
404       ;; TODO: This can be incorrect when the ratnum consists of bignums
405       [RAT (fp* x (%exact->inexact y))]
406       [COMP (%comp*comp (%make-complex x 0) y)]
407       [else (bad-number '* y)] ) ]
408    [BIG
409     (switchq (%check-number y)
410       [FIX (%fix*big y x)]
411       [FLO (fp* (%big->flo x) y)]
412       [BIG (%big*big x y)]
413       [RAT (%nonrat*rat x y)]
414       [COMP (%comp*comp (%make-complex x 0) y)]
415       [else (bad-number '* y)] ) ]
416    [RAT
417     (switchq (%check-number y)
418       ;; a/b * c/d = a*c / b*d  [generic]
419       ;;   = ((a / g1) * (c / g2)) / ((b / g2) * (d / g1))
420       ;; With   g1 = gcd(a, d)   and    g2 = gcd(b, c) [Knuth, 4.5.1]
421       [RAT (let* ((a (rat-numerator x)) (b (rat-denominator x))
422                   (c (rat-numerator y)) (d (rat-denominator y))
423                   (g1 (%gcd-0 '* a d)) (g2 (%gcd-0 '* b c)))
424              (ratnum (%* (%quotient '* a g1) (%quotient '* c g2))
425                      (%* (%quotient '* b g2) (%quotient '* d g1))))]
426       [COMP (%comp*comp (%make-complex x 0) y)]
427       ;; TODO: This can be incorrect when the ratnum consists of bignums
428       [FLO (fp* y (%exact->inexact x))]
429       [NONE (bad-number '* y)]
430       [else (%nonrat*rat y x)]) ]
431    [COMP
432     (switchq (%check-number y)
433       [COMP (%comp*comp x y)]
434       [NONE (bad-number '* y)] 
435       [else (%comp*comp x (%make-complex y 0))] ) ]
436    [else (bad-number '* x)] ) )
437
438(define (%comp*comp x y)
439  (let* ([a (complex-real x)]
440         [b (complex-imag x)]
441         [c (complex-real y)]
442         [d (complex-imag y)] 
443         [r (%- (%* a c) (%* b d))]
444         [i (%+ (%* a d) (%* b c))] )
445    (make-complex r i) ) )
446
447(define (/ arg1 . args)
448  (if (null? args) 
449      (%/ 1 arg1)
450      (let loop ([args (##sys#slot args 1)] [x (%/ arg1 (##sys#slot args 0))])
451        (if (null? args)
452            x
453            (loop (##sys#slot args 1) (%/ x (##sys#slot args 0))) ) ) ) )
454
455(define-inline (%nonrat/rat x y)
456  ;; a/b / c/d = a*d / b*c  [with b = 1]
457  ;;   = ((a / g1) * d * sign(a)) / abs(c / g1)
458  ;; With   g1 = gcd(a, c)   and   a = x  [Knuth, 4.5.1 ex. 4]
459  (let* ((c (rat-numerator y))
460         (g (%gcd-0 '/ x c)))
461    (%/ (%* (%quotient '/ x g) (rat-denominator y))
462        (%quotient '/ c g))))
463
464(define (%/ x y)
465  (switchq (%check-number x)
466    [FIX 
467     (switchq (%check-number y)
468       [FIX (let ((g (fxgcd x (fix-div/0 x y '/))))
469              (ratnum (fx/ x g) (fx/ y g)))]
470       [FLO (fp/ (%fix->flo x) y)]
471       ;; XXX TODO: Get rid of the gcd call here, which also divides!
472       [BIG (let ((g (%gcd-0 '/ x y)))
473              (ratnum (%quotient '/ x g) (%quotient '/ y g)))]
474       [RAT (%nonrat/rat x y)]
475       [COMP (%comp/comp (%make-complex x 0) y)]
476       [else (bad-number '/ y)] ) ]
477    [FLO 
478     (switchq (%check-number y)
479       [FIX (fp/ x (%fix->flo (fix-div/0 x y '/)))]
480       [FLO (fp/ x y)]
481       [BIG (fp/ x (%big->flo y))]
482       ;; TODO: This can be incorrect when the ratnum consists of bignums
483       [RAT (fp/ x (%exact->inexact y))]
484       [COMP (%comp/comp (%make-complex x 0) y)]
485       [else (bad-number '/ y)] ) ]
486    [BIG
487     (switchq (%check-number y)
488       [FIX (let ((g (%gcd-0 '/ x (fix-div/0 x y '/))))
489              (ratnum (%quotient '/ x g) (%quotient '/ y g)))]
490       [FLO (fp/ (%big->flo x) y)]
491       [BIG (let ((g (%gcd-0 '/ x y)))
492              (ratnum (%quotient '/ x g) (%quotient '/ y g)))]
493       [RAT (%nonrat/rat x y)]
494       [COMP (%comp/comp (%make-complex x 0) y)]
495       [else (bad-number '/ y)] ) ]
496    [RAT
497     (switchq (%check-number y)
498       ;; a/b / c/d = a*d / b*c  [generic]
499       ;;   = ((a / g1) * (d / g2) * sign(a)) / abs((b / g2) * (c / g1))
500       ;; With   g1 = gcd(a, c)   and    g2 = gcd(b, d) [Knuth, 4.5.1 ex. 4]
501       [RAT (let* ((a (rat-numerator x)) (b (rat-denominator x))
502                   (c (rat-numerator y)) (d (rat-denominator y))
503                   (g1 (%gcd-0 '/ a c)) (g2 (%gcd-0 '/ b d)))
504              (%/ (%* (%quotient '/ a g1) (%quotient '/ d g2))
505                  (%* (%quotient '/ b g2) (%quotient '/ c g1))))]
506       [COMP (%comp/comp (%make-complex x 0) y)]
507       ;; TODO: This can be incorrect when the ratnum consists of bignums
508       [FLO (fp/ (%exact->inexact x) y)]
509       [NONE (bad-number '/ y)]
510       ;; a/b / c/d = a*d / b*c  [with d = 1]
511       ;;   = ((a / g) * sign(a)) / abs(b * (c / g))
512       ;; With   g = gcd(a, c)   and  c = y  [Knuth, 4.5.1 ex. 4]
513       [else (let* ((a (rat-numerator x))
514                    (g (%gcd-0 '/ a y))) ;; TODO: Improve error message if /0
515               (%/ (%quotient '/ a g)
516                   (%* (rat-denominator x) (%quotient '/ y g))))] ) ]
517    [COMP
518     (switchq (%check-number y)
519       [COMP (%comp/comp x y)]
520       [NONE (bad-number '/ y)] 
521       [else (%comp/comp x (%make-complex y 0))] ) ]
522    [else (bad-number '/ x)] ) )
523
524(define (%comp/comp p q)
525  (let* ([a (complex-real p)]
526         [b (complex-imag p)]
527         [c (complex-real q)]
528         [d (complex-imag q)]
529         [r (%+ (%* c c) (%* d d))]
530         [x (%/ (%+ (%* a c) (%* b d)) r)]
531         [y (%/ (%- (%* b c) (%* a d)) r)] )
532    (make-complex x y) ) )
533
534
535;;; Comparisons:
536
537(define (%= x y)
538  (switchq (%check-number x)
539    [FIX 
540     (switchq (%check-number y)
541       [FIX (fx= x y)]
542       [FLO (and (finite? y)
543                 (%= x (%flo->rat '= y)))] ; Compare as ratnums (overflow)
544       [BIG #f] ;; Needs bignum representation?  Can't be equal to a fixnum!
545       [RAT #f] ;; Rats are never x/1, because those are normalised to just x
546       [COMP #f] ;; Comps are only ever equal to other comps
547       [else (bad-number '= y)] ) ]
548    [FLO 
549     (switchq (%check-number y)
550       [FIX (and (finite? x)
551                 (%= (%flo->rat '= x) y))] ; Compare as ratnums (overflow)
552       [FLO (fp= x y)]
553       [BIG (and (%flo-integer? x) (= (%flo->integer x) y))]
554       [RAT (and (not (or (fp= x +inf.0) (fp= x -inf.0)))
555                 (%= (%flo->rat '= x) y))] ; Compare as ratnums
556       [COMP #f] ;; Comps are only ever equal to other comps
557       [else (bad-number '= y)] ) ]
558    [BIG
559     (switchq (%check-number y)
560       [FIX #f]  ;; Needs bignum representation?  Can't be equal to a fixnum!
561       [FLO (and (%flo-integer? y) (= x (%flo->integer y)))]
562       [BIG (fx= (%big-cmp x y) 0)]
563       [RAT #f] ;; Rats are never x/1, because those are normalised to just x
564       [COMP #f] ;; Comps are only ever equal to other comps
565       [else (bad-number '= y)] ) ]
566    [RAT
567     (switchq (%check-number y)
568       [FIX #f] ;; Rats are never x/1, because those are normalised to just x
569       [FLO (and (not (or (fp= y +inf.0) (fp= y -inf.0)))
570                 (%= x (%flo->rat '= y)))] ; Compare as ratnums
571       [BIG #f] ;; Rats are never x/1, because those are normalised to just x
572       ;; TODO: Use integer= here, when we write it
573       [RAT (and (%= (rat-numerator x) (rat-numerator y))
574                 (%= (rat-denominator x) (rat-denominator y)))]
575       [COMP #f] ;; Comps are only ever equal to other comps
576       [else (bad-number '= y)] ) ]
577    [COMP
578     (switchq (%check-number y)
579       [COMP (and (%= (complex-real x) (complex-real y))
580                  (%= (complex-imag x) (complex-imag y)))]
581       [NONE (bad-number '= y)]
582       [else #f] ) ]
583    [else (bad-number '= x)] ))
584
585(define (= x1 x2 . xs)
586  (and (%= x1 x2)
587       (let loop ([x x2] [xs xs])
588         (or (null? xs)
589             (let ([h (##sys#slot xs 0)])
590               (and (%= x h)
591                    (loop h (##sys#slot xs 1)) ) ) ) ) ))
592
593(define (eqv? a b)
594  (let ((ta (%check-number a))
595        (tb (%check-number b)))
596    ;; If both are numbers of the same type and exactness, compare.
597    ;; Otherwise use eq? Characters are already compared correctly by eq?
598    (and (eq? ta tb)
599         (switchq ta
600           (NONE (eq? a b))
601           (FLO  (fp= a b))
602           (FIX  (fx= a b))
603           (BIG  (fx= (%big-cmp a b) 0))
604           ;; TODO: Use integer= here, when we write it
605           (RAT  (and (%= (rat-numerator a) (rat-numerator b))
606                      (%= (rat-denominator a) (rat-denominator b))))
607           ;; We use eqv? here because exactness of components needs to match
608           (COMP (and (eqv? (complex-real a) (complex-real b))
609                      (eqv? (complex-imag a) (complex-imag b))))
610           (else (error "This should not happen"))))))
611
612(define (> x1 x2 . xs)
613  (and (%> x1 x2 '>)
614       (let loop ([x x2] [xs xs])
615         (or (null? xs)
616             (let ([h (##sys#slot xs 0)])
617               (and (%> x h '>)
618                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )
619
620(define (%> x y loc)
621  (switchq (%check-number x)
622    (FIX 
623     (switchq (%check-number y)
624       (FIX (fx> x y))
625       ;; Compare as ratnum, to prevent overflows
626       (FLO (or (fp= y -inf.0)
627                (and (not (fp= y +inf.0)) (fp= y y)
628                     (%> x (%flo->rat loc y) loc))))
629       ;;   x neg?   y neg?   x > y?   reason
630       ;;  ---------------------------------------------------------------
631       ;;     no       no       no     abs(y) > abs(x), both positive
632       ;;     no      yes      yes     (a > b)  true if  (a > 0) & (b < 0)
633       ;;    yes       no       no     (a > b)  false if (a < 0) & (b > 0)
634       ;;    yes      yes      yes     abs(y) > abs(x), both negative
635       ;;
636       ;; It follows that x is only bigger than y when y is negative
637       (BIG (%big-negative? y))
638       ;; a/b > c/d  when  a*d > b*c  [with b = 1]
639       (RAT (%> (%* x (rat-denominator y))
640                (rat-numerator y) loc))
641       (COMP (bad-complex/o loc y))
642       (else (bad-number loc y)) ) )
643    (FLO
644     (switchq (%check-number y)
645       (FLO (fp> x y))
646       (COMP (bad-complex/o loc y))
647       (NONE (bad-number loc y))
648       ;; Compare as ratnums, to avoid errors when overflowing
649       ;; (this can happen for bignums, but also for fixnums on 64-bit)
650       (else (or (fp= x +inf.0)
651                 (and (not (fp= x -inf.0)) (fp= x x)
652                      (%> (%flo->rat loc x) y loc)))) ) )
653    (BIG 
654     (switchq (%check-number y)
655       ;;   x neg?   y neg?   x > y?   reason
656       ;;  ---------------------------------------------------------------
657       ;;     no       no      yes     abs(x) > abs(y), both positive
658       ;;     no      yes      yes     (a > b)  true if  (a > 0) & (b < 0)
659       ;;    yes       no       no     (a > b)  false if (a < 0) & (b > 0)
660       ;;    yes      yes       no     abs(x) > abs(y), both negative
661       ;;
662       ;; It follows that x is only bigger than y when x is not negative
663       (FIX (not (%big-negative? x)))
664       (FLO (or (fp= y -inf.0)
665                (and (not (fp= y +inf.0)) (fp= y y)
666                     (%> x (%flo->rat loc y) loc)))) ; Compare as ratnums
667       (BIG (fx> (%big-cmp x y) 0))
668       ;; a/b > c/d  when  a*d > b*c  [with b = 1]
669       (RAT (%> (%* x (rat-denominator y))
670                (rat-numerator y) loc))
671       (COMP (bad-complex/o loc y))
672       (else (bad-number loc y)) ) )
673    (RAT
674     (switchq (%check-number y)
675       ;; a/b > c/d  when  a*d > b*c  [generic]
676       (RAT (%> (%* (rat-numerator x) (rat-denominator y))
677                (%* (rat-denominator x) (rat-numerator y)) loc))
678       (FLO (or (fp= y -inf.0)
679                (and (not (fp= y +inf.0)) (fp= y y)
680                     (%> x (%flo->rat loc y) loc)))) ; Compare as ratnums
681       (COMP (bad-complex/o loc y))
682       (NONE (bad-number loc y))
683       ;; a/b > c/d  when  a*d > b*c  [with d = 1]
684       (else (%> (rat-numerator x)
685                 (%* (rat-denominator x) y) loc)) ) )
686    (COMP (bad-complex/o loc x))
687    (else (bad-number loc x)) ) )
688
689(define (< x1 x2 . xs)
690  (and (%< x1 x2 '<)
691       (let loop ([x x2] [xs xs])
692         (or (null? xs)
693             (let ([h (##sys#slot xs 0)])
694               (and (%< x h '<)
695                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )
696
697(define (%< x y loc)
698  (switchq (%check-number x)
699    (FIX 
700     (switchq (%check-number y)
701       (FIX (fx< x y))
702       ;; Compare as ratnum, to prevent overflows
703       (FLO (or (fp= y +inf.0)
704                (and (not (fp= y -inf.0)) (fp= y y)
705                     (%< x (%flo->rat loc y) loc))))
706       ;;   x neg?   y neg?   x < y?   reason
707       ;;  ---------------------------------------------------------------
708       ;;     no       no      yes     abs(x) < abs(y), both positive
709       ;;     no      yes       no     (a < b)  false if  (a > 0) & (b < 0)
710       ;;    yes       no      yes     (a < b)  true if (a < 0) & (b > 0)
711       ;;    yes      yes       no     abs(x) < abs(y), both negative
712       ;;
713       ;; It follows that x is only smaller than y when y is not negative
714       (BIG (not (%big-negative? y)))
715       ;; a/b < c/d  when  a*d < b*c  [with b = 1]
716       (RAT (%< (%* x (rat-denominator y))
717                (rat-numerator y) loc))
718       (COMP (bad-complex/o loc y))
719       (else (bad-number loc y)) ) )
720    (FLO
721     (switchq (%check-number y)
722       (FLO (fp< x y))
723       (COMP (bad-complex/o loc y))
724       (NONE (bad-number loc y))
725       ;; Compare as ratnums, to avoid errors when overflowing
726       ;; (this can happen for bignums, but also for fixnums on 64-bit)
727       (else (or (fp= x -inf.0)
728                (and (not (fp= x +inf.0)) (fp= x x)
729                     (%< (%flo->rat loc x) y loc))))) )
730    (BIG 
731     (switchq (%check-number y)
732       ;;   x neg?   y neg?   x < y?   reason
733       ;;  ---------------------------------------------------------------
734       ;;     no       no       no     abs(y) < abs(x), both positive
735       ;;     no      yes       no     (a < b)  false if  (a > 0) & (b < 0)
736       ;;    yes       no      yes     (a < b)  true if (a < 0) & (b > 0)
737       ;;    yes      yes      yes     abs(y) < abs(x), both negative
738       ;;
739       ;; It follows that x is only smaller than y when x is negative
740       (FIX (%big-negative? x))
741       (FLO (or (fp= y +inf.0)
742                (and (not (fp= y -inf.0)) (fp= y y)
743                     (%< x (%flo->rat loc y) loc)))) ; Compare as ratnums
744       (BIG (fx< (%big-cmp x y) 0))
745       ;; a/b < c/d  when  a*d < b*c  [with b = 1]
746       (RAT (%< (%* x (rat-denominator y))
747                (rat-numerator y) loc))
748       (COMP (bad-complex/o loc y))
749       (else (bad-number loc y)) ) )
750    (RAT
751     (switchq (%check-number y)
752       ;; a/b < c/d  when  a*d < b*c  [generic]
753       (RAT (%< (%* (rat-numerator x) (rat-denominator y))
754                (%* (rat-denominator x) (rat-numerator y)) loc))
755       (COMP (bad-complex/o loc y))
756       (FLO (or (fp= y +inf.0)
757                (and (not (fp= y -inf.0)) (fp= y y)
758                     (%< x (%flo->rat loc y) loc)))) ; Compare as ratnums
759       (NONE (bad-number loc y))
760       ;; a/b < c/d  when  a*d < b*c  [with d = 1]
761       (else (%< (rat-numerator x)
762                 (%* (rat-denominator x) y) loc)) ) )
763    (COMP (bad-complex/o loc x))
764    (else (bad-number loc x)) ) )
765
766(define (>= x1 x2 . xs)
767  (and (or (not (%flonum? x1)) (fp= x1 x1))
768       (or (not (%flonum? x2)) (fp= x2 x2))
769       (not (%< x1 x2 '>=))
770       (let loop ([x x2] [xs xs])
771         (or (null? xs)
772             (let ([h (##sys#slot xs 0)])
773               (and (or (not (%flonum? h)) (fp= h h)) ; not nan?
774                    (not (%< x h '>=))
775                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )
776
777(define (<= x1 x2 . xs)
778  (and (or (not (%flonum? x1)) (fp= x1 x1))
779       (or (not (%flonum? x2)) (fp= x2 x2))
780       (not (%> x1 x2 '<=))
781       (let loop ([x x2] [xs xs])
782         (or (null? xs)
783             (let ([h (##sys#slot xs 0)])
784               (and (or (not (%flonum? h)) (fp= h h)) ; not nan?
785                    (not (%> x h '<=))
786                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )
787
788
789;;; Complex numbers
790
791(define (make-complex r i)
792  (if (or (eq? i 0) (and (%flonum? i) (fp= i 0.0)))
793      r
794      (%make-complex (if (inexact? i) (exact->inexact r) r)
795                     (if (inexact? r) (exact->inexact i) i)) ) )
796
797(define (make-rectangular r i)
798  (switchq (%check-number r)
799    (COMP (bad-real 'make-rectangular r))
800    (NONE (bad-number 'make-rectangular r)) )
801  (switchq (%check-number i)
802    (COMP (bad-real 'make-rectangular i))
803    (NONE (bad-number 'make-rectangular i)) )
804  (make-complex r i) )
805
806(define (%make-polar r phi)
807  (switchq (%check-number r)
808    (COMP (bad-real 'make-polar r))
809    (NONE (bad-number 'make-polar r)) )
810  (switchq (%check-number phi)
811    (COMP (bad-real 'make-polar phi))
812    (NONE (bad-number 'make-polar phi)) )
813  (let ((fphi (exact->inexact phi)))
814    (make-complex (%* r (##core#inline_allocate ("C_a_i_cos" 4) fphi))
815                  (%* r (##core#inline_allocate ("C_a_i_sin" 4) fphi)))))
816
817(define make-polar %make-polar)
818
819(define (real-part x)
820  (switchq (%check-number x)
821    (COMP (complex-real x))
822    (NONE (bad-number 'real-part x))
823    (else x) ) )
824
825(define (imag-part x)
826  (switchq (%check-number x)
827    (NONE (bad-number 'imag-part x))
828    (COMP (complex-imag x))
829    (FLO 0.0)
830    (else 0) ) )
831
832(define (%magnitude x)
833  (switchq (%check-number x)
834    (COMP (let ((r (complex-real x))
835                (i (complex-imag x)) )
836            (%sqrt 'magnitude (%+ (%* r r) (%* i i))) ) )
837    (NONE (bad-number 'magnitude x))
838    (else (%abs x)) ) )
839
840(define magnitude %magnitude)
841
842(define (%angle x)
843  (switchq (%check-number x)
844    (NONE (bad-number 'angle x))
845    (COMP (##core#inline_allocate ("C_a_i_atan2" 4)
846                                  (%exact->inexact (complex-imag x))
847                                  (%exact->inexact (complex-real x))))
848    (else (##core#inline_allocate ("C_a_i_atan2" 4) 0.0 (%exact->inexact x))) ) )
849
850(define angle %angle)
851
852
853;;; Rationals
854
855(define (ratnum m n)
856  (cond
857   ((eq? n 1) m)
858   ((eq? n -1) (- m))
859   ((negative? n) (%make-rat (- m) (- n)))
860   (else (%make-rat m n))))
861
862;; Knuth, 4.5.1
863(define (rat+/- loc op x y)
864  (let ((a (rat-numerator x)) (b (rat-denominator x))
865        (c (rat-numerator y)) (d (rat-denominator y)))
866    (let ((g1 (%gcd-0 loc b d)))
867      (cond
868       ((eq? g1 1) (%make-rat (op (%* a d) (%* b c)) (%* b d)))
869       ;; Save a quotient and multiplication if the gcd is equal
870       ;; to one of the denominators since quotient of b or d and g1 = 1
871       ;; TODO: Check properties of the gcd to see if g2 and t are needed
872       ((%= g1 b) (let* ((t (op (%* a (%quotient loc d g1)) c))
873                         (g2 (%gcd-0 loc t g1)))
874                    (ratnum (%quotient loc t g2) (%quotient loc d g2))))
875       ((%= g1 d) (let* ((b/g1 (%quotient loc b g1))
876                         (t (op a (%* c b/g1))) ;; Is this worth it?
877                         (g2 (%gcd-0 loc t g1)))
878                    (ratnum (%quotient loc t g2)
879                            (%* b/g1 (%quotient loc d g2)))))
880       (else (let* ((b/g1 (%quotient loc b g1))
881                    (t (op (%* a (%quotient loc d g1))
882                           (%* c b/g1)))
883                    (g2 (%gcd-0 loc t g1)))
884               (%make-rat (%quotient loc t g2)
885                          (%* b/g1 (%quotient loc d g2)))))))))
886
887(define (numerator x)
888  (switchq (%check-number x)
889    (FIX x)
890    (FLO (cond
891          ((not (finite? x)) (bad-inexact 'numerator x))
892          ((%flo-integer? x) x)
893          (else (exact->inexact (numerator (%flo->rat 'numerator x))))))
894    (BIG x)
895    (RAT (rat-numerator x))
896    (COMP (bad-ratnum 'numerator x))
897    (else (bad-number 'numerator x)) ) )
898
899(define (denominator x)
900  (switchq (%check-number x)
901    (FIX 1)
902    (FLO (cond
903          ((not (finite? x)) (bad-inexact 'denominator x))
904          ((%flo-integer? x) 1.0)
905          (else (exact->inexact (denominator (%flo->rat 'denominator x))))))
906    (BIG 1)
907    (RAT (rat-denominator x))
908    (COMP (bad-ratnum 'denominator x))
909    (else (bad-number 'denominator x)) ) )
910
911
912;;; Enhanced versions of other standard procedures
913
914(define (%abs x)
915  (switchq (%check-number x)
916    (FIX (if (fx< x 0) (%fix-neg x) x))
917    (FLO (##core#inline_allocate ("C_a_i_abs" 4) x))
918    (BIG (%big-abs x))
919    (RAT (%make-rat (%abs (rat-numerator x)) (rat-denominator x)))
920    (COMP (##sys#signal-hook #:type-error 'abs "can not compute absolute value of complex number" x))
921    (NONE (bad-number 'abs x)) ) )
922
923(define abs %abs)
924
925(define (number? x)
926  (switchq (%check-number x)
927    (NONE #f)
928    (else #t) ) )
929
930(set! ##sys#number? number?)
931(define complex? number?)
932
933(define (real? x)
934  (switchq (%check-number x)
935    (COMP #f)
936    (NONE #f)
937    (else #t) ) )
938
939(define (rational? x) (and (real? x) (finite? x)))
940
941(define (%integer? x)
942  (switchq (%check-number x)
943    (FIX #t)
944    (FLO (%flo-integer? x))
945    (BIG #t)
946    (else #f) ) )
947
948(set! ##sys#integer? %integer?)
949(define integer? %integer?)
950
951(define (exact-integer? x)
952  (switchq (%check-number x)
953    (FIX #t)
954    (BIG #t)
955    (else #f) ) )
956
957(define (%exact? x)
958  (switchq (%check-number x)
959    (FLO #f)
960    (COMP (and (%exact? (complex-real x)) (%exact? (complex-imag x))))
961    (NONE (bad-number 'exact? x))
962    (else #t) ) )
963
964(define exact? %exact?)
965(define ##sys#exact? %exact?)
966
967(define (%inexact? x)
968  (switchq (%check-number x)
969    (FLO #t)
970    (COMP (and (%inexact? (complex-real x)) (%inexact? (complex-imag x))))
971    (NONE (bad-number 'inexact? x))
972    (else #f) ) )
973
974(define inexact? %inexact?)
975(define ##sys#inexact? %inexact?)
976
977(define (positive? x) (%> x 0 'positive?))
978(define (negative? x) (%< x 0 'negative?))
979
980(define (%zero? x)
981  (switchq (%check-number x)
982    (FIX (eq? x 0))
983    (FLO (fp= x 0.0))
984    (NONE (bad-number 'zero? x))
985    (else #f) ) )
986
987(define zero? %zero?)
988
989(define (odd? x)
990  (switchq (%check-number x)
991    (FIX (##core#inline "C_i_oddp" x))
992    (FLO (##core#inline "C_i_oddp" x))
993    (BIG (%big-odd? x))
994    (else (bad-integer 'odd? x)) ) )
995
996(define (even? x)
997  (switchq (%check-number x)
998    (FIX (##core#inline "C_i_evenp" x))
999    (FLO (##core#inline "C_i_evenp" x))
1000    (BIG (not (%big-odd? x)))
1001    (else (bad-integer 'even? x)) ) )
1002
1003(define (max x1 . xs)
1004  (let ((i (%flonum? x1)))
1005    (let loop ((m x1) (xs xs))
1006      (if (null? xs)
1007          (if i (%exact->inexact m) m)
1008          (let ((h (##sys#slot xs 0)))
1009            (switchq (%check-number h)
1010              (FLO (set! i #t))
1011              (COMP (bad-complex/o 'max h)) )
1012            (loop (if (%> h m 'max) h m) (##sys#slot xs 1)) ) ) ) ) )
1013
1014(define (min x1 . xs)
1015  (let ((i (%flonum? x1)))
1016    (let loop ((m x1) (xs xs))
1017      (if (null? xs)
1018          (if i (%exact->inexact m) m)
1019          (let ((h (##sys#slot xs 0)))
1020            (switchq (%check-number h)
1021              (FLO (set! i #t))
1022              (COMP (bad-complex/o 'min h)) )
1023            (loop (if (%< h m 'min) h m) (##sys#slot xs 1)) ) ) ) ) )
1024
1025;;; XXX TODO:
1026;; Normalise to match the floor/, truncate/, ceiling/, centering/,
1027;; round/, euclidian/ from Guile and Riastradh's division proposal:
1028;; http://people.csail.mit.edu/riastradh/tmp/division.txt
1029;;
1030;; The above should be the canonical names, the others should be
1031;; deprecated.  Unsure whether the full set ought to be supported
1032;; (floor/ and truncate/ for r7rs)
1033
1034(define (%quotient loc x y)
1035  (switchq (%check-number x)
1036    (FIX (switchq (%check-number y)
1037           ;; fix-quotient-big always returns 0 since abs(x) < abs(y)
1038           ;; But take care of MOST_NEGATIVE_FIXNUM (grrr!)
1039           (BIG (if (bignum? (- x)) -1 0))
1040           (RAT (bad-integer loc y)) ; Perhaps convert to flonum?
1041           (NONE (bad-number loc y))
1042           (else (%quotient-0 x y))))
1043    (BIG (switchq (%check-number y)
1044           (FIX (%big-quotient-fix x (fix-div/0 x y loc)))
1045           (BIG (%big-quotient-big x y))
1046           (FLO (if (not (%flo-integer? y))
1047                    (%quotient-0 (%big->flo x) y) ; Could overflow
1048                    (%exact->inexact
1049                     (%quotient loc x (%flo->integer y)))))
1050           (NONE (bad-number loc y))
1051           (else (bad-integer loc y))))
1052    (FLO (switchq (%check-number y)
1053           (BIG (if (%flo-integer? x)
1054                    (%exact->inexact (%quotient loc (%flo->integer x) y))
1055                    (%quotient-0 x (%big->flo y)))) ; Will probably overflow
1056           (RAT (bad-integer loc y))
1057           (NONE (bad-number loc y))
1058           (else (%quotient-0 x y))))
1059    (NONE (bad-number loc x))
1060    (else (bad-integer loc x))))
1061
1062(define (quotient x y) (%quotient 'quotient x y))
1063(define (truncate-quotient x y) (%quotient 'truncate-quotient x y))
1064
1065(define (%remainder loc x y)
1066  (switchq (%check-number x)
1067    [FIX (switchq (%check-number y)
1068           [FIX (fx- x (fx* (fx/ x y) y))]
1069           [FLO (let ((flx (%fix->flo x)))
1070                  (if (%flo-integer? y)
1071                      (fp- flx (fp* (##sys#truncate (fp/ flx y)) y))
1072                      (bad-integer loc y)))]
1073           ;; If abs(x) < abs(y), then remainder is always just x
1074           ;; But again, take care of MOST_NEGATIVE_FIXNUM
1075           [BIG (if (bignum? (- x)) 0 x)]
1076           [else (bad-integer loc y)])]
1077    [FLO (unless (%flo-integer? x)
1078           (bad-integer loc x))
1079         (switchq (%check-number y)
1080           [FLO (if (%flo-integer? y)
1081                    (fp- x (fp* (##sys#truncate (fp/ x y)) y))
1082                    (bad-integer loc y))]
1083           [FIX (let ((fly (%fix->flo (fix-div/0 x y loc))))
1084                  (fp- x (fp* (##sys#truncate (fp/ x fly)) fly)))]
1085           [BIG (%exact->inexact (%remainder loc (%flo->integer x) y))]
1086           [else (bad-integer loc y)])]
1087    [BIG (switchq (%check-number y)
1088           [FIX (%big-remainder-fix x (fix-div/0 x y loc))]
1089           [FLO (if (%flo-integer? y)
1090                    (%exact->inexact (%remainder loc x (%flo->integer y)))
1091                    (bad-integer loc y))]
1092           [BIG (%big-remainder-big x y)]
1093           [else (bad-integer loc y)])]
1094    [else (bad-integer loc x)]) )
1095
1096(define (remainder x y) (%remainder 'remainder x y))
1097(define (truncate-remainder x y) (%remainder 'truncate-remainder x y))
1098
1099;; Modulo's sign follows y  (whereas remainder's sign follows x)
1100(define (%modulo loc x y)
1101   (let ((r (%remainder loc x y)))
1102      (if (%> y 0 loc)
1103          (if (%< r 0 loc)
1104              (%+ r y)
1105              r)
1106          (if (%> r 0 loc)
1107              (%+ r y)
1108              r))))
1109
1110(define (modulo x y) (%modulo 'modulo x y))
1111(define (floor-remainder x y) (%modulo 'modulo x y))
1112
1113(define (%quotient&remainder loc x y)
1114  (switchq (%check-number x)
1115    [FIX (switchq (%check-number y)
1116           [FIX (values (fx/ x y) (%remainder loc x y))]
1117           [FLO (values (quotient x y) (%remainder loc x y))]
1118           ;; If abs(x) < abs(y), then remainder is always just x
1119           ;; But again, take care of MOST_NEGATIVE_FIXNUM
1120           [BIG (if (bignum? (- x)) (values -1 0) (values 0 x))]
1121           [else (bad-integer loc y)])]
1122    [FLO (switchq (%check-number y)
1123           [FLO (values (%quotient loc x y)
1124                        (%remainder loc x y))]
1125           [FIX (values (%quotient loc x y)
1126                        (%remainder loc x y))]
1127           [BIG (if (%flo-integer? x)
1128                    (receive (div rem)
1129                      (%quotient&remainder loc (%flo->integer x) y)
1130                      (values (%exact->inexact div) (%exact->inexact rem)))
1131                    (bad-integer loc x))]
1132           [else (bad-integer loc y)])]
1133    [BIG (switchq (%check-number y)
1134           [FIX (%big-divrem-fix x (fix-div/0 x y loc))]
1135           [FLO (if (%flo-integer? y)
1136                    (receive (div rem)
1137                      (%quotient&remainder loc x (%flo->integer y))
1138                      (values (%exact->inexact div) (%exact->inexact rem)))
1139                    (bad-integer loc y))]
1140           [BIG (%big-divrem-big x y)]
1141           [else (bad-integer loc y)])]
1142    [else (bad-integer loc x)]))
1143
1144(define (quotient&remainder x y) (%quotient&remainder 'quotient&remainder x y))
1145(define (truncate/ x y) (%quotient&remainder 'truncate/ x y))
1146
1147;; Modulo's sign follows y  (whereas remainder's sign follows x)
1148(define (quotient&modulo x y)
1149  (receive (div rem)
1150      (%quotient&remainder 'quotient&modulo x y)
1151    (if (%> y 0 'quotient&modulo)
1152        (if (%< rem 0 'quotient&modulo)
1153            (values div (%+ rem y))
1154            (values div rem))
1155        (if (%> rem 0 'quotient&modulo)
1156            (values div (%+ rem y))
1157            (values div rem)))))
1158
1159;; Same as above, but quotient gets adjusted along with the remainder
1160(define (floor/ x y)
1161  (receive (div rem)
1162      (%quotient&remainder 'floor/ x y)
1163    (if (%> y 0 'floor/)
1164        (if (%< rem 0 'floor/)
1165            (values (%- div 1) (%+ rem y))
1166            (values div rem))
1167        (if (%> rem 0 'floor/)
1168            (values (%- div 1) (%+ rem y))
1169            (values div rem)))))
1170
1171;; XXX This is a bad bad bad definition; very inefficient But to
1172;; improve it we would need to provide another implementation of the
1173;; quotient procedure which floors instead of truncates.
1174(define (floor-quotient x y)
1175  (receive (div rem) (floor/ x y) div))
1176
1177;; Try to multiply by two until we reach an integer
1178(define (%float-fraction-length x)
1179  (do ((x x (fp* x 2.0))
1180       (i 0 (+ i 1)))
1181      ((%flo-integer? x) i)
1182    ;; 3000 -> %float-precision?
1183    (if (> i 3000) (error "I'm bored." x))))
1184
1185(define (%flo->rat loc x)
1186  (define (deliver y d)
1187    (let ((q (expt 2 (%float-fraction-length y))))
1188      (if (%exact? q) ; XXX Ever untrue? float-fraction-length returns natnums
1189          (let ((e (%/ (%/ (%flo->integer
1190                            (%* y (%exact->inexact q)))
1191                           q)
1192                       d)))
1193             (if (%exact? e)
1194                e
1195                (bad-inexact loc x)))
1196          (bad-inexact loc x))))
1197  (if (and (< x (%fix->flo 1))    ; watch out for denormalized numbers
1198           (> x (%fix->flo -1)))
1199      (deliver (%* x (expt (%fix->flo 2) flonum-precision))
1200               ;; Can be bignum (is on 32-bit), so must wait until after init.
1201               ;; We shouldn't need to calculate this every single time, tho..
1202               (expt 2 flonum-precision))
1203      (deliver x 1)))
1204
1205(define (%inexact->exact x)
1206  (switchq (%check-number x)
1207    (FIX x)
1208    (FLO (cond
1209          ((not (finite? x)) (bad-inexact 'inexact->exact x))
1210          ((%flo-integer? x) (%flo->integer x))
1211          (else (%flo->rat 'inexact->exact x))))
1212    (BIG x)
1213    (RAT x)
1214    (COMP (make-complex (%inexact->exact (complex-real x))
1215                        (%inexact->exact (complex-imag x))))
1216    (NONE (bad-number 'inexact->exact x)) ) )
1217
1218(define inexact->exact %inexact->exact)
1219(define ##sys#inexact->exact %inexact->exact)
1220
1221;; Exponent of the lowest allowed flonum; if we get any lower we get zero.
1222;; In other words, this is the first (smallest) flonum after 0.
1223;; Equal to (expt 2.0 (- flonum-minimum-exponent flonum-precision))
1224(define minimum-denorm-flonum-expt (fx- flonum-minimum-exponent flonum-precision))
1225
1226;; This tries to keep the numbers within representable ranges and tries
1227;; to drop as few significant digits as possible by bringing the two numbers
1228;; to within the same powers of two.  See algorithms M & N in Knuth, 4.2.1
1229;;
1230;; TODO: Use (fp/ n d) if both are finite after conversion to flonums
1231(define (%rat->flo x)
1232  (let* ((n1 (rat-numerator x))
1233         (an (%abs n1))
1234         (d1 (rat-denominator x))
1235         ;; Approximate distance between the numbers in powers of 2
1236         ;; ie,  2^e-1 < n/d < 2^e+1  (e is the *un*biased value of e_w in M2)
1237         ;; XXX: What if b != 2 (ie, flonum-radix is not 2)
1238         (e (fx- (integer-length an) (integer-length d1)))
1239         (rnd (lambda (n d e)           ; Here, 1 <= n/d < 2  (normalized) [N5]
1240                ;; Cannot shift above the available precision, and can't have
1241                ;; an exponent that's below the minimum flonum exponent.
1242                (let* ((s (min (fx- flonum-precision 1)
1243                               (fx- e minimum-denorm-flonum-expt)))
1244                       (normalized (%/ (arithmetic-shift n s) d))
1245                       (r (round normalized))
1246                       (fraction (%exact->inexact r))
1247                       (exp (fx- e s)))
1248                  (let ((res (fp* fraction (expt 2.0 exp))))
1249                    (if (negative? n1) (%- 0 res) res)))))
1250         (scale (lambda (n d)                      ; Here, 1/2 <= n/d < 2   [N3]
1251                  (if (%< n d 'exact->inexact)     ; n/d < 1?
1252                      ;; Scale left [N3]; only needed once (see note in M3)
1253                      (rnd (arithmetic-shift n 1) d (fx- e 1))
1254                      ;; Already normalized
1255                      (rnd n d e)))))
1256    ;; After this step, which shifts the smaller number to align with the
1257    ;; larger, "f" in algorithm N is represented in the procedures above by n/d.
1258    (if (negative? e)
1259        (scale (arithmetic-shift an (%- 0 e)) d1)
1260        (scale an (arithmetic-shift d1 e)))))
1261
1262(define (%exact->inexact x)
1263  (switchq (%check-number x)
1264    (FIX (%fix->flo x))
1265    (FLO x)
1266    (BIG (%big->flo x))
1267    (RAT (%rat->flo x))
1268    (COMP (make-complex (%exact->inexact (complex-real x)) (%exact->inexact (complex-imag x))))
1269    (NONE (bad-number 'exact->inexact x)) ) )
1270
1271(define exact->inexact %exact->inexact)
1272(define inexact %exact->inexact)
1273(define ##sys#exact->inexact %exact->inexact)
1274(define exact %inexact->exact)
1275
1276(define (%gcd-0 loc x y)
1277  (switchq (%check-number x)
1278    [FIX (switchq (%check-number y)
1279           [FIX (fxgcd x y)]
1280           [FLO (if (%flo-integer? y)
1281                    (fpgcd (%fix->flo x) y)
1282                    (bad-integer loc y))]
1283           [BIG (if (eq? x 0) y (fxgcd x (%remainder loc y x)))]
1284           [else (bad-integer loc y)])]
1285    [FLO (switchq (%check-number y)
1286           [FIX (if (%flo-integer? x)
1287                    (fpgcd x (%fix->flo y))
1288                    (bad-integer loc x))]
1289           [FLO (if (%flo-integer? x)
1290                    (if (%flo-integer? y)
1291                        (fpgcd x y)
1292                        (bad-integer loc x))
1293                    (bad-integer loc x))]
1294           [BIG (if (fp= x 0.0) y (fpgcd x (%remainder loc y x)))]
1295           [else (bad-integer loc y)])]
1296    [BIG (switchq (%check-number y)
1297           [FIX (if (eq? y 0) x (fxgcd y (%remainder loc x y)))]
1298           [FLO (if (fp= y 0.0) x (fpgcd y (%remainder loc x y)))]
1299           [BIG (biggcd x y)]
1300           [else (bad-integer loc y)])]
1301    [else (bad-integer loc x)]) )
1302
1303(define (gcd . ns)
1304  (if (eq? ns '())
1305      0
1306      (let loop ([ns ns] [f #t])
1307        (let ([head (##sys#slot ns 0)]
1308              [next (##sys#slot ns 1)] )
1309          (if (null? next)
1310              (if f (%abs (%->integer 'gcd head)) (%abs head))
1311              (let ([n2 (##sys#slot next 0)])
1312                (loop (cons (%gcd-0 'gcd head n2) (##sys#slot next 1)) #f) ) ) ) ) ) )
1313
1314(define (%lcm-0 loc x y)
1315  (%quotient loc (%* x y) (%gcd-0 loc x y)) )
1316
1317(define (lcm . ns)
1318  (if (null? ns)
1319      1
1320      (let loop ([ns ns] [f #t])
1321        (let ([head (##sys#slot ns 0)]
1322              [next (##sys#slot ns 1)] )
1323          (if (null? next)
1324              (if f (%abs (%->integer 'lcm head)) (%abs head))
1325              (let ([n2 (##sys#slot next 0)])
1326                (loop (cons (%lcm-0 'lcm head (##sys#slot next 0)) (##sys#slot next 1)) #f) ) ) ) ) ) )
1327
1328(define (%floor x)
1329  (switchq (%check-number x)
1330    (FIX x)
1331    (FLO (##sys#floor x))
1332    (BIG x)
1333    ;; (floor x) = greatest integer <= x
1334    (RAT (let* ((n (rat-numerator x))
1335                (q (quotient n (rat-denominator x))))
1336           (if (>= n 0)
1337               q
1338               (%- q 1))))
1339    (else (bad-real 'floor x))) )
1340
1341(define floor %floor)
1342
1343(define (ceiling x)
1344  (switchq (%check-number x)
1345    (FIX x)
1346    (FLO (##sys#ceiling x))
1347    (BIG x)
1348    ;; (ceiling x) = smallest integer >= x
1349    (RAT (let* ((n (rat-numerator x))
1350                (q (quotient n (rat-denominator x))))
1351           (if (>= n 0)
1352               (%+ q 1)
1353               q)))
1354    (else (bad-real 'ceiling x))) )
1355
1356(define (truncate x)
1357  (switchq (%check-number x)
1358    (FIX x)
1359    (FLO (##sys#truncate x))
1360    (BIG x)
1361    ;; (rational-truncate x) = integer of largest magnitude <= (abs x)
1362    (RAT (%quotient 'truncate (rat-numerator x) (rat-denominator x)))
1363    (else (bad-real 'truncate x))) )
1364
1365(define (round x)
1366  (switchq (%check-number x)
1367    (FIX x)
1368    (FLO (##core#inline_allocate ("C_a_i_flonum_round_proper" 4) x))
1369    (BIG x)
1370    (RAT (let* ((x+1/2 (%+ x (%make-rat 1 2)))
1371                (r (%floor x+1/2)))
1372           (if (and (%= r x+1/2)
1373                    (odd? r))
1374               (%- r 1)
1375               r)))
1376    (else (bad-real 'round x)) ) )
1377
1378(define (find-ratio-between x y)
1379  (define (sr x y)
1380    (let ((fx (%inexact->exact (%floor x))) 
1381          (fy (%inexact->exact (%floor y))))
1382      (cond ((not (%< fx x 'rationalize)) (list fx 1))
1383            ((%= fx fy) 
1384             (let ((rat (sr (%/ 1 (%- y fy)) (%/ 1 (%- x fx)))))
1385               (list (%+ (cadr rat) (%* fx (car rat))) (car rat))))
1386            (else (list (%+ 1 fx) 1)))))
1387  (cond ((%< y x 'rationalize)
1388         (find-ratio-between y x))
1389        ((not (%< x y 'rationalize))
1390         (list x 1))
1391        ((%> x 0 'rationalize)
1392         (sr x y))
1393        ((%< y 0 'rationalize) 
1394         (let ((rat (sr (%- 0 y) (%- 0 x))))
1395           (list (%- 0 (car rat)) (cadr rat))))
1396        (else '(0 1))))
1397
1398(define (find-ratio x e) (find-ratio-between (%- x e) (%+ x e)))
1399
1400(define (rationalize x e)
1401  (let ((result (apply %/ (find-ratio x e))))
1402    (if (or (inexact? x) (inexact? e))
1403        (exact->inexact result)
1404        result)))
1405
1406(define (%exp n)
1407  (switchq (%check-number n)
1408    (NONE (bad-number 'exp n))
1409    (COMP (%* (##core#inline_allocate ("C_a_i_exp" 4) (complex-real n))
1410              (let ((p (complex-imag n)))
1411                (make-complex
1412                 (##core#inline_allocate ("C_a_i_cos" 4) p)
1413                 (##core#inline_allocate ("C_a_i_sin" 4) p) ) ) ) )
1414    (else (##core#inline_allocate ("C_a_i_exp" 4) (%exact->inexact n)) ) ))
1415
1416(define exp %exp)
1417
1418(define (%log x)
1419  (let ((type (%check-number x)))
1420    (cond
1421     ;; avoid calling inexact->exact on X here (to avoid overflow?)
1422     ((or (eq? type COMP) (%< x 0.0 'log)) ; General case
1423      (%+ (%log (%magnitude x)) (* (make-complex 0 1) (%angle x))))
1424     ((eq? x 0)                        ; Exact zero?  That's undefined
1425      (log0 'log x))
1426     ((eq? type NONE)
1427      (bad-number 'log x))
1428     (else                             ; Simple real number case
1429      (##core#inline_allocate ("C_a_i_log" 4) (%exact->inexact x))))))
1430
1431(define (log a #!optional b)
1432  (if b (%/ (%log a) (%log b)) (%log a)))
1433
1434(define %i (%make-complex 0 1))
1435(define %-i (%make-complex 0 -1))
1436(define %i2 (%make-complex 0 2))
1437
1438(define (%sin n)
1439  (switchq (%check-number n)
1440    (NONE (bad-number 'sin n))
1441    (COMP (let ((in (%* %i n)))
1442            (%/ (%- (%exp in) (%exp (%- 0 in))) %i2)))
1443    (else (##core#inline_allocate ("C_a_i_sin" 4) (%exact->inexact n)) ) ))
1444
1445(define sin %sin)
1446
1447(define (%cos n)
1448  (switchq (%check-number n)
1449    (NONE (bad-number 'cos n))
1450    (COMP (let ((in (%* %i n)))
1451            (%/ (%+ (%exp in) (%exp (%- 0 in))) 2) ) )
1452    (else (##core#inline_allocate ("C_a_i_cos" 4) (%exact->inexact n)) ) ) )
1453
1454(define cos %cos)
1455
1456(define (tan n)
1457  (switchq (%check-number n)
1458    (NONE (bad-number 'tan n))
1459    (COMP (%/ (%sin n) (%cos n)))
1460    (else (##core#inline_allocate ("C_a_i_tan" 4) (%exact->inexact n)) ) ))
1461
1462;; General case: sin^{-1}(z) = -i\ln(iz + \sqrt{1-z^2})
1463(define (%asin n)
1464  (let ((t (%check-number n)))
1465    (cond ((eq? t NONE) (bad-number 'asin n))
1466          ((and (eq? t FLO) (fp>= n -1.0) (fp<= n 1.0))
1467           (##core#inline_allocate ("C_a_i_asin" 4) n))
1468          ((and (eq? t FIX) (fx>= n -1) (fx<= n 1))
1469           (##core#inline_allocate ("C_a_i_asin" 4) (%fix->flo n)))
1470          ;; General definition can return compnums
1471          (else (%* %-i (%log (%+ (%* %i n) (%sqrt 'asin (%- 1 (%* n n))))))))))
1472
1473(define asin %asin)
1474
1475;; General case:
1476;; cos^{-1}(z) = 1/2\pi + i\ln(iz + \sqrt{1-z^2}) = 1/2\pi - sin^{-1}(z) = sin(1) - sin(z)
1477(define %acos
1478  (let ((asin1 (##core#inline_allocate ("C_a_i_asin" 4) 1)))
1479    (lambda (n)
1480      (let ((t (%check-number n)))
1481        (cond ((eq? t NONE) (bad-number 'acos n))
1482              ((and (eq? t FLO) (fp>= n -1.0) (fp<= n 1.0))
1483               (##core#inline_allocate ("C_a_i_acos" 4) n))
1484              ((and (eq? t FIX) (fx>= n -1) (fx<= n 1))
1485               (##core#inline_allocate ("C_a_i_acos" 4) (%fix->flo n)))
1486              ;; General definition can return compnums
1487              (else (%- asin1 (%asin n)))) ) ) ) )
1488
1489(define acos %acos)
1490
1491(define (atan n #!optional b)
1492  (switchq (%check-number n)
1493    (NONE (bad-number 'atan n))
1494    (COMP (if b
1495              (bad-real 'atan n)
1496              (let ((in (%* %i n)))
1497                (%/ (%- (%log (%+ 1 in)) (%log (%- 1 in))) %i2) ) ) )
1498    (else (if b
1499              (##core#inline_allocate ("C_a_i_atan2" 4) (%exact->inexact n) (%exact->inexact b))
1500              (##core#inline_allocate ("C_a_i_atan" 4) (%exact->inexact n))))) )
1501
1502(define (%exact-integer-sqrt loc k)
1503  (if (or (eq? 0 k) (eq? 1 k))
1504      (values k 0)
1505      ;; Hacker's Delight, figure 11-1 (Newton's method - see also SICP 1.1.7)
1506      ;; Initial guess is about the smallest x^2 which is bigger than sqrt(k)
1507      ;; We could subtract one from ilen(k) to avoid overshooting by 1, but
1508      ;; that's one more operation.
1509      (let* ((len (fx+ (quotient (integer-length k) 2) 1))
1510             (g0 (arithmetic-shift 1 len)))
1511        (let lp ((g0 g0)
1512                 (g1 (arithmetic-shift
1513                      (%+ g0 (arithmetic-shift k (fxneg len))) -1)))
1514          (if (%< g1 g0 loc)
1515              (lp g1 (arithmetic-shift (%+ g1 (quotient k g1)) -1))
1516              (values g0 (%- k (%* g0 g0))))))))
1517
1518(define (exact-integer-sqrt x)
1519  (switchq (%check-number x)
1520    (NONE (bad-number 'exact-integer-sqrt x))
1521    (FIX (if (fx< x 0)
1522             (bad-natural 'exact-integer-sqrt x)
1523             (%exact-integer-sqrt 'exact-integer-sqrt x)))
1524    (BIG (if (%big-negative? x)
1525             (bad-natural 'exact-integer-sqrt x)
1526             (%exact-integer-sqrt 'exact-integer-sqrt x)))
1527    (else (bad-natural 'exact-integer-sqrt x))))
1528
1529(define (%sqrt loc n)
1530  (switchq (%check-number n)
1531    (NONE (bad-number 'sqrt n))
1532    (COMP (let ((p (%/ (%angle n) 2))
1533                (m (##core#inline_allocate ("C_a_i_sqrt" 4) (%magnitude n))) )
1534            (make-complex (%* m (%cos p)) (%* m (%sin p)) ) ) )
1535    (RAT (let ((num (rat-numerator n))
1536               (den (rat-denominator n)))
1537           (if (and (>= num 0) (>= den 0))
1538               (receive (ns^2 nr)
1539                 (%exact-integer-sqrt loc num)
1540                 (if (eq? nr 0)
1541                     (receive (ds^2 dr)
1542                       (%exact-integer-sqrt loc den)
1543                       (if (eq? dr 0)
1544                           (%/ ns^2 ds^2)
1545                           (%sqrt loc (%exact->inexact n))))
1546                     (%sqrt loc (%exact->inexact n))))
1547               (%sqrt loc (%exact->inexact n)))))
1548    (else
1549     (cond
1550      ((negative? n)
1551       (make-complex
1552        0.0
1553        (##core#inline_allocate ("C_a_i_sqrt" 4) (%exact->inexact (- n)))))
1554      ((integer? n)
1555       (receive (s^2 r)
1556         (%exact-integer-sqrt loc (%->integer loc n))
1557         (if (eq? 0 r)
1558             (if (exact? n) s^2 (%exact->inexact s^2))
1559             (##core#inline_allocate ("C_a_i_sqrt" 4) (%exact->inexact n)))))
1560      (else (##core#inline_allocate ("C_a_i_sqrt" 4) (%exact->inexact n))) ) )))
1561
1562(define (sqrt x) (%sqrt 'sqrt x))
1563
1564(define (square x) (%* x x))
1565
1566;; Generalized Newton's algorithm for positive integers, with a little help
1567;; from Wikipedia ;)  https://en.wikipedia.org/wiki/Nth_root_algorithm
1568(define (%exact-integer-nth-root loc k n)
1569  (cond
1570   ((or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exact-integer-sqrt on n=2?
1571    (values k 0))
1572   ((< n 1)
1573    (bad-positive loc n))
1574   (else
1575    (let ((len (integer-length k)))
1576      (if (< len n)       ; Idea from Gambit: 2^{len-1} <= k < 2^{len}
1577          (values 1 (- k 1)) ; Since we know x >= 2, we know x^{n} can't exist
1578          ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n)
1579          (let* ((shift-amount (ceiling (/ (fx+ len 1) n)))
1580                 (g0 (arithmetic-shift 1 shift-amount))
1581                 (n-1 (%- n 1)))
1582            (let lp ((g0 g0)
1583                     (g1 (%quotient loc (%+ (%* n-1 g0) (%quotient loc k (%integer-power g0 n-1))) n)))
1584              (if (%< g1 g0 loc)
1585                  (lp g1 (%quotient loc (%+ (%* n-1 g1) (%quotient loc k (%integer-power g1 n-1))) n))
1586                  (values g0 (%- k (%integer-power g0 n)))))))))))
1587
1588(define (exact-integer-nth-root k n)
1589  (unless (exact-integer? n)
1590    (bad-natural 'exact-integer-nth-root n))
1591  (switchq (%check-number k)
1592    (NONE (bad-number 'exact-integer-nth-root k))
1593    (FIX (if (fx< k 0)
1594             (bad-natural 'exact-integer-nth-root k)
1595             (%exact-integer-nth-root 'exact-integer-nth-root k n)))
1596    (BIG (if (%big-negative? k)
1597             (bad-natural 'exact-integer-nth-root k)
1598             (%exact-integer-nth-root 'exact-integer-nth-root k n)))
1599    (else (bad-natural 'exact-integer-nth-root k))))
1600
1601(define (%integer-power base e)
1602  (if (negative? e)
1603      (%/ 1 (%integer-power base (- e)))
1604      (let lp ((res 1) (e2 e))
1605        (cond
1606          ((eq? e2 0) res)
1607          ((even? e2) ; recursion is faster than iteration here
1608           (%* res (square (lp 1 (arithmetic-shift e2 -1)))))
1609          (else
1610           (lp (%* res base) (%- e2 1)))))))
1611
1612(define (expt a b)
1613  (define (slow-expt a b)
1614    (if (eq? 0 a)
1615        (expt0 'expt a b)
1616        (%exp (%* b (%log a)))))
1617  (let ((ta (%check-number a))
1618        (tb (%check-number b)) )
1619    (cond ((eq? NONE ta) (bad-number 'expt a))
1620          ((eq? NONE tb) (bad-number 'expt b))
1621          ((and (eq? RAT ta) (not (inexact? b)))
1622           ;; (n*d)^b = n^b * d^b = n^b * x^{-b}  | x = 1/b
1623           ;; Hopefully faster than integer-power
1624           (%* (expt (rat-numerator a) b) (expt (rat-denominator a) (- b))))
1625          ;; x^{a/b} = (x^{1/b})^a
1626          ((eq? RAT tb)
1627           (switchq ta
1628             (FIX (if (fx< a 0)
1629                      (%expt-0 (%fix->flo a) (%rat->flo b))
1630                      (receive (ds^n r)
1631                        (%exact-integer-nth-root 'expt a (rat-denominator b))
1632                        (if (eq? r 0)
1633                            (expt ds^n (rat-numerator b))
1634                            (%expt-0 (%fix->flo a) (%rat->flo b))))))
1635             (BIG (if (%big-negative? a)
1636                      (%expt-0 (%big->flo a) (%rat->flo b))
1637                      (receive (ds^n r)
1638                        (%exact-integer-nth-root 'expt a (rat-denominator b))
1639                        (if (eq? r 0)
1640                            (expt ds^n (rat-numerator b))
1641                            (%expt-0 (%big->flo a) (%rat->flo b))))))
1642             (FLO (%expt-0 a (%rat->flo b)))
1643             (else (slow-expt a b))))
1644          ((or (eq? COMP tb) (and (eq? COMP ta) (not (integer? b))))
1645           (slow-expt a b))
1646          ((or (eq? FLO ta) (and (eq? FLO tb) (not (%flo-integer? b))))
1647           (%expt-0 (%exact->inexact a) (%exact->inexact b)))
1648          ;; this doesn't work that well, yet...
1649          ;; (XXX: What does this mean? why not? I do know this is ugly... :P)
1650          (else (if (or (inexact? a) (inexact? b))
1651                    (%exact->inexact (%integer-power a b))
1652                    (%integer-power a b))) ) ) )
1653
1654(define (conj n)
1655  (switchq (%check-number n)
1656    (NONE (bad-number 'conj n))
1657    (COMP (make-complex (complex-real n) (%- 0 (complex-imag n))))
1658    (else n) ) )
1659
1660(define (add1 n) (%+ n 1))
1661(define (sub1 n) (%- n 1))
1662
1663(define (signum n)
1664  (switchq (%check-number n)
1665    (FIX (cond ((eq? 0 n) 0)
1666               ((fx< n 0) -1)
1667               (else 1) ) )
1668    (FLO (cond ((fp= n 0.0) 0.0)
1669               ((fp< n 0.0) -1.0)
1670               (else 1.0) ) )
1671    (BIG (if (%big-negative? n) -1 1)) ; Can't be 0; it would be a fixnum then
1672    (RAT (signum (rat-numerator n)))
1673    (COMP (make-polar 1 (angle n)))     ; Definition from CLHS signum
1674    (else (bad-number 'signum n)) ) )
1675
1676(define (%->integer loc n)
1677  (switchq (%check-number n)
1678    (FIX n)
1679    (FLO (if (%integer? n)
1680             (%flo->integer n)
1681             (bad-integer loc n)))
1682    (BIG n)
1683    (else (bad-integer loc n)) ) )
1684
1685;; From SRFI-60 (oddly not in Chicken core)
1686(define (integer-length x)
1687  (%int-length (%->integer 'integer-length x)))
1688
1689(define (bitwise-and . xs)
1690  (let loop ((x -1) (xs xs))
1691    (if (null? xs)
1692        x
1693        (let ((xi (##sys#slot xs 0)))
1694          (loop
1695           (%int-bitwise-int 0 x (%->integer 'bitwise-and xi))
1696           (##sys#slot xs 1) ) ) ) ) )
1697
1698(define (bitwise-ior . xs)
1699  (let loop ((x 0) (xs xs))
1700    (if (null? xs)
1701        x
1702        (let ((xi (##sys#slot xs 0)))
1703          (loop
1704           (%int-bitwise-int 1 x (%->integer 'bitwise-ior xi))
1705           (##sys#slot xs 1) ) ) ) ) )
1706
1707(define (bitwise-xor . xs)
1708  (let loop ((x 0) (xs xs))
1709    (if (null? xs)
1710        x
1711        (let ((xi (##sys#slot xs 0)))
1712          (loop
1713           (%int-bitwise-int 2 x (%->integer 'bitwise-xor xi))
1714           (##sys#slot xs 1) ) ) ) ) )
1715
1716(define (bitwise-not n)
1717  (- -1 (%->integer 'bitwise-not n)) )
1718
1719(define (arithmetic-shift n m)
1720  (let ((n (%->integer 'arithmetic-shift n)))
1721    (switchq (%check-number m)
1722      (FIX (%int-shift-fix n m))
1723      (BIG (##sys#signal-hook #:type-error 'arithmetic-shift
1724                              "can not shift by bignum amounts" n m))
1725      (else (bad-exact 'arithmetic-shift m)))) )
1726
1727(define %number->string
1728  (let ((string-append string-append))
1729    (lambda (n #!optional (base 10))
1730      (unless (memq base '(2 8 10 16)) (bad-base 'number->string base))
1731      (let numstr ((n n))
1732        (switchq (%check-number n)
1733          (FIX (number->string-0 n base))
1734          (FLO (cond
1735                ((fp= n -inf.0) "-inf.0") ; Core does not handle these right
1736                ((fp= n +inf.0) "+inf.0")
1737                ((not (fp= n n)) "+nan.0")
1738                (else (number->string-0 n base))))
1739          (BIG (%big->string n base))
1740          (RAT (string-append (numstr (rat-numerator n))
1741                              "/"
1742                              (numstr (rat-denominator n))))
1743          (COMP (let ((r (complex-real n))
1744                      (i (complex-imag n)) )
1745                  (string-append
1746                   (numstr r)
1747                   ;; The infinities and NaN always print their sign
1748                   (if (and (finite? i) (%> i 0 'number->string)) "+" "")
1749                   (numstr i) "i") ) )
1750          (else (bad-number 'number->string n)) ) ) ) ) )
1751
1752(define number->string %number->string)
1753(define ##sys#number->string %number->string) ; for printer
1754
1755;; We try to prevent memory exhaustion attacks by limiting the
1756;; maximum exponent value.
1757;; TODO: Make this a parameter?  Would probably slow things down even more...
1758(define-constant +maximum-allowed-exponent+ 10000)
1759
1760(define (%string->compnum radix str offset exactness)
1761  (define (go-inexact!)
1762    ;; Go inexact unless exact was requested (with #e prefix)
1763    (unless (eq? exactness 'e) (set! exactness 'i)))
1764  (define (safe-exponent value e)
1765    (and e (if (not value)
1766               0
1767               (cond
1768                 ((> e +maximum-allowed-exponent+)
1769                  (and (eq? exactness 'i)
1770                       (cond ((zero? value) 0.0)
1771                             ((> value 0.0) +inf.0)
1772                             (else -inf.0))))
1773                 ((< e (- +maximum-allowed-exponent+))
1774                  (and (eq? exactness 'i) +0.0))
1775                 (else (%* value (expt 10 e)))))))
1776  (let* ((len (##sys#size str))
1777         (r..9 (integer->char (fx+ (char->integer #\0) radix)))
1778         (r..a (integer->char (fx+ (char->integer #\a) (fx- radix 10))))
1779         (r..A (integer->char (fx+ (char->integer #\A) (fx- radix 10))))
1780         ;; Ugly flag which we need (note that "exactness" is mutated too!)
1781         ;; Since there is (almost) no backtracking we can do this.
1782         (seen-hashes? #f)
1783         ;; All these procedures return #f or an object consed onto an end
1784         ;; position.  If the cdr is false, that's the end of the string.
1785         ;; If just #f is returned, the string contains invalid number syntax.
1786         (scan-digits
1787          (lambda (start)
1788            (let lp ((i start))
1789              (if (fx= i len)
1790                  (and (fx> i start) (cons i #f))
1791                  (let ((c (%subchar str i)))
1792                    (if (fx<= radix 10)
1793                        (if (and (char>=? c #\0) (char<=? c r..9))
1794                            (lp (fx+ i 1))
1795                            (and (fx> i start) (cons i i)))
1796                        (if (or (and (char>=? c #\0) (char<=? c #\9))
1797                                (and (char>=? c #\a) (char<=? c r..a))
1798                                (and (char>=? c #\A) (char<=? c r..A)))
1799                            (lp (fx+ i 1))
1800                            (and (fx> i start) (cons i i)))))))))
1801         (scan-hashes
1802          (lambda (start)
1803            (let lp ((i start))
1804              (if (fx= i len)
1805                  (and (fx> i start) (cons i #f))
1806                  (let ((c (%subchar str i)))
1807                    (if (eq? c #\#)
1808                        (lp (fx+ i 1))
1809                        (and (fx> i start) (cons i i))))))))
1810         (scan-digits+hashes
1811          (lambda (start neg? all-hashes-ok?)
1812            (let* ((digits (and (not seen-hashes?) (scan-digits start)))
1813                   (hashes (if digits
1814                               (and (cdr digits) (scan-hashes (cdr digits)))
1815                               (and all-hashes-ok? (scan-hashes start))))
1816                   (end (or hashes digits)))
1817              (and-let* ((end)
1818                         (num (%digits->integer str start (car end) radix neg?)))
1819                (when hashes            ; Eeewww. Feeling dirty yet?
1820                  (set! seen-hashes? #t)
1821                  (go-inexact!))
1822                (cons num (cdr end))))))
1823         (scan-exponent
1824          (lambda (start)
1825            (and (fx< start len)
1826                 (let ((sign (case (%subchar str start)
1827                               ((#\+) 'pos) ((#\-) 'neg) (else #f))))
1828                   (and-let* ((start (if sign (fx+ start 1) start))
1829                              (end (scan-digits start)))
1830                     (go-inexact!)
1831                     (cons (%digits->integer
1832                            str start (car end) radix (eq? sign 'neg))
1833                           (cdr end)))))))
1834         (scan-decimal-tail
1835          (lambda (start neg? decimal-head)
1836            (and (fx< start len)
1837                 (let* ((tail (scan-digits+hashes start neg? decimal-head))
1838                        (next (if tail (cdr tail) start)))
1839                   (and (or decimal-head (not next)
1840                            (fx> next start)) ; Don't allow empty "."
1841                        (case (and next (%subchar str next))
1842                          ((#\e #\s #\f #\d #\l
1843                            #\E #\S #\F #\D #\L)
1844                           (and-let* (((fx> len next))
1845                                      (ee (scan-exponent (fx+ next 1)))
1846                                      (e (car ee))
1847                                      (h (safe-exponent decimal-head e)))
1848                             (let* ((te (and tail (fx- e (fx- (cdr tail) start))))
1849                                    (num (and tail (car tail)))
1850                                    (t (safe-exponent num te)))
1851                               (cons (if t (%+ h t) h) (cdr ee)))))
1852                          (else (let* ((last (or next len))
1853                                       (te (and tail (fx- start last)))
1854                                       (num (and tail (car tail)))
1855                                       (t (safe-exponent num te))
1856                                       (h (or decimal-head 0)))
1857                                  (cons (if t (%+ h t) h) next)))))))))
1858         (scan-ureal
1859          (lambda (start neg?)
1860            (if (and (fx> len (fx+ start 1)) (eq? radix 10)
1861                     (eq? (%subchar str start) #\.))
1862                (begin
1863                  (go-inexact!)
1864                  (scan-decimal-tail (fx+ start 1) neg? #f))
1865                (and-let* ((end (scan-digits+hashes start neg? #f)))
1866                  (case (and (cdr end) (%subchar str (cdr end)))
1867                    ((#\.)
1868                     (go-inexact!)
1869                     (and (eq? radix 10)
1870                          (if (fx> len (fx+ (cdr end) 1))
1871                              (scan-decimal-tail (fx+ (cdr end) 1) neg? (car end))
1872                              (cons (car end) #f))))
1873                    ((#\e #\s #\f #\d #\l
1874                      #\E #\S #\F #\D #\L)
1875                     (and-let* (((eq? radix 10))
1876                                ((fx> len (cdr end)))
1877                                (ee (scan-exponent (fx+ (cdr end) 1)))
1878                                (num (car end))
1879                                (val (safe-exponent num (car ee))))
1880                       (cons val (cdr ee))))
1881                    ((#\/)
1882                     (set! seen-hashes? #f) ; Reset flag for denominator
1883                     (and-let* (((fx> len (cdr end)))
1884                                (d (scan-digits+hashes (fx+ (cdr end) 1) #f #f))
1885                                (num (car end))
1886                                (denom (car d)))
1887                       (if (not (eq? denom 0))
1888                           (cons (%/ num denom) (cdr d))
1889                           ;; Hacky: keep around an inexact until we decide we
1890                           ;; *really* need exact values, then fail at the end.
1891                           (and (not (eq? exactness 'e))
1892                            (case (signum num)
1893                              ((-1) (cons -inf.0 (cdr d)))
1894                              ((0)  (cons +nan.0 (cdr d)))
1895                              ((+1) (cons +inf.0 (cdr d))))))))
1896                    (else end))))))
1897         (scan-real
1898          (lambda (start)
1899            (and (fx< start len)
1900                 (let* ((sign (case (%subchar str start)
1901                                ((#\+) 'pos) ((#\-) 'neg) (else #f)))
1902                        (next (if sign (fx+ start 1) start)))
1903                   (and (fx< next len)
1904                        (case (%subchar str next)
1905                          ((#\i #\I)
1906                           (or (and sign
1907                                    (cond
1908                                     ((fx= (fx+ next 1) len) ; [+-]i
1909                                      (cons (if (eq? sign 'neg) -1 1) next))
1910                                     ((and (fx<= (fx+ next 5) len)
1911                                           (string-ci=? (substring str next (fx+ next 5)) "inf.0"))
1912                                      (go-inexact!)
1913                                      (cons (fp/ (if (eq? sign 'neg) -1.0 1.0) 0.0)
1914                                            (and (fx< (fx+ next 5) len)
1915                                                 (fx+ next 5))))
1916                                     (else #f)))
1917                               (scan-ureal next (eq? sign 'neg))))
1918                          ((#\n #\N)
1919                           (or (and sign
1920                                    (fx<= (fx+ next 5) len)
1921                                    (string-ci=? (substring str next (fx+ next 5)) "nan.0")
1922                                    (begin (go-inexact!)
1923                                           (cons (fp/ 0.0 0.0)
1924                                                 (and (fx< (fx+ next 5) len)
1925                                                      (fx+ next 5)))))
1926                               (scan-ureal next (eq? sign 'neg))))
1927                          (else (scan-ureal next (eq? sign 'neg)))))))))
1928         (number (and-let* ((r1 (scan-real offset)))
1929                   (case (and (cdr r1) (%subchar str (cdr r1)))
1930                     ((#f) (car r1))
1931                     ((#\i #\I) (and (fx= len (fx+ (cdr r1) 1))
1932                                     (or (eq? (%subchar str offset) #\+) ; ugh
1933                                         (eq? (%subchar str offset) #\-))
1934                                     (make-rectangular 0 (car r1))))
1935                     ((#\+ #\-)
1936                      (set! seen-hashes? #f) ; Reset flag for imaginary part
1937                      (and-let* ((r2 (scan-real (cdr r1)))
1938                                 ((cdr r2))
1939                                 ((fx= len (fx+ (cdr r2) 1)))
1940                                 ((or (eq? (%subchar str (cdr r2)) #\i)
1941                                      (eq? (%subchar str (cdr r2)) #\I))))
1942                        (make-rectangular (car r1) (car r2))))
1943                     ((#\@)
1944                      (set! seen-hashes? #f) ; Reset flag for angle
1945                      (and-let* ((r2 (scan-real (fx+ (cdr r1) 1)))
1946                                 ((not (cdr r2))))
1947                        (make-polar (car r1) (car r2))))
1948                     (else #f)))))
1949    (and number (if (eq? exactness 'i)
1950                    (exact->inexact number)
1951                    ;; Ensure we didn't encounter +inf.0 or +nan.0 with #e
1952                    (and (finite? number) number)))))
1953
1954(define (%string->number str #!optional (base 10))
1955  (##sys#check-string str 'string->number)
1956  (##sys#check-exact base 'string->number)
1957  (unless (< 1 base 37)           ; We only have 0-9 and the alphabet!
1958    (bad-base 'string->number base))
1959  (let scan-prefix ((i 0)
1960                    (exness #f)
1961                    (radix #f)
1962                    (len (##sys#size str)))
1963    (if (and (fx< (fx+ i 2) len) (eq? (%subchar str i) #\#))
1964        (case (%subchar str (fx+ i 1))
1965          ((#\i #\I) (and (not exness) (scan-prefix (fx+ i 2) 'i radix len)))
1966          ((#\e #\E) (and (not exness) (scan-prefix (fx+ i 2) 'e radix len)))
1967          ((#\b #\B) (and (not radix) (scan-prefix (fx+ i 2) exness 2 len)))
1968          ((#\o #\O) (and (not radix) (scan-prefix (fx+ i 2) exness 8 len)))
1969          ((#\d #\D) (and (not radix) (scan-prefix (fx+ i 2) exness 10 len)))
1970          ((#\x #\X) (and (not radix) (scan-prefix (fx+ i 2) exness 16 len)))
1971          (else #f))
1972        (%string->compnum (or radix base) str i exness))))
1973
1974(define (randomize #!optional (seed (##sys#fudge 2)))
1975  (switchq (%check-number seed)
1976    (FIX (%fix-randomize seed))
1977    (BIG (%big-randomize seed))
1978    (else (bad-integer 'randomize seed)) ) )
1979
1980(define (random n)
1981  (switchq (%check-number n)
1982    (FIX (%fix-random n))
1983    (BIG (%big-random n))
1984    (else (bad-integer 'random n)) ) )
1985
1986(define string->number %string->number)
1987
1988;;; Reader hook
1989(define (##sys#string->number str #!optional (radix 10) exactness)
1990  (%string->compnum radix str 0 exactness))
1991
1992
1993;;; Non-standard type procedures
1994
1995(define (bignum? x) (eq? (%check-number x) BIG)) ; big number
1996(define (ratnum? x) (eq? (%check-number x) RAT)) ; rational number
1997(define (cplxnum? x) (eq? (%check-number x) COMP)) ; complex number
1998
1999(define (nan? x)
2000  (switchq (%check-number x)
2001    (NONE (bad-number 'nan? x))
2002    (FLO (not (fp= x x)))
2003    (COMP (or (nan? (complex-real x)) (nan? (complex-imag x))))
2004    (else #f)))
2005(define (infinite? x)
2006  (switchq (%check-number x)
2007    (NONE (bad-number 'infinite? x))
2008    (FLO (or (fp= x +inf.0) (fp= x -inf.0)))
2009    (COMP (or (infinite? (complex-real x)) (infinite? (complex-imag x))))
2010    (else #f)))
2011(define (finite? x)
2012  (switchq (%check-number x)
2013    (NONE (bad-number 'finite? x))
2014    (FLO (and (fp= x x) (not (fp= x +inf.0)) (not (fp= x -inf.0))))
2015    (COMP (and (finite? (complex-real x)) (finite? (complex-imag x))))
2016    (else #t)))
2017
2018(define (rectnum? x)    ; "exact" complex number
2019  (and (eq? (%check-number x) COMP)
2020       (%integer? (complex-real x))
2021       (%integer? (complex-imag x))))
2022
2023(define (compnum? x)    ; inexact complex number
2024  (and (eq? (%check-number x) COMP)
2025       (%inexact? (complex-real x))
2026       (%inexact? (complex-imag x))))
2027
2028(define (cintnum? x)    ; integer number
2029  (switchq (%check-number x)
2030    (FIX #t)
2031    (BIG #t)
2032    (FLO (%flo-integer? x))
2033    (COMP (and (%integer? (complex-real x)) (%integer? (complex-imag x))))
2034    (else #f) ) )
2035
2036(define (cflonum? x)    ; floatingpoint number
2037  (switchq (%check-number x)
2038    (FLO #t)
2039    (COMP (and (%flonum? (complex-real x)) (%flonum? (complex-imag x))))
2040    (else #f) ) )
2041
2042;;; What we provide
2043
2044(register-feature! #:full-numeric-tower)
2045
2046)
Note: See TracBrowser for help on using the repository browser.