Raymond Toy pushed to branch rtoy-issue-78-unneded-code-code-in-complex-acos at cmucl / cmucl

Commits:

1 changed file:

Changes:

  • src/code/irrat.lisp
    ... ... @@ -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