# Diff of /src/code/irrat.lisp

revision 1.37.14.1 by rtoy, Fri May 14 16:26:07 2004 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))))))))

Legend:
 Removed from v.1.37.14.1 changed lines Added in v.1.37.14.2