# Diff of /src/code/irrat.lisp

revision 1.37 by toy, Wed Jan 29 18:51:48 2003 UTC revision 1.37.14.1 by rtoy, Fri May 14 16:26:07 2004 UTC
# 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))))