1;;; 5_3.ss
2;;; Copyright 1984-2017 Cisco Systems, Inc.
3;;;
4;;; Licensed under the Apache License, Version 2.0 (the "License");
5;;; you may not use this file except in compliance with the License.
6;;; You may obtain a copy of the License at
7;;;
8;;; http://www.apache.org/licenses/LICENSE-2.0
9;;;
10;;; Unless required by applicable law or agreed to in writing, software
11;;; distributed under the License is distributed on an "AS IS" BASIS,
12;;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13;;; See the License for the specific language governing permissions and
14;;; limitations under the License.
15
16;;; Care must be take with floating point constants to permit cross
17;;; compilation between machines with differing floating point styles.
18;;; Negative zero, infinities, large or small numbers, non-binary
19;;; fractions, and precise numbers are dangerous and should be calculated.
20;;; positive zero, NAN, small integers, and binary fractions with only a few
21;;; significant bits are safe on all current machines.
22;;; examples:
23;;; dangerous: -0.0, +inf.0, -inf.0, 1e100, 1e-100, 0.1
24;;; safe: 0.0, +nan.0, 1.0, 2.0, 0.5
25
26(begin
27(eval-when (compile)
28
29   (define-constant max-float-exponent
30      (float-type-case
31         [(ieee) 1023]))
32
33   (define-constant min-float-exponent
34      (float-type-case
35         [(ieee) -1023]))
36
37)
38
39(let ()
40; could use foreign-entry? primitive if foreign.ss were loaded first
41(define op-if-entry?
42   (let ()
43      (define lookup
44         (foreign-procedure "(cs)lookup_foreign_entry" (string)
45            void*))
46      (lambda (op name)
47         (and (not (eqv? (lookup name) 0))
48              (op name)))))
49
50(let ()
51
52(define cflop1
53   (lambda (x)
54      (foreign-procedure x (double-float) double-float)))
55
56(define cflop2
57   (lambda (x)
58      (foreign-procedure x (double-float double-float) double-float)))
59
60(define schemeop1
61   (lambda (x)
62      (foreign-procedure x (scheme-object) scheme-object)))
63
64(define schemeop2
65   (lambda (x)
66      (foreign-procedure x (scheme-object scheme-object) scheme-object)))
67
68(let ()
69
70(define biglength (schemeop1 "(cs)s_integer_length"))
71(define bigodd? (schemeop1 "(cs)s_bigoddp"))
72(define float (schemeop1 "(cs)s_float"))
73
74(define big=
75  (foreign-procedure "(cs)s_big_eq" (scheme-object scheme-object)
76    boolean))
77(define big<
78  (foreign-procedure "(cs)s_big_lt" (scheme-object scheme-object)
79    boolean))
80(define big-negate (schemeop1 "(cs)s_big_negate"))
81(define integer-ash (schemeop2 "(cs)s_ash"))
82(define integer+ (schemeop2 "(cs)add"))
83(define integer* (schemeop2 "(cs)mul"))
84(define integer- (schemeop2 "(cs)sub"))
85(define schoolbook-integer/ (schemeop2 "(cs)s_div"))
86(define schoolbook-intquotient (schemeop2 "(cs)ss_trunc"))
87(define schoolbook-intquotient-remainder (schemeop2 "(cs)ss_trunc_rem"))
88(define schoolbook-intremainder (schemeop2 "(cs)rem"))
89(define make-ratnum (schemeop2 "(cs)s_rational")) ; does not normalize, except detecting 1 as demoninator
90(define exgcd (schemeop2 "(cs)gcd"))
91
92(define $flsin (cflop1 "(cs)sin"))
93
94(define $flcos (cflop1 "(cs)cos"))
95
96(define $flasin (cflop1 "(cs)asin"))
97
98(define $flacos (cflop1 "(cs)acos"))
99(define $flfloor (cflop1 "(cs)floor"))
100(define $flceiling (cflop1 "(cs)ceil"))
101
102(let ()
103
104;; Burnikel-Ziegler division by Peter Bex from a series about CHICKEN's
105;; numeric tower:
106;;   https://www.more-magic.net/posts/numeric-tower-part-3.html
107;;   Licensed under the Creative Commons Attribution 3.0 License.
108;; The Scheme code here appears to be directly based on the C
109;; code in CHICKEN's BSD-licensed "runtime.c":
110;;   Copyright (c) 2008-2020, The CHICKEN Team
111;;   Copyright (c) 2000-2007, Felix L. Winkelmann
112;;   All rights reserved.
113
114(define DIV-LIMIT 100)
115
116(define (bigits->bits n) (fx* (constant bigit-bits) n))    ; Small helper
117
118(define (extract-bigits n start end)
119  (let ([s-bits (bigits->bits start)])
120    (bitwise-bit-field n s-bits (if end
121                                    (bigits->bits end)
122                                    (fxmax s-bits (integer-length n))))))
123
124;; Here and in 2n/1n we pass both b and [b1, b2] to avoid splitting
125;; up the number more than once.  This is a helper function for 2n/n.
126(define (burnikel-ziegler-3n/2n a12 a3 b b1 b2 n)
127  (let-values ([(q^ r1) (if (< (bitwise-arithmetic-shift-right a12 (bigits->bits n)) b1)
128                            (let* ((n/2 (fxsra n 1))                     ; (floor (/ n 2))
129                                   (b11 (extract-bigits b1 n/2 #f))      ; b1[n..n/2]
130                                   (b12 (extract-bigits b1 0 n/2)))      ; b1[n/2..0]
131                              (burnikel-ziegler-2n/1n a12 b1 b11 b12 n #t))
132                            ;; Don't bother dividing if a1 is a larger number than b1.
133	                    ;; We use a maximum guess instead (see paper for proof).
134                            (let ((base*n (bigits->bits n)))
135                              (values (- (bitwise-arithmetic-shift-left 1 base*n) 1)  ; B^n-1
136                                      (+ (- a12 (bitwise-arithmetic-shift-left b1 base*n)) b1))))])
137    (let ((r1a3 (+ (bitwise-arithmetic-shift-left r1 (bigits->bits n)) a3)))
138      (let lp ((r^ (- r1a3 (* q^ b2)))
139               (q^ q^))
140        (if (negative? r^)
141            (lp (+ r^ b) (- q^ 1))                     ; Adjust!
142            (values q^ r^))))))
143
144;; The main 2n/n algorithm which calls 3n/2n twice.  Here, a is the
145;; numerator, b the denominator, n is the length of b (in digits) and
146;; b1 and b2 are the two halves of b (these never change).
147(define (burnikel-ziegler-2n/1n a b b1 b2 n return-quot?)
148  (if (or (fxodd? n) (fx< n DIV-LIMIT))             ; Can't recur?
149      (let ([p (schoolbook-intquotient-remainder a b)]) ; Use school division
150        (values (car p) (cdr p)))
151      (let* ((n/2 (fxsra n 1))
152             ;; Split a and b into n-sized parts [a1, ..., a4] and [b1, b2]
153             (a12 (extract-bigits a n #f))             ; a[m..n]
154             (a3  (extract-bigits a n/2 n))            ; a[n..n/2]
155             (a4  (extract-bigits a 0 n/2)))           ; a[n..0]
156        ;; Calculate high quotient and intermediate remainder (first half)
157        (let-values ([(q1 r1) (burnikel-ziegler-3n/2n a12 a3 b b1 b2 n/2)])
158          ;; Calculate low quotient and final remainder (second half)
159          (let-values ([(q2 r) (burnikel-ziegler-3n/2n r1 a4 b b1 b2 n/2)])
160            ;; Recombine quotient parts as q = [q1,q2]
161            (values (and return-quot?
162                         (+ (bitwise-arithmetic-shift-left q1 (bigits->bits n/2)) q2))
163                    r))))))
164
165(define (quotient&remainder/burnikel-ziegler x y return-quot? return-rem?)
166  ;; Caller will have made sure that x and y are bignums
167  (let* ((q-neg? (if (negative? y) (not (negative? x)) (negative? x)))
168         (r-neg? (negative? x))
169         (abs-x (abs x))
170         (abs-y (abs y)))
171    (cond
172      [(> abs-x abs-y)
173       (let* ((x abs-x)
174              (y abs-y)
175              (s ($bignum-length y))
176              ;; Define m as min{2^k|(2^k)*DIV_LIMIT > s}.
177              ;; This ensures we shift as little as possible (less pressure
178              ;; on the GC) while maintaining a power of two until we drop
179              ;; below the threshold, so we can always split N in half.
180              (m (fxsll 1 (integer-length (fx/ s DIV-LIMIT))))
181              (j (fx/ (fx+ s (fx- m 1)) m))  ; j = s/m, rounded up
182              (n (fx* j m))
183              ;; Normalisation, just like with normal school division
184              (norm-shift (fx- (bigits->bits n) (integer-length y)))
185              (x (bitwise-arithmetic-shift-left x norm-shift))
186              (y (bitwise-arithmetic-shift-left y norm-shift))
187              ;; l needs to be the smallest value so that a < base^{l*n}/2
188              (l (fx/ (fx+ ($bignum-length x) n) n))
189              (l (if (fx= (bigits->bits l) (integer-length x)) (fx+ l 1) l))
190              (t (fxmax l 2))
191              (y-hi (extract-bigits y (fxsra n 1) #f))   ; y[n..n/2]
192              (y-lo (extract-bigits y 0 (fxsra n 1))))   ; y[n/2..0]
193         (let lp ((zi (bitwise-arithmetic-shift-right x (bigits->bits (fx* (fx- t 2) n))))
194                  (i (fx- t 2))
195                  (quot 0))
196           (let-values ([(qi ri) (burnikel-ziegler-2n/1n zi y y-hi y-lo n return-quot?)])
197             (let ((quot (and return-quot?
198                              (+ (bitwise-arithmetic-shift-left quot (bigits->bits n)) qi))))
199               (if (fx> i 0)
200                   (let ((zi-1 (let* ((base*n*i-1 (fx* n (fx- i 1)))
201                                      (base*n*i   (fx* n i))
202                                      ;; x_{i-1} = x[n*i..n*(i-1)]
203                                      (xi-1 (extract-bigits x base*n*i-1 base*n*i)))
204                                 (+ (bitwise-arithmetic-shift-left ri (bigits->bits n)) xi-1))))
205                     (lp zi-1 (fx- i 1) quot))
206                   ;; De-normalise remainder, but only if necessary
207                   (let ((rem (if (or (not return-rem?) (eq? 0 norm-shift))
208                                  ri
209                                  (bitwise-arithmetic-shift-right ri norm-shift))))
210                     ;; Return values (quot, rem or both) with correct sign:
211                     (cond ((and return-quot? return-rem?)
212                            (values (if q-neg? (- quot) quot)
213                                    (if r-neg? (- rem) rem)))
214                           (return-quot? (if q-neg? (- quot) quot))
215                           (else (if r-neg? (- rem) rem)))))))))]
216      [(< abs-x abs-y)
217       (cond
218         [(and return-quot? return-rem?) (values 0 x)]
219         [return-quot? 0]
220         [else x])]
221      [else
222       (cond
223         [(and return-quot? return-rem?) (values (if q-neg? -1 1) 0)]
224         [return-quot? (if q-neg? -1 1)]
225         [else 0])])))
226
227;; Only try to use Burnikel-Ziegler when we have large enough bignums:
228(define (big-divide-bignums? n d)
229  (and (bignum? n)
230       (bignum? d)
231       (fx>= ($bignum-length n) DIV-LIMIT)
232       (fx>= ($bignum-length d) DIV-LIMIT)))
233
234(define integer/
235  (lambda (n d)
236    (cond
237      [(big-divide-bignums? n d)
238       (let* ([g (exgcd n d)]
239              [g (if ($bigpositive? d)
240                     g
241                     (- g))])
242         (if (or (fixnum? g)
243                 (fx< ($bignum-length g) DIV-LIMIT))
244             (make-ratnum (schoolbook-intquotient n g)
245                          (schoolbook-intquotient d g))
246             (make-ratnum (quotient&remainder/burnikel-ziegler n g #t #f)
247                          (quotient&remainder/burnikel-ziegler d g #t #f))))]
248      [else (schoolbook-integer/ n d)])))
249
250(define intquotient
251  (lambda (n d)
252    (cond
253      [(big-divide-bignums? n d)
254       (quotient&remainder/burnikel-ziegler n d #t #f)]
255      [else
256       (schoolbook-intquotient n d)])))
257
258(define intquotient-remainder
259  (lambda (n d)
260    (cond
261      [(big-divide-bignums? n d)
262       (let-values ([(q r) (quotient&remainder/burnikel-ziegler n d #t #t)])
263         (cons q r))]
264      [else
265       (schoolbook-intquotient-remainder n d)])))
266
267(define intremainder
268  (lambda (n d)
269    (cond
270      [(big-divide-bignums? n d)
271       (quotient&remainder/burnikel-ziegler n d #f #t)]
272      [else
273       (schoolbook-intremainder n d)])))
274
275(let ()
276
277(define omega
278   (float-type-case
279      [(ieee) (float #e1.7976931348623157e308)]))
280
281(define $flexpt
282   (machine-case
283      [(i3nt ti3nt a6s2 ta6s2 i3s2 ti3s2 i3nb ti3nb a6nb ta6nb)
284       ; pow(nan,+0.0) => nan instead of +1.0
285       (let ([cexpt (cflop2 "(cs)pow")])
286         (lambda (x y)
287           (cond
288            [(fl= y 0.0) 1.0]
289            [else (cexpt x y)])))]
290      [else (cflop2 "(cs)pow")]))
291
292(define $fltan (cflop1 "(cs)tan"))
293
294(define flcosh (cflop1 "(cs)cosh"))
295
296(define fltanh
297   (machine-case
298      [(i3fb ti3fb)
299       ; broken for -0.0, +/-inf
300       (let ([ctanh (cflop1 "(cs)tanh")])
301        (lambda (x)
302          (cond
303           [(fl= x 0.0) x]
304           [(infinity? x) (if (negated-flonum? x) -1.0 1.0)]
305           [else (ctanh x)])))]
306      [(i3nb ti3nb a6nb ta6nb)
307       ; broken for -0.0
308       (let ([ctanh (cflop1 "(cs)tanh")])
309        (lambda (x)
310          (cond
311           [(fl= x 0.0) x]
312           [else (ctanh x)])))]
313      [else (cflop1 "(cs)tanh")]))
314
315(define $flexp (cflop1 "(cs)exp"))
316
317(define $fllog
318  (machine-case
319    [(a6s2 ta6s2 i3s2 ti3s2 i3ob ti3ob i3nb ti3nb a6nb ta6nb a6ob ta6ob)
320     ; broken for -inf.0
321     (let ([clog (cflop1 "(cs)log")])
322       (lambda (x) (if (and (infinity? x) (negated-flonum? x)) +nan.0 (clog x))))]
323    [else (cflop1 "(cs)log")]))
324
325(define $flsqrt (cflop1 "(cs)sqrt"))
326
327(define flatan2
328   (machine-case
329      [(i3nt ti3nt)
330       ; atan2(+inf.0,+inf.0) => pi/2 instead of pi/4
331       ; atan2(-inf.0,-inf.0) => -pi/2 instead of -3pi/4
332       ; atan2(+inf.0,-inf.0) => NAN instead of 3pi/4
333       ; atan2(-inf.0,+inf.0) => NAN instead of -pi/4
334       ; atan2(+0.0,-0.0) => +0.0 instead of +pi
335       ; atan2(-0.0,-0.0) => -0.0 instead of -pi
336       ; atan2(-0.0,-1.0) => pi instead of -pi
337       (let ([catan2 (cflop2 "(cs)atan2")])
338          (let ([pi (catan2 0.0 -1.0)])
339             (lambda (y x)
340                (cond
341                   [(and (infinity? y) (infinity? x))
342                    (let ([y (if (negated-flonum? y) -1.0 1.0)]
343                          [x (if (negated-flonum? x) -1.0 1.0)])
344                       (catan2 y x))]
345                   [(and (fl= y 0.0) (not ($nan? x)))
346                    (if (negated-flonum? y)
347                        (if (negated-flonum? x) (fl- pi) (fl- 0.0))
348                        (if (negated-flonum? x) pi 0.0))]
349                   [else (catan2 y x)]))))]
350      [(i3ob ti3ob a6ob ta6ob a6s2 ta6s2 i3s2 ti3s2 i3nb ti3nb a6nb ta6nb)
351       ; atan2(-0.0,+0.0) => +0.0 instead of -0.0
352       ; atan2(+0.0,-0.0) => +0.0 instead of +pi
353       ; atan2(-0.0,-0.0) => +0.0 instead of -pi
354       (let ([catan2 (cflop2 "(cs)atan2")])
355         (let ([pi (catan2 0.0 -1.0)])
356           (lambda (y x)
357             (cond
358              [(and (fl= y 0.0) (not ($nan? x)))
359               (if (negated-flonum? y)
360                   (if (negated-flonum? x) (fl- pi) (fl- 0.0))
361                   (if (negated-flonum? x) pi 0.0))]
362              [else (catan2 y x)]))))]
363       [else (cflop2 "(cs)atan2")]))
364
365(define $flatan (cflop1 "(cs)atan"))
366
367(define flsinh (cflop1 "(cs)sinh"))
368
369(define flatanh
370   (or (op-if-entry? cflop1 "(cs)atanh")
371       ; |x| <= 1
372       ; principle expression:
373       ; (log(1+x)-log(1-x))/2
374       ; should use "log1p" but it doesn't exist on the 88k
375       (let ([f (lambda (x)
376                   (fl* 0.5 (fl- ($fllog (fl+ 1.0 x)) ($fllog (fl- 1.0 x)))))])
377          (lambda (x)
378             (if (negated-flonum? x) (fl- (f (fl- x))) (f x))))))
379
380(define fllog1+
381   (or (op-if-entry? cflop1 "(cs)log1p")
382       (lambda (x) ($fllog (fl+ 1.0 x)))))
383
384(let ()
385
386(define log2 ($fllog 2.0))
387
388(define flhypot (cflop2 "(cs)hypot"))
389
390(define flasinh
391   ; scheme-coded version needs "log2"
392   (or (op-if-entry? cflop1 "(cs)asinh")
393       ; prinicple expression:
394       ; log(x + sqrt(xx + 1))
395       ; avoids spurious overflows
396       ; avoids underflow problems from negative x by using identity
397       ; asinh(-x) = -asinh(x)
398       ; should use "log1p" for small x but it doesn't exist on the 88k
399       (let ([f (lambda (x)
400                   (if (fl= (fl+ x 1.0) x)
401                       (fl+ ($fllog x) log2)
402                       ($fllog (fl+ x ($flsqrt (fl+ (fl* x x) 1.0))))))])
403          (lambda (x)
404             (if (negated-flonum? x) (fl- (f (fl- x))) (f x))))))
405
406(define flacosh
407   ; scheme-coded version needs "log2"
408   (or (op-if-entry? cflop1 "(cs)acosh")
409       ; x >= 1
410       ; prinicple expression:
411       ; log(x + sqrt(xx - 1))
412       ; avoids spurious overflows
413       (lambda (x)
414          (if (fl= (fl- x 1.0) x)
415              (fl+ ($fllog x) log2)
416              ($fllog (fl+ x ($flsqrt (fl- (fl* x x) 1.0))))))))
417
418(let ()
419
420(define pi (flatan2 0.0 -1.0))
421(define sqrt-omega ($flsqrt omega))
422(define log-omega ($fllog omega))
423(define acosh-omega (flacosh omega))
424
425(let ()
426
427(define-syntax define-trig-op
428  (syntax-rules ()
429    [(_ who flop cflop zero-value)
430     (set! who
431       (lambda (x)
432         (type-case x
433           [(flonum?) (flop x)]
434           [($inexactnum?) (cflop x)]
435           [(fixnum?) (if (fx= x 0) zero-value (who (fixnum->flonum x)))]
436           [(bignum? ratnum? $exactnum?) (who (inexact x))]
437           [else (nonnumber-error 'who x)])))]))
438
439(define $flinteger-or-inf?
440  (lambda (x)
441    (fl= ($flfloor x) x)))
442
443(define $flinteger?
444  (lambda (x)
445    (and ($flinteger-or-inf? x)
446         (not (exceptional-flonum? x)))))
447
448(define nonnumber-error
449   (lambda (who what)
450      ($oops who "~s is not a number" what)))
451
452(define noncomplex-error
453   (lambda (who what)
454      ($oops who "~s is not a complex number" what)))
455
456(define nonreal-error
457   (lambda (who what)
458      ($oops who "~s is not a real number" what)))
459
460(define nonrational-error
461   (lambda (who what)
462      ($oops who "~s is not a rational number" what)))
463
464(define noninteger-error
465   (lambda (who what)
466      ($oops who "~s is not an integer" what)))
467
468(define nonexact-integer-error
469   (lambda (who what)
470      ($oops who "~s is not an exact integer" what)))
471
472(define noncflonum-error
473   (lambda (who what)
474      ($oops who "~s is not a cflonum" what)))
475
476(define domain-error
477   (lambda (who what)
478      ($oops who "undefined for ~s" what)))
479
480(define domain-error2
481   (lambda (who x y)
482      ($oops who "undefined for values ~s and ~s" x y)))
483
484; note: (cfl*i z) =/= (* +i z) if RP(z) == -0.0
485(define cfl*i
486   (lambda (z)
487      (fl-make-rectangular (fl- (cfl-imag-part z)) (cfl-real-part z))))
488
489; note: (cfl/i z) =/= (/ z +i) or (* -i z) if IP(z) == -0.0
490(define cfl/i
491   (lambda (z)
492      (fl-make-rectangular (cfl-imag-part z) (fl- (cfl-real-part z)))))
493
494; Some of the following is based on
495; W. Kahan's "Branch Cuts for Complex Elementary Functions"
496; in "The State of the Art of Numerical Analysis"
497; (IMA/SIAM proceedings, 1986, pp 165-211)
498; ed. by A. Iserles and M.J.D. Powell
499
500; Kahan gives principal expressions and algorithms for several
501; complex functions.  The principal expressions are mathematically
502; correct, but not necessarily good computationally.  They
503; do, however, make good test expressions for ordinary inputs.
504
505; Steele's "Common Lisp: the Language" (second edition) was used
506; to determine valid domains for some of the functions.
507
508(define cflmagnitude
509   (lambda (z)
510      (flhypot (cfl-real-part z) (cfl-imag-part z))))
511
512(define cflangle
513   (lambda (z)
514      (flatan2 (cfl-imag-part z) (cfl-real-part z))))
515
516(define cfllog
517   ; principal expression from Kahan:
518   ; log(z) = log(|z|) + angle(z)i
519   ; Kahan uses a different algorithm to calculate the real part.
520   (let ([f (lambda (x y)
521               ; x >= y
522               (let ([r (fl/ y x)])
523                  (fl+ ($fllog x) (fl* .5 (fllog1+ (fl* r r))))))]
524         [k (fl* .5 log2)])
525      (lambda (z)
526         (let ([x (cfl-real-part z)] [y (cfl-imag-part z)])
527            (fl-make-rectangular
528               (let ([x (flabs x)] [y (flabs y)])
529                  (cond
530                     [(fl> x y) (f x y)]
531                     [(fl< x y) (f y x)]
532                     [(fl= x y) (fl+ ($fllog x) k)]
533                     [(infinity? x) x]
534                     [(infinity? y) y]
535                     [($nan? x) x]
536                     [else y]))
537               (flatan2 y x))))))
538
539(define cflsqrt
540   ; principal expression from Kahan:
541   ; sqrt(z) = expt(z,1/2)
542   ; Kahan's algorithm except for the calculation of "a"
543   (let ([f (let ([k ($flsqrt (fl* .5 (fl+ ($flsqrt 2.0) 1.0)))])
544               (lambda (x y)
545                  ; sqrt(|x+yi| + |x|)/2
546                  (cond
547                     [(fl> x y)
548                      (let ([r (fl/ y x)])
549                         (fl* ($flsqrt x)
550                              ($flsqrt (fl* .5 (fl+ ($flsqrt (fl+ 1.0 (fl* r r)))
551                                                   1.0)))))]
552                     [(fl< x y)
553                      (let ([r (fl/ x y)])
554                         (fl* ($flsqrt y)
555                              ($flsqrt (fl* .5 (fl+ ($flsqrt (fl+ (fl* r r) 1.0))
556                                                   r)))))]
557                     [(fl= x y) (fl* ($flsqrt x) k)]
558                     [(infinity? x) x]
559                     [(infinity? y) y]
560                     [($nan? x) x]
561                     [else y])))])
562      (lambda (z)
563         (let ([x (cfl-real-part z)] [y (cfl-imag-part z)])
564            (let ([a (f (flabs x) (flabs y))])
565               (if (fl= a 0.0)
566                   (fl-make-rectangular a y)
567                   (let ([b (if (infinity? y) y (fl* (fl/ y a) .5))])
568                      (if (fl< x 0.0)
569                          (fl-make-rectangular
570                             (flabs b)
571                             (if (negated-flonum? y) (fl- a) a))
572                          (fl-make-rectangular a b)))))))))
573
574(define cflexp
575   ; exp(a+bi) = exp(a)cos(b) + exp(a)sin(b)i
576   (lambda (z)
577      (let ([a (cfl-real-part z)] [b (cfl-imag-part z)])
578         (cond
579           ; perhaps misguidedly treat x+0.0i the same as x
580            [(fl= b 0.0) (fl-make-rectangular ($flexp a) b)]
581            [(fl<= a log-omega)
582             (let ([e^a ($flexp a)])
583                (fl-make-rectangular (fl* e^a ($flcos b)) (fl* e^a ($flsin b))))]
584            [else (fl-make-rectangular
585                    (let ([cosb ($flcos b)])
586                      (if (fl< cosb 0.0)
587                          (fl- ($flexp (fl+ a ($fllog (fl- cosb)))))
588                          ($flexp (fl+ a ($fllog cosb)))))
589                    (let ([sinb ($flsin b)])
590                      (if (fl< sinb 0.0)
591                          (fl- ($flexp (fl+ a ($fllog (fl- sinb)))))
592                          ($flexp (fl+ a ($fllog sinb))))))]))))
593
594(define cflslowsinh
595   ; probably not the best way to handle this
596   (let ([f (lambda (z -z)
597               (cfl- (cflexp (cfl- z log2)) (cfl* .5 (cflexp -z))))])
598      (lambda (z)
599         (if (fl< (cfl-real-part z) 0.0)
600             (cfl- (f (cfl- z) z))
601             (f z (cfl- z))))))
602
603(define cflslowcosh
604   ; probably not the best way to handle this
605   (let ([f (lambda (z -z)
606               (cfl+ (cflexp (cfl- z log2)) (cfl* .5 (cflexp -z))))])
607      (lambda (z)
608         (if (fl< (cfl-real-part z) 0.0)
609             (f (cfl- z) z)
610             (f z (cfl- z))))))
611
612(define cflsin
613   ; sin(a+bi) = sin(a)cosh(b)+cos(a)sinh(b)i
614   (lambda (z)
615      (let ([a (cfl-real-part z)] [b (cfl-imag-part z)])
616         (if (fl<= (flabs b) acosh-omega)
617             (fl-make-rectangular (fl* ($flsin a) (flcosh b))
618                                  (fl* ($flcos a) (flsinh b)))
619             (cfl/i (cflslowsinh (cfl*i z)))))))
620
621(define cflcos
622   ; cos(a+bi) = cos(a)cosh(b)-sin(a)sinh(b)i
623   (lambda (z)
624      (let ([a (cfl-real-part z)] [b (cfl-imag-part z)])
625         (if (fl<= (flabs b) acosh-omega)
626             (fl-make-rectangular (fl* ($flcos a) (flcosh b))
627                                  (fl- (fl* ($flsin a) (flsinh b))))
628             (cflslowcosh (cfl*i z))))))
629
630(define cfltan
631   ; from Kahan
632   (lambda (z)
633      (cfl/i (cfltanh (cfl*i z)))))
634
635(define cflacos
636   ; from Kahan
637   ; principal expression:
638   ; 2log(sqrt((1+z)/2) + sqrt((1-z)/2)i)/i = pi/2 - asin(z)
639   ; returns a+bi where
640   ; a = 2atan(RP(sqrt(1-z))/RP(sqrt(1+z)))
641   ; b = asinh(IP(conjugate(sqrt(1+z)))sqrt(1-z))
642   (lambda (z)
643      (let ([z- (cflsqrt (cfl- 1.0 z))]
644            [z+ (cflsqrt (cfl+ 1.0 z))])
645         (let ([a (cfl-real-part z-)] [b (cfl-imag-part z-)]
646               [c (cfl-real-part z+)] [d (cfl-imag-part z+)])
647            (fl-make-rectangular (fl* 2.0 (flatan2 a c))
648                                 (flasinh (fl- (fl* b c) (fl* a d))))))))
649
650(define cflasin
651   ; from Kahan
652   ; principal expression:
653   ; asinh(iz)/i
654   ; returns a+bi where
655   ; a = atan(RP(z)/RP(sqrt(1-z)sqrt(1+z)))
656   ; b = asinh(IP(conjugate(sqrt(1-z))sqrt(1+z)))
657   (lambda (z)
658      (let ([z- (cflsqrt (cfl- 1.0 z))]
659            [z+ (cflsqrt (cfl+ 1.0 z))])
660         (let ([a (cfl-real-part z-)] [b (cfl-imag-part z-)]
661               [c (cfl-real-part z+)] [d (cfl-imag-part z+)])
662            (fl-make-rectangular
663               (flatan2 (cfl-real-part z) (if (flonum? z)
664                                              0.0
665                                              (fl- (fl* a c) (fl* b d))))
666               (flasinh (fl- (fl* a d) (fl* b c))))))))
667
668(define cflasinh
669   ; from Kahan
670   ; principal expression:
671   ; log(z + sqrt(1 + zz))
672   (lambda (z)
673      (cfl/i (cflasin (cfl*i z)))))
674
675(define cflsinh
676   ; sinh(a+bi) = sinh(a)cos(b)+cosh(a)sin(b)i
677   (lambda (z)
678      (let ([a (cfl-real-part z)] [b (cfl-imag-part z)])
679         (if (fl<= a acosh-omega)
680             (fl-make-rectangular (fl* (flsinh a) ($flcos b))
681                                  (fl* (flcosh a) ($flsin b)))
682             (cflslowsinh z)))))
683
684(define cflcosh
685   ; cosh(a+bi) = cosh(a)cos(b)+sinh(a)sin(b)i
686   (lambda (z)
687      (let ([a (cfl-real-part z)] [b (cfl-imag-part z)])
688         (if (fl<= a acosh-omega)
689             (fl-make-rectangular (fl* (flcosh a) ($flcos b))
690                                  (fl* (flsinh a) ($flsin b)))
691             (cflslowcosh z)))))
692
693(define cfltanh
694   ; from Kahan
695   (let ([theta (fl/ acosh-omega 4.0)])
696      (lambda (z)
697         (let ([x (cfl-real-part z)] [y (cfl-imag-part z)])
698            (let ([ax (flabs x)])
699               (if (fl> ax theta)
700                   (fl-make-rectangular
701                      (if (negated-flonum? x) -1.0 1.0)
702                      (if (negated-flonum? y) (fl- 0.0) 0.0))
703                   (let ([t ($fltan y)]
704                         [s (flsinh x)])
705                      (let ([beta (fl+ 1.0 (fl* t t))]
706                            [ss (fl* s s)])
707                         (let ([rho ($flsqrt (fl+ 1.0 ss))])
708                            (if (infinity? t)
709                                (fl-make-rectangular (fl/ rho s) (/ t))
710                                (let ([k (/ (fl+ 1.0 (fl* beta ss)))])
711                                   (fl-make-rectangular (fl* beta rho s k)
712                                                        (fl* t k)))))))))))))
713
714(define cflacosh
715   ; from Kahan
716   ; principal expression:
717   ; 2log(sqrt((z+1)/2) + sqrt((z-1)/2))
718   ; returns a+bi where
719   ; a = (asinh (real-part (* (conjugate (sqrt (- z 1))) (sqrt (+ z 1)))))
720   ; b = (* 2 (atan (/ (imag-part (sqrt (- z 1))) (real-part (sqrt (+ z 1))))))
721   (lambda (z)
722      (let ([z- (cflsqrt (cfl- z 1.0))]
723            [z+ (cflsqrt (cfl+ z 1.0))])
724         (let ([a (cfl-real-part z-)] [b (cfl-imag-part z-)]
725               [c (cfl-real-part z+)] [d (cfl-imag-part z+)])
726            (fl-make-rectangular (flasinh (fl+ (fl* a c) (fl* b d)))
727                                 (fl* 2.0 ($flatan (fl/ b c))))))))
728
729(define cflatanh
730   ; principal expression from Kahan:
731   ; (log(1+z) - log(1-z))/2
732   (let ([f (let ([theta (fl/ sqrt-omega 4.0)] [pi/2 (flatan2 1.0 0.0)])
733               (let ([rho (fl/ theta)] [-pi/2 (fl- pi/2)])
734                  (lambda (x y)
735                     ; x is positive
736                     (let ([ay (abs y)])
737                        (cond
738                           [(or (fl> x theta) (fl> ay theta))
739                            ; RP(1/z) +/- (pi/2)i
740                            (fl-make-rectangular
741                               (cond
742                                  [(fl> x ay) (fl/ (fl+ x (fl* (fl/ y x) y)))]
743                                  [(fl< x ay) (let ([r (fl/ y x)])
744                                                 (fl/ r (fl+ (fl* x r) y)))]
745                                  [else (fl/ (fl+ x ay))])
746                               (if (negated-flonum? y) pi/2 -pi/2))]
747                           [(fl= x 1.0)
748                            (let ([k (fl+ ay rho)])
749                               (fl-make-rectangular
750                                  ($fllog (fl/ ($flsqrt ($flsqrt (fl+ 4.0
751                                                                   (* y y))))
752                                                      ($flsqrt k)))
753                                  (fl/ (fl+ pi/2 ($flatan (fl/ k 2.0)))
754                                       (if (negated-flonum? y) 2.0 -2.0))))]
755                           [else
756                            (let ([1-x (fl- 1.0 x)]
757                                  [k (let ([k (fl+ ay rho)]) (fl* k k))])
758                               (fl-make-rectangular
759                                  (fl/ (fllog1+ (fl/ (fl* 4.0 x)
760                                                     (fl+ (fl* 1-x 1-x) k)))
761                                       4.0)
762                                  (fl/ (flatan2 (fl* 2.0 y)
763                                                (fl- (fl* 1-x (fl+ 1.0 x)) k))
764                                       -2.0)))])))))])
765      (lambda (z)
766         (let ([x (cfl-real-part z)] [y (cfl-imag-part z)])
767            (if (negated-flonum? x)
768                (cfl- (f (fl- x) y))
769                (f x (fl- y)))))))
770
771(define cflatan
772   ; from Kahan
773   ; principal expression:
774   ; arctanh(zi)/i
775   (lambda (z)
776      (cfl/i (cflatanh (cfl*i z)))))
777
778(define exact-inexact+
779   (lambda (x y)
780      (cond
781         [(fixnum? x) (if (fx= x 0) y (fl+ (fixnum->flonum x) y))]
782         [(or (floatable? x) (fl= y 0.0)) (fl+ (inexact x) y)]
783         [(exceptional-flonum? y) y]
784         [else (inexact (+ x (exact y)))])))
785
786(define exact-inexact-
787   (lambda (x y)
788      (cond
789         [(fixnum? x) (if (fx= x 0) (fl- y) (fl- (fixnum->flonum x) y))]
790         [(or (floatable? x) (fl= y 0.0)) (fl- (inexact x) y)]
791         [(exceptional-flonum? y) (fl- y)]
792         [else (inexact (- x (exact y)))])))
793
794(define inexact-exact-
795   (lambda (x y)
796      (cond
797         [(fixnum? y) (fl- x (fixnum->flonum y))]
798         [(or (floatable? y) (fl= x 0.0)) (fl- x (inexact y))]
799         [(exceptional-flonum? x) x]
800         [else (inexact (- (exact x) y))])))
801
802(define exact-inexact*
803   (lambda (x y)
804      (cond
805         [(fixnum? x) (if (fx= x 0) 0 (fl* (fixnum->flonum x) y))]
806         [(floatable? x) (fl* (inexact x) y)]
807         [(or (fl= y 0.0) (exceptional-flonum? y)) (if (< x 0) (fl- y) y)]
808         [else (inexact (* x (exact y)))])))
809
810(define exact-inexact/
811   (lambda (x y)
812      (cond
813         [(fixnum? x) (if (fx= x 0) 0 (fl/ (fixnum->flonum x) y))]
814         [(floatable? x) (fl/ (inexact x) y)]
815         [(or (fl= y 0.0) (exceptional-flonum? y))
816          (if (< x 0) (fl/ -1.0 y) (fl/ y))]
817         [else (inexact (/ x (exact y)))])))
818
819(define inexact-exact/
820   (lambda (x y)
821      (cond
822         [(fixnum? y) (if (eq? y 0) (domain-error '/ y) (fl/ x (fixnum->flonum y)))]
823         [(floatable? y) (fl/ x (inexact y))]
824         [(or (fl= x 0.0) (exceptional-flonum? x)) (if (< y 0) (fl- x) x)]
825         [else (inexact (/ (exact x) y))])))
826
827(define floatable?
828   ; x is "floatable" if it can be made inexact without overflow or underflow
829   (lambda (x)
830      (type-case x
831         [(fixnum?) #t]
832         [(bignum?) (fx<= (integer-length x) (constant max-float-exponent))]
833         [(ratnum?) (fx<= (constant min-float-exponent)
834                          (fx- (integer-length (numerator x))
835                               (integer-length (denominator x)))
836                          (constant max-float-exponent))]
837         [($exactnum?) (and (floatable? (real-part x))
838                             (floatable? (imag-part x)))]
839         [else #t])))
840
841(define exact-integer-fits-float?
842  (lambda (x)
843    (<= (- (expt 2 53)) x (expt 2 53))))
844
845(define exact-inexact-compare?
846   ; e is an exact number, i is a flonum
847   ; Preserve transitivity by making i exact,
848   ; unless i is +/-infinity or a NAN, in which case any normal flonum
849   ; is a safe representation of e for comparative purposes.
850   (lambda (pred? e i)
851      (float-type-case
852         [(ieee)
853          (if (exceptional-flonum? i)
854              (pred? 0.0 i)
855              (pred? e (exact i)))]
856         [else (pred? e (exact i))])))
857
858(define exact-sqrt
859   ; x must be exact
860   ; returns the exact square root if it exists, otherwise an approximation
861   (lambda (x)
862      (type-case x
863         [(fixnum? bignum?)
864          (if (< x 0) (make-rectangular 0 (isqrt (abs x))) (isqrt x))]
865         [(ratnum?)
866          (/ (exact-sqrt (numerator x)) (exact-sqrt (denominator x)))]
867         [else
868          (let ([rp (exact-sqrt (/ (+ (exact-sqrt (magnitude-squared x))
869                                      (real-part x))
870                                   2))])
871             (make-rectangular rp (/ (imag-part x) (* 2 rp))))])))
872
873(define ($fldiv-and-mod x y)
874  (if (negated-flonum? y)
875      (let ([q ($flfloor (fl/ x (fl- y)))])
876        (values (fl- q) (fl+ x (fl* y q))))
877      (let ([q ($flfloor (fl/ x y))])
878        (values q (fl- x (fl* y q))))))
879
880(define ($fldiv x y)
881  (if (negated-flonum? y)
882      (fl- ($flfloor (fl/ x (fl- y))))
883      ($flfloor (fl/ x y))))
884
885(define ($flmod x y)
886  (if (negated-flonum? y)
887      (fl+ x (fl* y ($flfloor (fl/ x (fl- y)))))
888      (fl- x (fl* y ($flfloor (fl/ x y))))))
889
890(define ($fldiv0-and-mod0 x y)
891 ; there doesn't seem to be an easy way to do this...
892  (let-values ([(d m) ($fldiv-and-mod x y)])
893    (if (fl> y 0.0)
894        (if (fl< m (fl/ y 2.0))
895            (values d m)
896            (values (fl+ d 1.0) (fl- m y)))
897        (if (fl< m (fl/ y -2.0))
898            (values d m)
899            (values (fl- d 1.0) (fl+ m y))))))
900
901(define ($fldiv0 x y)
902  (let-values ([(d m) ($fldiv-and-mod x y)])
903    (if (fl> y 0.0)
904        (if (fl< m (fl/ y 2.0)) d (fl+ d 1.0))
905        (if (fl< m (fl/ y -2.0)) d (fl- d 1.0)))))
906
907(define ($flmod0 x y)
908  (let ([m ($flmod x y)])
909    (if (fl> y 0.0)
910        (if (fl< m (fl/ y 2.0)) m (fl- m y))
911        (if (fl< m (fl/ y -2.0)) m (fl+ m y)))))
912
913(define ($fxdiv-and-mod x y who) ; signal error on overflow if who != #f, otherwise return bignum
914  (if (fx< x 0)
915      (if (fx< y 0)
916          (if (fx> x y) ; |x| < |y| => q = 0, r = x != 0
917              (values 1 (fx- x y))
918              (if (and (fx= y -1) (fx= x (most-negative-fixnum)))
919                  (if who
920                      ($impoops who "fixnum overflow with arguments ~s and ~s" x y)
921                      (values (- (most-negative-fixnum)) 0))
922                  (let* ([q (fxquotient x y)] [r (fx- x (fx* y q))])
923                    (if (fx= r 0) (values q 0) (values (fx+ q 1) (fx- r y))))))
924          (if (fx> x (fx- y)) ; |x| < |y| => q = 0, r = x != 0
925              (values -1 (fx+ x y))
926              (let* ([q (fxquotient x y)] [r (fx- x (fx* y q))])
927                (if (fx= r 0) (values q 0) (values (fx- q 1) (fx+ r y))))))
928      (if (or (fx< x y) (fx> (fx- x) y)) ; |x| < |y| => q = 0, r = x
929          (values 0 x)
930          (let ([q (fxquotient x y)])
931            (values q (fx- x (fx* y q)))))))
932
933(define ($fxdiv x y who) ; signal error on overflow if who != #f, otherwise return bignum
934  (if (fx< x 0)
935      (if (fx< y 0)
936          (if (fx> x y) ; |x| < |y| => q = 0, r = x != 0
937              1
938              (if (and (fx= y -1) (fx= x (most-negative-fixnum)))
939                  (if who
940                      ($impoops who "fixnum overflow with arguments ~s and ~s" x y)
941                      (- (most-negative-fixnum)))
942                  (let ([q (fxquotient x y)])
943                    (if (fx= x (fx* y q)) q (fx+ q 1)))))
944          (if (fx> x (fx- y)) ; |x| < |y| => q = 0, r = x != 0
945              -1
946              (let ([q (fxquotient x y)])
947                (if (fx= x (fx* y q)) q (fx- q 1)))))
948      (if (or (fx< x y) (fx> (fx- x) y)) ; |x| < |y| => q = 0, r = x
949          0
950          (fxquotient x y))))
951
952(define ($fxmod x y) ; no overflow possible
953  (if (fx< x 0)
954      (if (fx< y 0)
955          (if (fx> x y) ; |x| < |y| => q = 0, r = x != 0
956              (fx- x y)
957              (if (and (fx= y -1) (fx= x (most-negative-fixnum)))
958                  0
959                  (let* ([q (fxquotient x y)] [r (fx- x (fx* y q))])
960                    (if (fx= r 0) 0 (fx- r y)))))
961          (if (fx> x (fx- y)) ; |x| < |y| => q = 0, r = x != 0
962              (fx+ x y)
963              (let* ([q (fxquotient x y)] [r (fx- x (fx* y q))])
964                (if (fx= r 0) 0 (fx+ r y)))))
965      (if (or (fx< x y) (fx> (fx- x) y)) ; |x| < |y| => q = 0, r = x
966          x
967          (fx- x (fx* y (fxquotient x y))))))
968
969(define ($fxdiv0-and-mod0 x y who)
970  (let-values ([(d m) ($fxdiv-and-mod x y who)])
971    (if (fx> y 0)
972        (if (fx< m (if (fx= y (most-positive-fixnum))
973                       (ash (+ (most-positive-fixnum) 1) -1)
974                       (fxsrl (fx+ y 1) 1)))
975            (values d m)
976            (values (fx+ d 1) (fx- m y)))
977        (if (fx< m (if (fx= y (most-negative-fixnum))
978                       (ash (- 1 (most-negative-fixnum)) -1)
979                       (fxsrl (fx- 1 y) 1)))
980            (values d m)
981            (values (fx- d 1) (fx+ m y))))))
982
983(define ($fxdiv0 x y who)
984  (let-values ([(d m) ($fxdiv-and-mod x y who)])
985    (if (fx> y 0)
986        (if (fx< m (if (fx= y (most-positive-fixnum))
987                       (ash (+ (most-positive-fixnum) 1) -1)
988                       (fxsrl (fx+ y 1) 1)))
989            d
990            (fx+ d 1))
991        (if (fx< m (if (fx= y (most-negative-fixnum))
992                       (ash (- 1 (most-negative-fixnum)) -1)
993                       (fxsrl (fx- 1 y) 1)))
994            d
995            (fx- d 1)))))
996
997(define ($fxmod0 x y)
998  (let ([m ($fxmod x y)])
999    (if (fx> y 0)
1000        (if (fx< m (if (fx= y (most-positive-fixnum))
1001                       (ash (+ (most-positive-fixnum) 1) -1)
1002                       (fxsrl (fx+ y 1) 1)))
1003            m
1004            (fx- m y))
1005        (if (fx< m (if (fx= y (most-negative-fixnum))
1006                       (ash (- 1 (most-negative-fixnum)) -1)
1007                       (fxsrl (fx- 1 y) 1)))
1008            m
1009            (fx+ m y)))))
1010
1011(define ($exdiv-and-mod x y) ; like $fldiv-and-mod
1012  (if (< y 0)
1013      (let ([q (floor (/ x (- y)))])
1014        (values (- q) (+ x (* y q))))
1015      (let ([q (floor (/ x y))])
1016        (values q (- x (* y q))))))
1017
1018(define ($exdiv0-and-mod0 x y)
1019  (let-values ([(d m) ($exdiv-and-mod x y)])
1020    (if (> y 0)
1021        (if (< m (/ y 2))
1022            (values d m)
1023            (values (+ d 1) (- m y)))
1024        (if (< m (/ y -2))
1025            (values d m)
1026            (values (- d 1) (+ m y))))))
1027
1028(define ($exdiv x y) ; like $fldiv
1029  (if (< y 0)
1030      (- (floor (/ x (- y))))
1031      (floor (/ x y))))
1032
1033(define ($exmod x y) ; like $flmod
1034  (if (< y 0)
1035      (+ x (* y (floor (/ x (- y)))))
1036      (- x (* y (floor (/ x y))))))
1037
1038(define $sll
1039  (lambda (who x n)
1040    (type-case n
1041      [(fixnum?)
1042       (unless (fx>= n 0) ($oops who "~s is not a nonnegative exact integer" n))
1043       (type-case x
1044         [(fixnum?)
1045          (let ([max-fx-shift (- (constant fixnum-bits) 1)])
1046            (if (fx> n max-fx-shift)
1047                (integer-ash x n)
1048                (let ([m (#3%fxsll x n)])
1049                  (if (fx= (fxsra m n) x)
1050                      m
1051                      (integer-ash x n)))))]
1052         [(bignum?) (integer-ash x n)]
1053         [else (nonexact-integer-error who x)])]
1054      [(bignum?)
1055       (unless ($bigpositive? n) ($oops who "~s is not a nonnegative exact integer" n))
1056       (type-case x
1057         [(fixnum? bignum?)
1058          (let ([k (most-positive-fixnum)])
1059            ($sll who ($sll who x k) (- n k)))]
1060         [else (nonexact-integer-error who x)])]
1061      [else (nonexact-integer-error who n)])))
1062
1063(define $sra
1064  (lambda (who x n)
1065    (type-case n
1066      [(fixnum?)
1067       (unless (fx>= n 0) ($oops who "~s is not a nonnegative exact integer" n))
1068       (type-case x
1069         [(fixnum?)
1070          (let ([max-fx-shift (- (constant fixnum-bits) 1)])
1071            (fxsra x (if (fx> n max-fx-shift) max-fx-shift n)))]
1072         [(bignum?) (integer-ash x (- n))]
1073         [else (nonexact-integer-error who x)])]
1074      [(bignum?)
1075       (unless ($bigpositive? n) ($oops who "~s is not a nonnegative exact integer" n))
1076       (type-case x
1077         [(fixnum? bignum?)
1078          (let ([k (most-positive-fixnum)])
1079            ($sra who ($sra who x k) (- n k)))]
1080         [else (nonexact-integer-error who x)])]
1081      [else (nonexact-integer-error who n)])))
1082
1083(define $negate
1084  (lambda (who x)
1085    (type-case x
1086      [(fixnum?)
1087       (if (fx= x (most-negative-fixnum))
1088           (let-syntax ([a (lambda (x) (- (constant most-negative-fixnum)))]) a)
1089           (fx- x))]
1090      [(bignum?) (big-negate x)]
1091      [(flonum?) (fl- x)]
1092      [(ratnum?) (make-ratnum (- ($ratio-numerator x)) ($ratio-denominator x))]
1093      [($exactnum? $inexactnum?) (make-rectangular (- (real-part x)) (- (imag-part x)))]
1094      [else (nonnumber-error who x)])))
1095
1096(set! integer?
1097  (lambda (x)
1098    (type-case x
1099      [(fixnum? bignum?) #t]
1100      [(flonum?) ($flinteger? x)]
1101      [else #f])))
1102
1103(set! integer-valued?
1104  (lambda (x)
1105    (type-case x
1106      [(fixnum? bignum?) #t]
1107      [(flonum?) ($flinteger? x)]
1108      [($inexactnum?)
1109       (and (fl= ($inexactnum-imag-part x) 0.0)
1110            ($flinteger? ($inexactnum-real-part x)))]
1111      [else #f])))
1112
1113(set! rational?
1114  (lambda (x)
1115    (type-case x
1116      [(fixnum? bignum? ratnum?) #t]
1117      [(flonum?) (not (exceptional-flonum? x))]
1118      [else #f])))
1119
1120(set! rational-valued?
1121  (lambda (x)
1122    (type-case x
1123      [(fixnum? bignum? ratnum?) #t]
1124      [(flonum?) (not (exceptional-flonum? x))]
1125      [($inexactnum?) (fl= ($inexactnum-imag-part x) 0.0)]
1126      [else #f])))
1127
1128(set! real?
1129  (lambda (x)
1130    (type-case x
1131      [(fixnum? flonum? bignum? ratnum?) #t]
1132      [else #f])))
1133
1134(set! real-valued?
1135  (lambda (x)
1136    (type-case x
1137      [(fixnum? flonum? bignum? ratnum?) #t]
1138      [($inexactnum?) (fl= ($inexactnum-imag-part x) 0.0)]
1139      [else #f])))
1140
1141(set! complex?
1142 ; same as number?
1143  (lambda (x)
1144    (type-case x
1145      [(fixnum? cflonum? bignum? ratnum? $exactnum?) #t]
1146      [else #f])))
1147
1148(set! number?
1149 ; same as complex?
1150  (lambda (x)
1151    (type-case x
1152      [(fixnum? cflonum? bignum? ratnum? $exactnum?) #t]
1153      [else #f])))
1154
1155(set! exact?
1156   (lambda (x)
1157      (type-case x
1158         [(fixnum?) #t]
1159         [(cflonum?) #f]
1160         [(bignum? ratnum? $exactnum?) #t]
1161         [else (nonnumber-error 'exact? x)])))
1162
1163(set! inexact?
1164   (lambda (x)
1165      (type-case x
1166         [(cflonum?) #t]
1167         [(fixnum? bignum? ratnum? $exactnum?) #f]
1168         [else (nonnumber-error 'inexact? x)])))
1169
1170(set-who! numerator
1171  (lambda (x)
1172    (type-case x
1173      [(ratnum?) ($ratio-numerator x)]
1174      [(fixnum? bignum?) x]
1175      [(flonum?)
1176       (cond
1177         [(exceptional-flonum? x) (nonrational-error who x)]
1178         [($flinteger-or-inf? x) x]
1179         [else (inexact (numerator (exact x)))])]
1180      [else (nonrational-error who x)])))
1181
1182(set-who! denominator
1183  (lambda (x)
1184    (type-case x
1185      [(ratnum?) ($ratio-denominator x)]
1186      [(fixnum? bignum?) 1]
1187      [(flonum?)
1188       (cond
1189         [(exceptional-flonum? x) (nonrational-error who x)]
1190         [($flinteger-or-inf? x) 1.0]
1191         [else (inexact (denominator (exact x)))])]
1192      [else (nonrational-error who x)])))
1193
1194(set! real-part
1195   (lambda (z)
1196      (type-case z
1197         [($inexactnum?) ($inexactnum-real-part z)]
1198         [($exactnum?) ($exactnum-real-part z)]
1199         [(flonum? fixnum? bignum? ratnum?) z]
1200         [else (noncomplex-error 'real-part z)])))
1201
1202(set! imag-part
1203   (lambda (z)
1204      (type-case z
1205         [($inexactnum?) ($inexactnum-imag-part z)]
1206         [($exactnum?) ($exactnum-imag-part z)]
1207         [(flonum?) 0]
1208         [(fixnum? bignum? ratnum?) 0]
1209         [else (noncomplex-error 'imag-part z)])))
1210
1211(set! modulo
1212   (lambda (x y)
1213      (unless (integer? x) (noninteger-error 'modulo x))
1214      (unless (integer? y) (noninteger-error 'modulo y))
1215      (when (= y 0) (domain-error 'modulo y))
1216      (let ([r (remainder x y)])
1217         (if (if (negative? y) (positive? r) (negative? r))
1218             (+ r y)
1219             r))))
1220
1221(set! expt-mod
1222   ; (modulo (expt x y) z)
1223   (lambda (x y z)
1224      (unless (integer? x)
1225         ($oops 'expt-mod "~s is not an integer" x))
1226      (unless (and (integer? y) (not (negative? y)))
1227         ($oops 'expt-mod "~s is not a nonnegative integer" y))
1228      (unless (and (integer? z) (not (zero? z)))
1229         ($oops 'expt-mod "~s is not a nonzero integer" z))
1230      (if (= y 0)
1231          (modulo 1 z)
1232          (do ([w 1 (if (even? y) w (remainder (* w b) z))]
1233               [y y (quotient y 2)]
1234               [b (remainder x z) (remainder (* b b) z)])
1235              ((= y 1) (modulo (* w b) z))))))
1236
1237(set-who! negative?
1238  (lambda (x)
1239    (type-case x
1240      [(fixnum?) (fx< x 0)]
1241      [(flonum?) (fl< x 0.0)]
1242      [(bignum?) (not ($bigpositive? x))]
1243      [(ratnum?) (< ($ratio-numerator x) 0)]
1244      [else (nonreal-error who x)])))
1245
1246(set-who! nonnegative?
1247  (lambda (x)
1248    (type-case x
1249      [(fixnum?) (fx>= x 0)]
1250      [(flonum?) (fl>= x 0.0)]
1251      [(bignum?) ($bigpositive? x)]
1252      [(ratnum?) (>= ($ratio-numerator x) 0)]
1253      [else (nonreal-error who x)])))
1254
1255(set-who! positive?
1256  (lambda (x)
1257    (type-case x
1258      [(fixnum?) (fx> x 0)]
1259      [(flonum?) (fl> x 0.0)]
1260      [(bignum?) ($bigpositive? x)]
1261      [(ratnum?) (> ($ratio-numerator x) 0)]
1262      [else (nonreal-error who x)])))
1263
1264(set-who! nonpositive?
1265  (lambda (x)
1266    (type-case x
1267      [(fixnum?) (fx<= x 0)]
1268      [(flonum?) (fl<= x 0.0)]
1269      [(bignum?) (not ($bigpositive? x))]
1270      [(ratnum?) (<= ($ratio-numerator x) 0)]
1271      [else (nonreal-error who x)])))
1272
1273(set-who! min
1274  (rec min
1275    (case-lambda
1276      [(x y)
1277       (type-case x
1278          [(flonum?)
1279           (type-case y
1280              [(flonum?) (if (or (fl< x y) ($nan? x)) x y)]
1281              [(fixnum? bignum? ratnum?) (min x (inexact y))]
1282              [else (nonreal-error who y)])]
1283          [(fixnum?)
1284           (type-case y
1285              [(fixnum?) (if (fx< x y) x y)]
1286              [(bignum? ratnum?) (if (< x y) x y)]
1287              [(flonum?) (min (inexact x) y)]
1288              [else (nonreal-error who y)])]
1289          [(bignum? ratnum?)
1290           (type-case y
1291              [(fixnum? bignum? ratnum?) (if (< x y) x y)]
1292              [(flonum?) (min (inexact x) y)]
1293              [else (nonreal-error who y)])]
1294          [else (nonreal-error who x)])]
1295      [(x) (if (real? x) x (nonreal-error who x))]
1296      [(x y . z)
1297       (let loop ([x (min x y)] [z z])
1298          (if (null? z)
1299              x
1300              (loop (min x (car z)) (cdr z))))])))
1301
1302(set-who! max
1303  (rec max
1304    (case-lambda
1305      [(x y)
1306       (type-case x
1307          [(flonum?)
1308           (type-case y
1309              [(flonum?) (if (or (fl> x y) ($nan? x)) x y)]
1310              [(fixnum? bignum? ratnum?) (max x (inexact y))]
1311              [else (nonreal-error who y)])]
1312          [(fixnum?)
1313           (type-case y
1314              [(fixnum?) (if (fx> x y) x y)]
1315              [(bignum? ratnum?) (if (> x y) x y)]
1316              [(flonum?) (max (inexact x) y)]
1317              [else (nonreal-error who y)])]
1318          [(bignum? ratnum?)
1319           (type-case y
1320              [(fixnum? bignum? ratnum?) (if (> x y) x y)]
1321              [(flonum?) (max (inexact x) y)]
1322              [else (nonreal-error who y)])]
1323          [else (nonreal-error who x)])]
1324      [(x) (if (real? x) x (nonreal-error who x))]
1325      [(x y . z)
1326       (let loop ([x (max x y)] [z z])
1327          (if (null? z)
1328              x
1329              (loop (max x (car z)) (cdr z))))])))
1330
1331(let ()
1332  (define (exlcm x1 x2)
1333    (if (or (eqv? x1 0) (eqv? x2 0))
1334        0
1335        (* (abs x1) (/ (abs x2) (exgcd x1 x2)))))
1336
1337  (set-who! gcd
1338    (rec gcd
1339      (case-lambda
1340        [() 0]
1341        [(x1) (gcd x1 x1)]
1342        [(x1 x2)
1343         (if (and (or (fixnum? x1) (bignum? x1))
1344                  (or (fixnum? x2) (bignum? x2)))
1345             (exgcd x1 x2)
1346             (begin
1347               (unless (integer? x1) (noninteger-error who x1))
1348               (unless (integer? x2) (noninteger-error who x2))
1349               (inexact (exgcd (exact x1) (exact x2)))))]
1350        [(x1 x2 . xr)
1351         (let f ([x1 x1] [x2 x2] [xr xr])
1352           (let ([x1 (gcd x1 x2)])
1353             (if (null? xr) x1 (f x1 (car xr) (cdr xr)))))])))
1354
1355  (set-who! lcm
1356    (rec lcm
1357      (case-lambda
1358        [() 1]
1359        [(x) (lcm x x)]
1360        [(x1 x2)
1361         (if (and (or (fixnum? x1) (bignum? x1))
1362                  (or (fixnum? x2) (bignum? x2)))
1363             (exlcm x1 x2)
1364             (begin
1365               (unless (integer? x1) (noninteger-error who x1))
1366               (unless (integer? x2) (noninteger-error who x2))
1367               (inexact (exlcm (exact x1) (exact x2)))))]
1368        [(x1 x2 . xr)
1369         (let f ([x1 x1] [x2 x2] [xr xr])
1370           (let ([x1 (lcm x1 x2)])
1371             (if (null? xr) x1 (f x1 (car xr) (cdr xr)))))]))))
1372
1373(let ()
1374  (define convert-to-inexact
1375    (lambda (z who)
1376      (type-case z
1377        [(fixnum?) (fixnum->flonum z)]
1378        [(bignum? ratnum?) (float z)]
1379        [($exactnum?)
1380         (fl-make-rectangular (inexact ($exactnum-real-part z))
1381                              (inexact ($exactnum-imag-part z)))]
1382        [(cflonum?) z]
1383        [else (nonnumber-error who z)])))
1384  (set-who! inexact (lambda (z) (convert-to-inexact z who)))
1385  (set-who! exact->inexact (lambda (z) (convert-to-inexact z who))))
1386
1387(let ()
1388  (define convert-to-exact
1389    (lambda (z who)
1390      (type-case z
1391        [(flonum?)
1392         (when (exceptional-flonum? z)
1393           ($oops 'exact "no exact representation for ~s" z))
1394         (let ([dx (decode-float z)])
1395           (let ([mantissa (* (vector-ref dx 0) (vector-ref dx 2))]
1396                 [exponent (vector-ref dx 1)])
1397             (if (fx< exponent 0)
1398                 (/ mantissa (ash 1 (fx- exponent)))
1399                 (* mantissa (ash 1 exponent)))))]
1400        [($inexactnum?)
1401         (make-rectangular
1402           (exact ($inexactnum-real-part z))
1403           (exact ($inexactnum-imag-part z)))]
1404        [(fixnum? bignum? ratnum? $exactnum?) z]
1405        [else (nonnumber-error who z)])))
1406  (set-who! exact (lambda (z) (convert-to-exact z who)))
1407  (set-who! inexact->exact (lambda (z) (convert-to-exact z who))))
1408
1409(set! rationalize
1410   ; Alan Bawden's algorithm
1411   (letrec
1412      ([rat1 ; x < y
1413        (lambda (x y)
1414           (cond
1415              [(> x 0) (rat2 x y)]
1416              [(< y 0) (- (rat2 (- y) (- x)))]
1417              [else (if (and (exact? x) (exact? y)) 0 0.0)]))]
1418       [rat2 ; 0 < x < y
1419        (lambda (x y)
1420           (let ([fx (floor x)] [fy (floor y)])
1421              (cond
1422                 [(= fx x) fx]
1423                 [(= fx fy) (+ fx (/ (rat2 (/ (- y fy)) (/ (- x fx)))))]
1424                 [else (+ fx 1)])))])
1425      (lambda (x e)
1426         (unless (real? x) (nonreal-error 'rationalize x))
1427         (unless (real? e) (nonreal-error 'rationalize e))
1428         (let ([x (- x e)] [y (+ x e)])
1429           (cond
1430              [(< x y) (rat1 x y)]
1431              [(< y x) (rat1 y x)]
1432              [else x])))))
1433
1434(set! abs
1435   (lambda (z)
1436      (type-case z
1437         [(fixnum?) (if (fx< z 0) (if (fx= z (most-negative-fixnum)) (- (most-negative-fixnum)) (fx- z)) z)]
1438         [(flonum?) (flabs z)]
1439         [(bignum?) (if ($bigpositive? z) z (- z))]
1440         [(ratnum?) (if (< z 0) (- z) z)]
1441         [else (nonreal-error 'abs z)])))
1442
1443(set! magnitude
1444   (lambda (z)
1445      (type-case z
1446         [(flonum?) (flabs z)]
1447         [(fixnum?) (if (fx< z 0) (- z) z)]
1448         [($inexactnum?) (cflmagnitude z)]
1449         [($exactnum?)
1450          (let ([x ($exactnum-real-part z)] [y ($exactnum-imag-part z)])
1451             (sqrt (+ (* x x) (* y y))))]
1452         [(bignum?) (if ($bigpositive? z) z (- z))]
1453         [(ratnum?) (if (< z 0) (- z) z)]
1454         [else (noncomplex-error 'magnitude z)])))
1455
1456(set! angle
1457   (lambda (z)
1458      (type-case z
1459         [(flonum?) (cond
1460                     [($nan? z) +nan.0]
1461                     [(negated-flonum? z) pi]
1462                     [else 0])]
1463         [($inexactnum?) (cflangle z)]
1464         [(fixnum? bignum? ratnum?)
1465          (cond
1466             [(< z 0) pi]
1467             [(> z 0) 0]
1468             [else (domain-error 'angle z)])]
1469         [($exactnum?)
1470          ; use single argument atan to avoid precision loss
1471          ; cases from Kahan
1472          (let ([x ($exactnum-real-part z)] [y ($exactnum-imag-part z)])
1473             (cond
1474                [(> (abs y) (abs x))
1475                 (- (fl* pi (if (< y 0) -.5 .5)) (atan (/ x y)))]
1476                [(< x 0)
1477                 (if (< y 0)
1478                     (- (atan (/ y x)) pi)
1479                     (+ (atan (/ y x)) pi))]
1480                [else (atan (/ y x))]))]
1481         [else (noncomplex-error 'angle z)])))
1482
1483(set-who! make-rectangular
1484  (lambda (x y)
1485    (type-case y
1486      [(flonum?)
1487       (fl-make-rectangular
1488          (type-case x
1489            [(flonum?) x]
1490            [(fixnum? bignum? ratnum?) (inexact x)]
1491            [else (nonreal-error who x)])
1492          y)]
1493      [(fixnum? bignum? ratnum?)
1494       (if (eq? y 0)
1495           (if (real? x) x (nonreal-error who x))
1496           (type-case x
1497             [(fixnum? bignum? ratnum?) ($make-exactnum x y)]
1498             [(flonum?) (fl-make-rectangular x (inexact y))]
1499             [else (nonreal-error who x)]))]
1500      [else (nonreal-error who y)])))
1501
1502(set-who! make-polar
1503  (lambda (x y)
1504    (unless (real? x) (nonreal-error 'make-polar x))
1505    (unless (real? y) (nonreal-error 'make-polar y))
1506    (cond
1507      [(eq? y 0) x]
1508      [(eq? x 0) 0]
1509      [else (make-rectangular (* x (cos y)) (* x (sin y)))])))
1510
1511(set! log
1512  (rec log
1513    (case-lambda
1514      [(x)
1515       (type-case x
1516         [(flonum?)
1517          (if (fl< x 0.0)
1518              (fl-make-rectangular ($fllog (fl- x)) pi)
1519              ($fllog x))]
1520         [($inexactnum?) (cfllog x)]
1521         [(fixnum?)
1522          (cond
1523             [(fx> x 1) ($fllog (fixnum->flonum x))]
1524             [(fx= x 1) 0]
1525             [(fx< x 0) (make-rectangular (log (- x)) pi)]
1526             [else (domain-error 'log x)])]
1527         [(bignum?)
1528          (let ([len (integer-length x)])
1529             (if (fx<= len (constant max-float-exponent))
1530                 (log (inexact x))
1531                 (+ (* len log2) (log (inexact (/ x (ash 1 len)))))))]
1532         [(ratnum?)
1533          (if (floatable? x)
1534              (log (inexact x))
1535              (- (log (numerator x)) (log (denominator x))))]
1536         [($exactnum?)
1537          (make-rectangular
1538            (/ (log (magnitude-squared x)) 2)
1539            (angle x))]
1540         [else (nonnumber-error 'log x)])]
1541      [(x y) (/ (log x) (log y))])))
1542
1543(define-trig-op exp $flexp cflexp 1)
1544(define-trig-op sin $flsin cflsin 0)
1545(define-trig-op cos $flcos cflcos 1)
1546(define-trig-op tan $fltan cfltan 0)
1547
1548(set! asin
1549   (lambda (x)
1550      (type-case x
1551         [(flonum?)
1552          ; make sure NANs go the "$flasin" route
1553          (if (or (fl< x -1.0) (fl> x 1.0))
1554              (cflasin x)
1555              ($flasin x))]
1556         [($inexactnum?) (cflasin x)]
1557         [(fixnum?) (if (fx= x 0) 0 (asin (fixnum->flonum x)))]
1558         [(bignum? ratnum? $exactnum?) (asin (inexact x))]
1559         [else (nonnumber-error 'asin x)])))
1560
1561(set! acos
1562   (lambda (x)
1563      (type-case x
1564         [(flonum?)
1565          ; make sure NANs go the "$flacos" route
1566          (if (or (fl< x -1.0) (fl> x 1.0))
1567              (cflacos x)
1568              ($flacos x))]
1569         [($inexactnum?) (cflacos x)]
1570         [(fixnum?) (if (fx= x 1) 0 (acos (fixnum->flonum x)))]
1571         [(bignum? ratnum? $exactnum?) (acos (inexact x))]
1572         [else (nonnumber-error 'acos x)])))
1573
1574(set! atan
1575   (case-lambda
1576      [(x)
1577       (type-case x
1578          [(flonum?) ($flatan x)]
1579          [($inexactnum?) (cflatan x)]
1580          [(fixnum?) (if (fx= x 0) 0 (atan (fixnum->flonum x)))]
1581          [(bignum? ratnum?) (atan (inexact x))]
1582          [($exactnum?)
1583           (when (or (= x +i) (= x -i)) (domain-error 'atan x))
1584           (atan (inexact x))]
1585          [else (nonnumber-error 'atan x)])]
1586      [(y x)
1587       (cond
1588          [(and (flonum? y) (flonum? x))
1589           (flatan2 y x)]
1590          [(and (fixnum? y) (fixnum? x))
1591           (if (fx= y 0)
1592               (cond
1593                  [(fx> x 0) 0]
1594                  [(fx< x 0) pi]
1595                  [else (domain-error2 'atan2 y x)])
1596               (flatan2 (fixnum->flonum y) (fixnum->flonum x)))]
1597          [else
1598           (unless (real? y) (nonreal-error 'atan y))
1599           (unless (real? x) (nonreal-error 'atan x))
1600           (angle (make-rectangular x y))])]))
1601
1602(define-trig-op sinh flsinh cflsinh 0)
1603(define-trig-op cosh flcosh cflcosh 1)
1604(define-trig-op tanh fltanh cfltanh 0)
1605(define-trig-op asinh flasinh cflasinh 0)
1606
1607(set! acosh
1608   (lambda (x)
1609      (type-case x
1610         [(flonum?)
1611          ; make sure NANs go the "flacosh" route
1612          (if (fl< x 1.0) (cflacosh x) (flacosh x))]
1613         [($inexactnum?) (cflacosh x)]
1614         [(fixnum?) (if (fx= x 1) 0 (acosh (fixnum->flonum x)))]
1615         [(bignum? ratnum? $exactnum?) (acosh (inexact x))]
1616         [else (nonnumber-error 'acosh x)])))
1617
1618(set! atanh
1619   (lambda (x)
1620      (type-case x
1621         [(flonum?)
1622          ; make sure NANs go the "flatanh" route
1623          (if (or (fl< x -1.0) (fl> x 1.0))
1624              (cflatanh x)
1625              (flatanh x))]
1626         [($inexactnum?) (cflatanh x)]
1627         [(fixnum?)
1628          (cond
1629             [(or (fx> x 1) (fx< x -1)) (atanh (fixnum->flonum x))]
1630             [(fx= x 0) 0]
1631             [else (domain-error 'atan x)])]
1632         [(bignum? ratnum? $exactnum?) (atanh (inexact x))]
1633         [else (nonnumber-error 'atanh x)])))
1634
1635; exceptional cases from Steele(CLtL), page 311
1636(set! expt
1637   (lambda (x y)
1638      (type-case y
1639         [(fixnum? bignum?)
1640          (cond
1641            [(and (eq? y 0) (number? x)) 1]
1642            [(eq? x 0)
1643             (if (< y 0)
1644                 ($impoops 'expt "undefined for values ~s and ~s" x y)
1645                 0)]
1646            [(eq? x 1) 1]
1647            [(eq? x 2) (if (< y 0) (/ (ash 1 (- y))) (ash 1 y))]
1648            [(and (flonum? x) (exact-integer-fits-float? y)) ($flexpt x (inexact y))]
1649            [(and ($inexactnum? x) (exact-integer-fits-float? y)) (exp (* y (log x)))]
1650            [(not (number? x)) (nonnumber-error 'expt x)]
1651            [(ratnum? x)
1652             (if (< y 0)
1653                 (let ([y (- y)])
1654                   (/ (expt (denominator x) y) (expt (numerator x) y)))
1655                 (/ (expt (numerator x) y) (expt (denominator x) y)))]
1656            [else
1657             (let ()
1658               (define (f x n)
1659                 (let loop ([i (integer-length n)] [a 1])
1660                   (let ([a (if (bitwise-bit-set? n i)
1661                                (* a x)
1662                                a)])
1663                     (if (fx= i 0)
1664                         a
1665                         (loop (fx- i 1)
1666                               (* a a))))))
1667               (if (< y 0)
1668                   (if (or (fixnum? x) (bignum? x))
1669                       (/ (f x (- y)))
1670                       (f (/ x) (- y)))
1671                   (f x y)))])]
1672         [(flonum?)
1673          (type-case x
1674             [(flonum?)
1675              (if (and (fl< x 0.0) (not ($flinteger? y)))
1676                  (exp (* y (log x)))
1677                  ($flexpt x y))]
1678             [($inexactnum? $exactnum?) (exp (* y (log x)))]
1679             [(fixnum? bignum? ratnum?)
1680              (cond
1681               [(eq? x 0)
1682                (cond
1683                 [(fl< y 0.0) ($impoops 'expt "undefined for values ~s and ~s" x y)]
1684                 [(fl= y 0.0) 1.0]
1685                 [($nan? y) +nan.0]
1686                 [else 0])]
1687               [(eq? x 1) 1]
1688               [else
1689                (if (floatable? x)
1690                    (expt (inexact x) y)
1691                    (exp (* y (log x))))])]
1692             [else (nonnumber-error 'expt x)])]
1693         [($inexactnum?)
1694          (cond
1695           [(eq? x 0)
1696            (let ([r ($inexactnum-real-part y)])
1697              (cond
1698               [(fl> r 0.0) 0]
1699               [else
1700                ($impoops 'expt "undefined for values ~s and ~s" x y)]))]
1701           [(eq? x 1) 1]
1702           [(and (flonum? x) (fl= x 0.0) (not (negated-flonum? ($inexactnum-real-part y))))
1703            0.0]
1704           [else
1705            (unless (number? x) (nonnumber-error 'expt x))
1706            (exp (* y (log x)))])]
1707         [(ratnum? $exactnum?)
1708          (unless (number? x) (nonnumber-error 'expt x))
1709          (cond
1710             [(eqv? y 1/2) (sqrt x)]
1711             [(eq? x 0)
1712              (if (> (real-part y) 0)
1713                  0
1714                  ($impoops 'expt "undefined for values ~s and ~s" x y))]
1715             [(eq? x 1) 1]
1716             [(and (floatable? y)
1717                   (let ([y (inexact y)])
1718                     ;; Don't use this case if `(inexact y)` loses
1719                     ;; precision and becomes an an integer, in which
1720                     ;; case the result would be real (but should be
1721                     ;; non-real complex)
1722                     (and (not (and (flonum? y)
1723                                    ($flinteger? y)))
1724                          y)))
1725              => (lambda (y) (expt x y))]
1726             [else (exp (* y (log x)))])]
1727         [else (nonnumber-error 'expt y)])))
1728
1729(set! sqrt
1730   (lambda (x)
1731      (type-case x
1732         [(flonum?)
1733          (if (fl< x 0.0)
1734              (fl-make-rectangular 0.0 ($flsqrt (flabs x)))
1735              ($flsqrt x))]
1736         [($inexactnum?) (cflsqrt x)]
1737         [(fixnum? bignum? ratnum? $exactnum?)
1738          (let ([y (exact-sqrt x)])
1739             (let ([yy (* y y)])
1740                (cond
1741                   [(= yy x) y]
1742                   [(floatable? x) (sqrt (inexact x))]
1743                   [else (* y (sqrt (inexact (/ x yy))))])))]
1744         [else (nonnumber-error 'sqrt x)])))
1745
1746(set! isqrt
1747   ; Based on code credited to "boyland@aspen.Berkeley.EDU" by
1748   ; Akira Kurihara (d34676@tansei.cc.u-tokyo.ac.jp)
1749   (lambda (n)
1750      (cond
1751         [(and (or (fixnum? n) (bignum? n)) (>= n 0))
1752          (let isqrt ([n n])
1753             (cond
1754                [(>= n 16) ; ensures k > 0
1755                 (let ([a1 (let ([k (ash (- (integer-length n) 1) -2)])
1756                              (ash (isqrt (ash n (- (ash k 1)))) k))])
1757                    (let ([q&r ($quotient-remainder n a1)])
1758                       (let ([a2 (car q&r)])
1759                          (let ([a3 (ash (+ a1 a2) -1)])
1760                             (if (odd? a2)
1761                                 a3
1762                                 (let ([d (- a3 a1)])
1763                                    (if (> (* d d) (cdr q&r))
1764                                        (- a3 1)
1765                                        a3)))))))]
1766                [(>= n  9) 3]
1767                [(>= n  4) 2]
1768                [(>= n  1) 1]
1769                [else 0]))]
1770         [(and (integer? n) (>= n 0)) (floor (sqrt n))]
1771         [else ($oops 'isqrt "~s is not a nonnegative integer" n)])))
1772
1773(set-who! floor
1774  (lambda (x)
1775    (type-case x
1776      [(fixnum? bignum?) x]
1777      [(flonum?) ($flfloor x)]
1778      [(ratnum?)
1779       (let ([y (quotient ($ratio-numerator x) ($ratio-denominator x))])
1780         (if (< x 0) (- y 1) y))]
1781      [else (nonreal-error who x)])))
1782
1783(set-who! ceiling
1784  (lambda (x)
1785    (type-case x
1786      [(fixnum? bignum?) x]
1787      [(flonum?) ($flceiling x)]
1788      [(ratnum?)
1789       (let ([y (quotient ($ratio-numerator x) ($ratio-denominator x))])
1790         (if (< x 0) y (+ y 1)))]
1791      [else (nonreal-error who x)])))
1792
1793(set-who! truncate
1794  (lambda (x)
1795    (type-case x
1796       [(fixnum? bignum?) x]
1797       [(flonum?) (if (negated-flonum? x) (fl- ($flfloor (flabs x))) ($flfloor x))]
1798       [(ratnum?) (quotient ($ratio-numerator x) ($ratio-denominator x))]
1799       [else (nonreal-error who x)])))
1800
1801(set-who! quotient
1802  (let ([f (lambda (x y) (truncate (/ x y)))])
1803    (lambda (x y)
1804      (type-case y
1805        [(fixnum?)
1806         (when (fx= y 0) (domain-error who y))
1807         (cond
1808           [(fx= y 1) (unless (integer? x) (noninteger-error who x)) x]
1809           [(fx= y -1) (unless (integer? x) (noninteger-error who x)) ($negate who x)]
1810           [else
1811             (type-case x
1812               [(fixnum?) (if (and (fx= y -1) (fx= x (most-negative-fixnum)))
1813                              (- (most-negative-fixnum))
1814                              (fxquotient x y))]
1815               [(bignum?) (intquotient x y)]
1816               [else
1817                 (unless (integer? x) (noninteger-error who x))
1818                 (f x y)])])]
1819        [(bignum?)
1820         (type-case x
1821           [(fixnum? bignum?) (intquotient x y)]
1822           [else
1823             (unless (integer? x) (noninteger-error who x))
1824             (f x y)])]
1825        [else
1826          (unless (integer? y) (noninteger-error who y))
1827          (unless (integer? x) (noninteger-error who x))
1828          (when (= y 0) (domain-error who y))
1829          (f x y)]))))
1830
1831(set-who! div-and-mod
1832  (lambda (x y)
1833    (type-case y
1834      [(fixnum?)
1835       (type-case x
1836         [(fixnum?)
1837          (when (fx= y 0) (domain-error who y))
1838          ($fxdiv-and-mod x y #f)]
1839         [(flonum?) ($fldiv-and-mod x (fixnum->flonum y))]
1840         [(bignum?)
1841          (cond
1842            [(fx= y 1) (values x 0)]
1843            [(fx= y -1) (values (big-negate x) 0)]
1844            [else
1845              (when (fx= y 0) (domain-error who y))
1846              (let ([q.r (intquotient-remainder x y)])
1847                (if ($bigpositive? x)
1848                    (values (car q.r) (cdr q.r))
1849                    (if (eq? (cdr q.r) 0)
1850                        (values (car q.r) 0)
1851                        (if (fx< y 0)
1852                            (values (+ (car q.r) 1) (fx- (cdr q.r) y))
1853                            (values (- (car q.r) 1) (fx+ (cdr q.r) y))))))])]
1854         [(ratnum?)
1855          (when (fx= y 0) (domain-error who y))
1856          ($exdiv-and-mod x y)]
1857         [else (domain-error who x)])]
1858      [(flonum?)
1859       (type-case x
1860         [(fixnum?) ($fldiv-and-mod (fixnum->flonum x) y)]
1861         [(flonum?) ($fldiv-and-mod x y)]
1862         [(bignum? ratnum?) ($fldiv-and-mod (real->flonum x) y)]
1863         [else (domain-error who x)])]
1864      [(bignum?)
1865       (type-case x
1866         [(fixnum?) ; know |x| < |y| => q = 0, r = x
1867          (if (fx< x 0)
1868              (if ($bigpositive? y) (values -1 (+ x y)) (values 1 (- x y)))
1869              (values 0 x))]
1870         [(flonum?) ($fldiv-and-mod x (real->flonum y))]
1871         [(bignum?)
1872          (let ([q.r (intquotient-remainder x y)])
1873            (if ($bigpositive? x)
1874                (values (car q.r) (cdr q.r))
1875                (if (eq? (cdr q.r) 0)
1876                    (values (car q.r) 0)
1877                    (if ($bigpositive? y)
1878                        (values (- (car q.r) 1) (+ (cdr q.r) y))
1879                        (values (+ (car q.r) 1) (- (cdr q.r) y))))))]
1880         [(ratnum?) ($exdiv-and-mod x y)]
1881         [else (domain-error who x)])]
1882      [(ratnum?)
1883       (type-case x
1884         [(fixnum? bignum? ratnum?) ($exdiv-and-mod x y)]
1885         [(flonum?) ($fldiv-and-mod x (real->flonum y))]
1886         [else (domain-error who x)])]
1887      [else (domain-error who y)])))
1888
1889(set-who! div
1890  (lambda (x y)
1891    (type-case y
1892      [(fixnum?)
1893       (type-case x
1894         [(fixnum?)
1895          (when (fx= y 0) (domain-error who y))
1896          ($fxdiv x y #f)]
1897         [(flonum?) ($fldiv x (fixnum->flonum y))]
1898         [(bignum?)
1899          (when (fx= y 0) (domain-error who y))
1900          (cond
1901            [(fx= y 1) x]
1902            [(fx= y -1) (big-negate x)]
1903            [else
1904             (if ($bigpositive? x)
1905                 (intquotient x y)
1906                 (let ([q.r (intquotient-remainder x y)])
1907                   (if (eq? (cdr q.r) 0)
1908                       (car q.r)
1909                       (if (fx< y 0)
1910                           (+ (car q.r) 1)
1911                           (- (car q.r) 1)))))])]
1912         [(ratnum?)
1913          (when (fx= y 0) (domain-error who y))
1914          ($exdiv x y)]
1915         [else (domain-error who x)])]
1916      [(flonum?)
1917       (type-case x
1918         [(fixnum?) ($fldiv (fixnum->flonum x) y)]
1919         [(flonum?) ($fldiv x y)]
1920         [(bignum? ratnum?) ($fldiv (real->flonum x) y)]
1921         [else (domain-error who x)])]
1922      [(bignum?)
1923       (type-case x
1924         [(fixnum?) ; know |x| < |y| => q = 0, r = x
1925          (if (fx< x 0) (if ($bigpositive? y) -1 1) 0)]
1926         [(flonum?) ($fldiv x (real->flonum y))]
1927         [(bignum?)
1928          (if ($bigpositive? x)
1929              (intquotient x y)
1930              (let ([q.r (intquotient-remainder x y)])
1931                (if (eq? (cdr q.r) 0)
1932                    (car q.r)
1933                    (if ($bigpositive? y)
1934                        (- (car q.r) 1)
1935                        (+ (car q.r) 1)))))]
1936         [(ratnum?) ($exdiv x y)]
1937         [else (domain-error who x)])]
1938      [(ratnum?)
1939       (type-case x
1940         [(fixnum? bignum? ratnum?) ($exdiv x y)]
1941         [(flonum?) ($fldiv x (real->flonum y))]
1942         [else (domain-error who x)])]
1943      [else (domain-error who y)])))
1944
1945(set-who! mod
1946  (lambda (x y)
1947    (type-case y
1948      [(fixnum?)
1949       (type-case x
1950         [(fixnum?)
1951          (when (fx= y 0) (domain-error who y))
1952          ($fxmod x y)]
1953         [(flonum?) ($flmod x (fixnum->flonum y))]
1954         [(bignum?)
1955          (when (fx= y 0) (domain-error who y))
1956          (cond
1957            [(or (fx= y 1) (fx= y -1)) 0]
1958            [else
1959             (if ($bigpositive? x)
1960                 (intremainder x y)
1961                 (let ([q.r (intquotient-remainder x y)])
1962                   (if (eq? (cdr q.r) 0)
1963                       0
1964                       (if (fx< y 0)
1965                           (fx- (cdr q.r) y)
1966                           (fx+ (cdr q.r) y)))))])]
1967         [(ratnum?)
1968          (when (fx= y 0) (domain-error who y))
1969          ($exmod x y)]
1970         [else (domain-error who x)])]
1971      [(flonum?)
1972       (type-case x
1973         [(fixnum?) ($flmod (fixnum->flonum x) y)]
1974         [(flonum?) ($flmod x y)]
1975         [(bignum? ratnum?) ($flmod (real->flonum x) y)]
1976         [else (domain-error who x)])]
1977      [(bignum?)
1978       (type-case x
1979         [(fixnum?) ; know |x| < |y| => q = 0, r = x
1980          (if (fx< x 0) (if ($bigpositive? y) (+ x y) (- x y)) x)]
1981         [(flonum?) ($flmod x (real->flonum y))]
1982         [(bignum?)
1983          (if ($bigpositive? x)
1984              (intremainder x y)
1985              (let ([q.r (intquotient-remainder x y)])
1986                (if (eq? (cdr q.r) 0)
1987                    0
1988                    (if ($bigpositive? y)
1989                        (+ (cdr q.r) y)
1990                        (- (cdr q.r) y)))))]
1991         [(ratnum?) ($exmod x y)]
1992         [else (domain-error who x)])]
1993      [(ratnum?)
1994       (type-case x
1995         [(fixnum? bignum? ratnum?) ($exmod x y)]
1996         [(flonum?) ($flmod x (real->flonum y))]
1997         [else (domain-error who x)])]
1998      [else (domain-error who y)])))
1999
2000(set-who! div0-and-mod0
2001  (lambda (x y)
2002    (type-case y
2003      [(fixnum?)
2004       (type-case x
2005         [(fixnum?)
2006          (when (fx= y 0) (domain-error who y))
2007          ($fxdiv0-and-mod0 x y #f)]
2008         [(flonum?) ($fldiv0-and-mod0 x (fixnum->flonum y))]
2009         [(bignum?)
2010          (cond
2011            [(fx= y 1) (values x 0)]
2012            [(fx= y -1) (values (big-negate x) 0)]
2013            [else
2014             (when (fx= y 0) (domain-error who y))
2015             ($exdiv0-and-mod0 x y)])]
2016         [(ratnum?)
2017          (when (fx= y 0) (domain-error who y))
2018          ($exdiv0-and-mod0 x y)]
2019         [else (domain-error who x)])]
2020      [(flonum?)
2021       (type-case x
2022         [(fixnum?) ($fldiv0-and-mod0 (fixnum->flonum x) y)]
2023         [(flonum?) ($fldiv0-and-mod0 x y)]
2024         [(bignum? ratnum?) ($fldiv0-and-mod0 (real->flonum x) y)]
2025         [else (domain-error who x)])]
2026      [(bignum?)
2027       (type-case x
2028         [(fixnum? bignum? ratnum?) ($exdiv0-and-mod0 x y)]
2029         [(flonum?) ($fldiv0-and-mod0 x (real->flonum y))]
2030         [else (domain-error who x)])]
2031      [(ratnum?)
2032       (type-case x
2033         [(fixnum? bignum? ratnum?) ($exdiv0-and-mod0 x y)]
2034         [(flonum?) ($fldiv0-and-mod0 x (real->flonum y))]
2035         [else (domain-error who x)])]
2036      [else (domain-error who y)])))
2037
2038(set-who! div0
2039  (lambda (x y)
2040    (define (exdiv0 x y)
2041      (let-values ([(d m) ($exdiv-and-mod x y)])
2042        (if (> y 0)
2043            (if (< m (/ y 2)) d (+ d 1))
2044            (if (< m (/ y -2)) d (- d 1)))))
2045    (type-case y
2046      [(fixnum?)
2047       (type-case x
2048         [(fixnum?)
2049          (when (fx= y 0) (domain-error who y))
2050          ($fxdiv0 x y #f)]
2051         [(flonum?) ($fldiv0 x (fixnum->flonum y))]
2052         [(bignum?)
2053          (cond
2054            [(fx= y 1) x]
2055            [(fx= y -1) (big-negate x)]
2056            [else
2057             (when (fx= y 0) (domain-error who y))
2058             (exdiv0 x y)])]
2059         [(ratnum?)
2060          (when (fx= y 0) (domain-error who y))
2061          (exdiv0 x y)]
2062         [else (domain-error who x)])]
2063      [(flonum?)
2064       (type-case x
2065         [(fixnum?) ($fldiv0 (fixnum->flonum x) y)]
2066         [(flonum?) ($fldiv0 x y)]
2067         [(bignum? ratnum?) ($fldiv0 (real->flonum x) y)]
2068         [else (domain-error who x)])]
2069      [(bignum?)
2070       (type-case x
2071         [(fixnum? bignum? ratnum?) (exdiv0 x y)]
2072         [(flonum?) ($fldiv0 x (real->flonum y))]
2073         [else (domain-error who x)])]
2074      [(ratnum?)
2075       (type-case x
2076         [(fixnum? bignum? ratnum?) (exdiv0 x y)]
2077         [(flonum?) ($fldiv0 x (real->flonum y))]
2078         [else (domain-error who x)])]
2079      [else (domain-error who y)])))
2080
2081(set-who! mod0
2082  (lambda (x y)
2083    (define (exmod0 x y)
2084      (let ([m ($exmod x y)])
2085        (if (> y 0)
2086            (if (< m (/ y 2)) m (- m y))
2087            (if (< m (/ y -2)) m (+ m y)))))
2088    (type-case y
2089      [(fixnum?)
2090       (type-case x
2091         [(fixnum?)
2092          (when (fx= y 0) (domain-error who y))
2093          ($fxmod0 x y)]
2094         [(flonum?) ($flmod0 x (fixnum->flonum y))]
2095         [(bignum?)
2096          (cond
2097            [(or (fx= y 1) (fx= y -1)) 0]
2098            [else
2099             (when (fx= y 0) (domain-error who y))
2100             (exmod0 x y)])]
2101         [(ratnum?)
2102          (when (fx= y 0) (domain-error who y))
2103          (exmod0 x y)]
2104         [else (domain-error who x)])]
2105      [(flonum?)
2106       (type-case x
2107         [(fixnum?) ($flmod0 (fixnum->flonum x) y)]
2108         [(flonum?) ($flmod0 x y)]
2109         [(bignum? ratnum?) ($flmod0 (real->flonum x) y)]
2110         [else (domain-error who x)])]
2111      [(bignum?)
2112       (type-case x
2113         [(fixnum? bignum? ratnum?) (exmod0 x y)]
2114         [(flonum?) ($flmod0 x (real->flonum y))]
2115         [else (domain-error who x)])]
2116      [(ratnum?)
2117       (type-case x
2118         [(fixnum? bignum? ratnum?) (exmod0 x y)]
2119         [(flonum?) ($flmod0 x (real->flonum y))]
2120         [else (domain-error who x)])]
2121      [else (domain-error who y)])))
2122
2123(set-who! remainder
2124  (let* ([fmod (cflop2 "(cs)mod")]
2125         [f (lambda (x y)
2126              (cond
2127                [(eqv? x 0) 0]
2128                [else
2129                 (let ([r (fmod (real->flonum x) (real->flonum y))])
2130                   (if (fl= r 0.0)
2131                       ;; Always return positive 0.0 --- not sure why,
2132                       ;; but Racket and other Schemes seem to agree
2133                       0.0
2134                       r))]))])
2135    (lambda (x y)
2136      (type-case y
2137        [(fixnum?)
2138         (when (fx= y 0) (domain-error who y))
2139         (cond
2140           [(or (fx= y 1) (fx= y -1)) (unless (integer? x) (noninteger-error who x)) 0]
2141           [else
2142             (type-case x
2143               [(fixnum?) (fxremainder x y)]
2144               [(bignum?) (intremainder x y)]
2145               [else
2146                 (unless (integer? x) (noninteger-error who x))
2147                 (f x y)])])]
2148        [(bignum?)
2149         (type-case x
2150           [(fixnum? bignum?) (intremainder x y)]
2151           [else
2152             (unless (integer? x) (noninteger-error who x))
2153             (f x y)])]
2154        [else
2155          (unless (integer? y) (noninteger-error who y))
2156          (unless (integer? x) (noninteger-error who x))
2157          (when (= y 0) (domain-error who y))
2158          (f x y)]))))
2159
2160(set-who! even?
2161   (lambda (x)
2162      (type-case x
2163         [(fixnum?) (fxeven? x)]
2164         [(bignum?) (not (bigodd? x))]
2165         [(flonum?)
2166          (when (exceptional-flonum? x) (noninteger-error who x))
2167          (let ([y (fl* ($flfloor (fl/ x 2.0)) 2.0)])
2168             (cond
2169                [(fl= x y) #t]
2170                [(fl= (fl+ y 1.0) x) #f]
2171                [else (noninteger-error who x)]))]
2172         [else
2173          (unless (integer? x) (noninteger-error who x))
2174          (even? (real-part x))])))
2175
2176(set-who! odd?
2177   (lambda (x)
2178      (type-case x
2179         [(fixnum?) (fxodd? x)]
2180         [(bignum?) (bigodd? x)]
2181         [(flonum?)
2182          (when (exceptional-flonum? x) (noninteger-error who x))
2183          (let ([y (fl* ($flfloor (fl/ x 2.0)) 2.0)])
2184             (cond
2185                [(fl= x y) #f]
2186                [(fl= (fl+ y 1.0) x) #t]
2187                [else (noninteger-error who x)]))]
2188         [else
2189          (unless (integer? x) (noninteger-error who x))
2190          (odd? (real-part x))])))
2191
2192(set-who! round
2193  (lambda (x)
2194    (type-case x
2195      [(flonum?) (flround x)]
2196      [(fixnum? bignum?) x]
2197      [(ratnum?)
2198       (let ([x1 (+ x 1/2)])
2199         (let ([x2 (floor x1)])
2200           (if (and (= x1 x2) (odd? x2))
2201               (- x2 1)
2202               x2)))]
2203      [else (nonreal-error who x)])))
2204
2205;;; help routines used by library entries
2206;;; they are fully generic, but the cases are organized to catch those
2207;;; the library routines don't check first
2208
2209(set! $=
2210   (lambda (who x y)
2211      (type-case x
2212         [(fixnum?)
2213          (type-case y
2214             [(fixnum?) (fx= x y)]
2215             [(bignum? ratnum? $exactnum?) #f]
2216             [(cflonum?) (cfl= (fixnum->flonum x) y)]
2217             [else (nonnumber-error who y)])]
2218         [(bignum?)
2219          (type-case y
2220             [(fixnum?) #f]
2221             [(bignum?) (big= x y)]
2222             [(ratnum? $exactnum?) #f]
2223             [(flonum?) (exact-inexact-compare? = x y)]
2224             [($inexactnum?) (and (fl= ($inexactnum-imag-part y) 0.0) (= x ($inexactnum-real-part y)))]
2225             [else (nonnumber-error who y)])]
2226         [(ratnum?)
2227          (type-case y
2228             [(fixnum? bignum? $exactnum?) #f]
2229             [(ratnum?)
2230              (and (= ($ratio-numerator x) ($ratio-numerator y))
2231                   (= ($ratio-denominator x) ($ratio-denominator y)))]
2232             [(flonum?) (exact-inexact-compare? = x y)]
2233             [($inexactnum?) (and (fl= ($inexactnum-imag-part y) 0.0) (= x ($inexactnum-real-part y)))]
2234             [else (nonnumber-error who y)])]
2235         [($exactnum? $inexactnum?)
2236          (unless (number? y) (nonnumber-error who y))
2237          (and (= (real-part x) (real-part y)) (= (imag-part x) (imag-part y)))]
2238         [(flonum?)
2239          (type-case y
2240             [(cflonum?) (cfl= x y)]
2241             [(fixnum?) (fl= x (fixnum->flonum y))]
2242             [(bignum? ratnum?) (exact-inexact-compare? = y x)]
2243             [($exactnum?) #f]
2244             [else (nonnumber-error who y)])]
2245         [else (nonnumber-error who x)])))
2246
2247(set! $<
2248   (lambda (who x y)
2249      (type-case x
2250         [(fixnum?)
2251          (type-case y
2252             [(fixnum?) (fx< x y)]
2253             [(bignum?) ($bigpositive? y)]
2254             [(ratnum?) (< (* ($ratio-denominator y) x) ($ratio-numerator y))]
2255             [(flonum?) (< (fixnum->flonum x) y)]
2256             [else (nonreal-error who y)])]
2257         [(bignum?)
2258          (type-case y
2259             [(bignum?) (big< x y)]
2260             [(fixnum?) (not ($bigpositive? x))]
2261             [(ratnum?) (< (* ($ratio-denominator y) x) ($ratio-numerator y))]
2262             [(flonum?) (exact-inexact-compare? < x y)]
2263             [else (nonreal-error who y)])]
2264         [(ratnum?)
2265          (type-case y
2266             [(fixnum? bignum?)
2267              (< ($ratio-numerator x) (* ($ratio-denominator x) y))]
2268             [(ratnum?)
2269              (< (* ($ratio-numerator x) ($ratio-denominator y))
2270                 (* ($ratio-numerator y) ($ratio-denominator x)))]
2271             [(flonum?) (exact-inexact-compare? < x y)]
2272             [else (nonreal-error who y)])]
2273         [(flonum?)
2274          (type-case y
2275             [(flonum?) (fl< x y)]
2276             [(fixnum?) (fl< x (fixnum->flonum y))]
2277             [(bignum? ratnum?) (exact-inexact-compare? > y x)]
2278             [else (nonreal-error who y)])]
2279         [else (nonreal-error who x)])))
2280
2281(set! $<=
2282   (lambda (who x y)
2283      (type-case x
2284         [(fixnum?)
2285          (type-case y
2286             [(fixnum?) (fx<= x y)]
2287             [(bignum?) ($bigpositive? y)]
2288             [(ratnum?)
2289              (<= (* ($ratio-denominator y) x) ($ratio-numerator y))]
2290             [(flonum?) (<= (fixnum->flonum x) y)]
2291             [else (nonreal-error who y)])]
2292         [(bignum?)
2293          (type-case y
2294             [(bignum?) (not (big< y x))]
2295             [(fixnum?) (not ($bigpositive? x))]
2296             [(ratnum?)
2297              (<= (* ($ratio-denominator y) x) ($ratio-numerator y))]
2298             [(flonum?) (exact-inexact-compare? <= x y)]
2299             [else (nonreal-error who y)])]
2300         [(ratnum?)
2301          (type-case y
2302             [(fixnum? bignum?)
2303              (<= ($ratio-numerator x) (* ($ratio-denominator x) y))]
2304             [(ratnum?)
2305              (<= (* ($ratio-numerator x) ($ratio-denominator y))
2306                 (* ($ratio-numerator y) ($ratio-denominator x)))]
2307             [(flonum?) (exact-inexact-compare? <= x y)]
2308             [else (nonreal-error who y)])]
2309         [(flonum?)
2310          (type-case y
2311             [(flonum?) (fl<= x y)]
2312             [(fixnum?) (fl<= x (fixnum->flonum y))]
2313             [(bignum? ratnum?) (exact-inexact-compare? >= y x)]
2314             [else (nonreal-error who y)])]
2315         [else (nonreal-error who x)])))
2316
2317(set! $+
2318  (lambda (who x y)
2319    (define (exint-unknown+ who x y)
2320      (type-case y
2321        [(fixnum? bignum?) (integer+ x y)]
2322        [(ratnum?)
2323         (let ([d ($ratio-denominator y)])
2324           (make-ratnum (+ (* x d) ($ratio-numerator y)) d))]
2325        [(flonum?) (exact-inexact+ x y)]
2326        [($exactnum? $inexactnum?)
2327         (make-rectangular (+ x (real-part y)) (imag-part y))]
2328        [else (nonnumber-error who y)]))
2329    (cond
2330      [(eqv? y 0) (unless (number? x) (nonnumber-error who x)) x]
2331      [else
2332        (type-case x
2333          [(fixnum?)
2334           (cond
2335             [(fx= x 0) (unless (number? y) (nonnumber-error who y)) y]
2336             [else (exint-unknown+ who x y)])]
2337          [(bignum?) (exint-unknown+ who x y)]
2338          [(ratnum?)
2339           (type-case y
2340             [(fixnum? bignum?)
2341              (let ([d ($ratio-denominator x)])
2342                (make-ratnum (+ (* y d) ($ratio-numerator x)) d))]
2343             [(ratnum?)
2344	      ;; adapted from Gambit, see gambit/lib/_num.scm
2345	      (let ((p ($ratio-numerator x))
2346		    (q ($ratio-denominator x))
2347		    (r ($ratio-numerator y))
2348		    (s ($ratio-denominator y)))
2349		(let ((d1 (exgcd q s)))
2350		  (if (eqv? d1 1)
2351		      (make-ratnum (+ (* p s)
2352				      (* r q))
2353				   (* q s))
2354		      (let* ((s-prime (intquotient s d1))
2355			     (t (+ (* p s-prime)
2356				   (* r (intquotient q d1))))
2357			     (d2 (exgcd d1 t))
2358			     (num (intquotient t d2))
2359			     (den (* (intquotient q d2)
2360				     s-prime)))
2361			(if (eqv? den 1)
2362			    num
2363			    (make-ratnum num den))))))]
2364             [($exactnum? $inexactnum?)
2365              (make-rectangular (+ x (real-part y)) (imag-part y))]
2366             [(flonum?) (exact-inexact+ x y)]
2367             [else (nonnumber-error who y)])]
2368          [(flonum?)
2369           (type-case y
2370             [(cflonum?) (cfl+ x y)]
2371             [(fixnum? bignum? ratnum?) (exact-inexact+ y x)]
2372             [($exactnum?)
2373              (make-rectangular (+ x (real-part y)) (imag-part y))]
2374             [else (nonnumber-error who y)])]
2375          [($exactnum? $inexactnum?)
2376           (type-case y
2377             [(fixnum? bignum? ratnum? flonum?)
2378              (make-rectangular (+ (real-part x) y) (imag-part x))]
2379             [($exactnum? $inexactnum?)
2380              (make-rectangular (+ (real-part x) (real-part y))
2381                (+ (imag-part x) (imag-part y)))]
2382             [else (nonnumber-error who y)])]
2383          [else (nonnumber-error who x)])])))
2384
2385(set! $*
2386  (let ([$bignum-trailing-zero-bits (foreign-procedure "(cs)s_big_trailing_zero_bits" (ptr) ptr)])
2387   (lambda (who x y)
2388    (cond
2389      [(and (fixnum? y) ($fxu< (#3%fx+ y 1) 3))
2390       (cond
2391         [(fx= y 0) (unless (number? x) (nonnumber-error who x)) 0]
2392         [(fx= y 1) (unless (number? x) (nonnumber-error who x)) x]
2393         [else ($negate who x)])]
2394      [else
2395       (type-case x
2396         [(fixnum? bignum?)
2397          (type-case y
2398             [(fixnum?) (integer* x y)]
2399            [(bignum?) (if (fixnum? x)
2400                           (cond
2401                            [($fxu< (#3%fx+ x 1) 3)
2402                             (cond
2403                              [(fx= x 0) (unless (number? y) (nonnumber-error who y)) 0]
2404                              [(fx= x 1) (unless (number? y) (nonnumber-error who y)) y]
2405                              [else ($negate who y)])]
2406                            [else (integer* x y)])
2407                           (let ([slim 32]
2408				 [klim 100]
2409				 [t3lim 512])
2410                             ; both of the following functions were adapted from
2411                             ; https://github.com/casevh/DecInt/blob/master/DecInt.py#L451
2412                             ; under the BSD license
2413                             (define (toom3 x y)
2414                               (define xl (if (bignum? x) ($bignum-length x) 0))
2415                               (define yl (if (bignum? y) ($bignum-length y) 0))
2416                               (cond
2417                                 [(and (fx< xl slim) (fx< yl slim))
2418                                  (integer* x y)]
2419                                 [(and (fx< xl klim) (fx< yl klim))
2420                                  (karatsuba x y)]
2421                                 [else
2422                                  (let* ([k (fx* (fxquotient (fxmax xl yl) 3) (constant bigit-bits))]
2423                                         [x-hi (ash x (fx* -2 k))]
2424                                         [y-hi (ash y (fx* -2 k))]
2425                                         [x-mid (bitwise-bit-field x k (fx* 2 k))]
2426                                         [y-mid (bitwise-bit-field y k (fx* 2 k))]
2427                                         [x-lo (bitwise-bit-field x 0 k)]
2428                                         [y-lo (bitwise-bit-field y 0 k)]
2429                                         [z0 (toom3 x-hi y-hi)]
2430                                         [z4 (toom3 x-lo y-lo)]
2431                                         [t1 (toom3 (+ x-hi x-mid x-lo) (+ y-hi y-mid y-lo))]
2432                                         [t2 (toom3 (+ (- x-hi x-mid) x-lo) (+ (- y-hi y-mid) y-lo))]
2433                                         [t3 (* (+ x-hi (ash x-mid 1) (ash x-lo 2))
2434                                                (+ y-hi (ash y-mid 1) (ash y-lo 2)))]
2435                                         [z2 (- (ash (+ t1 t2) -1) z0 z4)]
2436                                         [t4 (- t3 z0 (ash z2 2) (ash z4 4))]
2437                                         [z3 (quotient (+ (- t4  t1) t2) 6)]
2438                                         [z1 (- (ash (- t1 t2) -1) z3)])
2439                                    (+ (ash z0 (* k 4))
2440                                       (ash z1 (* k 3))
2441                                       (ash z2 (* k 2))
2442                                       (ash z3 (* k 1))
2443                                       (ash z4 (* k 0))))]))
2444
2445                             (define (toom4 x y)
2446                               (define xl (if (bignum? x) ($bignum-length x) 0))
2447                               (define yl (if (bignum? y) ($bignum-length y) 0))
2448                               (cond
2449                                 [(and (fx< xl slim) (fx< yl slim))
2450                                  (integer* x y)]
2451                                 [(and (fx< xl klim) (fx< yl klim))
2452                                  (karatsuba x y)]
2453                                 [(and (fx< xl t3lim) (fx< yl t3lim))
2454                                  (toom3 x y)]
2455                                 [else
2456                                  (let* ((k (fx* (fxquotient (fxmax xl yl) 4) (constant bigit-bits)))
2457                                         (x0 (ash x (fx* -3 k)))
2458                                         (y0 (ash y (fx* -3 k)))
2459                                         (x1 (bitwise-bit-field x (fx* 2 k) (fx* 3 k)))
2460                                         (y1 (bitwise-bit-field y (fx* 2 k) (fx* 3 k)))
2461                                         (x2 (bitwise-bit-field x (fx* 1 k) (fx* 2 k)))
2462                                         (y2 (bitwise-bit-field y (fx* 1 k) (fx* 2 k)))
2463                                         (x3 (bitwise-bit-field x 0 k))
2464                                         (y3 (bitwise-bit-field y 0 k))
2465                                         (z0 (toom4 x0 y0))
2466                                         (z6 (toom4 x3 y3))
2467                                         (t0 (+ z0 z6))
2468                                         (xeven (+ x0 x2))
2469                                         (xodd (+ x1 x3))
2470                                         (yeven (+ y0 y2))
2471                                         (yodd (+ y1 y3))
2472                                         (t1 (- (toom4 (+ xeven xodd) (+ yeven yodd)) t0))
2473                                         (t2 (- (toom4 (- xeven xodd) (- yeven yodd)) t0))
2474                                         (xeven (+ x0 (ash x2 2)))
2475                                         (xodd (+ (ash x1 1) (ash x3 3)))
2476                                         (yeven (+ y0 (ash y2 2)))
2477                                         (yodd (+ (ash y1 1) (ash y3 3)))
2478                                         (t0 (+ z0 (ash z6 6)))
2479                                         (t3 (- (toom4 (+ xeven xodd) (+ yeven yodd)) t0))
2480                                         (t4 (- (toom4 (- xeven xodd) (- yeven yodd)) t0))
2481                                         (t5 (- (* (+ x0 (* 3 x1) (* 9 x2) (* 27 x3))
2482                                                   (+ y0 (* 3 y1) (* 9 y2) (* 27 y3)))
2483                                                (+ z0 (* 729 z6))))
2484                                         (t6 (+ t1 t2))
2485                                         (t7 (+ t3 t4))
2486                                         (z4 (quotient (- t7 (ash t6 2)) 24))
2487                                         (z2 (- (ash t6 -1) z4))
2488                                         (t8 (- t1 z2 z4))
2489                                         (t9 (- t3 (ash z2 2) (ash z4 4)))
2490                                         (t10 (- t5 (* 9 z2) (* 81 z4)))
2491                                         (t11 (- t10 (* 3 t8)))
2492                                         (t12 (- t9 (ash t8 1)))
2493                                         (z5 (quotient (- t11 (ash t12 2)) 120))
2494                                         (z3 (quotient (- (ash t12 3) t11) 24))
2495                                         (z1 (- t8 z3 z5)))
2496                                    (+ (ash z0 (* k 6))
2497                                       (ash z1 (* k 5))
2498                                       (ash z2 (* k 4))
2499                                       (ash z3 (* k 3))
2500                                       (ash z4 (* k 2))
2501                                       (ash z5 (* k 1))
2502                                       (ash z6 (* k 0))))]))
2503
2504                              ;; _Modern Computer Arithmetic_, Brent and Zimmermann
2505                              (define (karatsuba x y)
2506                                (define xl (if (bignum? x) ($bignum-length x) 0))
2507                                (define yl (if (bignum? y) ($bignum-length y) 0))
2508                                (cond
2509                                 [(and (fx< xl 30) (fx< yl 30))
2510                                  (integer* x y)]
2511                                 [else
2512                                  (let* ([k (fx* (fxquotient (fxmax xl yl) 2) (constant bigit-bits))]
2513                                         [x-hi (ash x (fx- k))]
2514                                         [y-hi (ash y (fx- k))]
2515					 [x-lo (bitwise-bit-field x 0 k)]
2516					 [y-lo (bitwise-bit-field y 0 k)]
2517                                         [c0 (karatsuba x-lo y-lo)]
2518                                         [c1 (karatsuba x-hi y-hi)]
2519                                         [c1-c2 (cond
2520                                                 [(< x-lo x-hi)
2521                                                  (cond
2522                                                   [(< y-lo y-hi)
2523                                                    (- c1 (karatsuba (- x-hi x-lo) (- y-hi y-lo)))]
2524                                                   [else
2525                                                    (+ c1 (karatsuba (- x-hi x-lo) (- y-lo y-hi)))])]
2526                                                 [else
2527                                                  (cond
2528                                                   [(< y-lo y-hi)
2529                                                    (+ c1 (karatsuba (- x-lo x-hi) (- y-hi y-lo)))]
2530                                                   [else
2531                                                    (- c1 (karatsuba (- x-lo x-hi) (- y-lo y-hi)))])])])
2532                                    (+ c0 (integer-ash (+ c0 c1-c2) k) (integer-ash c1 (fx* 2 k))))]))
2533                              ;; Multiplying numbers with trailing 0s is common, so
2534                              ;; check for that case:
2535                              (let ([xz ($bignum-trailing-zero-bits x)]
2536                                    [yz (if (bignum? y) ($bignum-trailing-zero-bits y) 0)])
2537                                (let ([z (fx+ xz yz)])
2538                                  (if (fx= z 0)
2539                                      (toom4 x y)
2540                                      (bitwise-arithmetic-shift-left
2541                                       (toom4 (bitwise-arithmetic-shift-right x xz)
2542					      (bitwise-arithmetic-shift-right y yz))
2543                                       z))))))]
2544             [(ratnum?) (/ (* x ($ratio-numerator y)) ($ratio-denominator y))]
2545             [($exactnum? $inexactnum?)
2546              (make-rectangular (* x (real-part y)) (* x (imag-part y)))]
2547             [(flonum?) (exact-inexact* x y)]
2548             [else (nonnumber-error who y)])]
2549         [(ratnum?)
2550          (type-case y
2551             [(fixnum? bignum?)
2552              (integer/ (* y ($ratio-numerator x)) ($ratio-denominator x))]
2553             [(ratnum?)
2554	      ;; adapted from Gambit, see gambit/lib/_num.scm
2555	      (let ((p ($ratio-numerator x))
2556		    (q ($ratio-denominator x))
2557		    (r ($ratio-numerator y))
2558		    (s ($ratio-denominator y)))
2559		(if (eq? x y)
2560		    (make-ratnum (magnitude-squared p) (magnitude-squared q))     ;; already in lowest form
2561		    (let* ((gcd-ps (exgcd p s))
2562			   (gcd-rq (exgcd r q))
2563			   (num (* (intquotient p gcd-ps) (intquotient r gcd-rq)))
2564			   (den (* (intquotient q gcd-rq) (intquotient s gcd-ps))))
2565		      (if (eqv? den 1)
2566			  num
2567			  (make-ratnum num den)))))]
2568             [($exactnum? $inexactnum?)
2569              (make-rectangular (* x (real-part y)) (* x (imag-part y)))]
2570             [(flonum?) (exact-inexact* x y)]
2571             [else (nonnumber-error who y)])]
2572          [(flonum?)
2573           (type-case y
2574             [(cflonum?) (cfl* x y)]
2575             [(fixnum? bignum? ratnum?) (exact-inexact* y x)]
2576             [($exactnum?)
2577              (make-rectangular (* x (real-part y)) (* x (imag-part y)))]
2578             [else (nonnumber-error who y)])]
2579          [($exactnum? $inexactnum?)
2580           (type-case y
2581             [(fixnum? bignum? ratnum? flonum?)
2582              (make-rectangular (* (real-part x) y) (* (imag-part x) y))]
2583             [($exactnum? $inexactnum?)
2584              (let ([a (real-part x)] [b (imag-part x)]
2585                    [c (real-part y)] [d (imag-part y)])
2586                (make-rectangular (- (* a c) (* b d)) (+ (* a d) (* b c))))]
2587             [else (nonnumber-error who y)])]
2588          [else (nonnumber-error who x)])]))))
2589
2590(set! $-
2591  (lambda (who x y)
2592    (define (exint-unknown- who x y)
2593      (type-case y
2594        [(fixnum? bignum?) (integer- x y)]
2595        [(ratnum?)
2596         (let ([d ($ratio-denominator y)])
2597           (make-ratnum (- (* x d) ($ratio-numerator y)) d))]
2598        [($exactnum? $inexactnum?)
2599         (make-rectangular (- x (real-part y)) (- (imag-part y)))]
2600        [(flonum?) (exact-inexact- x y)]
2601        [else (nonnumber-error who y)]))
2602    (cond
2603      [(eqv? y 0) (unless (number? x) (nonnumber-error who x)) x]
2604      [else
2605        (type-case x
2606          [(fixnum?)
2607           (cond
2608             [(eqv? x 0) ($negate who y)]
2609             [else (exint-unknown- who x y)])]
2610          [(bignum?) (exint-unknown- who x y)]
2611          [(ratnum?)
2612           (type-case y
2613             [(fixnum? bignum?)
2614              (let ([d ($ratio-denominator x)])
2615                (make-ratnum (- ($ratio-numerator x) (* y d)) d))]
2616             [(ratnum?)
2617	      ;; adapted from Gambit, see gambit/lib/_num.scm
2618	      (let ((p ($ratio-numerator x))
2619		    (q ($ratio-denominator x))
2620		    (r ($ratio-numerator y))
2621		    (s ($ratio-denominator y)))
2622		(let ((d1 (gcd q s)))
2623		  (if (eqv? d1 1)
2624		      (make-ratnum (- (* p s)
2625				      (* r q))
2626				   (* q s))
2627		      (let* ((s-prime (intquotient s d1))
2628			     (t (- (* p s-prime)
2629				   (* r (intquotient q d1))))
2630			     (d2 (exgcd d1 t))
2631			     (num (intquotient t d2))
2632			     (den (* (intquotient q d2)
2633				     s-prime)))
2634			(if (eqv? den 1)
2635			    num
2636			    (make-ratnum num den))))))]
2637             [($exactnum? $inexactnum?)
2638              (make-rectangular (- x (real-part y)) (- (imag-part y)))]
2639             [(flonum?) (exact-inexact- x y)]
2640             [else (nonnumber-error who y)])]
2641          [(flonum?)
2642           (type-case y
2643             [(cflonum?) (cfl- x y)]
2644             [(fixnum? bignum? ratnum?) (inexact-exact- x y)]
2645             [($exactnum?)
2646              (make-rectangular (- x (real-part y)) (- (imag-part y)))]
2647             [else (nonnumber-error who y)])]
2648          [($exactnum? $inexactnum?)
2649           (type-case y
2650             [(fixnum? bignum? ratnum? flonum?)
2651              (make-rectangular (- (real-part x) y) (imag-part x))]
2652             [($exactnum? $inexactnum?)
2653              (make-rectangular (- (real-part x) (real-part y))
2654                (- (imag-part x) (imag-part y)))]
2655             [else (nonnumber-error who y)])]
2656          [else (nonnumber-error who x)])])))
2657
2658(set! $/
2659   (lambda (who x y)
2660      (type-case y
2661         [(fixnum?)
2662          (cond
2663           [(fx= y 1) (unless (number? x) (nonnumber-error who x)) x]
2664           [(fx= y -1) (unless (number? x) (nonnumber-error who x)) ($negate who x)]
2665           [else
2666            (type-case x
2667             [(fixnum?)
2668              ;; Trying `fxquotient` followed by a `fx*` check
2669              ;; is so much faster (in the case that it works)
2670              ;; that it's worth a try
2671              (when (eq? y 0) (domain-error who y))
2672              (let ([q (fxquotient x y)])
2673                (if (fx= x (fx* y q))
2674                    q
2675                    (integer/ x y)))]
2676             [(bignum?)
2677              (when (eq? y 0) (domain-error who y))
2678              (integer/ x y)]
2679             [(ratnum?)
2680              (when (eq? y 0) (domain-error who y))
2681              (/ ($ratio-numerator x) (* y ($ratio-denominator x)))]
2682             [($exactnum?)
2683              (when (eq? y 0) (domain-error who y))
2684              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2685             [($inexactnum?)
2686              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2687             [(flonum?) (inexact-exact/ x y)]
2688             [else (nonnumber-error who x)])])]
2689         [(bignum?)
2690          (type-case x
2691             [(fixnum? bignum?)
2692              (integer/ x y)]
2693             [(ratnum?)
2694              (/ ($ratio-numerator x) (* y ($ratio-denominator x)))]
2695             [($exactnum?)
2696              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2697             [($inexactnum?)
2698              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2699             [(flonum?) (inexact-exact/ x y)]
2700             [else (nonnumber-error who x)])]
2701         [(ratnum?)
2702          (type-case x
2703	     [(fixnum? bignum?)
2704              (cond
2705                [(eq? x 1) (if (negative? ($ratio-numerator y))
2706                               (make-ratnum ($negate who ($ratio-denominator y)) ($negate who ($ratio-numerator y)))
2707                               (make-ratnum ($ratio-denominator y) ($ratio-numerator y)))]
2708                [(eq? x -1) (if (negative? ($ratio-numerator y))
2709                                (make-ratnum ($ratio-denominator y) ($negate who ($ratio-numerator y)))
2710                                (make-ratnum ($negate who ($ratio-denominator y)) ($ratio-numerator y)))]
2711                [else
2712                 (integer/ (* x ($ratio-denominator y)) ($ratio-numerator y))])]
2713             [(ratnum?)
2714	      ;; adapted from Gambit, see gambit/lib/_num.scm
2715	      (let ((p ($ratio-numerator x))
2716		    (q ($ratio-denominator x))
2717		    (r ($ratio-denominator y))
2718		    (s ($ratio-numerator y)))
2719		(if (eq? x y)
2720		    1
2721		    (let* ((gcd-ps (exgcd p s))
2722			   (gcd-rq (exgcd r q))
2723			   (num (* (intquotient p gcd-ps) (intquotient r gcd-rq)))
2724			   (den (* (intquotient q gcd-rq) (intquotient s gcd-ps))))
2725		      (if (negative? den)
2726			  (if (eqv? den -1)
2727			      (- num)
2728			      (make-ratnum (- num) (- den)))
2729			  (if (eqv? den 1)
2730			      num
2731			      (make-ratnum num den))))))]
2732             [($exactnum? $inexactnum?)
2733              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2734             [(flonum?) (inexact-exact/ x y)]
2735             [else (nonnumber-error who x)])]
2736         [(flonum?)
2737          (type-case x
2738             [(cflonum?) (cfl/ x y)]
2739             [(fixnum? bignum? ratnum?) (exact-inexact/ x y)]
2740             [($exactnum?)
2741              (make-rectangular (/ (real-part x) y) (/ (imag-part x) y))]
2742             [else (nonnumber-error who x)])]
2743         [($exactnum? $inexactnum?)
2744          (type-case x
2745             [(fixnum? bignum? ratnum? flonum?)
2746              ;; a / c+di => c(a/(cc+dd)) + (-d(a/cc+dd))i
2747              (if (eq? x 0)
2748                  0
2749                  (let ([c (real-part y)] [d (imag-part y)])
2750                    (let ([t (/ x (+ (* c c) (* d d)))])
2751                      (make-rectangular (* c t) (- (* d t))))))]
2752             [($exactnum? $inexactnum?)
2753              (let ([a (real-part x)] [b (imag-part x)]
2754                    [c (real-part y)] [d (imag-part y)])
2755                ;; a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
2756                (define (simpler-divide a b c d)
2757                  ;; Direct calculuation does not work as well for complex numbers with
2758                  ;; large parts, such as `(/ 1e+300+1e+300i 4e+300+4e+300i)`, but it
2759                  ;; works better for small parts, as in `(/ 0.0+0.0i 1+1e-320i)`
2760                  (let ([t (+ (* c c) (* d d))])
2761                    (make-rectangular (/ (+ (* a c) (* b d)) t)
2762                                      (/ (- (* b c) (* a d)) t))))
2763                ;; Let r = c/d or d/c, depending on which is larger
2764                (cond
2765                 [(or (eq? c 0) (and ($exactnum? x) ($exactnum? y)))
2766                  (simpler-divide a b c d)]
2767                 [(< (abs c) (abs d))
2768                  (let ([r (/ d c)])
2769                    (if (infinite? r)
2770                        ;; Too large; try form that works better with small c or d
2771                        (simpler-divide a b c d)
2772                        ;; a+bi / c+di =>
2773                        (let ([x (+ c (* d r))]) ; x = c+dd/c = (cc+dd)/c
2774                          ;; (a+br)/x + ((b-ar)/x)i = (a+bd/c)c/(cc+dd) + ((b-ad/c)c/(cc+dd))i
2775                          ;; = (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
2776                          (make-rectangular (/ (+ a (* b r)) x)
2777                                            (/ (- b (* a r)) x)))))]
2778                 [else
2779                  (let ([r (/ c d)])
2780                    (if (infinite? r)
2781                        ;; Too large; try form that works better with small c or d
2782                        (simpler-divide a b c d)
2783                        (let ([x (+ d (* c r))]) ; x = d+cc/d = (cc+dd)/d
2784                          ;; (b+ar)/x + ((br-a)/x)i = (b+ac/d)d/(cc+dd) + ((bc/d-a)d/(cc+dd))i
2785                          ;; = (bd+ac)/(cc+dd) + ((bc-ad)/(cc+dd))i
2786                          (make-rectangular (/ (+ b (* a r)) x)
2787                                            (/ (- (* b r) a) x)))))]))]
2788             [else (nonnumber-error who x)])]
2789         [else (nonnumber-error who y)])))
2790
2791(set! conjugate
2792   (lambda (x)
2793      (type-case x
2794         [(flonum? fixnum? ratnum? bignum?) x]
2795         [($inexactnum?)
2796          (fl-make-rectangular ($inexactnum-real-part x)
2797                               (fl- ($inexactnum-imag-part x)))]
2798         [($exactnum?)
2799          ($make-exactnum ($exactnum-real-part x)
2800                           (- ($exactnum-imag-part x)))]
2801         [else (nonnumber-error 'conjugate x)])))
2802
2803(set! magnitude-squared
2804   (lambda (x)
2805      (type-case x
2806         [(flonum?) (fl* x x)]
2807         [($inexactnum?)
2808          (let ([a ($inexactnum-real-part x)] [b ($inexactnum-imag-part x)])
2809             (fl+ (fl* a a) (fl* b b)))]
2810         [(fixnum? ratnum? bignum?) (* x x)]
2811         [($exactnum?)
2812          (let ([a ($exactnum-real-part x)] [b ($exactnum-imag-part x)])
2813             (+ (* a a) (* b b)))]
2814         [else (nonnumber-error 'magnitude-squared x)])))
2815
2816(set! cfl-magnitude-squared
2817   (lambda (x)
2818      (type-case x
2819         [(flonum?) (fl* x x)]
2820         [($inexactnum?)
2821          (let ([a ($inexactnum-real-part x)] [b ($inexactnum-imag-part x)])
2822             (fl+ (fl* a a) (fl* b b)))]
2823         [else (noncflonum-error 'cfl-magnitude-squared x)])))
2824
2825(set! zero?
2826   (lambda (z)
2827      (type-case z
2828         [(fixnum?) (fx= z 0)]
2829         [(flonum?) (fl= z 0.0)]
2830         [($inexactnum?) (cfl= z 0.0)]
2831         [(bignum? ratnum? $exactnum?) #f]
2832         [else (nonnumber-error 'zero? z)])))
2833
2834(set-who! nan?
2835  (lambda (x)
2836    (type-case x
2837      [(flonum?) ($nan? x)]
2838      [(fixnum? bignum? ratnum?) #f]
2839      [else (nonreal-error who x)])))
2840
2841(set-who! infinite?
2842  (lambda (x)
2843    (type-case x
2844      [(flonum?) (infinity? x)]
2845      [(fixnum? bignum? ratnum?) #f]
2846      [else (nonreal-error who x)])))
2847
2848(set-who! finite?
2849  (lambda (x)
2850    (type-case x
2851      [(flonum?) (not (exceptional-flonum? x))]
2852      [(fixnum? bignum? ratnum?) #t]
2853      [else (nonreal-error who x)])))
2854
2855(let ()
2856  (define $ash
2857    (lambda (who x n)
2858      (type-case n
2859        [(fixnum?)
2860         (type-case x
2861           [(fixnum?)
2862            (let ([max-fx-shift (- (constant fixnum-bits) 1)])
2863              (if (fx< n 0)
2864                 ; can't just go for it since (- n) may not be representable
2865                  (if (fx< n (- max-fx-shift))
2866                      (fxsra x max-fx-shift)
2867                      (fxsra x (fx- n)))
2868                  (if (fx> n max-fx-shift)
2869                      (integer-ash x n)
2870                      (let ([m (#3%fxsll x n)])
2871                        (if (fx= (fxsra m n) x)
2872                            m
2873                            (integer-ash x n))))))]
2874           [(bignum?) (integer-ash x n)]
2875           [else (nonexact-integer-error who x)])]
2876        [(bignum?)
2877         (type-case x
2878           [(fixnum?)
2879            (cond
2880             [(fx= x 0) 0]
2881             [($bigpositive? n) ($oops who "out of memory")]
2882             [(fxpositive? x) 0]
2883             [else -1])]
2884           [(bignum?)
2885            (cond
2886             [($bigpositive? n) ($oops who "out of memory")]
2887             [($bigpositive? x) 0]
2888             [else -1])]
2889           [else (nonexact-integer-error who x)])]
2890        [else (nonexact-integer-error who n)])))
2891
2892  (set-who! ash (lambda (x n) ($ash who x n)))
2893
2894  (set-who! bitwise-arithmetic-shift (lambda (x n) ($ash who x n))))
2895
2896(set-who! bitwise-arithmetic-shift-left (lambda (x n) ($sll who x n)))
2897
2898(set-who! bitwise-arithmetic-shift-right (lambda (x n) ($sra who x n)))
2899
2900(set-who! integer-length
2901  (lambda (x)
2902    (type-case x
2903      [(fixnum?) (fxlength x)]
2904      [(bignum?) (biglength x)]
2905      [else (nonexact-integer-error who x)])))
2906
2907(set-who! bitwise-length ; same as integer-length
2908  (lambda (x)
2909    (type-case x
2910      [(fixnum?) (fxlength x)]
2911      [(bignum?) (biglength x)]
2912      [else (nonexact-integer-error who x)])))
2913
2914(set-who! bitwise-if
2915  (lambda (x y z)
2916    (define big-if
2917      (lambda (ei1 ei2 ei3)
2918        (bitwise-ior (bitwise-and ei1 ei2)
2919                     (bitwise-and (bitwise-not ei1) ei3))))
2920    (type-case x
2921      [(fixnum?)
2922       (type-case y
2923         [(fixnum?)
2924          (type-case z
2925            [(fixnum?) (fxif x y z)]
2926            [(bignum?) (big-if x y z)]
2927            [else (nonexact-integer-error who z)])]
2928         [(bignum?)
2929          (type-case z
2930            [(fixnum? bignum?) (big-if x y z)]
2931            [else (nonexact-integer-error who z)])]
2932         [else (nonexact-integer-error who y)])]
2933      [(bignum?)
2934       (type-case y
2935         [(fixnum? bignum?)
2936          (type-case z
2937            [(fixnum? bignum?) (big-if x y z)]
2938            [else (nonexact-integer-error who z)])]
2939         [else (nonexact-integer-error who y)])]
2940      [else (nonexact-integer-error who x)])))
2941
2942(set-who! bitwise-copy-bit
2943  (lambda (x y b)
2944    (unless (and (integer? x) (exact? x))
2945      ($oops who "~s is not an exact integer" x))
2946    (unless (or (and (fixnum? y) (fxnonnegative? y))
2947                (and (bignum? y) ($bigpositive? y)))
2948      ($oops who "~s is not a nonnegative exact integer" y))
2949    (cond
2950      [(eq? b 0) (logbit0 y x)]
2951      [(eq? b 1) (logbit1 y x)]
2952      [else ($oops who "bit argument ~s is not 0 or 1" b)])))
2953
2954(let ()
2955  (define count-table
2956    (let ()
2957      (define-syntax make-count-table
2958        (lambda (x)
2959          #`'#,(let ([bv (make-bytevector 256)])
2960                 (define slow-bit-count
2961                   (lambda (x)
2962                     (do ([x x (fxsrl x 1)] [cnt 0 (if (fxodd? x) (fx+ cnt 1) cnt)])
2963                         ((fx= x 0) cnt))))
2964                 (do ([i 0 (fx+ i 1)])
2965                     ((fx= i 256))
2966                   (bytevector-u8-set! bv i (slow-bit-count i)))
2967                 bv)))
2968      (make-count-table)))
2969  (define $fxbit-count
2970    (lambda (n)
2971      (if (fx< n 0)
2972          (fxnot (fxpopcount (fxnot n)))
2973          (fxpopcount n))))
2974  (define $big-bit-count
2975    (lambda (n)
2976      (let ([end (fx+ (fx* ($bignum-length n) (constant bigit-bytes)) (constant bignum-data-disp))])
2977        (do ([i (constant bignum-data-disp) (fx+ i 1)]
2978             [cnt 0 (+ cnt (bytevector-u8-ref count-table ($object-ref 'unsigned-8 n i)))])
2979            ((fx= i end) cnt)))))
2980  (set-who! fxbit-count
2981    (lambda (n)
2982      (unless (fixnum? n) ($oops who "~s is not a fixnum" n))
2983      ($fxbit-count n)))
2984  (set-who! bitwise-bit-count
2985    (lambda (n)
2986      (cond
2987        [(fixnum? n)
2988         (if (fx< n 0)
2989             (fxnot ($fxbit-count (fxnot n)))
2990             ($fxbit-count n))]
2991        [(bignum? n)
2992         (if ($bigpositive? n)
2993             ($big-bit-count n)
2994             (fxnot ($big-bit-count (bitwise-not n))))]
2995        [else ($oops who "~s is not an exact integer" n)]))))
2996
2997(set-who! bitwise-first-bit-set
2998  (let ()
2999    (define $big-first-bit-set
3000      (foreign-procedure "(cs)s_big_first_bit_set" (ptr) ptr))
3001    (lambda (n)
3002      (cond
3003        [(fixnum? n) (fxfirst-bit-set n)]
3004        [(bignum? n) ($big-first-bit-set n)]
3005        [else ($oops who "~s is not an exact integer" n)]))))
3006
3007(set-who! bitwise-bit-field
3008  (let ()
3009   ; big-positive-bit-field assumes n is a positive bignum, start and
3010   ; end are nonnegative fixnums, and end > start
3011    (define big-positive-bit-field
3012      (foreign-procedure "(cs)s_big_positive_bit_field" (ptr ptr ptr) ptr))
3013    (define (generic-bit-field n start end)
3014      (bitwise-and
3015        ($sra who n start)
3016        (- ($sll who 1 (- end start)) 1)))
3017    (lambda (n start end)
3018      (unless (or (fixnum? n) (bignum? n))
3019        ($oops who "~s is not an exact integer" n))
3020      (cond
3021        [(and (fixnum? start) (fixnum? end))
3022         (unless (fx>= start 0) ($oops who "invalid start index ~s" start))
3023         (unless (fx>= end start) ($oops who "invalid end index ~s" end))
3024         (cond
3025           [(fx= end start) 0]
3026           [(and (fixnum? n) (fx< end (fx- (fixnum-width) 1)))
3027            (fxsra (fxand n (fxnot (fxsll -1 end))) start)]
3028           [(and (bignum? n) ($bigpositive? n))
3029            (big-positive-bit-field n start end)]
3030           [else (generic-bit-field n start end)])]
3031        [else
3032         (unless (or (and (fixnum? start) (fx>= start 0))
3033                     (and (bignum? start) ($bigpositive? start)))
3034           ($oops who "invalid start index ~s" start))
3035         (unless (or (and (fixnum? end) (>= end start))
3036                     (and (bignum? end) (>= end start)))
3037           ($oops who "invalid end index ~s" end))
3038         (generic-bit-field n start end)]))))
3039
3040(set-who! exact-integer-sqrt
3041  (lambda (n)
3042    (define (big-integer-sqrt n)
3043     ; adapted from SRFI 77 mail-archive posting by Brad Lucier, who derived
3044     ; it from "Karatsuba Square Root" by Paul Zimmermann, INRIA technical report
3045     ; RR-3805, 1999.
3046      (if (and (fixnum? n) (or (not (fixnum? (expt 2 52))) (< n (expt 2 52))))
3047          (let ([q (flonum->fixnum (flsqrt (fixnum->flonum n)))])
3048            (values q (fx- n (fx* q q))))
3049          (let ([b ($sra who (+ (integer-length n) 1) 2)])
3050            (let-values ([(s^ r^) (big-integer-sqrt ($sra who n (+ b b)))])
3051              (let* ([q&u (intquotient-remainder
3052                            (+ ($sll who r^ b)
3053                               (bitwise-bit-field n b (+ b b)))
3054                            ($sll who s^ 1))]
3055                     [q (car q&u)]
3056                     [u (cdr q&u)])
3057                (let ([s (+ ($sll who s^ b) q)]
3058                      [r (- (+ ($sll who u b)
3059                               (bitwise-bit-field n 0 b))
3060                            (* q q))])
3061                  (if (negative? r)
3062                      (values
3063                        (- s 1)
3064                        (+ r (- ($sll who s 1) 1)))
3065                      (values s r))))))))
3066    (cond
3067      [(and (fixnum? n) (fx>= n 0))
3068       (if (or (not (fixnum? (expt 2 52)))
3069               (fx< n (expt 2 52)))
3070           (let ([q (flonum->fixnum (flsqrt (fixnum->flonum n)))])
3071             (values q (fx- n (fx* q q))))
3072           (big-integer-sqrt n))]
3073      [(and (bignum? n) (#%$bigpositive? n)) (big-integer-sqrt n)]
3074      [else ($oops who "~s is not a nonnegative exact integer" n)])))
3075
3076(set-who! $quotient-remainder
3077  (lambda (x y)
3078    (type-case y
3079      [(fixnum? bignum?)
3080       (when (eq? y 0) (domain-error who y))
3081       (type-case x
3082         [(fixnum? bignum?) (intquotient-remainder x y)]
3083         [else (nonexact-integer-error who x)])]
3084      [else (nonexact-integer-error who y)])))
3085
3086(let ()
3087  (define-record pseudo-random-generator
3088    ((mutable double x10)
3089     (mutable double x11)
3090     (mutable double x12)
3091     (mutable double x20)
3092     (mutable double x21)
3093     (mutable double x22))
3094    ()
3095    ((constructor create-pseudo-random-generator)
3096     (predicate is-pseudo-random-generator?)))
3097
3098  (set! pseudo-random-generator?
3099        (lambda (x) (is-pseudo-random-generator? x)))
3100
3101  (let ([init! (foreign-procedure "(cs)s_random_state_init" (scheme-object unsigned) void)])
3102    (set! make-pseudo-random-generator
3103          (lambda ()
3104            (let ([s (create-pseudo-random-generator 0.0 0.0 0.0 0.0 0.0 0.0)]
3105                  [t (current-time 'time-utc)])
3106              (init! s (bitwise-and (+ (time-second t) (time-nanosecond t))
3107                                    #xFFFFFFFF))
3108              s)))
3109    (set-who! pseudo-random-generator-seed!
3110      (lambda (s n)
3111        (unless (is-pseudo-random-generator? s) ($oops who "not a pseudo-random generator ~s" s))
3112        (unless (or (and (fixnum? n) (fx>= n 0))
3113                    (and (bignum? n)  ($bigpositive? n)))
3114          ($oops who "not a nonnegative exact integer ~s" n))
3115        (init! s (bitwise-and n #xFFFFFFFF)))))
3116
3117  (set-who! pseudo-random-generator-next!
3118     (let ([random-double (foreign-procedure "(cs)s_random_state_next_double"
3119                            (scheme-object) double)]
3120           [random-int (foreign-procedure "(cs)s_random_state_next_integer"
3121                            (scheme-object uptr) uptr)])
3122       (case-lambda
3123        [(s)
3124         (unless (is-pseudo-random-generator? s) ($oops who "not a pseudo-random generator ~s" s))
3125         (random-double s)]
3126        [(s x)
3127         (define (random-integer s x)
3128           (let ([bits (integer-length x)])
3129             (let loop ([shift 0])
3130               (cond
3131                 [(<= bits shift) 0]
3132                 [else
3133                  ;; Assuming that a `uptr` is at least 32 bits:
3134                  (bitwise-ior (loop (+ shift 32))
3135                               (let ([n (bitwise-bit-field x shift (+ shift 32))])
3136                                 (if (zero? n)
3137                                     0
3138                                     (bitwise-arithmetic-shift-left
3139                                      (random-int s n)
3140                                      shift))))]))))
3141         (unless (is-pseudo-random-generator? s) ($oops who "not a pseudo-random generator ~s" s))
3142         (cond
3143          [(fixnum? x)
3144           (unless (fxpositive? x) ($oops who "not a positive exact integer ~s" x))
3145           (meta-cond
3146            [(<= (constant most-negative-fixnum) 4294967087 (constant most-positive-fixnum))
3147             (if (fx<= x 4294967087)
3148                 (random-int s x)
3149                 (random-integer s x))]
3150            [else
3151             (random-int s x)])]
3152          [(bignum? x)
3153           (unless ($bigpositive? x) ($oops who "not a positive exact integer ~s" x))
3154           (random-integer s x)]
3155          [else
3156           ($oops who "not a positive exact integer ~s" x)])])))
3157
3158  (set-who! pseudo-random-generator->vector
3159    (lambda (s)
3160      (unless (is-pseudo-random-generator? s) ($oops who "not a pseudo-random generator ~s" s))
3161      (vector (inexact->exact (pseudo-random-generator-x10 s))
3162              (inexact->exact (pseudo-random-generator-x11 s))
3163              (inexact->exact (pseudo-random-generator-x12 s))
3164              (inexact->exact (pseudo-random-generator-x20 s))
3165              (inexact->exact (pseudo-random-generator-x21 s))
3166              (inexact->exact (pseudo-random-generator-x22 s)))))
3167
3168  (let ([vector->prgen
3169         (let ([ok? (foreign-procedure "(cs)s_random_state_check" (double double double double double double) boolean)])
3170           (lambda (who s v)
3171             (define (bad-vector)
3172               ($oops who "not a valid pseudo-random generator state vector ~s" v))
3173             (define (int->double i)
3174               (unless (and (exact? i) (integer? i)) (bad-vector))
3175               (exact->inexact i))
3176             (unless (and (vector? v) (= 6 (vector-length v))) (bad-vector))
3177             (let ([x10 (int->double (vector-ref v 0))]
3178                   [x11 (int->double (vector-ref v 1))]
3179                   [x12 (int->double (vector-ref v 2))]
3180                   [x20 (int->double (vector-ref v 3))]
3181                   [x21 (int->double (vector-ref v 4))]
3182                   [x22 (int->double (vector-ref v 5))])
3183               (unless (ok? x10 x11 x12 x20 x21 x22) (bad-vector))
3184               (cond
3185                [s
3186                 (set-pseudo-random-generator-x10! s x10)
3187                 (set-pseudo-random-generator-x11! s x11)
3188                 (set-pseudo-random-generator-x12! s x12)
3189                 (set-pseudo-random-generator-x20! s x20)
3190                 (set-pseudo-random-generator-x21! s x21)
3191                 (set-pseudo-random-generator-x22! s x22)]
3192                [else
3193                 (create-pseudo-random-generator x10 x11 x12 x20 x21 x22)]))))])
3194
3195    (set-who! vector->pseudo-random-generator
3196      (lambda (vec) (vector->prgen who #f vec)))
3197    (set-who! vector->pseudo-random-generator!
3198      (lambda (s vec)
3199        (unless (is-pseudo-random-generator? s) ($oops who "not a pseudo-random generator ~s" s))
3200        (vector->prgen who s vec)))))
3201
3202(set! random
3203   (let ([fxrandom (foreign-procedure "(cs)s_fxrandom"
3204                      (scheme-object) scheme-object)]
3205         [flrandom (foreign-procedure "(cs)s_flrandom"
3206                      (scheme-object) scheme-object)])
3207      (lambda (x)
3208         (cond
3209            [(and (fixnum? x) (fx> x 0)) (fxrandom x)]
3210            [(and (flonum? x) (fl> x 0.0)) (flrandom x)]
3211            [(and (bignum? x) (> x 0))
3212             (let ([radix (most-positive-fixnum)])
3213                (do ([i x (quotient i radix)]
3214                     [a (fxrandom radix) (+ (* a radix) (fxrandom radix))])
3215                    ((<= i radix) (remainder a x))))]
3216            [else ($oops 'random "invalid argument ~s" x)]))))
3217
3218(set! random-seed ; must follow \#-
3219   (let ([limit #xFFFFFFFF]
3220         [get-seed (foreign-procedure "(cs)s_random_seed"
3221                      () unsigned-32)]
3222         [set-seed (foreign-procedure "(cs)s_set_random_seed"
3223                      (unsigned-32) void)])
3224      (case-lambda
3225         [() (get-seed)]
3226         [(n)
3227          (unless (and (or (fixnum? n) (bignum? n)) (<= 1 n limit))
3228             ($oops 'random-seed "invalid argument ~s" n))
3229          (set-seed n)])))
3230
3231(let ()
3232  (define-syntax fl-op
3233    (syntax-rules ()
3234      [(_ name $name x ...)
3235       (set-who! name
3236         (lambda (x ...)
3237           (unless (flonum? x) ($oops who "~s is not a flonum" x))
3238           ...
3239           ($name x ...)))]))
3240  (fl-op flexp $flexp x)
3241  (fl-op flsin $flsin x)
3242  (fl-op flcos $flcos x)
3243  (fl-op fltan $fltan x)
3244  (fl-op flasin $flasin x)
3245  (fl-op flacos $flacos x)
3246  (fl-op flsqrt $flsqrt x)
3247  (fl-op flexpt $flexpt x y)
3248  (fl-op flfloor $flfloor x)
3249  (fl-op flceiling $flceiling x))
3250
3251(set-who! flinteger?
3252  (lambda (x)
3253    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3254    ($flinteger? x)))
3255
3256(set-who! fllog
3257  (rec fllog
3258    (case-lambda
3259      [(x)
3260       (unless (flonum? x) ($oops who "~s is not a flonum" x))
3261       ($fllog x)]
3262      [(x y)
3263       (unless (flonum? x) ($oops who "~s is not a flonum" x))
3264       (unless (flonum? y) ($oops who "~s is not a flonum" y))
3265       (/ ($fllog x) ($fllog y))])))
3266
3267(set-who! flatan
3268  (rec flatan
3269    (case-lambda
3270      [(x)
3271       (unless (flonum? x) ($oops who "~s is not a flonum" x))
3272       ($flatan x)]
3273      [(x y)
3274       (unless (flonum? x) ($oops who "~s is not a flonum" x))
3275       (unless (flonum? y) ($oops who "~s is not a flonum" y))
3276       (flatan2 x y)])))
3277
3278(set-who! fltruncate
3279  (lambda (x)
3280    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3281    (if (negated-flonum? x) (fl- ($flfloor (flabs x))) ($flfloor x))))
3282
3283(set-who! flnan?
3284  (lambda (x)
3285    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3286    ($nan? x)))
3287
3288(set-who! flinfinite?
3289  (lambda (x)
3290    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3291    (infinity? x)))
3292
3293(set-who! flfinite?
3294  (lambda (x)
3295    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3296    (not (exceptional-flonum? x))))
3297
3298(set-who! flzero?
3299  (lambda (x)
3300    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3301    (fl= x 0.0)))
3302
3303(set-who! flpositive?
3304  (lambda (x)
3305    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3306    (fl> x 0.0)))
3307
3308(set-who! flnegative?
3309  (lambda (x)
3310    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3311    (fl< x 0.0)))
3312
3313(set-who! flnonpositive?
3314  (lambda (x)
3315    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3316    (fl<= x 0.0)))
3317
3318(set-who! flnonnegative?
3319  (lambda (x)
3320    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3321    (fl>= x 0.0)))
3322
3323(set-who! fleven?
3324  (lambda (x)
3325    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3326    (when (exceptional-flonum? x) (noninteger-error who x))
3327    (let ([y (fl* ($flfloor (fl/ x 2.0)) 2.0)])
3328      (cond
3329        [(fl= x y) #t]
3330        [(fl= (fl+ y 1.0) x) #f]
3331        [else (noninteger-error who x)]))))
3332
3333(set-who! flodd?
3334  (lambda (x)
3335    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3336    (when (exceptional-flonum? x) (noninteger-error who x))
3337    (let ([y (fl* ($flfloor (fl/ x 2.0)) 2.0)])
3338      (cond
3339        [(fl= x y) #f]
3340        [(fl= (fl+ y 1.0) x) #t]
3341        [else (noninteger-error who x)]))))
3342
3343(set-who! flmin
3344  (let ([$flmin (lambda (x y) (if (or (fl< x y) ($nan? x)) x y))])
3345    (case-lambda
3346       [(x y)
3347        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3348        (unless (flonum? y) ($oops who "~s is not a flonum" y))
3349        ($flmin x y)]
3350       [(x)
3351        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3352        x]
3353       [(x y . r)
3354        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3355        (unless (flonum? y) ($oops who "~s is not a flonum" y))
3356        (let loop ([x ($flmin x y)] [r r])
3357          (if (null? r)
3358              x
3359              (let ([y (car r)])
3360                (unless (flonum? y) ($oops who "~s is not a flonum" y))
3361                (loop ($flmin x y) (cdr r)))))])))
3362
3363(set-who! flmax
3364  (let ([$flmax (lambda (x y) (if (or (fl> x y) ($nan? x)) x y))])
3365    (case-lambda
3366       [(x y)
3367        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3368        (unless (flonum? y) ($oops who "~s is not a flonum" y))
3369        ($flmax x y)]
3370       [(x)
3371        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3372        x]
3373       [(x y . r)
3374        (unless (flonum? x) ($oops who "~s is not a flonum" x))
3375        (unless (flonum? y) ($oops who "~s is not a flonum" y))
3376        (let loop ([x ($flmax x y)] [r r])
3377          (if (null? r)
3378              x
3379              (let ([y (car r)])
3380                (unless (flonum? y) ($oops who "~s is not a flonum" y))
3381                (loop ($flmax x y) (cdr r)))))])))
3382
3383(set-who! flnumerator
3384  (lambda (x)
3385    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3386    (cond
3387      [($flinteger-or-inf? x) x]
3388      [($nan? x) x]
3389      [else (inexact (numerator (exact x)))])))
3390
3391(set-who! fldenominator
3392  (lambda (x)
3393    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3394    (cond
3395      [($flinteger-or-inf? x) 1.0]
3396      [($nan? x) x]
3397      [else (inexact (denominator (exact x)))])))
3398
3399(set-who! fldiv-and-mod
3400  (lambda (x y)
3401    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3402    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3403    ($fldiv-and-mod x y)))
3404
3405(set-who! fldiv
3406  (lambda (x y)
3407    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3408    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3409    ($fldiv x y)))
3410
3411(set-who! flmod
3412  (lambda (x y)
3413    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3414    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3415    ($flmod x y)))
3416
3417(set-who! fldiv0-and-mod0
3418  (lambda (x y)
3419    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3420    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3421    ($fldiv0-and-mod0 x y)))
3422
3423(set-who! fldiv0
3424  (lambda (x y)
3425    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3426    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3427    ($fldiv0 x y)))
3428
3429(set-who! flmod0
3430  (lambda (x y)
3431    (unless (flonum? x) ($oops who "~s is not a flonum" x))
3432    (unless (flonum? y) ($oops who "~s is not a flonum" y))
3433    ($flmod0 x y)))
3434
3435(set-who! fxdiv-and-mod
3436  (lambda (x y)
3437    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3438    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3439    (when (fx= y 0) (domain-error who y))
3440    ($fxdiv-and-mod x y who)))
3441
3442(set-who! fxdiv
3443  (lambda (x y)
3444    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3445    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3446    (when (fx= y 0) (domain-error who y))
3447    ($fxdiv x y who)))
3448
3449(set-who! fxmod
3450  (lambda (x y)
3451    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3452    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3453    (when (fx= y 0) (domain-error who y))
3454    ($fxmod x y)))
3455
3456(set-who! fxdiv0-and-mod0
3457  (lambda (x y)
3458    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3459    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3460    (when (fx= y 0) (domain-error who y))
3461    ($fxdiv0-and-mod0 x y who)))
3462
3463(set-who! fxdiv0
3464  (lambda (x y)
3465    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3466    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3467    (when (fx= y 0) (domain-error who y))
3468    ($fxdiv0 x y who)))
3469
3470(set-who! fxmod0
3471  (lambda (x y)
3472    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3473    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3474    (when (fx= y 0) (domain-error who y))
3475    ($fxmod0 x y)))
3476
3477(let ()
3478  (define (return n)
3479    (if (fixnum? n)
3480        (values n 0)
3481        (if ($bigpositive? n)
3482            (values (- n (expt 2 (fixnum-width))) 1)
3483            (values (+ n (expt 2 (fixnum-width))) -1))))
3484
3485  (set-who! fx+/carry
3486    (lambda (x y z)
3487      (cond
3488        [($fx+? ($fx+? x y) z) => (lambda (n) (values n 0))]
3489        [else
3490         (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3491         (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3492         (unless (fixnum? z) ($oops who "~s is not a fixnum" z))
3493         (return (+ x y z))])))
3494
3495  (set-who! fx-/carry
3496    (lambda (x y z)
3497      (cond
3498        [($fx-? ($fx-? x y) z) => (lambda (n) (values n 0))]
3499        [else
3500         (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3501         (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3502         (unless (fixnum? z) ($oops who "~s is not a fixnum" z))
3503         (return (- x y z))]))))
3504
3505(set-who! fx*/carry
3506  (lambda (x y z)
3507    (unless (fixnum? x) ($oops who "~s is not a fixnum" x))
3508    (unless (fixnum? y) ($oops who "~s is not a fixnum" y))
3509    (let ([t (* x y)])
3510      (cond
3511        [($fx+? t z) => (lambda (n) (values n 0))]
3512        [else
3513         (unless (fixnum? z) ($oops who "~s is not a fixnum" z))
3514         (let-values ([(q r) ($exdiv0-and-mod0 (+ (* x y) z) (expt 2 (fixnum-width)))])
3515           (values r q))]))))
3516
3517(set-who! bitwise-copy-bit-field
3518  (lambda (n start end m)
3519    (unless (or (fixnum? n) (bignum? n))
3520      ($oops who "~s is not an exact integer" n))
3521    (unless (or (and (fixnum? start) (fx>= start 0))
3522                (and (bignum? start) ($bigpositive? start)))
3523      ($oops who "invalid start index ~s" start))
3524    (unless (or (and (fixnum? end) (fixnum? start) (fx>= end start))
3525                (and (bignum? end) (>= end start)))
3526      ($oops who "invalid end index ~s" end))
3527    (unless (or (fixnum? m) (bignum? m))
3528      ($oops who "~s is not an exact integer" m))
3529    (let ([mask (- ($sll who 1 (- end start)) 1)])
3530      (logor
3531        (logand n (lognot ($sll who mask start)))
3532        ($sll who (logand m mask) start)))))
3533
3534(set-who! bitwise-rotate-bit-field
3535  (lambda (n start end count)
3536    (unless (or (fixnum? n) (bignum? n))
3537      ($oops who "~s is not an exact integer" n))
3538    (unless (or (and (fixnum? start) (fx>= start 0))
3539                (and (bignum? start) ($bigpositive? start)))
3540      ($oops who "invalid start index ~s" start))
3541    (unless (or (and (fixnum? end) (fixnum? start) (fx>= end start))
3542                (and (bignum? end) (>= end start)))
3543      ($oops who "invalid end index ~s" end))
3544    (unless (or (and (fixnum? count) (fx>= count 0))
3545                (and (bignum? count) ($bigpositive? count)))
3546      ($oops who "invalid count ~s" count))
3547    (let ([width (- end start)])
3548      (if (positive? width)
3549          (let ([count (mod count width)]
3550                [mask ($sll who (- ($sll who 1 width) 1) start)])
3551            (let ([field (logand n mask)])
3552              (logxor
3553                (logxor
3554                  (logand
3555                    (logor ($sll who field count)
3556                           ($sra who field (- width count)))
3557                    mask)
3558                  field)
3559                n)))
3560          n))))
3561
3562(set-who! fxrotate-bit-field
3563  (lambda (n start end count)
3564    (unless (fixnum? n) ($oops who "~s is not a fixnum" n))
3565    (unless (and (fixnum? end) ($fxu< end (fixnum-width)))
3566      ($oops who "invalid end index ~s" end))
3567    (unless (and (fixnum? start) (not ($fxu< end start)))
3568      (if (and (fixnum? start) ($fxu< start (fixnum-width)))
3569          ($oops who "start index ~s is greater than end index ~s" start end)
3570          ($oops who "invalid start index ~s" start)))
3571    (let ([width (fx- end start)])
3572      (unless (and (fixnum? count) (not ($fxu< width count)))
3573        (if (and (fixnum? count) ($fxu< count (fixnum-width)))
3574            ($oops who "count ~s is greater than difference between end index ~s and start index ~s" count end start)
3575            ($oops who "invalid count ~s" count)))
3576      (let ([mask (fxsll (fxsrl -1 (fx- (fixnum-width) width)) start)])
3577        (let ([field (fxlogand n mask)])
3578          (fxlogor
3579            (fxlogxor n field)
3580            (fxlogand
3581              (fxlogor
3582                (fxsll (fxlogand field (fxsrl mask count)) count)
3583                (fxsrl field (fx- width count)))
3584              mask)))))))
3585
3586(let ()
3587  (define rev-table
3588    (let ()
3589      (define-syntax make-rev-table
3590        (lambda (x)
3591          #`'#,(let ([bv (make-bytevector 256)])
3592                 (for-each
3593                   (lambda (m)
3594                     (bytevector-u8-set! bv m
3595                       (do ([m m (fxsrl m 1)]
3596                            [m^ 0 (fxior (fxsll m^ 1) (fxand m 1))]
3597                            [k 8 (fx- k 1)])
3598                           ((fx= k 0) m^))))
3599                   (iota 256))
3600                 bv)))
3601      (make-rev-table)))
3602
3603  (define $fxreverse
3604    (lambda (m k)
3605      (do ([m m (fxsrl m 8)]
3606           [m^ 0 (fxior (fxsll m^ 8) (bytevector-u8-ref rev-table (fxand m #xff)))]
3607           [k k (fx- k 8)])
3608          ((fx< k 8)
3609           (fxior
3610             (fxsll m^ k)
3611             (fxsrl (bytevector-u8-ref rev-table m) (fx- 8 k)))))))
3612
3613  (set-who! fxreverse-bit-field
3614    (lambda (n start end)
3615      (unless (fixnum? n) ($oops who "~s is not a fixnum" n))
3616      (unless (and (fixnum? start) ($fxu< start (fixnum-width)))
3617        ($oops who "invalid start index ~s" start))
3618      (unless (and (fixnum? end) ($fxu< end (fixnum-width)))
3619        ($oops who "invalid end index ~s" end))
3620      (unless (fx<= start end)
3621        ($oops who "start index ~s is greater than end index ~s" start end))
3622      (fxcopy-bit-field n start end
3623        ($fxreverse (fxbit-field n start end) (fx- end start)))))
3624
3625  (set-who! bitwise-reverse-bit-field
3626    (lambda (n start end)
3627      (define sra bitwise-arithmetic-shift-right)
3628      (define sll bitwise-arithmetic-shift-left)
3629      (define w-1 (fx- (fixnum-width) 1))
3630      (define mask (- (sll 1 w-1) 1))
3631      (unless (or (fixnum? n) (bignum? n))
3632        ($oops who "~s is not an exact integer" n))
3633      (unless (or (and (fixnum? start) (fx>= start 0))
3634                  (and (bignum? start) ($bigpositive? start)))
3635        ($oops who "invalid start index ~s" start))
3636      (unless (or (and (fixnum? end) (fx>= end 0))
3637                  (and (bignum? end) ($bigpositive? end)))
3638        ($oops who "invalid end index ~s" end))
3639      (unless (<= start end)
3640        ($oops who "start index ~s is greater than end index ~s" start end))
3641      (bitwise-copy-bit-field n start end
3642        (do ([m (bitwise-bit-field n start end) (sra m w-1)]
3643             [m^ 0 (logor (sll m^ w-1) ($fxreverse (logand m mask) w-1))]
3644             [k (- end start) (- k w-1)])
3645            ((<= k w-1) (logor (sll m^ k) ($fxreverse m k))))))))
3646))))))))
3647)
3648