home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / runtime / arith.scm < prev    next >
Text File  |  1999-01-02  |  60KB  |  1,938 lines
  1. #| -*-Scheme-*-
  2.  
  3. $Id: arith.scm,v 1.45 1999/01/02 06:11:34 cph Exp $
  4.  
  5. Copyright (c) 1989-1999 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. |#
  21.  
  22. ;;;; Scheme Arithmetic
  23. ;;; package: (runtime number)
  24.  
  25. (declare (usual-integrations))
  26.  
  27. ;;;; Utilities
  28.  
  29. (define-macro (copy x)
  30.   `(LOCAL-DECLARE ((INTEGRATE ,x)) ,x))
  31.  
  32. ;;;; Primitives
  33.  
  34. (define-primitives
  35.   (listify-bignum 2)
  36.   (integer->flonum 2)
  37.   (fixnum->flonum 1)
  38.   (flo:denormalize flonum-denormalize 2)
  39.   (integer-length-in-bits 1)
  40.   (integer-shift-left 2))
  41.  
  42. (define-integrable (int:bignum? object)
  43.   (object-type? (ucode-type big-fixnum) object))
  44.  
  45. (define-integrable (make-ratnum n d)
  46.   (system-pair-cons (ucode-type ratnum) n d))
  47.  
  48. (define-integrable (ratnum? object)
  49.   (object-type? (ucode-type ratnum) object))
  50.  
  51. (define-integrable (ratnum-numerator ratnum)
  52.   (system-pair-car ratnum))
  53.  
  54. (define-integrable (ratnum-denominator ratnum)
  55.   (system-pair-cdr ratnum))
  56.  
  57. (define-integrable (flonum? object)
  58.   (object-type? (ucode-type big-flonum) object))
  59.  
  60. (define (flo:normalize x)
  61.   (let ((r ((ucode-primitive flonum-normalize 1) x)))
  62.     (values (car r) (cdr r))))
  63.  
  64. (define-integrable flo:->integer
  65.   flo:truncate->exact)
  66.  
  67. (define-integrable (recnum? object)
  68.   (object-type? (ucode-type recnum) object))
  69.  
  70. (define-integrable (make-recnum real imag)
  71.   (system-pair-cons (ucode-type recnum) real imag))
  72.  
  73. (define-integrable (rec:real-part recnum)
  74.   (system-pair-car recnum))
  75.  
  76. (define-integrable (rec:imag-part recnum)
  77.   (system-pair-cdr recnum))
  78.  
  79. ;;;; Constants
  80.  
  81. (define-integrable flo:0 0.)
  82. (define-integrable flo:1 1.)
  83. (define rec:pi/2 (flo:* 2. (flo:atan2 1. 1.)))
  84. (define rec:pi (flo:* 2. rec:pi/2))
  85.  
  86. (define flo:significand-digits-base-2)
  87. (define flo:significand-digits-base-10)
  88. (define int:flonum-integer-limit)
  89.  
  90. (define (initialize-microcode-dependencies!)
  91.   (let ((p microcode-id/floating-mantissa-bits))
  92.     (set! flo:significand-digits-base-2 p)
  93.     ;; Add two here because first and last digits may be
  94.     ;; "partial" in the sense that each represents less than the
  95.     ;; `flo:log10/log2' bits.  This is a kludge, but doing the
  96.     ;; "right thing" seems hard.  See Steele&White for a discussion of
  97.     ;; this phenomenon.
  98.     (set! flo:significand-digits-base-10
  99.       (int:+ 2
  100.          (flo:floor->exact
  101.           (flo:/ (int:->flonum p)
  102.              (flo:/ (flo:log 10.) (flo:log 2.))))))
  103.     (set! int:flonum-integer-limit (int:expt 2 p)))
  104.   unspecific)
  105.  
  106. (define (initialize-package!)
  107.   (initialize-microcode-dependencies!)
  108.   (add-event-receiver! event:after-restore initialize-microcode-dependencies!)
  109.   (initialize-*maximum-fixnum-radix-powers*!)
  110.   (let ((fixed-objects-vector (get-fixed-objects-vector)))
  111.     (let ((set-trampoline!
  112.        (lambda (slot operator)
  113.          (vector-set! fixed-objects-vector
  114.               (fixed-objects-vector-slot slot)
  115.               operator))))
  116.       (set-trampoline! 'GENERIC-TRAMPOLINE-ZERO? complex:zero?)
  117.       (set-trampoline! 'GENERIC-TRAMPOLINE-POSITIVE? complex:positive?)
  118.       (set-trampoline! 'GENERIC-TRAMPOLINE-NEGATIVE? complex:negative?)
  119.       (set-trampoline! 'GENERIC-TRAMPOLINE-ADD-1 complex:1+)
  120.       (set-trampoline! 'GENERIC-TRAMPOLINE-SUBTRACT-1 complex:-1+)
  121.       (set-trampoline! 'GENERIC-TRAMPOLINE-EQUAL? complex:=)
  122.       (set-trampoline! 'GENERIC-TRAMPOLINE-LESS? complex:<)
  123.       (set-trampoline! 'GENERIC-TRAMPOLINE-GREATER? complex:>)
  124.       (set-trampoline! 'GENERIC-TRAMPOLINE-ADD complex:+)
  125.       (set-trampoline! 'GENERIC-TRAMPOLINE-SUBTRACT complex:-)
  126.       (set-trampoline! 'GENERIC-TRAMPOLINE-MULTIPLY complex:*)
  127.       (set-trampoline! 'GENERIC-TRAMPOLINE-DIVIDE complex:/)
  128.       (set-trampoline! 'GENERIC-TRAMPOLINE-QUOTIENT complex:quotient)
  129.       (set-trampoline! 'GENERIC-TRAMPOLINE-REMAINDER complex:remainder)
  130.       (set-trampoline! 'GENERIC-TRAMPOLINE-MODULO complex:modulo)))
  131.  
  132.   ;; The binary cases for the following operators rely on the fact that the
  133.   ;; &<mumble> operators, either interpreted or open-coded by the
  134.   ;; compiler, calls the GENERIC-TRAMPOLINE version above, are set to
  135.   ;; the appropriate binary procedures when this package is
  136.   ;; initialized.  We could have just replaced (ucode-primitive &+)
  137.   ;; with + etc and relied on + being integrated, but that is not
  138.   ;; very clear.
  139.  
  140.   (let-syntax
  141.       ((commutative
  142.     (macro (name generic-binary identity primitive-binary)
  143.       `(SET! ,name
  144.          (MAKE-ENTITY
  145.           (NAMED-LAMBDA (,name SELF . ZS)
  146.             SELF        ; ignored
  147.             (REDUCE ,generic-binary ,identity ZS))
  148.           (VECTOR (FIXED-OBJECTS-ITEM 'arity-dispatcher-tag)
  149.               (NAMED-LAMBDA (,(symbol-append 'NULLARY- name))
  150.                 ,identity)
  151.               (NAMED-LAMBDA (,(symbol-append 'UNARY- name) Z)
  152.                 (IF (NOT (COMPLEX:COMPLEX? Z))
  153.                 (ERROR:WRONG-TYPE-ARGUMENT Z FALSE ',name))
  154.                 Z)
  155.               (NAMED-LAMBDA (,(symbol-append 'BINARY- name) Z1 Z2)
  156.                 ((UCODE-PRIMITIVE ,primitive-binary) Z1 Z2))))))))
  157.     (commutative + complex:+ 0 &+)
  158.     (commutative * complex:* 1 &*))
  159.  
  160.   (let-syntax
  161.       ((non-commutative
  162.     (macro (name generic-unary generic-binary
  163.              generic-inverse inverse-identity primitive-binary)
  164.       `(SET! ,name
  165.          (MAKE-ENTITY
  166.           (NAMED-LAMBDA (,name SELF Z1 . ZS)
  167.             SELF        ; ignored
  168.             (,generic-binary
  169.              Z1
  170.              (REDUCE ,generic-inverse ,inverse-identity ZS)))
  171.           (VECTOR (FIXED-OBJECTS-ITEM 'arity-dispatcher-tag)
  172.               #F
  173.               ,generic-unary
  174.               (NAMED-LAMBDA (,(symbol-append 'BINARY- name) Z1 Z2)
  175.                 ((UCODE-PRIMITIVE ,primitive-binary) Z1 Z2))))))))
  176.     (non-commutative -  complex:negate  complex:-  complex:+  0  &-)
  177.     (non-commutative /  complex:invert  complex:/  complex:*  1  &/))
  178.  
  179.   (let-syntax
  180.       ((relational
  181.     (macro (name generic-binary primitive-binary correct-type? negated?)
  182.       `(SET! ,name
  183.          (MAKE-ENTITY
  184.           (NAMED-LAMBDA (,name SELF . ZS)
  185.             SELF        ; ignored
  186.             (REDUCE-COMPARATOR ,generic-binary ZS ',name))
  187.           (VECTOR
  188.            (FIXED-OBJECTS-ITEM 'arity-dispatcher-tag)
  189.            (NAMED-LAMBDA (,(symbol-append 'NULLARY- name)) #T)
  190.            (NAMED-LAMBDA (,(symbol-append 'UNARY- name) Z)
  191.              (IF (NOT (,correct-type? Z))
  192.              (ERROR:WRONG-TYPE-ARGUMENT Z FALSE ',name))
  193.              #T)
  194.            ,(if negated?
  195.             `(NAMED-LAMBDA (,(symbol-append 'BINARY- name) Z1 Z2)
  196.                (NOT ((UCODE-PRIMITIVE ,primitive-binary) Z1 Z2)))
  197.             `(NAMED-LAMBDA (,(symbol-append 'BINARY- name) Z1 Z2)
  198.                ((UCODE-PRIMITIVE ,primitive-binary) Z1 Z2)))))))))
  199.     (relational =  complex:=  &=  complex:complex? #F)
  200.     (relational <  complex:<  &<  complex:real?    #F)
  201.     (relational >  complex:>  &>  complex:real?    #F)
  202.     (relational <= (lambda (x y) (not (complex:< y x)))  &>  complex:real? #T)
  203.     (relational >= (lambda (x y) (not (complex:< x y)))  &<  complex:real? #T))
  204.  
  205.   (let-syntax
  206.       ((max/min
  207.     (macro (name generic-binary)
  208.       `(SET! ,name
  209.          (MAKE-ENTITY
  210.           (NAMED-LAMBDA (,name SELF X . XS)
  211.             SELF        ; ignored
  212.             (REDUCE-MAX/MIN ,generic-binary X XS ',name))
  213.           (VECTOR
  214.            (FIXED-OBJECTS-ITEM 'arity-dispatcher-tag)
  215.            #F
  216.            (NAMED-LAMBDA (,(symbol-append 'UNARY- name) X)
  217.              (IF (NOT (COMPLEX:REAL? X))
  218.              (ERROR:WRONG-TYPE-ARGUMENT X FALSE ',name))
  219.              X)
  220.            ,generic-binary))))))
  221.     (max/min max complex:max)
  222.     (max/min min complex:min))
  223.  
  224.   unspecific)
  225.  
  226. (define (int:max n m)
  227.   (if (int:< n m) m n))
  228.  
  229. (define (int:min n m)
  230.   (if (int:< n m) n m))
  231.  
  232. (define (int:abs n)
  233.   (if (int:negative? n) (int:negate n) n))
  234.  
  235. (define (int:even? n)
  236.   (int:zero? (int:remainder n 2)))
  237.  
  238. (define (int:modulo n d)
  239.   (let ((r (int:remainder n d)))
  240.     (if (or (int:zero? r)
  241.         (if (int:negative? n)
  242.         (int:negative? d)
  243.         (not (int:negative? d))))
  244.     r
  245.     (int:+ r d))))
  246.  
  247. (define (int:gcd n m)
  248.   (let loop ((n n) (m m))
  249.     (cond ((not (int:zero? m)) (loop m (int:remainder n m)))
  250.       ((int:negative? n) (int:negate n))
  251.       (else n))))
  252.  
  253. (define (int:lcm n m)
  254.   (if (or (int:zero? n) (int:zero? m))
  255.       0
  256.       (int:quotient (let ((n (int:* n m)))
  257.               (if (int:negative? n)
  258.               (int:negate n)
  259.               n))
  260.             (int:gcd n m))))
  261.  
  262. (define (int:floor n d)
  263.   (let ((qr (int:divide n d)))
  264.     (let ((q (integer-divide-quotient qr)))
  265.       (if (or (int:zero? (integer-divide-remainder qr))
  266.           (if (int:negative? n)
  267.           (int:negative? d)
  268.           (not (int:negative? d))))
  269.       q
  270.       (int:-1+ q)))))
  271.  
  272. (define (int:ceiling n d)
  273.   (let ((qr (int:divide n d)))
  274.     (let ((q (integer-divide-quotient qr)))
  275.       (if (or (int:zero? (integer-divide-remainder qr))
  276.           (if (int:negative? n)
  277.           (not (int:negative? d))
  278.           (int:negative? d)))
  279.       q
  280.       (int:1+ q)))))
  281.  
  282. (define (int:round n d)
  283.   (let ((positive-case
  284.      (lambda (n d)
  285.        (let ((qr (int:divide n d)))
  286.          (let ((q (integer-divide-quotient qr))
  287.            (2r (int:* 2 (integer-divide-remainder qr))))
  288.            (if (or (int:> 2r d)
  289.                (and (int:= 2r d)
  290.                 (not (fix:zero? (int:remainder q 2)))))
  291.            (int:1+ q)
  292.            q))))))
  293.     (if (int:negative? n)
  294.     (if (int:negative? d)
  295.         (positive-case (int:negate n) (int:negate d))
  296.         (int:negate (positive-case (int:negate n) d)))
  297.     (if (int:negative? d)
  298.         (int:negate (positive-case n (int:negate d)))
  299.         (positive-case n d)))))
  300.  
  301. (define (int:expt b e)
  302.   (cond ((int:positive? e)
  303.      (cond ((or (int:= 1 e)
  304.             (int:zero? b)
  305.             (int:= 1 b))
  306.         b)
  307.            ((int:= 2 b)
  308.         (integer-shift-left 1 e))
  309.            (else
  310.         (let loop ((b b) (e e) (answer 1))
  311.           (let ((qr (int:divide e 2)))
  312.             (let ((b (int:* b b))
  313.               (e (integer-divide-quotient qr))
  314.               (answer
  315.                (if (fix:= 0 (integer-divide-remainder qr))
  316.                    answer
  317.                    (int:* answer b))))
  318.               (if (int:= 1 e)
  319.               (int:* answer b)
  320.               (loop b e answer))))))))
  321.     ((int:zero? e) 1)
  322.     (else (error:bad-range-argument e 'EXPT))))
  323.  
  324. ;; A vector indexed by radix of pairs of the form (N . (expt RADIX N))
  325. ;; where N is the maximum value for which the cdr is a fixnum.  Used
  326. ;; to quickly determine how many digits to process at a time to
  327. ;; optimize the use of fixnum arithmetic.
  328.  
  329. (define *maximum-fixnum-radix-powers*)
  330.  
  331. (define (initialize-*maximum-fixnum-radix-powers*!)
  332.   (set! *maximum-fixnum-radix-powers*
  333.     (make-initialized-vector 37
  334.       (lambda (radix)
  335.         (and (fix:>= radix 2)
  336.          (let loop ((digits 0) (factor 1))
  337.            (let ((nf (int:* factor radix)))
  338.              (if (fix:fixnum? nf)
  339.              (loop (int:+ digits 1) nf)
  340.              (cons digits factor)))))))))
  341.  
  342. ;; INT:->STRING chooses between 3 strategies for generating the digits:
  343. ;;
  344. ;;  PRINT-FIXNUM exploits fast fixnum arithmetic
  345. ;;  PRINT-MEDIUM chops off groups of digits that can be printed by PRINT-FIXNUM
  346. ;;  PRINT-LARGE works by dividing the problem into approximately equal sizes,
  347. ;;   which is asympotically faster but requires more operations for moderate
  348. ;;   values.
  349.  
  350. (define (int:->string number radix)
  351.   ;; Pre: (and (exact-integer? NUMBER) (fixnum? radix) (<= 2 radix 36))
  352.  
  353.   (define-integrable (digit->char digit radix)
  354.     radix ; ignored
  355.     (string-ref "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" digit))
  356.  
  357.   (define (print-fixnum n min-digits tail)
  358.     (let loop ((n n) (n-digits 0) (tail tail))
  359.       (cond ((not (fix:zero? n))
  360.          (loop (fix:quotient n radix)
  361.            (fix:+ n-digits 1)
  362.            (cons (digit->char (fix:remainder n radix)
  363.                       radix)
  364.              tail)))
  365.         ((fix:< n-digits min-digits)
  366.          (loop n (fix:+ n-digits 1) (cons #\0 tail)))
  367.         (else
  368.          tail))))
  369.  
  370.   (define (print-medium value split-factor split-digits)
  371.     (let loop ((n value) (tail '()))
  372.       (if (fix:fixnum? n)
  373.       (print-fixnum n 0 tail)
  374.       (let ((qr (integer-divide n split-factor)))
  375.         (loop (integer-divide-quotient qr)
  376.           (print-fixnum (integer-divide-remainder qr)
  377.                 split-digits
  378.                 tail))))))
  379.  
  380.   (define (fast-test-to-avoid-ultimate-multiply quantum value)
  381.     ;; Uses the number of bignum `digits' or words to test if
  382.     ;; QUANTUM^2>VALUE (the `-1' skips the bignum internal header).
  383.     ;; Since an N digit multiply takes time O(N^2), the benefit of
  384.     ;; avoiding the last squaring is detectable for VALUE>10^100,
  385.     ;; increases with VALUE, but is limited to about 40-50% overall
  386.     ;; improvement by the division operations.
  387.     (define-integrable (bignum-digits n) (fix:+ -1 (system-vector-length n)))
  388.     (and (not (fixnum? quantum))    ; i.e. bignum
  389.      (fix:> (fix:- (fix:* (bignum-digits quantum) 2) 1)
  390.         (bignum-digits value))))
  391.  
  392.   (define (make-power-stack value quantum stack n-digits)
  393.     (cond ((> quantum value)
  394.        (use-power-stack value stack n-digits))
  395.       ((fast-test-to-avoid-ultimate-multiply quantum value)
  396.        (use-power-stack value (cons quantum stack) (* 2 n-digits)))
  397.       (else
  398.        (make-power-stack value
  399.                  (* quantum quantum)
  400.                  (cons quantum stack)
  401.                  (* 2 n-digits)))))
  402.  
  403.   (define (use-power-stack value stack digits)
  404.     ;; Test at [1] could be (null? stack), but (fixnum? value) is true
  405.     ;; in this case by construction and accelerates printing of numbers
  406.     ;; with a large number of zero digits.
  407.     (define (separate leftmost? value stack n-digits tail)
  408.       (if (fix:fixnum? value)        ; [1]
  409.       (print-fixnum value (if leftmost? 0 n-digits) tail)
  410.       (let ((split (integer-divide value (car stack)))
  411.         (rest (cdr stack)))
  412.         (let ((next-left (integer-divide-quotient split))
  413.           (n-digits/2 (fix:quotient n-digits 2)))
  414.           (if (and leftmost? (zero? next-left))
  415.           (separate true
  416.                 (integer-divide-remainder split)
  417.                 rest
  418.                 n-digits/2
  419.                 tail)
  420.           (separate leftmost?
  421.                 next-left
  422.                 rest
  423.                 n-digits/2
  424.                 (separate false
  425.                       (integer-divide-remainder split)
  426.                       rest
  427.                       n-digits/2
  428.                       tail)))))))
  429.  
  430.     (separate true value stack digits '()))
  431.  
  432.   (define (n>0 value)
  433.     (if (fix:fixnum? value)
  434.     (print-fixnum value 1 '())
  435.     (let* ((split-info (vector-ref *maximum-fixnum-radix-powers* radix))
  436.            (split-digits (car split-info))
  437.            (split-factor (cdr split-info))
  438.            (sl (system-vector-length value)))
  439.       (if (< sl 10)
  440.           (print-medium value split-factor split-digits)
  441.           (make-power-stack value split-factor '() split-digits)))))
  442.   
  443.   (cond ((not (int:integer? number))
  444.      (error:wrong-type-argument number false 'NUMBER->STRING))
  445.     ((int:negative? number)
  446.      (list->string (cons #\- (n>0 (int:negate number)))))
  447.     (else
  448.      (list->string (n>0 number)))))
  449.  
  450. (declare (integrate-operator rat:rational?))
  451. (define (rat:rational? object)
  452.   (or (ratnum? object)
  453.       (int:integer? object)))
  454.  
  455. (define (rat:integer? object)
  456.   (and (not (ratnum? object))
  457.        (int:integer? object)))
  458.  
  459. (define (rat:= q r)
  460.   (if (ratnum? q)
  461.       (if (ratnum? r)
  462.       (and (int:= (ratnum-numerator q) (ratnum-numerator r))
  463.            (int:= (ratnum-denominator q) (ratnum-denominator r)))
  464.       (if (int:integer? r)
  465.           #f
  466.           (error:wrong-type-argument r false '=)))
  467.       (if (ratnum? r)
  468.       (if (int:integer? q)
  469.           #f
  470.           (error:wrong-type-argument q false '=))
  471.       (int:= q r))))
  472.  
  473. (define (rat:< q r)
  474.   (if (ratnum? q)
  475.       (if (ratnum? r)
  476.       (int:< (int:* (ratnum-numerator q) (ratnum-denominator r))
  477.          (int:* (ratnum-numerator r) (ratnum-denominator q)))
  478.       (int:< (ratnum-numerator q) (int:* r (ratnum-denominator q))))
  479.       (if (ratnum? r)
  480.       (int:< (int:* q (ratnum-denominator r)) (ratnum-numerator r))
  481.       (int:< q r))))
  482.  
  483. (define (rat:zero? q)
  484.   (and (not (ratnum? q))
  485.        (int:zero? q)))
  486.  
  487. (define (rat:negative? q)
  488.   (if (ratnum? q)
  489.       (int:negative? (ratnum-numerator q))
  490.       (int:negative? q)))
  491.  
  492. (define (rat:positive? q)
  493.   (if (ratnum? q)
  494.       (int:positive? (ratnum-numerator q))
  495.       (int:positive? q)))
  496.  
  497. (define (rat:max m n)
  498.   (if (rat:< m n) n m))
  499.  
  500. (define (rat:min m n)
  501.   (if (rat:< m n) m n))
  502.  
  503. ;;; The notation here is from Knuth (p. 291).
  504. ;;; In various places we take the gcd of two numbers and then call
  505. ;;; quotient to reduce those numbers.  We could check for 1 here, but
  506. ;;; this is generally important only for bignums, and the bignum
  507. ;;; quotient already performs that check.
  508.  
  509. (let-syntax
  510.     ((define-addition-operator
  511.        (macro (name int:op)
  512.      `(define (,name u/u* v/v*)
  513.         (rat:binary-operator u/u* v/v*
  514.           ,int:op
  515.           (lambda (u v v*)
  516.         (make-rational (,int:op (int:* u v*) v) v*))
  517.           (lambda (u u* v)
  518.         (make-rational (,int:op u (int:* v u*)) u*))
  519.           (lambda (u u* v v*)
  520.         (let ((d1 (int:gcd u* v*)))
  521.           (if (int:= d1 1)
  522.               (make-rational (,int:op (int:* u v*) (int:* v u*))
  523.                      (int:* u* v*))
  524.               (let* ((u*/d1 (int:quotient u* d1))
  525.                  (t
  526.                   (,int:op (int:* u (int:quotient v* d1))
  527.                        (int:* v u*/d1))))
  528.             (if (int:zero? t)
  529.                 0 ;(make-rational 0 1)
  530.                 (let ((d2 (int:gcd t d1)))
  531.                   (make-rational
  532.                    (int:quotient t d2)
  533.                    (int:* u*/d1 (int:quotient v* d2))))))))))))))
  534.   (define-addition-operator rat:+ int:+)
  535.   (define-addition-operator rat:- int:-))
  536.  
  537. (define (rat:1+ v/v*)
  538.   (if (ratnum? v/v*)
  539.       (let ((v* (ratnum-denominator v/v*)))
  540.     (make-ratnum (int:+ (ratnum-numerator v/v*) v*) v*))
  541.       (int:1+ v/v*)))
  542.  
  543. (define (rat:-1+ v/v*)
  544.   (if (ratnum? v/v*)
  545.       (let ((v* (ratnum-denominator v/v*)))
  546.     (make-ratnum (int:- (ratnum-numerator v/v*) v*) v*))
  547.       (int:-1+ v/v*)))
  548.  
  549. (define (rat:negate v/v*)
  550.   (if (ratnum? v/v*)
  551.       (make-ratnum (int:negate (ratnum-numerator v/v*))
  552.            (ratnum-denominator v/v*))
  553.       (int:negate v/v*)))
  554.  
  555. (define (rat:* u/u* v/v*)
  556.   (rat:binary-operator u/u* v/v*
  557.     int:*
  558.     (lambda (u v v*)
  559.       (let ((d (int:gcd u v*)))
  560.     (make-rational (int:* (int:quotient u d) v)
  561.                (int:quotient v* d))))
  562.     (lambda (u u* v)
  563.       (let ((d (int:gcd v u*)))
  564.     (make-rational (int:* u (int:quotient v d))
  565.                (int:quotient u* d))))
  566.     (lambda (u u* v v*)
  567.       (let ((d1 (int:gcd u v*))
  568.         (d2 (int:gcd v u*)))
  569.     (make-rational (int:* (int:quotient u d1) (int:quotient v d2))
  570.                (int:* (int:quotient u* d2) (int:quotient v* d1)))))))
  571.  
  572. (define (rat:square q)
  573.   (if (ratnum? q)
  574.       (make-ratnum (let ((n (ratnum-numerator q))) (int:* n n))
  575.            (let ((d (ratnum-denominator q))) (int:* d d)))
  576.       (int:* q q)))
  577.  
  578. (define (rat:/ u/u* v/v*)
  579.   (declare (integrate-operator rat:sign-correction))
  580.   (define (rat:sign-correction u v cont)
  581.     (declare (integrate u v))
  582.     (if (int:negative? v)
  583.     (cont (int:negate u) (int:negate v))
  584.     (cont u v)))
  585.   (rat:binary-operator u/u* v/v*
  586.     (lambda (u v)
  587.       (if (int:zero? v)
  588.       (error:divide-by-zero '/ (list u/u* v/v*))
  589.       (rat:sign-correction u v
  590.         (lambda (u v)
  591.           (let ((d (int:gcd u v)))
  592.         (make-rational (int:quotient u d)
  593.                    (int:quotient v d)))))))
  594.     (lambda (u v v*)
  595.       (rat:sign-correction u v
  596.     (lambda (u v)
  597.       (let ((d (int:gcd u v)))
  598.         (make-rational (int:* (int:quotient u d) v*)
  599.                (int:quotient v d))))))
  600.     (lambda (u u* v)
  601.       (if (int:zero? v)
  602.       (error:divide-by-zero '/ (list u/u* v/v*))
  603.       (rat:sign-correction u v
  604.         (lambda (u v)
  605.           (let ((d (int:gcd u v)))
  606.         (make-rational (int:quotient u d)
  607.                    (int:* u* (int:quotient v d))))))))
  608.     (lambda (u u* v v*)
  609.       (let ((d1 (int:gcd u v))
  610.         (d2
  611.          (let ((d2 (int:gcd v* u*)))
  612.            (if (int:negative? v)
  613.            (int:negate d2)
  614.            d2))))
  615.     (make-rational (int:* (int:quotient u d1) (int:quotient v* d2))
  616.                (int:* (int:quotient u* d2) (int:quotient v d1)))))))
  617.  
  618. (define (rat:invert v/v*)
  619.   (if (ratnum? v/v*)
  620.       (let ((v (ratnum-numerator v/v*))
  621.         (v* (ratnum-denominator v/v*)))
  622.     (cond ((int:positive? v)
  623.            (make-rational v* v))
  624.           ((int:negative? v)
  625.            (make-rational (int:negate v*) (int:negate v)))
  626.           (else
  627.            (error:divide-by-zero '/ (list 1 v/v*)))))
  628.       (cond ((int:positive? v/v*) (make-rational 1 v/v*))
  629.         ((int:negative? v/v*) (make-rational -1 (int:negate v/v*)))
  630.         (else (error:divide-by-zero '/ (list 1 v/v*))))))
  631.  
  632. (define-integrable (rat:binary-operator u/u* v/v*
  633.                     int*int int*rat rat*int rat*rat)
  634.   (if (ratnum? u/u*)
  635.       (if (ratnum? v/v*)
  636.       (rat*rat (ratnum-numerator u/u*)
  637.            (ratnum-denominator u/u*)
  638.            (ratnum-numerator v/v*)
  639.            (ratnum-denominator v/v*))
  640.       (rat*int (ratnum-numerator u/u*)
  641.            (ratnum-denominator u/u*)
  642.            v/v*))
  643.       (if (ratnum? v/v*)
  644.       (int*rat u/u*
  645.           (ratnum-numerator v/v*)
  646.           (ratnum-denominator v/v*))
  647.       (int*int u/u* v/v*))))
  648.  
  649. (define (rat:abs q)
  650.   (cond ((ratnum? q)
  651.      (let ((numerator (ratnum-numerator q)))
  652.        (if (int:negative? numerator)
  653.            (make-ratnum (int:negate numerator) (ratnum-denominator q))
  654.            q)))
  655.     ((int:negative? q) (int:negate q))
  656.     (else q)))
  657.  
  658. (define (rat:numerator q)
  659.   (cond ((ratnum? q) (ratnum-numerator q))
  660.     ((int:integer? q) q)
  661.     (else (error:wrong-type-argument q false 'NUMERATOR))))
  662.  
  663. (define (rat:denominator q)
  664.   (cond ((ratnum? q) (ratnum-denominator q))
  665.     ((int:integer? q) 1)
  666.     (else (error:wrong-type-argument q false 'DENOMINATOR))))
  667.  
  668. (let-syntax
  669.     ((define-integer-coercion
  670.        (macro (name operation-name coercion)
  671.      `(DEFINE (,name Q)
  672.         (COND ((RATNUM? Q)
  673.            (,coercion (RATNUM-NUMERATOR Q) (RATNUM-DENOMINATOR Q)))
  674.           ((INT:INTEGER? Q) Q)
  675.           (ELSE
  676.            (ERROR:WRONG-TYPE-ARGUMENT Q FALSE ',operation-name)))))))
  677.   (define-integer-coercion rat:floor floor int:floor)
  678.   (define-integer-coercion rat:ceiling ceiling int:ceiling)
  679.   (define-integer-coercion rat:truncate truncate int:quotient)
  680.   (define-integer-coercion rat:round round int:round))
  681.  
  682. (define (rat:rationalize q e)
  683.   (rat:simplest-rational (rat:- q e) (rat:+ q e)))
  684.  
  685. (define (rat:simplest-rational x y)
  686.   ;; Courtesy of Alan Bawden.
  687.   ;; Produces the simplest rational between X and Y inclusive.
  688.   ;; (In the comments that follow, [x] means (rat:floor x).)
  689.   (let ((x<y
  690.      (lambda (x y)
  691.        (define (loop x y)
  692.          (if (int:integer? x)
  693.          x
  694.          (let ((fx (rat:floor x)) ; [X] <= X < [X]+1
  695.                (fy (rat:floor y))) ; [Y] <= Y < [Y]+1, also [X] <= [Y]
  696.            (if (rat:= fx fy)
  697.                ;; [Y] = [X] < X < Y so expand the next term in
  698.                ;; the continued fraction:
  699.                (rat:+ fx
  700.                   (rat:invert (loop (rat:invert (rat:- y fy))
  701.                         (rat:invert (rat:- x fx)))))
  702.                ;; [X] < X < [X]+1 <= [Y] <= Y so [X]+1 is the answer:
  703.                (rat:1+ fx)))))
  704.        (cond ((rat:positive? x)
  705.           ;; 0 < X < Y
  706.           (loop x y))
  707.          ((rat:negative? y)
  708.           ;; X < Y < 0 so 0 < -Y < -X and we negate the answer:
  709.           (rat:negate (loop (rat:negate y) (rat:negate x))))
  710.          (else
  711.           ;; X <= 0 <= Y so zero is the answer:
  712.           0)))))
  713.     (cond ((rat:< x y) (x<y x y))
  714.       ((rat:< y x) (x<y y x))
  715.       (else x))))
  716.  
  717. (define (rat:expt b e)
  718.   (if (int:integer? e)
  719.       (if (int:integer? b)
  720.       (if (int:negative? e)
  721.           (rat:invert (int:expt b (int:negate e)))
  722.           (int:expt b e))
  723.       (let ((exact-method
  724.          (lambda (e)
  725.            (if (int:= 1 e)
  726.                b
  727.                (let loop ((b b) (e e) (answer 1))
  728.              (let ((qr (int:divide e 2)))
  729.                (let ((b (rat:* b b))
  730.                  (e (integer-divide-quotient qr))
  731.                  (answer
  732.                   (if (int:zero? (integer-divide-remainder qr))
  733.                       answer
  734.                       (rat:* answer b))))
  735.                  (if (int:= 1 e)
  736.                  (rat:* answer b)
  737.                  (loop b e answer)))))))))
  738.         (cond ((int:negative? e)
  739.            (rat:invert (exact-method (int:negate e))))
  740.           ((int:positive? e)
  741.            (exact-method e))
  742.           (else 1))))
  743.       (error:bad-range-argument e 'EXPT)))
  744.  
  745. (define (rat:->string q radix)
  746.   (if (ratnum? q)
  747.       (string-append (int:->string (ratnum-numerator q) radix)
  748.              "/"
  749.              (int:->string (ratnum-denominator q) radix))
  750.       (int:->string q radix)))
  751.  
  752. (define (make-rational n d)
  753.   (if (or (int:zero? n) (int:= 1 d))
  754.       n
  755.       (make-ratnum n d)))
  756.  
  757. (define (rat:->inexact q)
  758.   (if (ratnum? q)
  759.       (ratio->flonum (ratnum-numerator q) (ratnum-denominator q))
  760.       (int:->inexact q)))
  761.  
  762. (define (ratio->flonum n d)
  763.   (define (n>0 n d)
  764.     (let ((k (int:- (integer-length-in-bits n)
  765.             (integer-length-in-bits d)))
  766.       (p flo:significand-digits-base-2))
  767.       (letrec
  768.       ((step1
  769.         (lambda (n d)
  770.           ;; (assert (< (expt 2 (- k 1)) (/ n d) (expt 2 (+ k 1))))
  771.           (if (int:negative? k)
  772.           (step2 (integer-shift-left n (int:negate k)) d)
  773.           (step2 n (integer-shift-left d k)))))
  774.        (step2
  775.         (lambda (n d)
  776.           ;; (assert (< 1/2 (/ n d) 2))
  777.           (if (int:< n d)
  778.           (step3 n d (int:- k p))
  779.           (step3 n (int:* 2 d) (int:- (int:1+ k) p)))))
  780.        (step3
  781.         (lambda (n d e)
  782.           ;; (assert (and (<= 1/2 (/ n d)) (< (/ n d) 1)))
  783.           (let ((n (int:round (integer-shift-left n p) d)))
  784.         (if (int:= n int:flonum-integer-limit)
  785.             (step4 (int:quotient n 2) (int:1+ e))
  786.             (step4 n e)))))
  787.        (step4
  788.         (lambda (n e)
  789.           (flo:denormalize (integer->flonum n #b11) e))))
  790.     (step1 n d))))
  791.   
  792.   (define (slow-method n d)
  793.     (if (int:positive? n)
  794.     (n>0 n d)
  795.     (flo:negate (n>0 (int:negate n) d))))
  796.  
  797.   (cond ((eq? n 0) flo:0)
  798.     ((integer->flonum n #b01)
  799.      => (lambda (n-exact-flonum)
  800.           (cond ((integer->flonum d #b01)
  801.              => (lambda (d-exact-flonum)
  802.               (flo:/ n-exact-flonum d-exact-flonum)))
  803.             (else (slow-method n d)))))
  804.     (else (slow-method n d))))
  805.  
  806. (define (int:->inexact n)
  807.   (if (fixnum? n)
  808.       (fixnum->flonum n) ;; 8.0 compiler open-codes when is N fixnum (by test)
  809.       (integer->flonum n #b10)))
  810.  
  811. (define (flo:significand-digits radix)
  812.   (cond ((int:= radix 10)
  813.      flo:significand-digits-base-10)
  814.     ((int:= radix 2)
  815.      flo:significand-digits-base-2)
  816.     (else
  817.      (int:+ 2
  818.         (flo:floor->exact
  819.          (flo:/ (int:->flonum flo:significand-digits-base-2)
  820.             (flo:/ (flo:log (int:->flonum radix))
  821.                    (flo:log 2.))))))))
  822.  
  823. (declare (integrate flo:integer?))
  824. (define (flo:integer? x)
  825.   (flo:= x (flo:round x)))
  826.  
  827. (define (flo:rationalize x e)
  828.   (flo:simplest-rational (flo:- x e) (flo:+ x e)))
  829.  
  830. (define (flo:simplest-rational x y)
  831.   ;; See comments at `rat:simplest-rational'.
  832.   (let ((x<y
  833.      (lambda (x y)
  834.        (define (loop x y)
  835.          (let ((fx (flo:floor x))
  836.            (fy (flo:floor y)))
  837.            (cond ((not (flo:< fx x)) fx)
  838.              ((flo:= fx fy)
  839.               (flo:+ fx
  840.                  (flo:/ flo:1
  841.                     (loop (flo:/ flo:1 (flo:- y fy))
  842.                       (flo:/ flo:1 (flo:- x fx))))))
  843.              (else (flo:+ fx flo:1)))))
  844.        (cond ((flo:positive? x) (loop x y))
  845.          ((flo:negative? y)
  846.           (flo:negate (loop (flo:negate y) (flo:negate x))))
  847.          (else flo:0)))))
  848.     (cond ((flo:< x y) (x<y x y))
  849.       ((flo:< y x) (x<y y x))
  850.       (else x))))
  851.  
  852. (define (flo:rationalize->exact x e)
  853.   (flo:simplest-exact-rational (flo:- x e) (flo:+ x e)))
  854.  
  855. (define (flo:simplest-exact-rational x y)
  856.   ;; See comments at `rat:simplest-rational'.
  857.   (let ((x<y
  858.      (lambda (x y)
  859.        (define (loop x y)
  860.          (let ((fx (flo:floor x))
  861.            (fy (flo:floor y)))
  862.            (cond ((not (flo:< fx x))
  863.               (flo:->integer fx))
  864.              ((flo:= fx fy)
  865.               (rat:+ (flo:->integer fx)
  866.                  (rat:invert (loop (flo:/ flo:1 (flo:- y fy))
  867.                            (flo:/ flo:1 (flo:- x fx))))))
  868.              (else
  869.               (rat:1+ (flo:->integer fx))))))
  870.        (cond ((flo:positive? x) (loop x y))
  871.          ((flo:negative? y)
  872.           (rat:negate (loop (flo:negate y) (flo:negate x))))
  873.          (else 0)))))
  874.     (cond ((flo:< x y) (x<y x y))
  875.       ((flo:< y x) (x<y y x))
  876.       (else (flo:->rational x)))))
  877.  
  878. (define (flo:->rational x)
  879.   (with-values (lambda () (flo:normalize x))
  880.     (lambda (f e-p)
  881.       (let ((p flo:significand-digits-base-2))
  882.     (rat:* (flo:->integer (flo:denormalize f p))
  883.            (rat:expt 2 (int:- e-p p)))))))
  884.  
  885. (define (real:real? object)
  886.   (or (flonum? object)
  887.       (rat:rational? object)))
  888.  
  889. (define-integrable (real:0 exact?)
  890.   (if exact? 0 0.0))
  891.  
  892. (define (real:exact1= x)
  893.   (and (real:exact? x)
  894.        (real:= 1 x)))
  895.  
  896. (define (real:rational? x)
  897.   (or (flonum? x) (rat:rational? x)))
  898.  
  899. (define (real:integer? x)
  900.   (if (flonum? x) (flo:integer? x) ((copy rat:integer?) x)))
  901.  
  902. (define (real:exact? x)
  903.   (and (not (flonum? x))
  904.        (or (rat:rational? x)
  905.        (error:wrong-type-argument x false 'EXACT?))))
  906.  
  907. (define (real:zero? x)
  908.   (if (flonum? x) (flo:zero? x) ((copy rat:zero?) x)))
  909.  
  910. (define (real:exact0= x)
  911.   (and (not (flonum? x)) ((copy rat:zero?) x)))
  912.  
  913. (define (real:negative? x)
  914.   (if (flonum? x) (flo:negative? x) ((copy rat:negative?) x)))
  915.  
  916. (define (real:positive? x)
  917.   (if (flonum? x) (flo:positive? x) ((copy rat:positive?) x)))
  918.  
  919. (let-syntax
  920.     ((define-standard-unary
  921.        (macro (name flo:op rat:op)
  922.      `(DEFINE (,name X)
  923.         (IF (FLONUM? X)
  924.         (,flo:op X)
  925.         (,rat:op X))))))
  926.   (define-standard-unary real:1+ (lambda (x) (flo:+ x flo:1)) (copy rat:1+))
  927.   (define-standard-unary real:-1+ (lambda (x) (flo:- x flo:1)) (copy rat:-1+))
  928.   (define-standard-unary real:negate flo:negate (copy rat:negate))
  929.   (define-standard-unary real:invert (lambda (x) (flo:/ flo:1 x)) rat:invert)
  930.   (define-standard-unary real:abs flo:abs rat:abs)
  931.   (define-standard-unary real:square (lambda (x) (flo:* x x)) rat:square)
  932.   (define-standard-unary real:floor flo:floor rat:floor)
  933.   (define-standard-unary real:ceiling flo:ceiling rat:ceiling)
  934.   (define-standard-unary real:truncate flo:truncate rat:truncate)
  935.   (define-standard-unary real:round flo:round rat:round)
  936.   (define-standard-unary real:floor->exact flo:floor->exact rat:floor)
  937.   (define-standard-unary real:ceiling->exact flo:ceiling->exact rat:ceiling)
  938.   (define-standard-unary real:truncate->exact flo:truncate->exact rat:truncate)
  939.   (define-standard-unary real:round->exact flo:round->exact rat:round)
  940.   (define-standard-unary real:exact->inexact (lambda (x) x) rat:->inexact)
  941.   (define-standard-unary real:inexact->exact flo:->rational
  942.     (lambda (q)
  943.       (if (rat:rational? q)
  944.       q
  945.       (error:wrong-type-argument q false 'INEXACT->EXACT)))))
  946.  
  947. (let-syntax
  948.     ((define-standard-binary
  949.        (macro (name flo:op rat:op)
  950.      `(DEFINE (,name X Y)
  951.         (IF (FLONUM? X)
  952.         (IF (FLONUM? Y)
  953.             (,flo:op X Y)
  954.             (,flo:op X (RAT:->INEXACT Y)))
  955.         (IF (FLONUM? Y)
  956.             (,flo:op (RAT:->INEXACT X) Y)
  957.             (,rat:op X Y)))))))
  958.   (define-standard-binary real:+ flo:+ (copy rat:+))
  959.   (define-standard-binary real:- flo:- (copy rat:-))
  960.   (define-standard-binary real:rationalize
  961.     flo:rationalize
  962.     rat:rationalize)
  963.   (define-standard-binary real:rationalize->exact
  964.     flo:rationalize->exact
  965.     rat:rationalize)
  966.   (define-standard-binary real:simplest-rational
  967.     flo:simplest-rational
  968.     rat:simplest-rational)
  969.   (define-standard-binary real:simplest-exact-rational
  970.     flo:simplest-exact-rational
  971.     rat:simplest-rational))
  972.  
  973. (define (real:= x y)
  974.   (if (flonum? x)
  975.       (if (flonum? y)
  976.       (flo:= x y)
  977.       (rat:= (flo:->rational x) y))
  978.       (if (flonum? y)
  979.       (rat:= x (flo:->rational y))
  980.       ((copy rat:=) x y))))
  981.  
  982. (define (real:< x y)
  983.   (if (flonum? x)
  984.       (if (flonum? y)
  985.       (flo:< x y)
  986.       (rat:< (flo:->rational x) y))
  987.       (if (flonum? y)
  988.       (rat:< x (flo:->rational y))
  989.       ((copy rat:<) x y))))
  990.  
  991. (define (real:max x y)
  992.   (if (flonum? x)
  993.       (if (flonum? y)
  994.       (if (flo:< x y) y x)
  995.       (if (rat:< (flo:->rational x) y) (rat:->inexact y) x))
  996.       (if (flonum? y)
  997.       (if (rat:< x (flo:->rational y)) y (rat:->inexact x))
  998.       (if (rat:< x y) y x))))
  999.  
  1000. (define (real:min x y)
  1001.   (if (flonum? x)
  1002.       (if (flonum? y)
  1003.       (if (flo:< x y) x y)
  1004.       (if (rat:< (flo:->rational x) y) x (rat:->inexact y)))
  1005.       (if (flonum? y)
  1006.       (if (rat:< x (flo:->rational y)) (rat:->inexact x) y)
  1007.       (if (rat:< x y) x y))))
  1008.  
  1009. (define (real:* x y)
  1010.   (cond ((flonum? x)
  1011.      (cond ((flonum? y) (flo:* x y))
  1012.            ((rat:zero? y) y)
  1013.            (else (flo:* x (rat:->inexact y)))))
  1014.     ((rat:zero? x) x)
  1015.     ((flonum? y) (flo:* (rat:->inexact x) y))
  1016.     (else ((copy rat:*) x y))))
  1017.  
  1018. (define (real:/ x y)
  1019.   (cond ((flonum? x) (flo:/ x (if (flonum? y) y (rat:->inexact y))))
  1020.     ((flonum? y) (if (rat:zero? x) x (flo:/ (rat:->inexact x) y)))
  1021.     (else ((copy rat:/) x y))))
  1022.  
  1023. (define (real:even? n)
  1024.   ((copy int:even?)
  1025.    (if (flonum? n)
  1026.        (if (flo:integer? n)
  1027.        (flo:->integer n)
  1028.        (error:wrong-type-argument n false 'EVEN?))
  1029.        n)))
  1030.  
  1031. (let-syntax
  1032.     ((define-integer-binary
  1033.        (macro (name operator-name operator)
  1034.      (let ((flo->int
  1035.         (lambda (n)
  1036.           `(IF (FLO:INTEGER? ,n)
  1037.                (FLO:->INTEGER ,n)
  1038.                (ERROR:WRONG-TYPE-ARGUMENT ,n FALSE ',operator-name)))))
  1039.        `(DEFINE (,name N M)
  1040.           (IF (FLONUM? N)
  1041.           (INT:->INEXACT
  1042.            (,operator ,(flo->int 'N)
  1043.                   (IF (FLONUM? M)
  1044.                   ,(flo->int 'M)
  1045.                   M)))
  1046.           (IF (FLONUM? M)
  1047.               (INT:->INEXACT (,operator N ,(flo->int 'M)))
  1048.               (,operator N M))))))))
  1049.   (define-integer-binary real:quotient quotient int:quotient)
  1050.   (define-integer-binary real:remainder remainder int:remainder)
  1051.   (define-integer-binary real:modulo modulo int:modulo)
  1052.   (define-integer-binary real:integer-floor integer-floor int:floor)
  1053.   (define-integer-binary real:integer-ceiling integer-ceiling int:ceiling)
  1054.   (define-integer-binary real:integer-round integer-round int:round)
  1055.   (define-integer-binary real:divide integer-divide int:divide)
  1056.   (define-integer-binary real:gcd gcd int:gcd)
  1057.   (define-integer-binary real:lcm lcm int:lcm))
  1058.  
  1059. (let-syntax
  1060.     ((define-rational-unary
  1061.        (macro (name operator)
  1062.      `(DEFINE (,name Q)
  1063.         (IF (FLONUM? Q)
  1064.         (RAT:->INEXACT (,operator (FLO:->RATIONAL Q)))
  1065.         (,operator Q))))))
  1066.   (define-rational-unary real:numerator rat:numerator)
  1067.   (define-rational-unary real:denominator rat:denominator))
  1068.  
  1069. (let-syntax
  1070.     ((define-transcendental-unary
  1071.        (macro (name hole? hole-value function)
  1072.      `(DEFINE (,name X)
  1073.         (IF (,hole? X)
  1074.         ,hole-value
  1075.         (,function (REAL:->INEXACT X)))))))
  1076.   (define-transcendental-unary real:exp real:exact0= 1 flo:exp)
  1077.   (define-transcendental-unary real:log real:exact1= 0 flo:log)
  1078.   (define-transcendental-unary real:sin real:exact0= 0 flo:sin)
  1079.   (define-transcendental-unary real:cos real:exact0= 1 flo:cos)
  1080.   (define-transcendental-unary real:tan real:exact0= 0 flo:tan)
  1081.   (define-transcendental-unary real:asin real:exact0= 0 flo:asin)
  1082.   (define-transcendental-unary real:acos real:exact1= 0 flo:acos)
  1083.   (define-transcendental-unary real:atan real:exact0= 0 flo:atan))
  1084.  
  1085. (define (real:atan2 y x)
  1086.   (if (and (real:exact0= y)
  1087.        (real:exact? x))
  1088.       (if (real:negative? x) rec:pi 0)
  1089.       (flo:atan2 (real:->inexact y) (real:->inexact x))))
  1090.  
  1091. (define (rat:sqrt x)
  1092.   (let ((guess (flo:sqrt (rat:->inexact x))))
  1093.     (if (int:integer? x)
  1094.     (let ((n (flo:round->exact guess)))
  1095.       (if (int:= x (int:* n n))
  1096.           n
  1097.           guess))
  1098.     (let ((q (flo:->rational guess)))
  1099.       (if (rat:= x (rat:square q))
  1100.           q
  1101.           guess)))))
  1102.  
  1103. (define (real:sqrt x)
  1104.   (if (flonum? x) (flo:sqrt x) (rat:sqrt x)))
  1105.  
  1106. (define (real:->inexact x)
  1107.   (if (flonum? x)
  1108.       x
  1109.       (rat:->inexact x)))
  1110.  
  1111. (define (real:->string x radix)
  1112.   (if (flonum? x)
  1113.       (flo:->string x radix)
  1114.       (rat:->string x radix)))
  1115.  
  1116. (define (real:expt x y)
  1117.   (let ((general-case
  1118.      (lambda (x y)
  1119.        (cond ((flo:zero? y) flo:1)
  1120.          ((flo:zero? x)
  1121.           (if (flo:positive? y)
  1122.               x
  1123.               (error:divide-by-zero 'EXPT (list x y))))
  1124.          ((and (flo:negative? x)
  1125.                (not (flo:integer? y)))
  1126.           (error:bad-range-argument x 'EXPT))
  1127.          (else
  1128.           (flo:expt x y))))))
  1129.     (if (flonum? x)
  1130.     (cond ((flonum? y)
  1131.            (general-case x y))
  1132.           ((int:integer? y)
  1133.            (let ((exact-method
  1134.               (lambda (y)
  1135.             (if (int:= 1 y)
  1136.                 x
  1137.                 (let loop ((x x) (y y) (answer flo:1))
  1138.                   (let ((qr (int:divide y 2)))
  1139.                 (let ((x (flo:* x x))
  1140.                       (y (integer-divide-quotient qr))
  1141.                       (answer
  1142.                        (if (int:zero?
  1143.                         (integer-divide-remainder qr))
  1144.                        answer
  1145.                        (flo:* answer x))))
  1146.                   (if (int:= 1 y)
  1147.                       (flo:* answer x)
  1148.                       (loop x y answer)))))))))
  1149.          (cond ((int:positive? y) (exact-method y))
  1150.                ((int:negative? y)
  1151.             (flo:/ flo:1 (exact-method (int:negate y))))
  1152.                (else flo:1))))
  1153.           (else
  1154.            (general-case x (rat:->inexact y))))
  1155.     (cond ((flonum? y)
  1156.            (general-case (rat:->inexact x) y))
  1157.           ((int:integer? y)
  1158.            (rat:expt x y))
  1159.           ((and (rat:positive? x)
  1160.             (int:= 1 (rat:numerator y)))
  1161.            (let ((d (rat:denominator y)))
  1162.          (if (int:= 2 d)
  1163.              (rat:sqrt x)
  1164.              (let ((guess
  1165.                 (flo:expt (rat:->inexact x) (rat:->inexact y))))
  1166.                (let ((q
  1167.                   (if (int:integer? x)
  1168.                   (flo:round->exact guess)
  1169.                   (flo:->rational guess))))
  1170.              (if (rat:= x (rat:expt q d))
  1171.                  q
  1172.                  guess))))))
  1173.           (else
  1174.            (general-case (rat:->inexact x) (rat:->inexact y)))))))
  1175.  
  1176. (define (complex:complex? object)
  1177.   (or (recnum? object) ((copy real:real?) object)))
  1178.  
  1179. (define (complex:real? object)
  1180.   (if (recnum? object)
  1181.       (real:zero? (rec:imag-part object))
  1182.       ((copy real:real?) object)))
  1183.  
  1184. (define (complex:rational? object)
  1185.   (if (recnum? object)
  1186.       (and (real:zero? (rec:imag-part object))
  1187.        (real:rational? (rec:real-part object)))
  1188.       ((copy real:rational?) object)))
  1189.  
  1190. (define (complex:integer? object)
  1191.   (if (recnum? object)
  1192.       (and (real:zero? (rec:imag-part object))
  1193.        (real:integer? (rec:real-part object)))
  1194.       ((copy real:integer?) object)))
  1195.  
  1196. (define (complex:exact? z)
  1197.   (if (recnum? z)
  1198.       ((copy rec:exact?) z)
  1199.       ((copy real:exact?) z)))
  1200.  
  1201. (define (rec:exact? z)
  1202.   (and (real:exact? (rec:real-part z))
  1203.        (real:exact? (rec:imag-part z))))
  1204.  
  1205. (define (complex:real-arg name x)
  1206.   (if (recnum? x) (rec:real-arg name x) x))
  1207.  
  1208. (define (rec:real-arg name x)
  1209.   (if (not (real:zero? (rec:imag-part x)))
  1210.       (error:wrong-type-argument x false name))
  1211.   (rec:real-part x))
  1212.  
  1213. (define (complex:= z1 z2)
  1214.   (if (recnum? z1)
  1215.       (if (recnum? z2)
  1216.       (and (real:= (rec:real-part z1) (rec:real-part z2))
  1217.            (real:= (rec:imag-part z1) (rec:imag-part z2)))
  1218.       (and (real:zero? (rec:imag-part z1))
  1219.            (real:= (rec:real-part z1) z2)))
  1220.       (if (recnum? z2)
  1221.       (and (real:zero? (rec:imag-part z2))
  1222.            (real:= z1 (rec:real-part z2)))
  1223.       ((copy real:=) z1 z2))))
  1224.  
  1225. (define (complex:< x y)
  1226.   (if (recnum? x)
  1227.       (if (recnum? y)
  1228.       (real:< (rec:real-arg '< x) (rec:real-arg '< y))
  1229.       (real:< (rec:real-arg '< x) y))
  1230.       (if (recnum? y)
  1231.       (real:< x (rec:real-arg '< y))
  1232.       ((copy real:<) x y))))
  1233.  
  1234. (define (complex:> x y)
  1235.   (complex:< y x))
  1236.  
  1237. (define (complex:zero? z)
  1238.   (if (recnum? z)
  1239.       (and (real:zero? (rec:real-part z))
  1240.        (real:zero? (rec:imag-part z)))
  1241.       ((copy real:zero?) z)))
  1242.  
  1243. (define (complex:positive? x)
  1244.   (if (recnum? x)
  1245.       (real:positive? (rec:real-arg 'POSITIVE? x))
  1246.       ((copy real:positive?) x)))
  1247.  
  1248. (define (complex:negative? x)
  1249.   (if (recnum? x)
  1250.       (real:negative? (rec:real-arg 'NEGATIVE? x))
  1251.       ((copy real:negative?) x)))
  1252.  
  1253. (define (complex:even? x)
  1254.   (if (recnum? x) (real:even? (rec:real-arg 'EVEN? x)) ((copy real:even?) x)))
  1255.  
  1256. (define (complex:max x y)
  1257.   (if (recnum? x)
  1258.       (if (recnum? y)
  1259.       (real:max (rec:real-arg 'MAX x) (rec:real-arg 'MAX y))
  1260.       (real:max (rec:real-arg 'MAX x) y))
  1261.       (if (recnum? y)
  1262.       (real:max x (rec:real-arg 'MAX y))
  1263.       ((copy real:max) x y))))
  1264.  
  1265. (define (complex:min x y)
  1266.   (if (recnum? x)
  1267.       (if (recnum? y)
  1268.       (real:min (rec:real-arg 'MIN x) (rec:real-arg 'MIN y))
  1269.       (real:min (rec:real-arg 'MIN x) y))
  1270.       (if (recnum? y)
  1271.       (real:min x (rec:real-arg 'MIN y))
  1272.       ((copy real:min) x y))))
  1273.  
  1274. (define (complex:+ z1 z2)
  1275.   (if (recnum? z1)
  1276.       (if (recnum? z2)
  1277.       (complex:%make-rectangular
  1278.        (real:+ (rec:real-part z1) (rec:real-part z2))
  1279.        (real:+ (rec:imag-part z1) (rec:imag-part z2)))
  1280.       (make-recnum (real:+ (rec:real-part z1) z2)
  1281.                (rec:imag-part z1)))
  1282.       (if (recnum? z2)
  1283.       (make-recnum (real:+ z1 (rec:real-part z2))
  1284.                (rec:imag-part z2))
  1285.       ((copy real:+) z1 z2))))
  1286.  
  1287. (define (complex:1+ z)
  1288.   (if (recnum? z)
  1289.       (make-recnum (real:1+ (rec:real-part z)) (rec:imag-part z))
  1290.       ((copy real:1+) z)))
  1291.  
  1292. (define (complex:-1+ z)
  1293.   (if (recnum? z)
  1294.       (make-recnum (real:-1+ (rec:real-part z)) (rec:imag-part z))
  1295.       ((copy real:-1+) z)))
  1296.  
  1297. (define (complex:* z1 z2)
  1298.   (if (recnum? z1)
  1299.       (if (recnum? z2)
  1300.       (let ((z1r (rec:real-part z1))
  1301.         (z1i (rec:imag-part z1))
  1302.         (z2r (rec:real-part z2))
  1303.         (z2i (rec:imag-part z2)))
  1304.         (complex:%make-rectangular
  1305.          (real:- (real:* z1r z2r) (real:* z1i z2i))
  1306.          (real:+ (real:* z1r z2i) (real:* z1i z2r))))
  1307.       (complex:%make-rectangular (real:* (rec:real-part z1) z2)
  1308.                      (real:* (rec:imag-part z1) z2)))
  1309.       (if (recnum? z2)
  1310.       (complex:%make-rectangular (real:* z1 (rec:real-part z2))
  1311.                      (real:* z1 (rec:imag-part z2)))
  1312.       ((copy real:*) z1 z2))))
  1313.  
  1314. (define (complex:+i* z)
  1315.   (if (recnum? z)
  1316.       (complex:%make-rectangular (real:negate (rec:imag-part z))
  1317.                  (rec:real-part z))
  1318.       (complex:%make-rectangular 0 z)))
  1319.  
  1320. (define (complex:-i* z)
  1321.   (if (recnum? z)
  1322.       (complex:%make-rectangular (rec:imag-part z)
  1323.                  (real:negate (rec:real-part z)))
  1324.       (complex:%make-rectangular 0 (real:negate z))))
  1325.  
  1326. (define (complex:- z1 z2)
  1327.   (if (recnum? z1)
  1328.       (if (recnum? z2)
  1329.       (complex:%make-rectangular
  1330.        (real:- (rec:real-part z1) (rec:real-part z2))
  1331.        (real:- (rec:imag-part z1) (rec:imag-part z2)))
  1332.       (make-recnum (real:- (rec:real-part z1) z2)
  1333.                (rec:imag-part z1)))
  1334.       (if (recnum? z2)
  1335.       (make-recnum (real:- z1 (rec:real-part z2))
  1336.                (real:negate (rec:imag-part z2)))
  1337.       ((copy real:-) z1 z2))))
  1338.  
  1339. (define (complex:negate z)
  1340.   (if (recnum? z)
  1341.       (make-recnum (real:negate (rec:real-part z))
  1342.            (real:negate (rec:imag-part z)))
  1343.       ((copy real:negate) z)))
  1344.  
  1345. (define (complex:conjugate z)
  1346.   (cond ((recnum? z)
  1347.      (make-recnum (rec:real-part z)
  1348.               (real:negate (rec:imag-part z))))
  1349.     ((real:real? z)
  1350.      z)
  1351.     (else
  1352.      (error:wrong-type-argument z false 'CONJUGATE))))
  1353.  
  1354. (define (complex:/ z1 z2)
  1355.   (if (recnum? z1)
  1356.       (if (recnum? z2)
  1357.       (let ((z1r (rec:real-part z1))
  1358.         (z1i (rec:imag-part z1))
  1359.         (z2r (rec:real-part z2))
  1360.         (z2i (rec:imag-part z2)))
  1361.         (let ((d (real:+ (real:square z2r) (real:square z2i))))
  1362.           (complex:%make-rectangular
  1363.            (real:/ (real:+ (real:* z1r z2r) (real:* z1i z2i)) d)
  1364.            (real:/ (real:- (real:* z1i z2r) (real:* z1r z2i)) d))))
  1365.       (make-recnum (real:/ (rec:real-part z1) z2)
  1366.                (real:/ (rec:imag-part z1) z2)))
  1367.       (if (recnum? z2)
  1368.       (let ((z2r (rec:real-part z2))
  1369.         (z2i (rec:imag-part z2)))
  1370.         (let ((d (real:+ (real:square z2r) (real:square z2i))))
  1371.           (complex:%make-rectangular
  1372.            (real:/ (real:* z1 z2r) d)
  1373.            (real:/ (real:negate (real:* z1 z2i)) d))))
  1374.       ((copy real:/) z1 z2))))
  1375.  
  1376. (define (complex:invert z)
  1377.   (if (recnum? z)
  1378.       (let ((zr (rec:real-part z))
  1379.         (zi (rec:imag-part z)))
  1380.     (let ((d (real:+ (real:square zr) (real:square zi))))
  1381.       (make-recnum (real:/ zr d)
  1382.                (real:/ (real:negate zi) d))))
  1383.       ((copy real:invert) z)))
  1384.  
  1385. (define (complex:abs x)
  1386.   (if (recnum? x) (real:abs (rec:real-arg 'ABS x)) ((copy real:abs) x)))
  1387.  
  1388. (define (complex:quotient n d)
  1389.   (real:quotient (complex:real-arg 'QUOTIENT n)
  1390.          (complex:real-arg 'QUOTIENT d)))
  1391.  
  1392. (define (complex:remainder n d)
  1393.   (real:remainder (complex:real-arg 'REMAINDER n)
  1394.           (complex:real-arg 'REMAINDER d)))
  1395.  
  1396. (define (complex:modulo n d)
  1397.   (real:modulo (complex:real-arg 'MODULO n)
  1398.            (complex:real-arg 'MODULO d)))
  1399.  
  1400. (define (complex:integer-floor n d)
  1401.   (real:integer-floor (complex:real-arg 'INTEGER-FLOOR n)
  1402.               (complex:real-arg 'INTEGER-FLOOR d)))
  1403.  
  1404. (define (complex:integer-ceiling n d)
  1405.   (real:integer-ceiling (complex:real-arg 'INTEGER-CEILING n)
  1406.             (complex:real-arg 'INTEGER-CEILING d)))
  1407.  
  1408. (define (complex:integer-round n d)
  1409.   (real:integer-round (complex:real-arg 'INTEGER-ROUND n)
  1410.               (complex:real-arg 'INTEGER-ROUND d)))
  1411.  
  1412. (define (complex:divide n d)
  1413.   (real:divide (complex:real-arg 'DIVIDE n)
  1414.            (complex:real-arg 'DIVIDE d)))
  1415.  
  1416. (define (complex:gcd n m)
  1417.   (real:gcd (complex:real-arg 'GCD n)
  1418.         (complex:real-arg 'GCD m)))
  1419.  
  1420. (define (complex:lcm n m)
  1421.   (real:lcm (complex:real-arg 'LCM n)
  1422.         (complex:real-arg 'LCM m)))
  1423.  
  1424. (define (complex:numerator q)
  1425.   (real:numerator (complex:real-arg 'NUMERATOR q)))
  1426.  
  1427. (define (complex:denominator q)
  1428.   (real:denominator (complex:real-arg 'DENOMINATOR q)))
  1429.  
  1430. (define (complex:floor x)
  1431.   (if (recnum? x)
  1432.       (real:floor (rec:real-arg 'FLOOR x))
  1433.       ((copy real:floor) x)))
  1434.  
  1435. (define (complex:ceiling x)
  1436.   (if (recnum? x)
  1437.       (real:ceiling (rec:real-arg 'CEILING x))
  1438.       ((copy real:ceiling) x)))
  1439.  
  1440. (define (complex:truncate x)
  1441.   (if (recnum? x)
  1442.       (real:truncate (rec:real-arg 'TRUNCATE x))
  1443.       ((copy real:truncate) x)))
  1444.  
  1445. (define (complex:round x)
  1446.   (if (recnum? x)
  1447.       (real:round (rec:real-arg 'ROUND x))
  1448.       ((copy real:round) x)))
  1449.  
  1450. (define (complex:floor->exact x)
  1451.   (if (recnum? x)
  1452.       (real:floor->exact (rec:real-arg 'FLOOR->EXACT x))
  1453.       ((copy real:floor->exact) x)))
  1454.  
  1455. (define (complex:ceiling->exact x)
  1456.   (if (recnum? x)
  1457.       (real:ceiling->exact (rec:real-arg 'CEILING->EXACT x))
  1458.       ((copy real:ceiling->exact) x)))
  1459.  
  1460. (define (complex:truncate->exact x)
  1461.   (if (recnum? x)
  1462.       (real:truncate->exact (rec:real-arg 'TRUNCATE->EXACT x))
  1463.       ((copy real:truncate->exact) x)))
  1464.  
  1465. (define (complex:round->exact x)
  1466.   (if (recnum? x)
  1467.       (real:round->exact (rec:real-arg 'ROUND->EXACT x))
  1468.       ((copy real:round->exact) x)))
  1469.  
  1470. (define (complex:rationalize x e)
  1471.   (real:rationalize (complex:real-arg 'RATIONALIZE x)
  1472.             (complex:real-arg 'RATIONALIZE e)))
  1473.  
  1474. (define (complex:rationalize->exact x e)
  1475.   (real:rationalize->exact (complex:real-arg 'RATIONALIZE x)
  1476.                (complex:real-arg 'RATIONALIZE e)))
  1477.  
  1478. (define (complex:simplest-rational x y)
  1479.   (real:simplest-rational (complex:real-arg 'SIMPLEST-RATIONAL x)
  1480.               (complex:real-arg 'SIMPLEST-RATIONAL y)))
  1481.  
  1482. (define (complex:simplest-exact-rational x y)
  1483.   (real:simplest-exact-rational (complex:real-arg 'SIMPLEST-RATIONAL x)
  1484.                 (complex:real-arg 'SIMPLEST-RATIONAL y)))
  1485.  
  1486. (define (complex:exp z)
  1487.   (if (recnum? z)
  1488.       (complex:%make-polar (real:exp (rec:real-part z))
  1489.                (rec:imag-part z))
  1490.       ((copy real:exp) z)))
  1491.  
  1492. (define (complex:log z)
  1493.   (cond ((recnum? z)
  1494.      (complex:%make-rectangular (real:log (complex:magnitude z))
  1495.                     (complex:angle z)))
  1496.     ((real:negative? z)
  1497.      (make-recnum (real:log (real:negate z)) rec:pi))
  1498.     (else
  1499.      ((copy real:log) z))))
  1500.  
  1501. (define (complex:sin z)
  1502.   (if (recnum? z)
  1503.       (complex:/ (let ((iz (complex:+i* z)))
  1504.            (complex:- (complex:exp iz)
  1505.                   (complex:exp (complex:negate iz))))
  1506.          +2i)
  1507.       ((copy real:sin) z)))
  1508.  
  1509. (define (complex:cos z)
  1510.   (if (recnum? z)
  1511.       (complex:/ (let ((iz (complex:+i* z)))
  1512.            (complex:+ (complex:exp iz)
  1513.                   (complex:exp (complex:negate iz))))
  1514.          2)
  1515.       ((copy real:cos) z)))
  1516.  
  1517. (define (complex:tan z)
  1518.   (if (recnum? z)
  1519.       (complex:-i*
  1520.        (let ((iz (complex:+i* z)))
  1521.      (let ((e+iz (complex:exp iz))
  1522.            (e-iz (complex:exp (complex:negate iz))))
  1523.        (complex:/ (complex:- e+iz e-iz)
  1524.               (complex:+ e+iz e-iz)))))
  1525.       ((copy real:tan) z)))
  1526.  
  1527. ;;; Complex arguments -- ASIN
  1528. ;;;   The danger in the complex case happens for large y when
  1529. ;;;     z = iy.  In this case iz + sqrt(1-z^2) --> -y + y.
  1530. ;;;   A clever way out of this difficulty uses symmetry to always
  1531. ;;;     take the benevolent branch of the square root.
  1532. ;;;   That is, make iz and sqrt(1-z^2) always end up in the same
  1533. ;;;     quadrant so catastrophic cancellation cannot occur.
  1534. ;;;  This is ensured if z is in quadrants III or IV.
  1535.  
  1536. (define (complex:asin z)
  1537.   (let ((safe-case
  1538.      (lambda (z)
  1539.        (complex:-i*
  1540.         (complex:log
  1541.          (complex:+ (complex:+i* z)
  1542.             (complex:sqrt (complex:- 1 (complex:* z z)))))))))
  1543.     (let ((unsafe-case
  1544.        (lambda (z)
  1545.          (complex:negate (safe-case (complex:negate z))))))
  1546.       (cond ((recnum? z)
  1547.          (if (let ((imag (rec:imag-part z)))
  1548.            (or (real:positive? imag)    ;get out of Q I and II
  1549.                (and (real:zero? imag)    ;and stay off negative reals
  1550.                 (real:negative? (rec:real-part z)))))
  1551.          (unsafe-case z)
  1552.          (safe-case z)))
  1553.         ((real:< z -1)
  1554.          (unsafe-case z))
  1555.         ((real:< 1 z)
  1556.          (safe-case z))
  1557.         (else
  1558.          ((copy real:asin) z))))))
  1559.  
  1560. (define (complex:acos z)
  1561.   (if (or (recnum? z)
  1562.       (real:< z -1)
  1563.       (real:< 1 z))
  1564.       (complex:-i*
  1565.        (complex:log
  1566.     (complex:+ z
  1567.            (complex:+i*
  1568.             (complex:sqrt (complex:- 1 (complex:* z z)))))))
  1569.       ((copy real:acos) z)))
  1570.  
  1571. (define (complex:atan z)
  1572.   (if (recnum? z)
  1573.       (rec:atan z)
  1574.       ((copy real:atan) z)))
  1575.  
  1576. (define (complex:atan2 y x)
  1577.   (let ((rec-case
  1578.      (lambda (y x)
  1579.        (rec:atan (make-recnum (real:exact->inexact x)
  1580.                   (real:exact->inexact y))))))
  1581.     (cond ((recnum? y)
  1582.        (rec-case (rec:real-arg 'ATAN y) (complex:real-arg 'ATAN x)))
  1583.       ((recnum? x)
  1584.        (rec-case y (rec:real-arg 'ATAN x)))
  1585.       (else
  1586.        ((copy real:atan2) y x)))))
  1587.  
  1588. (define (rec:atan z)
  1589.   (complex:/ (let ((iz (complex:+i* z)))
  1590.            (complex:- (complex:log (complex:1+ iz))
  1591.               (complex:log (complex:- 1 iz))))
  1592.          +2i))
  1593.  
  1594. (define (complex:angle z)
  1595.   (cond ((recnum? z)
  1596.      (if (and (real:zero? (rec:real-part z))
  1597.           (real:zero? (rec:imag-part z)))
  1598.          (real:0 (complex:exact? z))
  1599.          (real:atan2 (rec:imag-part z) (rec:real-part z))))
  1600.     ((real:negative? z) rec:pi)
  1601.     (else (real:0 (real:exact? z)))))
  1602.  
  1603. (define (complex:magnitude z)
  1604.   (if (recnum? z)
  1605.       (let ((ar (real:abs (rec:real-part z)))
  1606.         (ai (real:abs (rec:imag-part z))))
  1607.     (let ((v (real:max ar ai))
  1608.           (w (real:min ar ai)))
  1609.       (if (real:zero? v)
  1610.           v
  1611.           (real:* v (real:sqrt (real:1+ (real:square (real:/ w v))))))))
  1612.       (real:abs z)))
  1613.  
  1614. (define (complex:sqrt z)
  1615.   (cond ((recnum? z)
  1616.      (complex:%make-polar (real:sqrt (complex:magnitude z))
  1617.                   (real:/ (complex:angle z) 2)))
  1618.     ((real:negative? z)
  1619.      (complex:%make-rectangular 0 (real:sqrt (real:negate z))))
  1620.     (else
  1621.      ((copy real:sqrt) z))))
  1622.  
  1623. (define (complex:expt z1 z2)
  1624.   (let ((general-case
  1625.      (lambda ()
  1626.        (complex:exp (complex:* (complex:log z1) z2)))))
  1627.     (cond ((recnum? z1)
  1628.        (if (and (rec:exact? z1)
  1629.             (int:integer? z2))
  1630.            (let ((exact-method
  1631.               (lambda (z2)
  1632.             (if (int:= 1 z2)
  1633.                 z1
  1634.                 (let loop ((z1 z1) (z2 z2) (answer 1))
  1635.                   (let ((qr (int:divide z2 2)))
  1636.                 (let ((z1 (complex:* z1 z1))
  1637.                       (z2 (integer-divide-quotient qr))
  1638.                       (answer
  1639.                        (if (int:zero?
  1640.                         (integer-divide-remainder qr))
  1641.                        answer
  1642.                        (complex:* answer z1))))
  1643.                   (if (int:= 1 z2)
  1644.                       (complex:* answer z1)
  1645.                       (loop z1 z2 answer)))))))))
  1646.          (cond ((int:positive? z2) (exact-method z2))
  1647.                ((int:negative? z2)
  1648.             (complex:/ 1 (exact-method (int:negate z2))))
  1649.                (else 1)))
  1650.            (general-case)))
  1651.       ((or (recnum? z2)
  1652.            (and (real:negative? z1)
  1653.             (not (real:integer? z2))))
  1654.        (general-case))
  1655.       (else
  1656.        (real:expt z1 z2)))))
  1657.  
  1658. (define (complex:make-rectangular real imag)
  1659.   (let ((check-arg
  1660.      (lambda (x)
  1661.        (if (recnum? x)
  1662.            (rec:real-arg 'MAKE-RECTANGULAR x)
  1663.            (begin
  1664.          (if (not (real:real? x))
  1665.              (error:wrong-type-argument x false 'MAKE-RECTANGULAR))
  1666.          x)))))
  1667.     ((copy complex:%make-rectangular) (check-arg real) (check-arg imag))))
  1668.  
  1669. (define (complex:make-polar real imag)
  1670.   ((copy complex:%make-polar) (complex:real-arg 'MAKE-POLAR real)
  1671.                   (complex:real-arg 'MAKE-POLAR imag)))
  1672.  
  1673. (define (complex:%make-rectangular real imag)
  1674.   (if (real:exact0= imag)
  1675.       real
  1676.       (make-recnum real imag)))
  1677.  
  1678. (define (complex:%make-polar magnitude angle)
  1679.   (complex:%make-rectangular (real:* magnitude (real:cos angle))
  1680.                  (real:* magnitude (real:sin angle))))
  1681.  
  1682. (define (complex:real-part z)
  1683.   (cond ((recnum? z) (rec:real-part z))
  1684.     ((real:real? z) z)
  1685.     (else (error:wrong-type-argument z false 'REAL-PART))))
  1686.  
  1687. (define (complex:imag-part z)
  1688.   (cond ((recnum? z) (rec:imag-part z))
  1689.     ((real:real? z) 0)
  1690.     (else (error:wrong-type-argument z false 'IMAG-PART))))
  1691.  
  1692. (define (complex:exact->inexact z)
  1693.   (if (recnum? z)
  1694.       (complex:%make-rectangular (real:exact->inexact (rec:real-part z))
  1695.                  (real:exact->inexact (rec:imag-part z)))
  1696.       ((copy real:exact->inexact) z)))
  1697.  
  1698. (define (complex:inexact->exact z)
  1699.   (if (recnum? z)
  1700.       (complex:%make-rectangular (real:inexact->exact (rec:real-part z))
  1701.                  (real:inexact->exact (rec:imag-part z)))
  1702.       ((copy real:inexact->exact) z)))
  1703.  
  1704. (define (complex:->string z radix)
  1705.   (if (recnum? z)
  1706.       (string-append
  1707.        (let ((r (rec:real-part z)))
  1708.      (if (real:exact0= r)
  1709.          ""
  1710.          (real:->string r radix)))
  1711.        (let ((i (rec:imag-part z))
  1712.          (positive-case
  1713.           (lambda (i)
  1714.         (if (real:exact1= i)
  1715.             ""
  1716.             (real:->string i radix)))))
  1717.      (if (real:negative? i)
  1718.          (string-append "-" (positive-case (real:negate i)))
  1719.          (string-append "+" (positive-case i))))
  1720.        (if imaginary-unit-j? "j" "i"))
  1721.       (real:->string z radix)))
  1722.  
  1723. (define imaginary-unit-j? #f)
  1724.  
  1725. (define number? complex:complex?)
  1726. (define complex? complex:complex?)
  1727. (define real? complex:real?)
  1728. (define rational? complex:rational?)
  1729. (define integer? complex:integer?)
  1730. (define exact? complex:exact?)
  1731. (define exact-rational? rat:rational?)
  1732. (define exact-integer? int:integer?)
  1733.  
  1734. (define (exact-nonnegative-integer? object)
  1735.   (and (int:integer? object)
  1736.        (not (int:negative? object))))
  1737.  
  1738. (define (inexact? z)
  1739.   (not (complex:exact? z)))
  1740.  
  1741. ;; Replaced with arity-dispatched version in INITIALIZE-PACKAGE!
  1742.  
  1743. (define =)
  1744. (define <)
  1745. (define >)
  1746. (define <=)
  1747. (define >=)
  1748.  
  1749. (define (reduce-comparator binary-comparator numbers procedure)
  1750.   (cond ((null? numbers)
  1751.      true)
  1752.     ((null? (cdr numbers))
  1753.      (if (not (complex:complex? (car numbers)))
  1754.          (error:wrong-type-argument (car numbers) false procedure))
  1755.      true)
  1756.     (else
  1757.      (let loop ((x (car numbers)) (rest (cdr numbers)))
  1758.        (or (null? rest)
  1759.            (let ((y (car rest)))
  1760.          (and (binary-comparator x y)
  1761.               (loop y (cdr rest)))))))))
  1762.  
  1763. (define zero? complex:zero?)
  1764. (define positive? complex:positive?)
  1765. (define negative? complex:negative?)
  1766.  
  1767. (define (odd? n)
  1768.   (not (complex:even? n)))
  1769.  
  1770. (define even? complex:even?)
  1771.  
  1772.  
  1773. ;; Replaced with arity-dispatched version in INITIALIZE-PACKAGE!
  1774.  
  1775. (define +)
  1776. (define *)
  1777. (define -)
  1778. (define /)
  1779.  
  1780. (define max)
  1781. (define min)
  1782.  
  1783. (define (reduce-max/min max/min x1 xs procedure)
  1784.   (if (null? xs)
  1785.       (begin
  1786.     (if (not (complex:complex? x1))
  1787.         (error:wrong-type-argument x1 false procedure))
  1788.     x1)
  1789.       (let loop ((x1 x1) (xs xs))
  1790.     (let ((x1 (max/min x1 (car xs)))
  1791.           (xs (cdr xs)))
  1792.       (if (null? xs)
  1793.           x1
  1794.           (loop x1 xs))))))
  1795.  
  1796. (define 1+ complex:1+)
  1797. (define -1+ complex:-1+)
  1798.  
  1799. (define conjugate complex:conjugate)
  1800.  
  1801. (define abs complex:abs)
  1802.  
  1803. ;;; The following three procedures were originally just renamings of
  1804. ;;; their COMPLEX: equivalents.  They have been rewritten this way to
  1805. ;;; cause the compiler to generate better code for them.
  1806.  
  1807. (define (quotient n d)
  1808.   ((ucode-primitive quotient 2) n d))
  1809.  
  1810. (define (remainder n d)
  1811.   ((ucode-primitive remainder 2) n d))
  1812.  
  1813. (define (modulo n d)
  1814.   (let ((r ((ucode-primitive remainder 2) n d)))
  1815.     (if (or (zero? r)
  1816.         (if (negative? n)
  1817.         (negative? d)
  1818.         (not (negative? d))))
  1819.     r
  1820.     (+ r d))))
  1821.  
  1822. (define integer-floor complex:integer-floor)
  1823. (define integer-ceiling complex:integer-ceiling)
  1824. (define integer-truncate complex:quotient)
  1825. (define integer-round complex:integer-round)
  1826. (define integer-divide complex:divide)
  1827. (define-integrable integer-divide-quotient car)
  1828. (define-integrable integer-divide-remainder cdr)
  1829.  
  1830. (define (gcd . integers)
  1831.   (fold-left complex:gcd 0 integers))
  1832.  
  1833. (define (lcm . integers)
  1834.   (fold-left complex:lcm 1 integers))
  1835.  
  1836. (define numerator complex:numerator)
  1837. (define denominator complex:denominator)
  1838. (define floor complex:floor)
  1839. (define ceiling complex:ceiling)
  1840. (define truncate complex:truncate)
  1841. (define round complex:round)
  1842. (define floor->exact complex:floor->exact)
  1843. (define ceiling->exact complex:ceiling->exact)
  1844. (define truncate->exact complex:truncate->exact)
  1845. (define round->exact complex:round->exact)
  1846. (define rationalize complex:rationalize)
  1847. (define rationalize->exact complex:rationalize->exact)
  1848. (define simplest-rational complex:simplest-rational)
  1849. (define simplest-exact-rational complex:simplest-exact-rational)
  1850. (define exp complex:exp)
  1851. (define log complex:log)
  1852. (define sin complex:sin)
  1853. (define cos complex:cos)
  1854. (define tan complex:tan)
  1855. (define asin complex:asin)
  1856. (define acos complex:acos)
  1857.  
  1858. (define (atan z #!optional x)
  1859.   (if (default-object? x)
  1860.       (complex:atan z)
  1861.       (complex:atan2 z x)))
  1862.  
  1863. (define sqrt complex:sqrt)
  1864. (define expt complex:expt)
  1865. (define make-rectangular complex:make-rectangular)
  1866. (define make-polar complex:make-polar)
  1867. (define real-part complex:real-part)
  1868. (define imag-part complex:imag-part)
  1869. (define magnitude complex:magnitude)
  1870. (define angle complex:angle)
  1871. (define exact->inexact complex:exact->inexact)
  1872. (define inexact->exact complex:inexact->exact)
  1873.  
  1874. (define (square z)
  1875.   (complex:* z z))
  1876.  
  1877. (define (number->string z #!optional radix)
  1878.   (complex:->string
  1879.    z
  1880.    (cond ((default-object? radix)
  1881.       10)
  1882.      ((and (exact-integer? radix)
  1883.            (<= 2 radix 36))
  1884.       radix)
  1885.      ((and (pair? radix)
  1886.            (eq? (car radix) 'HEUR)
  1887.            (list? radix))
  1888.       (parse-format-tail (cdr radix)))
  1889.      (else
  1890.       (error:bad-range-argument radix 'NUMBER->STRING)))))
  1891.  
  1892. (define (parse-format-tail tail)
  1893.   (let loop
  1894.       ((tail tail)
  1895.        (exactness-expressed false)
  1896.        (radix false)
  1897.        (radix-expressed false))
  1898.     (if (null? tail)
  1899.     (case radix ((B) 2) ((O) 8) ((#F D) 10) ((X) 16))
  1900.     (let ((modifier (car tail))
  1901.           (tail (cdr tail)))
  1902.       (let ((specify-modifier
  1903.          (lambda (old)
  1904.            (if old
  1905.                (error "Respecification of format modifier"
  1906.                   (cadr modifier)))
  1907.            (cadr modifier))))
  1908.         (cond ((and (pair? modifier)
  1909.             (eq? (car modifier) 'EXACTNESS)
  1910.             (pair? (cdr modifier))
  1911.             (memq (cadr modifier) '(E S))
  1912.             (null? (cddr modifier)))
  1913.            (if (eq? (cadr modifier) 'E)
  1914.                (warn "NUMBER->STRING: ignoring exactness modifier"
  1915.                  modifier))
  1916.            (loop tail
  1917.              (specify-modifier exactness-expressed)
  1918.              radix
  1919.              radix-expressed))
  1920.           ((and (pair? modifier)
  1921.             (eq? (car modifier) 'RADIX)
  1922.             (pair? (cdr modifier))
  1923.             (memq (cadr modifier) '(B O D X))
  1924.             (or (null? (cddr modifier))
  1925.                 (and (pair? (cddr modifier))
  1926.                  (memq (caddr modifier) '(E S))
  1927.                  (null? (cdddr modifier)))))
  1928.            (if (and (pair? (cddr modifier))
  1929.                 (eq? (caddr modifier) 'E))
  1930.                (warn
  1931.             "NUMBER->STRING: ignoring radix expression modifier"
  1932.             modifier))
  1933.            (loop tail
  1934.              exactness-expressed
  1935.              (specify-modifier radix)
  1936.              (if (null? (cddr modifier)) 'E (caddr modifier))))
  1937.           (else
  1938.            (error "Illegal format modifier" modifier))))))))