| ... |
... |
@@ -1510,7 +1510,7 @@ Z may be any number, but the result is always a complex." |
|
1510
|
1510
|
;; NOTE: this differs from what the CLHS says for the continuity.
|
|
1511
|
1511
|
;; Instead of the text in the CLHS, we choose to use the definition
|
|
1512
|
1512
|
;; to derive the correct values.
|
|
1513
|
|
- (if (and nil (realp z))
|
|
|
1513
|
+ (if (realp z)
|
|
1514
|
1514
|
(complex-atanh (complex (float z) (- (* 0 (float z)))))
|
|
1515
|
1515
|
(let* ( ;; Constants
|
|
1516
|
1516
|
(theta (/ (sqrt most-positive-double-float) 4.0d0))
|
| ... |
... |
@@ -1518,9 +1518,10 @@ Z may be any number, but the result is always a complex." |
|
1518
|
1518
|
(half-pi (/ pi 2.0d0))
|
|
1519
|
1519
|
(rp (float (realpart z) 1.0d0))
|
|
1520
|
1520
|
(beta (float-sign rp 1.0d0))
|
|
1521
|
|
- (z* (conjugate z))
|
|
1522
|
|
- (x (* beta (realpart z*)))
|
|
1523
|
|
- (y (* beta (imagpart z*)))
|
|
|
1521
|
+ ;; x+iy = beta*conjugate(z), but being careful to produce
|
|
|
1522
|
+ ;; a signed-zero if z is rational.
|
|
|
1523
|
+ (x (* beta (realpart z)))
|
|
|
1524
|
+ (y (* beta (- (float (imagpart z) 1d0))))
|
|
1524
|
1525
|
(eta 0.0d0)
|
|
1525
|
1526
|
(nu 0.0d0))
|
|
1526
|
1527
|
;; Shouldn't need this declare.
|
| ... |
... |
@@ -1562,8 +1563,7 @@ Z may be any number, but the result is always a complex." |
|
1562
|
1563
|
(atan (* 2.0d0 y)
|
|
1563
|
1564
|
(- (* (- 1.0d0 x)
|
|
1564
|
1565
|
(+ 1.0d0 x))
|
|
1565
|
|
- (square t1)))))
|
|
1566
|
|
- (format t "eta = ~A nu ~A~%" eta nu))))
|
|
|
1566
|
+ (square t1))))))))
|
|
1567
|
1567
|
(coerce-to-complex-type (* beta eta)
|
|
1568
|
1568
|
(- (* beta nu))
|
|
1569
|
1569
|
z)))))
|
| ... |
... |
@@ -1715,16 +1715,13 @@ Z may be any number, but the result is always a complex." |
|
1715
|
1715
|
#+double-double
|
|
1716
|
1716
|
(when (typep z '(or double-double-float (complex double-double-float)))
|
|
1717
|
1717
|
(return-from complex-acos (dd-complex-acos z)))
|
|
1718
|
|
- (if (and nil (realp z) (> z 1))
|
|
1719
|
|
- ;; acos is continuous in quadrant IV in this case.
|
|
1720
|
|
- (complex-acos (complex z -0f0))
|
|
1721
|
|
- (let ((sqrt-1+z (complex-sqrt (1+z z)))
|
|
1722
|
|
- (sqrt-1-z (complex-sqrt (1-z z))))
|
|
1723
|
|
- (with-float-traps-masked (:divide-by-zero)
|
|
1724
|
|
- (complex (* 2 (atan (/ (realpart sqrt-1-z)
|
|
1725
|
|
- (realpart sqrt-1+z))))
|
|
1726
|
|
- (asinh (imagpart (* (conjugate sqrt-1+z)
|
|
1727
|
|
- sqrt-1-z))))))))
|
|
|
1718
|
+ (let ((sqrt-1+z (complex-sqrt (1+z z)))
|
|
|
1719
|
+ (sqrt-1-z (complex-sqrt (1-z z))))
|
|
|
1720
|
+ (with-float-traps-masked (:divide-by-zero)
|
|
|
1721
|
+ (complex (* 2 (atan (/ (realpart sqrt-1-z)
|
|
|
1722
|
+ (realpart sqrt-1+z))))
|
|
|
1723
|
+ (asinh (imagpart (* (conjugate sqrt-1+z)
|
|
|
1724
|
+ sqrt-1-z)))))))
|
|
1728
|
1725
|
|
|
1729
|
1726
|
;; acosh(z) = 2*log(sqrt((x+1)/2) + sqrt((x-1)/2))
|
|
1730
|
1727
|
;;
|
| ... |
... |
@@ -1823,16 +1820,13 @@ Z may be any number, but the result is always a complex." |
|
1823
|
1820
|
#+double-double
|
|
1824
|
1821
|
(when (typep z '(or double-double-float (complex double-double-float)))
|
|
1825
|
1822
|
(return-from complex-asin (dd-complex-asin z)))
|
|
1826
|
|
- (if (and nil (realp z) (> z 1))
|
|
1827
|
|
- ;; asin is continuous in quadrant IV in this case.
|
|
1828
|
|
- (complex-asin (complex z -0f0))
|
|
1829
|
|
- (let ((sqrt-1-z (complex-sqrt (1-z z)))
|
|
1830
|
|
- (sqrt-1+z (complex-sqrt (1+z z))))
|
|
1831
|
|
- (with-float-traps-masked (:divide-by-zero)
|
|
1832
|
|
- (complex (atan (/ (realpart z)
|
|
1833
|
|
- (realpart (* sqrt-1-z sqrt-1+z))))
|
|
1834
|
|
- (asinh (imagpart (* (conjugate sqrt-1-z)
|
|
1835
|
|
- sqrt-1+z))))))))
|
|
|
1823
|
+ (let ((sqrt-1-z (complex-sqrt (1-z z)))
|
|
|
1824
|
+ (sqrt-1+z (complex-sqrt (1+z z))))
|
|
|
1825
|
+ (with-float-traps-masked (:divide-by-zero)
|
|
|
1826
|
+ (complex (atan (/ (realpart z)
|
|
|
1827
|
+ (realpart (* sqrt-1-z sqrt-1+z))))
|
|
|
1828
|
+ (asinh (imagpart (* (conjugate sqrt-1-z)
|
|
|
1829
|
+ sqrt-1+z)))))))
|
|
1836
|
1830
|
|
|
1837
|
1831
|
(defun complex-asinh (z)
|
|
1838
|
1832
|
"Compute asinh z = log(z + sqrt(1 + z*z))
|