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

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

numbers: Convert "remainder" to new style

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