1261
1272
(return (cond ((null r) (list '(mabs) (car l)))
1262
1273
(t (bcons (fpabs (cdr r))))))))
1276
;;;; Bigfloat implementations of special functions.
1279
;;; This is still a bit messy. Some functions here take bigfloat
1280
;;; numbers, represented by ((bigfloat) <mant> <exp>), but others want
1281
;;; just the FP number, represented by (<mant> <exp>). Likewise, some
1282
;;; return a bigfloat, some return just the FP.
1284
;;; This needs to be systemized somehow. It isn't helped by the fact
1285
;;; that some of the routines above also do the samething.
1287
;;; The implementation for the special functions for a complex
1288
;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex
1289
;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in
1290
;;; Iserles and Powell (eds.) "The State of the Art in Numerical
1291
;;; Analysis", pp 165-211, Clarendon Press, 1987
1293
;; Compute exp(x) - 1, but do it carefully to preserve precision when
1294
;; |x| is small. X is a FP number, and a FP number is returned. That
1295
;; is, no bigfloat stuff.
1297
;; What is the right breakpoint here? Is 1 ok? Perhaps 1/e is better?
1298
(cond ((fpgreaterp (fpabs x) (fpone))
1300
(fpdifference (fpexp x) (fpone)))
1302
;; Use Taylor series for exp(x) - 1
1308
(setf term (fpquotient (fptimes* x term) (intofp n)))
1310
(setf ans (fpplus ans term)))
1313
;; log(1+x) for small x. X is FP number, and a FP number is returned.
1315
;; Use the same series as given above for fplog. For small x we use
1316
;; the series, otherwise fplog is accurate enough.
1317
(cond ((fpgreaterp (fpabs x) (fpone))
1318
(fplog (fpplus x (fpone))))
1320
(let* ((sum (intofp 0))
1321
(term (fpquotient x (fpplus x (intofp 2))))
1322
(f (fptimes* term term))
1325
((equal sum oldans))
1327
(setq sum (fpplus sum
1328
(fpquotient term (intofp n))))
1329
(setq term (fptimes* term f)))
1330
(fptimes* sum (intofp 2))))))
1332
;; sinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
1334
;; X must be a maxima bigfloat
1336
;; See, for example, Hart et al., Computer Approximations, 6.2.27:
1338
;; sinh(x) = 1/2*(D(x) + D(x)/(1+D(x)))
1340
;; where D(x) = exp(x) - 1.
1341
(let ((d (fpexpm1 (cdr (bigfloatp x)))))
1343
(fpquotient (fpplus d
1345
(fpplus d (fpone))))
1348
(defun big-float-sinh (x &optional y)
1349
;; The rectform for sinh for complex args should be numerically
1350
;; accurate, so return nil in that case.
1354
;; asinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
1356
;; asinh(x) = sign(x) * log(|x| + sqrt(1+x*x))
1360
;; asinh(x) = x, if 1+x*x = 1
1361
;; = sign(x) * (log(2) + log(x)), large |x|
1362
;; = sign(x) * log(2*|x| + 1/(|x|+sqrt(1+x*x))), if |x| > 2
1363
;; = sign(x) * log1p(|x|+x^2/(1+sqrt(1+x*x))), otherwise.
1365
;; But I'm lazy right now and we only implement the last 2 cases.
1366
;; We should implement all cases.
1367
(let* ((fp-x (cdr (bigfloatp x)))
1371
(minus (minusp (car fp-x)))
1373
;; We only use two formulas here. |x| <= 2 and |x| > 2. Should
1374
;; we add one for very big x and one for very small x, as given above.
1375
(cond ((fpgreaterp absx two)
1378
;; log(2*|x| + 1/(|x|+sqrt(1+x^2)))
1380
(fplog (fpplus (fptimes* absx two)
1383
(fproot (bcons (fpplus one
1384
(fptimes* absx absx)))
1389
;; log1p(|x|+x^2/(1+sqrt(1+x^2)))
1390
(let ((x*x (fptimes* absx absx)))
1392
(fplog1p (fpplus absx
1395
(fproot (bcons (fpplus one x*x))
1398
(bcons (fpminus result))
1401
(defun complex-asinh (x y)
1402
;; asinh(z) = -%i * asin(%i*z)
1403
(multiple-value-bind (u v)
1404
(complex-asin (mul -1 y) x)
1405
(values v (bcons (fpminus (cdr u))))))
1407
(defun big-float-asinh (x &optional y)
1409
(multiple-value-bind (u v)
1415
(defun fpasin-core (x)
1416
;; asin(x) = atan(x/(sqrt(1-x^2))
1417
;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
1419
;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
1421
;; If |x| > 1, we need to do something else.
1423
;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
1424
;; = -%i*log(%i*x + %i*sqrt(x^2-1))
1425
;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
1426
;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
1428
(let ((fp-x (cdr (bigfloatp x))))
1429
(cond ((minusp (car fp-x))
1430
;; asin(-x) = -asin(x);
1431
(mul -1 (fpasin (bcons (fpminus fp-x)))))
1432
((fplessp fp-x (cdr bfhalf))
1434
;; asin(x) = atan(x/sqrt(1-x^2))
1436
(fpatan (fpquotient fp-x
1438
(fptimes* (fpdifference (fpone) fp-x)
1439
(fpplus (fpone) fp-x)))
1441
((fpgreaterp fp-x (fpone))
1443
;; asin(x) = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
1445
;; Should we try to do something a little fancier with the
1446
;; argument to log and use log1p for better accuracy?
1447
(let ((arg (fpplus fp-x
1448
(fproot (bcons (fptimes* (fpdifference fp-x (fpone))
1449
(fpplus fp-x (fpone))))
1453
(bcons (fplog arg))))))
1457
;; asin(x) = %pi/2 - atan(sqrt(1-x^2)/x)
1458
(let ((piby2 (fpquotient (fppi) (intofp 2))))
1464
(bcons (fptimes* (fpdifference (fpone)
1466
(fpplus (fpone) fp-x)))
1470
;; asin(x) for real x. X is a bigfloat, and a maxima number (real or
1471
;; complex) is returned.
1473
;; asin(x) = atan(x/(sqrt(1-x^2))
1474
;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
1476
;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
1478
;; If |x| > 1, we need to do something else.
1480
;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
1481
;; = -%i*log(%i*x + %i*sqrt(x^2-1))
1482
;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
1483
;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
1485
($bfloat (fpasin-core x)))
1487
;; Square root of a complex number (xx, yy). Both are bigfloats. FP
1488
;; (non-bigfloat) numbers are returned.
1489
(defun complex-sqrt (xx yy)
1490
(let* ((x (cdr (bigfloatp xx)))
1491
(y (cdr (bigfloatp yy)))
1492
(rho (fpplus (fptimes* x x)
1494
(setf rho (fpplus (fpabs x) (fproot (bcons rho) 2)))
1495
(setf rho (fpplus rho rho))
1496
(setf rho (fpquotient (fproot (bcons rho) 2) (intofp 2)))
1500
(when (fpgreaterp rho (intofp 0))
1501
(setf nu (fpquotient (fpquotient nu rho) (intofp 2)))
1502
(when (fplessp x (intofp 0))
1503
(setf eta (fpabs nu))
1504
(setf nu (if (minusp (car y))
1509
;; asin(z) for complex z = x + %i*y. X and Y are bigfloats. The real
1510
;; and imaginary parts are returned as bigfloat numbers.
1511
(defun complex-asin (x y)
1512
(let ((x (cdr (bigfloatp x)))
1513
(y (cdr (bigfloatp y))))
1514
(multiple-value-bind (re-sqrt-1-z im-sqrt-1-z)
1515
(complex-sqrt (bcons (fpdifference (intofp 1) x))
1516
(bcons (fpminus y)))
1517
(multiple-value-bind (re-sqrt-1+z im-sqrt-1+z)
1518
(complex-sqrt (bcons (fpplus (intofp 1) x))
1520
;; Realpart is atan(x/Re(sqrt(1-z)*sqrt(1+z)))
1521
;; Imagpart is asinh(Im(conj(sqrt(1-z))*sqrt(1+z)))
1523
(fpatan (fpquotient x
1524
(fpdifference (fptimes* re-sqrt-1-z
1526
(fptimes* im-sqrt-1-z
1529
(fpdifference (fptimes* re-sqrt-1-z
1531
(fptimes* im-sqrt-1-z
1532
re-sqrt-1+z)))))))))
1534
(defun big-float-asin (x &optional y)
1536
(multiple-value-bind (u v)
1543
;; tanh(x) for real x. X is a bigfloat, and a bigfloat is returned.
1545
;; X is Maxima bigfloat
1546
;; tanh(x) = D(2*x)/(2+D(2*x))
1547
(let* ((two (intofp 2))
1548
(fp (cdr (bigfloatp x)))
1549
(d (fpexpm1 (fptimes* fp two))))
1551
(fpquotient d (fpplus d two)))))
1553
;; tanh(z), z = x + %i*y. X, Y are bigfloats, and a maxima number is
1555
(defun complex-tanh (x y)
1556
(let* ((fpx (cdr (bigfloatp x)))
1557
(fpy (cdr (bigfloatp y)))
1558
(tv (cdr (tanbigfloat (list y))))
1559
(beta (fpplus (fpone) (fptimes* tv tv)))
1560
(s (cdr (fpsinh x)))
1561
(s^2 (fptimes* s s))
1562
(rho (fproot (bcons (fpplus (fpone) s^2))
1564
(den (fpplus (fpone) (fptimes* beta s^2))))
1565
(values (bcons (fpquotient (fptimes* beta (fptimes* rho s))
1567
(bcons (fpquotient tv den)))))
1569
(defun big-float-tanh (x &optional y)
1571
(multiple-value-bind (u v)
1573
(add u (mul '$%i v)))
1576
;; atanh(x) for real x, |x| <= 1. X is a bigfloat, and a bigfloat is
1579
;; atanh(x) = -atanh(-x)
1580
;; = 1/2*log1p(2*x/(1-x)), x >= 0.5
1581
;; = 1/2*log1p(2*x+2*x*x/(1-x)), x <= 0.5
1583
(let* ((fp-x (cdr (bigfloatp x))))
1584
(cond ((fplessp fp-x (intofp 0))
1585
;; atanh(x) = -atanh(-x)
1586
(mul -1 (fpatanh (bcons (fpminus fp-x)))))
1587
((fpgreaterp fp-x (fpone))
1588
;; x > 1, so use complex version.
1589
(multiple-value-bind (u v)
1590
(complex-atanh x (bcons (intofp 0)))
1591
(add u (mul '$%i v))))
1592
((fpgreaterp fp-x (cdr bfhalf))
1593
;; atanh(x) = 1/2*log1p(2*x/(1-x))
1596
(fplog1p (fpquotient (fptimes* (intofp 2) fp-x)
1597
(fpdifference (fpone) fp-x))))))
1599
;; atanh(x) = 1/2*log1p(2*x + 2*x*x/(1-x))
1600
(let ((2x (fptimes* (intofp 2) fp-x)))
1602
(fptimes* (cdr bfhalf)
1604
(fpquotient (fptimes* 2x fp-x)
1605
(fpdifference (fpone) fp-x)))))))))))
1607
;; Stuff which follows is derived from atanh z = (log(1 + z) - log(1 - z))/2
1608
;; which apparently originates with Kahan's "Much ado" paper.
1610
;; The formulas for eta and nu below can be easily derived from
1611
;; rectform(atanh(x+%i*y)) =
1613
;; 1/4*log(((1+x)^2+y^2)/((1-x)^2+y^2)) + %i/2*(arg(1+x+%i*y)-arg(1-x+%i*(-y)))
1615
;; Expand the argument of log out and divide it out and we get
1617
;; log(((1+x)^2+y^2)/((1-x)^2+y^2)) = log(1+4*x/((1-x)^2+y^2))
1619
;; When y = 0, Im atanh z = 1/2 (arg(1 + x) - arg(1 - x))
1620
;; = if x < -1 then %pi/2 else if x > 1 then -%pi/2 else <whatever>
1622
;; Otherwise, arg(1 - x + %i*(-y)) = - arg(1 - x + %i*y),
1623
;; and Im atanh z = 1/2 (arg(1 + x + %i*y) + arg(1 - x + %i*y)).
1624
;; Since arg(x)+arg(y) = arg(x*y) (almost), we can simplify the
1625
;; imaginary part to
1627
;; arg((1+x+%i*y)*(1-x+%i*y)) = arg((1-x)*(1+x)-y^2+2*y*%i)
1628
;; = atan2(2*y,((1-x)*(1+x)-y^2))
1630
;; These are the eta and nu forms below.
1631
(defun complex-atanh (x y)
1632
(let* ((fpx (cdr (bigfloatp x)))
1633
(fpy (cdr (bigfloatp y)))
1634
(beta (if (minusp (car fpx))
1637
(x-lt-minus-1 (mevalp `((mlessp) ,x -1)))
1638
(x-gt-plus-1 (mevalp `((mgreaterp) ,x 1)))
1639
(y-equals-0 (like y '((bigfloat) 0 0)))
1640
(x (fptimes* beta fpx))
1641
(y (fptimes* beta (fpminus fpy)))
1642
;; Kahan has rho = 4/most-positive-float. What should we do
1643
;; here about that? Our big floats don't really have a
1644
;; most-positive float value.
1646
(t1 (fpplus (fpabs y) rho))
1647
(t1^2 (fptimes* t1 t1))
1648
(1-x (fpdifference (fpone) x))
1649
;; eta = log(1+4*x/((1-x)^2+y^2))/4
1651
(fplog1p (fpquotient (fptimes* (intofp 4) x)
1652
(fpplus (fptimes* 1-x 1-x)
1655
;; If y = 0, then Im atanh z = %pi/2 or -%pi/2.
1656
;; Otherwise nu = 1/2*atan2(2*y,(1-x)*(1+x)-y^2)
1658
;; EXTRA FPMINUS HERE TO COUNTERACT FPMINUS IN RETURN VALUE
1659
(fpminus (if x-lt-minus-1 (cdr ($bfloat '((mquotient) $%pi 2))) (if x-gt-plus-1 (cdr ($bfloat '((mminus) ((mquotient) $%pi 2)))) (merror "COMPLEX-ATANH: HOW DID I GET HERE?"))))
1660
(fptimes* (cdr bfhalf)
1662
(fptimes* (intofp 2) y)
1663
(fpdifference (fptimes* 1-x
1666
(values (bcons (fptimes* beta eta))
1667
;; WTF IS FPMINUS DOING HERE ??
1668
(bcons (fpminus (fptimes* beta nu))))))
1671
(defun big-float-atanh (x &optional y)
1673
(multiple-value-bind (u v)
1675
(add u (mul '$%i v)))
1678
;; acos(x) for real x. X is a bigfloat, and a maxima number is returned.
1680
;; acos(x) = %pi/2 - asin(x)
1681
($bfloat (add (div '$%pi 2)
1682
(mul -1 (fpasin-core x)))))
1684
(defun complex-acos (x y)
1685
(let ((x (cdr (bigfloatp x)))
1686
(y (cdr (bigfloatp y))))
1687
(multiple-value-bind (re-sqrt-1-z im-sqrt-1-z)
1688
(complex-sqrt (bcons (fpdifference (intofp 1) x))
1689
(bcons (fpminus y)))
1690
(multiple-value-bind (re-sqrt-1+z im-sqrt-1+z)
1691
(complex-sqrt (bcons (fpplus (intofp 1) x))
1694
(fptimes* (intofp 2)
1695
(fpatan (fpquotient re-sqrt-1-z
1699
(fptimes* re-sqrt-1+z im-sqrt-1-z)
1700
(fptimes* im-sqrt-1+z re-sqrt-1-z)))))))))
1703
(defun big-float-acos (x &optional y)
1705
(multiple-value-bind (u v)
1707
(add u (mul '$%i v)))
1710
(defun complex-log (x y)
1711
(let* ((x (cdr (bigfloatp x)))
1712
(y (cdr (bigfloatp y)))
1713
(t1 (let (($float2bf t))
1714
;; No warning message, please.
1717
(rho (fpplus (fptimes* x x)
1721
(beta (fpmax abs-x abs-y))
1722
(theta (fpmin abs-x abs-y)))
1723
(values (if (or (fpgreaterp t1 beta)
1725
(fpquotient (fplog1p (fpplus (fptimes* (fpdifference beta (fpone))
1726
(fpplus beta (fpone)))
1727
(fptimes* theta theta)))
1729
(fpquotient (fplog rho) (intofp 2)))
1732
(defun big-float-log (x &optional y)
1734
(multiple-value-bind (u v)
1736
(add (bcons u) (mul '$%i (bcons v))))
1737
(let ((fp-x (cdr (bigfloatp x))))
1738
(if (fplessp fp-x (intofp 0))
1739
;; ??? Do we want to return an exact %i*%pi or a float
1741
(add (bcons (fplog (fpminus fp-x)))
1742
(mul '$%i (bcons (fppi))))
1743
(bcons (fplog fp-x))))))
1745
(defun big-float-sqrt (x &optional y)
1747
(multiple-value-bind (u v)
1749
(add (bcons u) (mul '$%i (bcons v))))
1750
(let ((fp-x (cdr (bigfloatp x))))
1751
(if (fplessp fp-x (intofp 0))
1752
(mul '$%i (bcons (fproot (bcons (fpminus fp-x)) 2)))
1753
(bcons (fproot x 2))))))
1266
1757
#-gcl (:load-toplevel)
added
removed
src/float.lisp
