Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: