... |
... |
@@ -1334,6 +1334,27 @@ and Y are coerced to single-float." |
1334
|
1334
|
(t
|
1335
|
1335
|
(values rho 0)))))))
|
1336
|
1336
|
|
|
1337
|
+;; sqrt(z) = z^(1/2)
|
|
1338
|
+;; = exp(1/2*log(z))
|
|
1339
|
+;;
|
|
1340
|
+;; Some special values, with x > 0
|
|
1341
|
+;; sqrt(x + i*0) = exp(1/2*log(x+i*0))
|
|
1342
|
+;; = exp(1/2*(log(x)+i*0))
|
|
1343
|
+;; = exp(1/2*log(x)+i*0)
|
|
1344
|
+;; = exp(1/2*log(x))*exp(i*0)
|
|
1345
|
+;; = sqrt(x)*(cos 0 + i*sin(0))
|
|
1346
|
+;; = sqrt(x) + i*0
|
|
1347
|
+;; sqrt(x - i*0) = exp(1/2*log(x-i*0))
|
|
1348
|
+;; = exp(1/2*(log(x) - i*0))
|
|
1349
|
+;; = exp(1/2*log(x) - i*0)
|
|
1350
|
+;; = exp(1/2*log(x))*exp(-i*0)
|
|
1351
|
+;; = sqrt(x)*(cos(-0) + i*sin(-0))
|
|
1352
|
+;; = sqrt(x) - i*0
|
|
1353
|
+;; Kahan also gives the following for b > 0
|
|
1354
|
+;; sqrt(-b+/-i*0) = +0 +/- i*sqrt(b)
|
|
1355
|
+;; sqrt(x +/- i*inf) = +inf +/- i*inf for all finite x.
|
|
1356
|
+;; sqrt(inf +/- i*b) = +inf +/- i*0
|
|
1357
|
+;; sqrt(-inf +/- i*b) = +0 +/- i*inf
|
1337
|
1358
|
(defun complex-sqrt (z)
|
1338
|
1359
|
"Principle square root of Z
|
1339
|
1360
|
|
... |
... |
@@ -1604,6 +1625,51 @@ Z may be any number, but the result is always a complex." |
1604
|
1625
|
(+ z 1)
|
1605
|
1626
|
(complex (+ (realpart z) 1) (imagpart z))))
|
1606
|
1627
|
|
|
1628
|
+;; acos(z) = 2*log(sqrt((1+z)/2) + i*sqrt((1-z)/2))/i
|
|
1629
|
+;; = pi/2 - asin(z)
|
|
1630
|
+;;
|
|
1631
|
+;; In particular for z = x > 1:
|
|
1632
|
+;; acos(x) = 2/i*log(sqrt((x+1)/2) + i*sqrt((1-x)/2))
|
|
1633
|
+;; = 2/i*log(sqrt((x+1)/2) + i*i*sqrt((x-1)/2))
|
|
1634
|
+;; = 2/i*log(sqrt((x+1)/2) - sqrt((x-1)/2))
|
|
1635
|
+;; = -2*i*log(sqrt((x+1)/2) - sqrt((x-1)/2))
|
|
1636
|
+;; = -i*log(x-sqrt(x-1)*sqrt(x+1))
|
|
1637
|
+;; Thus, acos(2) = -i*log(2-sqrt(3) = -i*1.31695...
|
|
1638
|
+;;
|
|
1639
|
+;; Similarly for z = -x, x > 1:
|
|
1640
|
+;; acos(-x) = 2/i*log(sqrt((1-x)/2) + i*sart((1+x)/2))
|
|
1641
|
+;; = 2/i*log(i*sqrt((x-1)/2) + i*sqrt((1+x)/2))
|
|
1642
|
+;; = 2/i*log(i*(sqrt((x+1)/2)+sqrt((x-1)/2)))
|
|
1643
|
+;; = 2/i*(log(sqrt((x+1)/2)+sqrt((x-1)/2)) + i*pi/2)
|
|
1644
|
+;; = -i*2*log(sqrt((x+1)/2)+sqrt((x-1)/2)) - pi
|
|
1645
|
+;; = -pi - i*log(x+sqrt(x+1)*(sqrt(x-1)))
|
|
1646
|
+;;
|
|
1647
|
+;; Thus acos(-2) = -pi - i*log(2+sqrt(3)) = -pi -i*1.31695...
|
|
1648
|
+;;
|
|
1649
|
+;; Now see what aacos(x-i0) for x > 1 is:
|
|
1650
|
+;; acos(x-i*0) = 2/i*log(sqrt((x+1)/2-i*0) + i*sqrt((1-x)/2+i*0))
|
|
1651
|
+;; = 2/i*log(sqrt((x+1)/2) -i*0+ i*(0+i*sqrt((x-1)/2)))
|
|
1652
|
+;; = 2/i*log(sqrt((x+1)/2) - sqrt((x-1)/2) - i*0)
|
|
1653
|
+;; = 2/i*[log(sqrt((x+1)/2) - sqrt((x-1)/2)) - i*0]
|
|
1654
|
+;; = -i*[2*log(sqrt((x+1)/2) - sqrt((x-1)/2)) - i*0]
|
|
1655
|
+;; = -i*log(x-sqrt(x-1)*sqrt(x+1)) + 0
|
|
1656
|
+;;
|
|
1657
|
+;; Thus acos(2 - i*0) is the same as acos(2). That is, acos is
|
|
1658
|
+;; continuous with quadrant IV for x > 1.
|
|
1659
|
+;;
|
|
1660
|
+;; For acos(-x+i0), x > 1:
|
|
1661
|
+;; acos(-x+i0) = 2/i*log(sqrt((1-x+i0)/2) + i*sqrt((1+x-i0)/2))
|
|
1662
|
+;; = 2/i*log(sqrt((1-x)/2+i0) + i*sqrt((1+x)/2-i0))
|
|
1663
|
+;; = 2/i*log((i*sqrt((x-1)/2) + 0) + i*(sqrt((1+x)/2)-i0))
|
|
1664
|
+;; = 2/i*log(i*sqrt((x-1)/2) + i*sqrt((1+x)/2) + 0)
|
|
1665
|
+;; = 2/i*log(i*(sqrt((x-1)/2 + sqrt(1+x)/2)) + 0)
|
|
1666
|
+;; = 2/i*(log(sqrt(x+sqrt(x-1)*sqrt(x+1))) + i*pi/2)
|
|
1667
|
+;; = -i*log(x+sqrt(x-1)*sqrt(x+1)) + pi
|
|
1668
|
+;; = pi - i*log(x+sqrt(x-1)*sqrt(x+1))
|
|
1669
|
+;;
|
|
1670
|
+;; Thus acos(-2+i0) = pi - i*log(sqrt(3)+2) = pi -i*1.31695, and this
|
|
1671
|
+;; equal acos(-2). This means acos is continuous with quadrant II.
|
|
1672
|
+;;
|
1607
|
1673
|
(defun complex-acos (z)
|
1608
|
1674
|
"Compute acos z = pi/2 - asin z
|
1609
|
1675
|
|