home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / garith.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  21.5 KB  |  758 lines  |  [TEXT/????]
  1. ;;; $Header: garith.scm,v 1.1 87/08/25 17:40:01 GMT gjs Exp $
  2. ;;;; Generic Complex and Rational Arithmetic with Operators
  3.  
  4. ;;; Requires macro file GENERAL\SYN.SCM
  5.  
  6. ;;; The definitions follow the order in
  7. ;;; Revised^3 Report on the Algorithmic Language Scheme.
  8.  
  9. ;;; Primitives are inherited from underlying Scheme.
  10.  
  11. (move-definition number? s-number?)
  12. (move-definition real? s-real?)
  13.  
  14. (move-definition zero? s-zero?)
  15.  
  16. (move-definition = s-=)
  17.  
  18. (move-definition 1+ s-1+)
  19. (move-definition -1+ s--1+)
  20.  
  21. (move-definition + s-+)
  22. (move-definition - s--)
  23. (move-definition * s-*)
  24. (move-definition / s-/)
  25.  
  26. (move-definition sin s-sin)
  27. (move-definition cos s-cos)
  28. (move-definition atan s-atan)
  29. (move-definition acos s-acos)
  30. (move-definition asin s-asin)
  31.  
  32. (move-definition sqrt s-sqrt)
  33. (move-definition exp s-exp)
  34. (move-definition log s-log)
  35. (move-definition expt s-expt)
  36.  
  37. (move-definition abs s-abs)
  38. (move-definition positive? s-positive?)
  39. (move-definition negative? s-negative?)
  40. (move-definition < s-<)
  41. (move-definition > s->)
  42. (move-definition <= s-<=)
  43. (move-definition >= s->=)
  44.  
  45. ;;; Representation of rational numbers
  46.  
  47. (define-integrable construct-rational
  48.   (lambda (n d) (list* '*rat* n d)))
  49.  
  50. (define-integrable &numer cadr)
  51. (define-integrable &denom cddr)
  52.  
  53. (define-integrable &rational?
  54.   (lambda (q)
  55.     (and (pair? q)
  56.          (eq? (car q) '*rat*))))
  57.  
  58. (define-integrable rational?
  59.   (lambda (q)
  60.     (or (integer? q)
  61.         (&rational? q))))
  62.  
  63. (define (make-rational n d)        ;gcd already done
  64.   (cond ((s-= d 0) (error "zero denominator -- rational"))
  65.     ((s-= n 0) 0)
  66.     ((s-= d 1) n)
  67.     ((s-< d 0) (make-rational (s-- n) (s-- d)))
  68.     (else (construct-rational n d))))
  69.  
  70. ;;; The following procedures try to remove common factors as early as
  71. ;;;  possible to reduce the cost of manipulation of big integers.
  72.  
  73. (define (+_rat x y)
  74.   (let ((dx (&denom x)) (dy (&denom y)))
  75.     (let ((d1 (gcd dx dy)))
  76.       (let ((dx/d1 (quotient dx d1)) (dy/d1 (quotient dy d1)))
  77.     (if (s-= 1 d1)
  78.     (make-rational (s-+ (s-* (&numer x) dy)
  79.                     (s-* dx (&numer y)))
  80.                (s-* dx dy))
  81.         (let ((tem (s-+ (s-* dy/d1 (&numer x))
  82.                 (s-* dx/d1 (&numer y)))))
  83.           (if (s-= tem 0)
  84.           0
  85.           (let ((d2 (gcd tem d1)))
  86.             (make-rational (quotient tem d2)
  87.                    (s-* dx/d1
  88.                         (quotient dy d2)))))))))))
  89.  
  90. (define (-_rat x y)
  91.   (+_rat x (construct-rational (s-- (&numer y)) (&denom y))))
  92.  
  93. (define (*_rat x y)
  94.   (let ((nx (&numer x)) (ny (&numer y))
  95.         (dx (&denom x)) (dy (&denom y)))
  96.     (let ((gnxdy (gcd nx dy)) (gnydx (gcd ny dx)))
  97.       (let ((nx/gcd (quotient nx gnxdy)) (dy/gcd (quotient dy gnxdy))
  98.         (ny/gcd (quotient ny gnydx)) (dx/gcd (quotient dx gnydx)))
  99.     (make-rational (s-* nx/gcd ny/gcd)
  100.                (s-* dx/gcd dy/gcd))))))
  101.  
  102. (define (/_rat x y)
  103.   (*_rat x (construct-rational (&denom y) (&numer y))))
  104.  
  105. (define =_rat
  106.   (lambda (x y)
  107.     (and (s-= (&numer x) (&numer y)) (s-= (&denom x) (&denom y)))))
  108.  
  109. (define (make-rel_rat rel)
  110.   (lambda (x y)
  111.     (rel (s-* (&numer x) (&denom y)) (s-* (&denom x) (&numer y)))))
  112.  
  113. (define <_rat (make-rel_rat s-<))
  114. (define >_rat (make-rel_rat s->))
  115. (define <=_rat (make-rel_rat s-<=))
  116. (define >=_rat (make-rel_rat s->=))
  117.  
  118. (define (negate_rat x)
  119.   (make-rational (s-- (&numer x)) (&denom x)))
  120.  
  121. (define (recip_rat x) (make-rational (&denom x) (&numer x)))
  122.  
  123. (define negative_rat? (lambda (p) (s-negative? (&numer p))))
  124.  
  125. (define positive_rat? (lambda (p) (s-positive? (&numer p))))
  126.  
  127. (define (abs_rat p) (if (negative_rat? p) (negate_rat p) p))
  128.  
  129. (define-integrable 1_rat (construct-rational 1 1))
  130. (define-integrable -1_rat (construct-rational -1 1))
  131.  
  132. (define -1+_rat (lambda (p) (+_rat p -1_rat)))
  133. (define 1+_rat (lambda (p) (+_rat p 1_rat)))
  134.  
  135.  
  136. (define (numerator q)
  137.   (cond ((integer? q) q)
  138.     ((eq? '*rat* (car q)) (cadr q))
  139.     (else (error "bad rational"))))
  140.  
  141. (define (denominator q)
  142.   (cond ((integer? q) 1)
  143.     ((eq? '*rat* (car q)) (caddr q))
  144.     (else (error "bad rational"))))
  145.  
  146. (define (rational->float q)
  147.   (s-/ (&numer q) (&denom q)))
  148.  
  149. (define (integer->rational z)
  150.   (construct-rational z 1))
  151.  
  152. ;;; Real arithmetic with rationals
  153.  
  154. (define (make-unary-real-operation s-op rat-op)
  155.   (lambda (p)
  156.     (cond ((integer? p) (s-op p))
  157.           ((&rational? p) (rat-op p))
  158.           (else (s-op p)))))
  159.  
  160. (define (make-binary-real-operation s-op rat-op)
  161.   (lambda (p q)
  162.     (cond ((integer? p)
  163.            (cond ((integer? q) (s-op p q))
  164.                  ((&rational? q) (rat-op (integer->rational p) q))
  165.                  (else (s-op p q))))
  166.           ((&rational? p)
  167.            (cond ((integer? q) (rat-op p (integer->rational q)))
  168.                  ((&rational? q) (rat-op p q))
  169.                  (else (s-op (rational->float p) q))))
  170.           (else
  171.            (cond ((integer? q) (s-op p q))
  172.                  ((&rational? q) (s-op p (rational->float q)))
  173.                  (else (s-op p q)))))))
  174.  
  175. (define (make-binary-special-operation s-op rat-op)
  176.   (lambda (p q)
  177.     (cond ((integer? p)
  178.            (cond ((integer? q)
  179.                   (rat-op (integer->rational p)
  180.                           (integer->rational q)))
  181.                  ((&rational? q) (rat-op (integer->rational p) q))
  182.                  (else (s-op p q))))
  183.           ((&rational? p)
  184.            (cond ((integer? q) (rat-op p (integer->rational q)))
  185.                  ((&rational? q) (rat-op p q))
  186.                  (else (s-op (rational->float p) q))))
  187.           (else                                 ; p is a float
  188.            (cond ((integer? q) (s-op p q))
  189.                  ((&rational? q) (s-op p (rational->float q)))
  190.                  (else (s-op p q)))))))
  191.  
  192. (define (make-irrational-operation s-op)
  193.   (lambda ps
  194.     (apply s-op
  195.            (map (lambda (p)
  196.                   (cond ((integer? p) p)
  197.                         ((&rational? p) (rational->float p))
  198.                         (else p)))
  199.                 ps))))
  200.  
  201. (define-integrable real?
  202.   (lambda (x) (or (s-real? x) (&rational? x))))
  203.  
  204. (define-integrable r-zero?
  205.   (lambda (z) (and (s-number? z) (s-zero? z))))
  206.  
  207. (define-integrable r-= (make-binary-real-operation s-= =_rat))
  208.  
  209. (define-integrable r-1+ (make-unary-real-operation s-1+ 1+_rat))
  210. (define-integrable r--1+ (make-unary-real-operation s--1+ -1+_rat))
  211.  
  212. (define-integrable r-+ (make-binary-real-operation s-+ +_rat))
  213. (define-integrable r-- (make-binary-real-operation s-- -_rat))
  214. (define-integrable r-* (make-binary-real-operation s-* *_rat))
  215. (define-integrable r-/ (make-binary-special-operation s-/ /_rat))
  216.  
  217. (define-integrable r-abs (make-unary-real-operation s-abs abs_rat))
  218. (define-integrable positive?
  219.   (make-unary-real-operation s-positive? positive_rat?))
  220. (define-integrable negative?
  221.   (make-unary-real-operation s-negative? negative_rat?))
  222. (define-integrable < (make-binary-real-operation s-< <_rat))
  223. (define-integrable > (make-binary-real-operation s-> >_rat))
  224. (define-integrable <= (make-binary-real-operation s-<= <=_rat))
  225. (define-integrable >= (make-binary-real-operation s->= >=_rat))
  226.  
  227. (define-integrable r-sin (make-irrational-operation s-sin))
  228. (define-integrable r-cos (make-irrational-operation s-cos))
  229. (define-integrable r-atan (make-irrational-operation s-atan))
  230. (define-integrable r-acos (make-irrational-operation s-acos))
  231. (define-integrable r-asin (make-irrational-operation s-asin))
  232. (define-integrable r-sqrt (make-irrational-operation s-sqrt))
  233. (define-integrable r-exp (make-irrational-operation s-exp))
  234. (define-integrable r-log (make-irrational-operation s-log))
  235. (define-integrable r-expt (make-irrational-operation s-expt))
  236.  
  237. ;;; Operator stuff
  238.  
  239. (define-integrable operator-type "*operator*")
  240.  
  241. (define &procedure->operator
  242.   (lambda (procedure)
  243.     (cons operator-type procedure)))
  244.  
  245. (define &operator?
  246.   (lambda (object)
  247.     (and (pair? object)
  248.          (eq? operator-type (car object)))))
  249.  
  250. (define &operator->procedure cdr)
  251.  
  252. (define identity-operator
  253.   (&procedure->operator (lambda (object) object)))
  254.  
  255. (define (operation*operator->operator operation operator)
  256.   (&procedure->operator
  257.    (lambda (object)
  258.      (operation ((&operator->procedure operator) object)))))
  259.  
  260. (define (operation*constant*operator->operator operation constant operator)
  261.   (&procedure->operator
  262.    (lambda (object)
  263.      (operation constant ((&operator->procedure operator) object)))))
  264.  
  265. (define (operation*operator*constant->operator operation operator constant)
  266.   (&procedure->operator
  267.    (lambda (object)
  268.      (operation ((&operator->procedure operator) object) constant))))
  269.  
  270. (define (operation*operator*operator->operator operation operator1 operator2)
  271.   (&procedure->operator
  272.    (lambda (object)
  273.      (operation ((&operator->procedure operator1) object)
  274.         ((&operator->procedure operator2) object)))))
  275.  
  276. (define (procedure->operator procedure)
  277.   (if (procedure? procedure)
  278.       (&procedure->operator procedure)
  279.       (error "Illegal operator definition" procedure)))
  280.  
  281. (define (operator->procedure operator)
  282.   (if (&operator? operator)
  283.       (&operator->procedure operator)
  284.       (lambda (object) operator)))
  285.  
  286. (define (operator-compose operator1 operator2)
  287.   (&procedure->operator
  288.    (lambda (object)
  289.      ((operator->procedure operator1)
  290.       ((operator->procedure operator2)
  291.        object)))))
  292.  
  293.  
  294. (define (operator-value operator object)
  295.   (if (&operator? operator)
  296.       (with-reasonable-object object (&operator->procedure operator))
  297.       operator))
  298.  
  299. (define (with-reasonable-object object procedure)
  300.   (fluid-let ((operator-value-retry
  301.                 (call-with-current-continuation
  302.                   (lambda (k)
  303.                     (lambda ()
  304.                       (set! object (perturb-object object))
  305.                       (k k))))))
  306.     (procedure object)))
  307.  
  308. ;;; Because of the following definitions, these operators are
  309. ;;; restricted to numbers.
  310.  
  311. (define operator-perturbation-epsilon 1e-11)
  312.  
  313. (define (perturb-object object)
  314.   (+ object operator-perturbation-epsilon))
  315.  
  316.  
  317. (define (operation*operator*object->operator operation operator object)
  318.   (cond ((or (real? object) (&complex? object))
  319.      (operation*operator*constant->operator operation operator object))
  320.     ((&operator? object)
  321.      (operation*operator*operator->operator operation operator object))
  322.     (else
  323.      (error "Bad number" object))))
  324.  
  325. (define (operation*object*operator->operator operation object operator)
  326.   (cond ((or (real? object) (&complex? object))
  327.      (operation*constant*operator->operator operation object operator))
  328.     ((&operator? object)
  329.      (operation*operator*operator->operator operation object operator))
  330.     (else
  331.      (error "Bad number" object))))
  332.  
  333. ;;;; Complex number data abstraction
  334.  
  335. (define-integrable complex-type '*rect*)
  336.  
  337. (define-integrable &complex?
  338.   (lambda (z)
  339.     (and (pair? z) (eq? (car z) complex-type))))
  340.  
  341. (define-integrable &complex
  342.   (lambda (re im)
  343.     (list* complex-type re im)))
  344.  
  345. (define-integrable &real-part cadr)
  346.  
  347. (define-integrable &imag-part cddr)
  348.  
  349. (define (&conjugate z)
  350.   (&complex (&real-part z)
  351.         (r-- 0 (&imag-part z))))
  352.  
  353. (define (&sqr-magnitude z)
  354.   (define r-sqr (lambda (x) (r-* x x)))
  355.   (r-+ (r-sqr (&real-part z))
  356.        (r-sqr (&imag-part z))))
  357.  
  358. (define (&magnitude z)
  359.   (r-sqrt (&sqr-magnitude z)))
  360.  
  361. (define (&angle z m)
  362.   (if (zero? m)
  363.       (if (fluid-bound? operator-value-retry)
  364.           ((fluid operator-value-retry))
  365.           (error "angle: magnitude is zero"))
  366.       (r-atan (&imag-part z) (&real-part z))))
  367.  
  368.  
  369.  
  370. ;;; Complex stuff
  371.  
  372. (define (&make-rectangular re im)
  373.   (if (r-zero? im)
  374.       re
  375.       (&complex re im)))
  376.  
  377. (define (&make-polar r theta)
  378.   (if (r-zero? r)
  379.       0
  380.       (&make-rectangular (r-* r (r-cos theta))
  381.              (r-* r (r-sin theta)))))
  382.  
  383. (define i (&make-rectangular 0 1))
  384. (define -i (&make-rectangular 0 -1))
  385.  
  386. ;;;; Bottom layer
  387.  
  388. (define (make-binary-componentwise-operation r-opr accum error-string)
  389.   (define (operation z1 z2)
  390.     (cond ((real? z1)
  391.        (cond ((real? z2) (r-opr z1 z2))
  392.          ((&complex? z2)
  393.           (accum (r-opr z1 (&real-part z2))
  394.              (r-opr 0 (&imag-part z2))))
  395.          ((&operator? z2)
  396.           (operation*constant*operator->operator operation z1 z2))
  397.          (else (error error-string z2))))
  398.       ((&complex? z1)
  399.        (cond ((real? z2)
  400.           (accum (r-opr (&real-part z1) z2)
  401.              (r-opr (&imag-part z1) 0)))
  402.          ((&complex? z2)
  403.           (accum (r-opr (&real-part z1) (&real-part z2))
  404.              (r-opr (&imag-part z1) (&imag-part z2))))
  405.          ((&operator? z2)
  406.           (operation*constant*operator->operator operation z1 z2))
  407.          (else (error error-string z2))))
  408.       ((&operator? z1)
  409.        (operation*operator*object->operator operation z1 z2))
  410.       (else
  411.        (error error-string z1))))
  412.   operation)
  413.  
  414. (define &+
  415.   (make-binary-componentwise-operation r-+
  416.      &make-rectangular
  417.      "+: bad number"))
  418.  
  419. (define &-
  420.   (make-binary-componentwise-operation r--
  421.      &make-rectangular
  422.      "-: bad number"))
  423.  
  424. (define &=
  425.   (make-binary-componentwise-operation r-=
  426.      (lambda (p q) (and p q))
  427.      "=: bad number"))
  428.  
  429. (define (&* z1 z2)
  430.   (cond ((real? z1)
  431.      (cond ((real? z2) (r-* z1 z2))
  432.            ((&complex? z2)
  433.         (&make-rectangular (r-* z1 (&real-part z2))
  434.                    (r-* z1 (&imag-part z2))))
  435.            ((&operator? z2)
  436.         (operation*constant*operator->operator &* z1 z2))
  437.            (else (error "*: bad number" z2))))
  438.     ((&complex? z1)
  439.      (cond ((real? z2)
  440.         (&make-rectangular (r-* (&real-part z1) z2)
  441.                    (r-* (&imag-part z1) z2)))
  442.            ((&complex? z2)
  443.         (&make-rectangular
  444.          (r-- (r-* (&real-part z1) (&real-part z2))
  445.               (r-* (&imag-part z1) (&imag-part z2)))
  446.          (r-+ (r-* (&real-part z1) (&imag-part z2))
  447.               (r-* (&imag-part z1) (&real-part z2)))))
  448.            ((&operator? z2)
  449.         (operation*constant*operator->operator &* z1 z2))
  450.            (else (error "*: bad number" z2))))
  451.     ((&operator? z1)
  452.      (operation*operator*object->operator &* z1 z2))
  453.     (else
  454.      (error "*: bad number" z1))))
  455.  
  456. (define (square z) (&* z z))
  457.  
  458. (define (&/ z1 z2)
  459.   (cond ((real? z2)
  460.      (/real z1 z2))
  461.     ((&complex? z2)
  462.      (/real (&* z1 (&conjugate z2)) (&sqr-magnitude z2)))
  463.     ((&operator? z2)
  464.      (operation*object*operator->operator &/ z1 z2))
  465.     (else
  466.      (error "/: bad number" z2))))
  467.  
  468. (define (/real z x)
  469.   (cond ((r-zero? x)
  470.      (if (fluid-bound? operator-value-retry)
  471.          ((fluid operator-value-retry))
  472.          (error "/: division by zero" z x)))
  473.     ((real? z)
  474.      (r-/ z x))
  475.     ((&complex? z)
  476.      (&make-rectangular (r-/ (&real-part z) x)
  477.                             (r-/ (&imag-part z) x)))
  478.     ((&operator? z)
  479.      (operation*operator*constant->operator /real z x))
  480.     (else
  481.      (error "/: bad number" z))))
  482.  
  483. ;;;; User selector, constructors, and type predicates
  484.  
  485. (define-integrable complex?
  486.   (lambda (z)
  487.     (or (real? z) (&complex? z))))
  488.  
  489. (define (make-rectangular x y)
  490.   (if (and (real? x) (real? y))
  491.       (&make-rectangular x y)
  492.       (error "Make-rectangular: bad arguments" x y)))
  493.  
  494. (define (make-polar r theta)
  495.   (if (and (real? r) (real? theta))
  496.       (&make-polar r theta)
  497.       (error "Make-polar: bad arguments" r theta)))
  498.  
  499. (define (real-part z)
  500.   (cond ((real? z) z)
  501.     ((&complex? z) (&real-part z))
  502.     ((&operator? z) (operation*operator->operator real-part z))
  503.     (else (error "Real-part: bad number" z))))
  504.  
  505. (define (imag-part z)
  506.   (cond ((real? z) 0)
  507.     ((&complex? z) (&imag-part z))
  508.     ((&operator? z) (operation*operator->operator imag-part z))
  509.     (else (error "Imag-part: bad number" z))))
  510.  
  511. (define (magnitude z)
  512.   (cond ((real? z) (r-abs z))
  513.     ((&complex? z) (&magnitude z))
  514.     ((&operator? z) (operation*operator->operator magnitude z))
  515.     (else (error "Magnitude: bad number" z))))
  516.  
  517. (define-integrable abs magnitude)
  518.  
  519. (define (angle z)
  520.   (cond ((real? z) (if (negative? z) pi 0))
  521.     ((&complex? z) (&angle z (&magnitude z)))
  522.     ((&operator? z) (operation*operator->operator angle z))
  523.     (else (error "Angle: bad number" z))))
  524.  
  525. (define (conjugate z)
  526.   (cond ((real? z) z)
  527.     ((&complex? z) (&conjugate z))
  528.     ((&operator? z) (operation*operator->operator conjugate z))
  529.     (else (error "Conjugate: bad number" z))))
  530.  
  531. ;;;; Unary arithmetic and predicates
  532.  
  533. (define zero?
  534.   (lambda (z)
  535.     (cond ((real? z)
  536.        (r-zero? z))
  537.           ((&complex? z)
  538.        (and (r-zero? (&real-part z))
  539.             (r-zero? (&imag-part z))))
  540.           (else
  541.        (error "Zero?: Bad complex number" z)))))
  542.  
  543. (define (make-unary-componentwise-operation r-op error-string)
  544.   (define (operation z)
  545.     (cond ((real? z)
  546.        (r-op z))
  547.       ((&complex? z)
  548.        (&make-rectangular (r-op (&real-part z))
  549.                   (&imag-part z)))
  550.       ((&operator? z)
  551.        (operation*operator->operator operation z))
  552.       (else
  553.        (error error-string z))))
  554.   operation)
  555.  
  556. (define-integrable 1+
  557.   (make-unary-componentwise-operation r-1+ "1+: bad number"))
  558.  
  559. (define-integrable -1+
  560.   (make-unary-componentwise-operation r--1+ "-1+: bad number"))
  561.  
  562.  
  563. ;;;; N-ary predicates
  564.  
  565. ;;; Predicates are extensions of common Lisp because they can be
  566. ;;; given no arguments, in which case they return `#t'.
  567.  
  568. (define (make-pairwise-test &pred)
  569.   (lambda args
  570.     (define (loop x y rem)
  571.       (and (&pred x y)
  572.        (or (null? rem)
  573.            (loop y (car rem) (cdr rem)))))
  574.     (or (null? args)
  575.     (null? (cdr args))
  576.     (loop (car args) (cadr args) (cddr args)))))
  577.  
  578. (define-integrable = (make-pairwise-test &=))
  579.  
  580. ;;;; N-ary arithmetic
  581.  
  582. (define (accumulation operation identity)
  583.   (lambda rest
  584.     (define (loop accum rem)
  585.       (if (null? rem)
  586.       accum
  587.       (loop (operation accum (car rem)) (cdr rem))))
  588.     (cond ((null? rest) identity)
  589.       ((null? (cdr rest)) (car rest))
  590.       (else (operation (car rest) (loop (cadr rest) (cddr rest)))))))
  591.  
  592. (define-integrable + (accumulation &+ 0))
  593. (define-integrable * (accumulation &* 1))
  594.  
  595. (define (inverse-accumulation operation1 operation2 identity)
  596.   (lambda rest
  597.     (define (loop accum rem)
  598.       (if (null? rem)
  599.       accum
  600.       (loop (operation2 accum (car rem)) (cdr rem))))
  601.     (cond ((null? rest) identity)
  602.       ((null? (cdr rest)) (operation1 identity (car rest)))
  603.       ((null? (cddr rest)) (operation1 (car rest) (cadr rest)))
  604.       (else (operation1 (car rest) (loop (cadr rest) (cddr rest)))))))
  605.  
  606. (define-integrable - (inverse-accumulation &- &+ 0))
  607. (define-integrable / (inverse-accumulation &/ &* 1))
  608.  
  609. ;;;; Transcendental functions
  610.  
  611. (define (exp z)
  612.   (cond ((real? z) (r-exp z))
  613.     ((&complex? z)
  614.      (&make-polar (r-exp (&real-part z)) (&imag-part z)))
  615.     ((&operator? z)
  616.      (operation*operator->operator exp z))
  617.     (else (error "Exp: bad number" z))))
  618.  
  619. (define (log z)
  620.   (cond ((real? z)
  621.      (if (negative? z)
  622.          (&make-rectangular (error-log (r-- 0 z)) pi)
  623.          (error-log z)))
  624.     ((&complex? z)
  625.      (let ((m (&magnitude z)))
  626.        (&make-rectangular (error-log m) (&angle z m))))
  627.     ((&operator? z)
  628.      (operation*operator->operator log z))
  629.     (else
  630.      (error "log: bad number" z))))
  631.  
  632. (define (error-log x)
  633.   (if (r-zero? x)
  634.       (if (fluid-bound? operator-value-retry)
  635.           ((fluid operator-value-retry))
  636.           (error "log: Log of zero"))
  637.       (r-log x)))
  638.  
  639. (define (expt z1 z2)
  640.   (if (and (real? z1) (real? z2))
  641.       (r-expt z1 z2)
  642.       (exp (&* z2 (log z1)))))
  643.  
  644. (define log10
  645.   (let ((e^base (log 10)))
  646.     (lambda (z)
  647.       (/ (log z) e^base))))
  648.  
  649. (define (sqrt z)
  650.   (cond ((real? z)
  651.      (if (negative? z)
  652.          (&make-rectangular 0 (r-sqrt (r-- 0 z)))
  653.          (r-sqrt z)))
  654.     ((&complex? z)
  655.      (let ((m (&magnitude z)))
  656.        (&make-polar (r-sqrt m)
  657.             (r-/ (&angle z m) 2))))
  658.     ((&operator? z)
  659.      (operation*operator->operator sqrt z))
  660.     (else
  661.      (error "sqrt: bad number" z))))
  662.  
  663. (define (sin z)
  664.   (define 2i (&make-rectangular 0 2))
  665.   (cond ((real? z) (r-sin z))
  666.     ((&complex? z)
  667.      (let ((iz (&* i z)))
  668.        (&/ (&- (exp iz) (exp (&- 0 iz)))
  669.         2i)))
  670.     ((&operator? z)
  671.      (operation*operator->operator sin z))
  672.     (else
  673.      (error "sin: bad number" z))))
  674.  
  675. (define (cos z)
  676.   (cond ((real? z) (r-cos z))
  677.     ((&complex? z)
  678.      (let ((iz (&* i z)))
  679.        (&/ (&+ (exp iz) (exp (&- 0 iz)))
  680.         2)))
  681.     ((&operator? z)
  682.      (operation*operator->operator cos z))
  683.     (else
  684.      (error "cos: bad number" z))))
  685.  
  686. (define (tan z)
  687.   (cond ((real? z)
  688.      (let ((den (r-cos z)))
  689.        (if (r-zero? den)
  690.                (if (fluid-bound? operator-value-retry)
  691.                    ((fluid operator-value-retry))
  692.                (error "tan: Overflow" z))
  693.            (r-/ (r-sin z) den))))
  694.     ((&complex? z)
  695.      (let* ((iz (&* i z))
  696.         (den (&+ (exp iz) (exp (&- 0 iz)))))
  697.        (if (zero? den)
  698.            (if (fluid-bound? operator-value-retry)
  699.                    ((fluid operator-value-retry))
  700.                (error "tan: Overflow" z))
  701.             (&* -i (&/ (&- (exp iz) (exp (&- 0 iz))) den)))))
  702.     ((&operator? z)
  703.      (operation*operator->operator tan z))
  704.     (else
  705.      (error "tan: bad number" z))))
  706.  
  707. ;;; There's a problem trying to decide what meaning to attach to (atan y x)
  708. ;;;  when y and x can both be complex. The solution taken below is to
  709. ;;;  handle the one- and two-argument cases separately, as follows. With
  710. ;;;  one argument, the analytic function (atan z) is returned. With two
  711. ;;;  arguments -- (atan z1 z2) -- the returned value is (angle z), which is
  712. ;;;  to say, (r-atan (imaginary-part z) (real-part z)), r-atan = old atan, and
  713. ;;;  z is (+ (* i z1) z2); this has the advantage of doing the right thing
  714. ;;;  when z1 and z2 are real. Notice that this is not an analytic function;
  715. ;;;  e.g., it is purely real-valued. Perhaps the only virtue to not simply
  716. ;;;  flagging an error in the two-argument case when both arguments are
  717. ;;;  not real is that the user may have planned for them to be real and
  718. ;;;  missed by a small amount due to roundoff; in this case the answer will
  719. ;;;  be essentially what's wanted. -- MH.
  720.  
  721. (define atan 
  722.   (let ((2i (&* 2 i)))
  723.     (lambda (y . optionals)
  724.       (let ((x (if (not (null? optionals)) (car optionals) 1)))
  725.         (if (and (real? y) (real? x))
  726.             (r-atan y x)
  727.             (let ((iy (&* i y)))
  728.               (&/ (&- (log (&+ x iy)) (log (&- x iy))) 2i)))))))
  729.  
  730. (define (acos z)
  731.   (if (and (real? z) (<= z 1) (>= z -1))
  732.       (r-acos z)
  733.       (&* -i (log (&- z (&* -i (sqrt (&- 1 (&* z z)))))))))
  734.  
  735. (define (asin z)
  736.   (if (and (real? z) (<= z 1) (>= z -1))
  737.       (r-asin z)
  738.       (&* -i (log (&+ (&* i z) (sqrt (&- 1 (&* z z))))))))
  739.  
  740. (define-integrable number? complex?)
  741.  
  742. (define-integrable operator?
  743.   (lambda (object)
  744.     (or (number? object) (&operator? object))))
  745.  
  746.  
  747.  
  748. ;;; 6.003 useful operators
  749.  
  750. (define s identity-operator)
  751. (define t identity-operator)
  752.  
  753.  
  754. ;;; For electrical engineers
  755.  
  756. (define j  i)
  757. (define -j -i)
  758.