Diff of /src/code/irrat-dd.lisp

revision 1.7 by rtoy, Wed Mar 21 18:08:25 2007 UTC revision 1.8 by rtoy, Wed May 23 13:16:33 2007 UTC
# Line 1872  Z may be any number, but the result is a Line 1872  Z may be any number, but the result is a
1872
1873  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1874    (declare (number z))    (declare (number z))
1875    (if (and (realp z) (> z 1))    (cond ((realp z)
1876        ;; asin is continuous in quadrant IV in this case.           ;; Z is real, and |Z| > 1.
1877        (dd-complex-asin (complex z -0f0))           ;;
1878        (let* ((sqrt-1-z (complex-sqrt (1-z z)))           ;; Note that
1879               (sqrt-1+z (complex-sqrt (1+z z)))           ;;
1880               (den (realpart (* sqrt-1-z sqrt-1+z))))           ;;   sin(pi/2+i*y) = cosh(y)
1881          (cond ((zerop den)           ;;
1882                 ;; Like below but we handle atan part ourselves.           ;; and
1883                 (complex (if (minusp (float-sign den))           ;;
1884                              (- dd-pi/2)           ;;   sin(-pi/2+i*y) = -cosh(y).
1885                              dd-pi/2)           ;;
1886                     (asinh (imagpart (* (conjugate sqrt-1-z)           ;; Since cosh(y) >= 1 for all real y, we see that
1887                                         sqrt-1+z)))))           ;;
1888                (t           ;;   asin(z) = pi/2*sign(z) + i*acosh(abs(z))
1889                 (with-float-traps-masked (:divide-by-zero)           ;;
1890                   ;; We get a invalid operation here when z is real and |z| > 1.           (complex (if (minusp z)
1891                   (complex (atan (/ (realpart z)                        (- dd-pi/2)
1892                                     (realpart (* sqrt-1-z sqrt-1+z))))                        dd-pi/2)
1893                            (asinh (imagpart (* (conjugate sqrt-1-z)                    (dd-%acosh (abs z))))
1894                                                sqrt-1+z))))))))))          (t
1895             (let* ((sqrt-1-z (complex-sqrt (1-z z)))
1896                    (sqrt-1+z (complex-sqrt (1+z z)))
1897                    (den (realpart (* sqrt-1-z sqrt-1+z))))
1898               (cond ((zerop den)
1899                      ;; Like below but we handle atan part ourselves.
1900                      (complex (if (minusp (float-sign den))
1901                                   (- dd-pi/2)
1902                                   dd-pi/2)
1903                               (asinh (imagpart (* (conjugate sqrt-1-z)
1904                                                   sqrt-1+z)))))
1905                     (t