Diff of /src/code/irrat.lisp

revision 1.40 by rtoy, Wed Jun 9 14:48:15 2004 UTC revision 1.41 by rtoy, Tue Oct 19 15:07:30 2004 UTC
# Line 916  Z may be any number, but the result is a Line 916  Z may be any number, but the result is a
916  (defun complex-atanh (z)  (defun complex-atanh (z)
917    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
918    (declare (number z))    (declare (number z))
919    (let* (;; Constants    (if (and (realp z) (< z -1))
920           (theta (/ (sqrt most-positive-double-float) 4.0d0))        ;; atanh is continuous in quadrant III in this case.
921           (rho (/ 4.0d0 (sqrt most-positive-double-float)))        (complex-atanh (complex z -0f0))
922           (half-pi (/ pi 2.0d0))        (let* ( ;; Constants
923           (rp (float (realpart z) 1.0d0))               (theta (/ (sqrt most-positive-double-float) 4.0d0))
924           (beta (float-sign rp 1.0d0))               (rho (/ 4.0d0 (sqrt most-positive-double-float)))
925           (x (* beta rp))               (half-pi (/ pi 2.0d0))
926           (y (* beta (- (float (imagpart z) 1.0d0))))               (rp (float (realpart z) 1.0d0))
927           (eta 0.0d0)               (beta (float-sign rp 1.0d0))
928           (nu 0.0d0))               (x (* beta rp))
929      ;; Shouldn't need this declare.               (y (* beta (- (float (imagpart z) 1.0d0))))
930      (declare (double-float x y))               (eta 0.0d0)
931      (locally               (nu 0.0d0))
932          (declare (optimize (speed 3)))          ;; Shouldn't need this declare.
933        (cond ((or (> x theta)          (declare (double-float x y))
934                   (> (abs y) theta))          (locally
935               ;; To avoid overflow...              (declare (optimize (speed 3)))
936               (setf nu (float-sign y half-pi))            (cond ((or (> x theta)
937               ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),                       (> (abs y) theta))
938               ;; which can cause overflow.  Arrange this computation so                   ;; To avoid overflow...
939               ;; that it won't overflow.                   (setf nu (float-sign y half-pi))
940               (setf eta (let* ((x-bigger (> x (abs y)))                   ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
941                                (r (if x-bigger (/ y x) (/ x y)))                   ;; which can cause overflow.  Arrange this computation so
942                                (d (+ 1.0d0 (* r r))))                   ;; that it won't overflow.
943                           (if x-bigger                   (setf eta (let* ((x-bigger (> x (abs y)))
944                               (/ (/ x) d)                                    (r (if x-bigger (/ y x) (/ x y)))
945                               (/ (/ r y) d)))))                                    (d (+ 1.0d0 (* r r))))
946              ((= x 1.0d0)                               (if x-bigger
947               ;; Should this be changed so that if y is zero, eta is set                                   (/ (/ x) d)
948               ;; to +infinity instead of approx 176?  In any case                                   (/ (/ r y) d)))))
949               ;; tanh(176) is 1.0d0 within working precision.                  ((= x 1.0d0)
950               (let ((t1 (+ 4d0 (square y)))                   ;; Should this be changed so that if y is zero, eta is set
951                     (t2 (+ (abs y) rho)))                   ;; to +infinity instead of approx 176?  In any case
952                 (setf eta (log (/ (sqrt (sqrt t1))                   ;; tanh(176) is 1.0d0 within working precision.
953                                   (sqrt t2))))                   (let ((t1 (+ 4d0 (square y)))
954                 (setf nu (* 0.5d0                         (t2 (+ (abs y) rho)))
955                             (float-sign y                     (setf eta (log (/ (sqrt (sqrt t1))
956                                         (+ half-pi (atan (* 0.5d0 t2))))))))                                       (sqrt t2))))
957              (t                     (setf nu (* 0.5d0
958               (let ((t1 (+ (abs y) rho)))                                 (float-sign y
959                 ;; Normal case using log1p(x) = log(1 + x)                                             (+ half-pi (atan (* 0.5d0 t2))))))))
960                 (setf eta (* 0.25d0                  (t
961                              (%log1p (/ (* 4.0d0 x)                   (let ((t1 (+ (abs y) rho)))
962                                         (+ (square (- 1.0d0 x))                     ;; Normal case using log1p(x) = log(1 + x)
963                                            (square t1))))))                     (setf eta (* 0.25d0
964                 (setf nu (* 0.5d0                                  (%log1p (/ (* 4.0d0 x)
965                             (atan (* 2.0d0 y)                                             (+ (square (- 1.0d0 x))
966                                   (- (* (- 1.0d0 x)                                                (square t1))))))
967                                         (+ 1.0d0 x))                     (setf nu (* 0.5d0
968                                      (square t1))))))))                                 (atan (* 2.0d0 y)
969        (coerce-to-complex-type (* beta eta)                                       (- (* (- 1.0d0 x)
970                                (- (* beta nu))                                             (+ 1.0d0 x))
971                                z))))                                          (square t1))))))))
972              (coerce-to-complex-type (* beta eta)
973                                      (- (* beta nu))
974                                      z)))))
975
976  (defun complex-tanh (z)  (defun complex-tanh (z)
977    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"
# Line 1056  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-1+z (complex-sqrt (1+z z)))    (if (and (realp z) (> z 1))
1063          (sqrt-1-z (complex-sqrt (1-z z))))        ;; acos is continuous in quadrant IV in this case.
1064      (with-float-traps-masked (:divide-by-zero)        (complex-acos (complex z -0f0))
1065        (complex (* 2 (atan (/ (realpart sqrt-1-z)        (let ((sqrt-1+z (complex-sqrt (1+z z)))
1066                               (realpart sqrt-1+z))))              (sqrt-1-z (complex-sqrt (1-z z))))
1067                 (asinh (imagpart (* (conjugate sqrt-1+z)          (with-float-traps-masked (:divide-by-zero)
1068                                     sqrt-1-z)))))))            (complex (* 2 (atan (/ (realpart sqrt-1-z)
1069                                     (realpart sqrt-1+z))))
1070                       (asinh (imagpart (* (conjugate sqrt-1+z)
1071                                           sqrt-1-z))))))))
1072
1073  (defun complex-acosh (z)  (defun complex-acosh (z)
1074    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
# Line 1083  Z may be any number, but the result is a Line 1089  Z may be any number, but the result is a
1089
1090  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1091    (declare (number z))    (declare (number z))
1092    (let ((sqrt-1-z (complex-sqrt (1-z z)))    (if (and (realp z) (> z 1))
1093          (sqrt-1+z (complex-sqrt (1+z z))))        ;; asin is continuous in quadrant IV in this case.
1094      (with-float-traps-masked (:divide-by-zero)        (complex-asin (complex z -0f0))
1095        (complex (atan (/ (realpart z)        (let ((sqrt-1-z (complex-sqrt (1-z z)))
1096                          (realpart (* sqrt-1-z sqrt-1+z))))              (sqrt-1+z (complex-sqrt (1+z z))))
1097                 (asinh (imagpart (* (conjugate sqrt-1-z)          (with-float-traps-masked (:divide-by-zero)
1098                                     sqrt-1+z)))))))            (complex (atan (/ (realpart z)
1099                                (realpart (* sqrt-1-z sqrt-1+z))))
1100                       (asinh (imagpart (* (conjugate sqrt-1-z)
1101                                           sqrt-1+z))))))))
1102
1103  (defun complex-asinh (z)  (defun complex-asinh (z)
1104    "Compute asinh z = log(z + sqrt(1 + z*z))    "Compute asinh z = log(z + sqrt(1 + z*z))

Legend:
 Removed from v.1.40 changed lines Added in v.1.41

 ViewVC Help Powered by ViewVC 1.1.5