# Diff of /src/code/irrat.lisp

revision 1.37 by toy, Wed Jan 29 18:51:48 2003 UTC revision 1.37.14.2 by rtoy, Thu Jun 3 13:37:02 2004 UTC
# Line 923  Z may be any number, but the result is a Line 923  Z may be any number, but the result is a
923        (cond ((or (> x theta)        (cond ((or (> x theta)
924                   (> (abs y) theta))                   (> (abs y) theta))
925               ;; To avoid overflow...               ;; To avoid overflow...
926               (setf eta (float-sign y half-pi))               (setf nu (float-sign y half-pi))
927               ;; nu is real part of 1/(x + iy).  This is x/(x^2+y^2),               ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
928               ;; which can cause overflow.  Arrange this computation so               ;; which can cause overflow.  Arrange this computation so
929               ;; that it won't overflow.               ;; that it won't overflow.
930               (setf nu (let* ((x-bigger (> x (abs y)))               (setf eta (let* ((x-bigger (> x (abs y)))
931                               (r (if x-bigger (/ y x) (/ x y)))                                (r (if x-bigger (/ y x) (/ x y)))
932                               (d (+ 1.0d0 (* r r))))                                (d (+ 1.0d0 (* r r))))
933                          (if x-bigger                           (if x-bigger
934                              (/ (/ x) d)                               (/ (/ x) d)
935                              (/ (/ r y) d)))))                               (/ (/ r y) d)))))
936              ((= x 1.0d0)              ((= x 1.0d0)
937               ;; Should this be changed so that if y is zero, eta is set               ;; Should this be changed so that if y is zero, eta is set
938               ;; to +infinity instead of approx 176?  In any case               ;; to +infinity instead of approx 176?  In any case
939               ;; tanh(176) is 1.0d0 within working precision.               ;; tanh(176) is 1.0d0 within working precision.
940               (let ((t1 (+ 4d0 (square y)))               (let ((t1 (+ 4d0 (square y)))
941                     (t2 (+ (abs y) rho)))                     (t2 (+ (abs y) rho)))
942                 (setf eta (log (/ (sqrt (sqrt t1)))                 (setf eta (log (/ (sqrt (sqrt t1))
943                                (sqrt t2)))                                   (sqrt t2))))
944                 (setf nu (* 0.5d0                 (setf nu (* 0.5d0
945                             (float-sign y                             (float-sign y
946                                         (+ half-pi (atan (* 0.5d0 t2))))))))                                         (+ half-pi (atan (* 0.5d0 t2))))))))
# Line 1012  Z may be any number, but the result is a Line 1012  Z may be any number, but the result is a
1012  ;; and these two expressions are equal if and only if arg conj z =  ;; and these two expressions are equal if and only if arg conj z =
1013  ;; -arg z, which is clearly true for all z.  ;; -arg z, which is clearly true for all z.
1014
1015    ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1016    ;; complex, the real is converted to a complex before performing the
1017    ;; operation.  However, Kahan says in this paper (pg 176):
1018    ;;
1019    ;; (iii) Careless handling can turn infinity or the sign of zero into
1020    ;;       misinformation that subsequently disappears leaving behind
1021    ;;       only a plausible but incorrect result.  That is why compilers
1022    ;;       must not transform z-1 into z-(1+i*0), as we have seen above,
1023    ;;       nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1024    ;;       subsequent logarithm or square root produce a non-zero
1025    ;;       imaginary part whose sign is opposite to what was intended.
1026    ;;
1027    ;; The interesting examples are too long and complicated to reproduce
1028    ;; here.  We refer the reader to his paper.
1029    ;;
1030    ;; The functions below are intended to handle the cases where a real
1031    ;; is mixed with a complex and we don't want CL complex contagion to
1032    ;; occur..
1033
1034    (declaim (inline 1+z 1-z z-1 z+1))
1035    (defun 1+z (z)
1036      (complex (+ 1 (realpart z)) (imagpart z)))
1037    (defun 1-z (z)
1038      (complex (- 1 (realpart z)) (- (imagpart z))))
1039    (defun z-1 (z)
1040      (complex (- (realpart z) 1) (imagpart z)))
1041    (defun z+1 (z)
1042      (complex (+ (realpart z) 1) (imagpart z)))
1043
1044  (defun complex-acos (z)  (defun complex-acos (z)
1045    "Compute acos z = pi/2 - asin z    "Compute acos z = pi/2 - asin z
1046
1047  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1048    (declare (number z))    (declare (number z))
1049    (let ((sqrt-1+z (complex-sqrt (+ 1 z)))    (let ((sqrt-1+z (complex-sqrt (1+z z)))
1050          (sqrt-1-z (complex-sqrt (- 1 z))))          (sqrt-1-z (complex-sqrt (1-z z))))
1052        (complex (* 2 (atan (/ (realpart sqrt-1-z)        (complex (* 2 (atan (/ (realpart sqrt-1-z)
1053                               (realpart sqrt-1+z))))                               (realpart sqrt-1+z))))
# Line 1030  Z may be any number, but the result is a Line 1059  Z may be any number, but the result is a
1059
1060  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1061    (declare (number z))    (declare (number z))
1062    (let ((sqrt-z-1 (complex-sqrt (- z 1)))    (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1063          (sqrt-z+1 (complex-sqrt (+ z 1))))          (sqrt-z+1 (complex-sqrt (z+1 z))))
1065        (complex (asinh (realpart (* (conjugate sqrt-z-1)        (complex (asinh (realpart (* (conjugate sqrt-z-1)
1066                                     sqrt-z+1)))                                     sqrt-z+1)))
# Line 1044  Z may be any number, but the result is a Line 1073  Z may be any number, but the result is a
1073
1074  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1075    (declare (number z))    (declare (number z))
1076    (let ((sqrt-1-z (complex-sqrt (- 1 z)))    (let ((sqrt-1-z (complex-sqrt (1-z z)))
1077          (sqrt-1+z (complex-sqrt (+ 1 z))))          (sqrt-1+z (complex-sqrt (1+z z))))