| blob | 02bf7b09dae1fbdfa967e9eac32e85aeb2a06cf5 |
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 ;; EXPERIMENTAL BIGFLOAT PACKAGE VERSION 2- USING BINARY MANTISSA
16 ;; AND POWER-OF-2 EXPONENT.
17 ;; EXPONENTS MAY BE BIG NUMBERS NOW (AUG. 1975 --RJF)
18 ;; Modified: July 1979 by CWH to run on the Lisp Machine and to comment
19 ;; the code.
20 ;; August 1980 by CWH to run on Multics and to install
21 ;; new FIXFLOAT.
22 ;; December 1980 by JIM to fix BIGLSH not to pass LSH a second
23 ;; argument with magnitude greater than MACHINE-FIXNUM-PRECISION.
25 ;; Number of bits of precision in a fixnum and in the fields of a flonum for
26 ;; a particular machine. These variables should only be around at eval
27 ;; and compile time. These variables should probably be set up in a prelude
28 ;; file so they can be accessible to all Macsyma files.
35 ;; External variables
38 "If TRUE, no MAXIMA-ERROR message is printed when a floating point number is
42 "Controls the conversion of bigfloat numbers to rational numbers. If
43 FALSE, RATEPSILON will be used to control the conversion (this results in
44 relatively small rational numbers). If TRUE, the rational number generated
48 "If TRUE, printing of bigfloat numbers will truncate trailing zeroes.
52 "Controls the number of significant digits printed for floats. If
53 0, then full precision is used."
54 fixnum)
60 "Number of decimal digits of precision to use when creating new bigfloats.
87 ;; Internal specials
89 ;; Number of bits of precision in the mantissa of newly created bigfloats.
90 ;; FPPREC = ($FPPREC+1)*(Log base 2 of 10)
94 ;; FPROUND uses this to return a second value, i.e. it sets it before
95 ;; returning. This number represents the number of binary digits its input
96 ;; bignum had to be shifted right to be aligned into the mantissa. For
97 ;; example, aligning 1 would mean shifting it FPPREC-1 places left, and
98 ;; aligning 7 would mean shifting FPPREC-3 places left.
102 ;; *DECFP = T if the computation is being done in decimal radix. NIL implies
103 ;; base 2. Decimal radix is used only during output.
115 ;; Representation of a Bigfloat: ((BIGFLOAT SIMP precision) mantissa exponent)
116 ;; precision -- number of bits of precision in the mantissa.
117 ;; precision = (integer-length mantissa)
118 ;; mantissa -- a signed integer representing a fractional portion computed by
119 ;; fraction = (// mantissa (^ 2 precision)).
120 ;; exponent -- a signed integer representing the scale of the number.
121 ;; The actual number represented is (* fraction (^ 2 exponent)).
137 q)
139 ;; FPSCAN is called by lexical scan when a
140 ;; bigfloat is encountered. For example, 12.01B-3
141 ;; would be the result of (FPSCAN '(/1 /2) '(/0 /1) '(/- /3))
142 ;; Arguments to FPSCAN are a list of characters to the left of the
143 ;; decimal point, to the right of the decimal point, and in the exponent.
150 $float temp)
161 ;; Converts the bigfloat L to list of digits including |.| and the
162 ;; exponent marker |b|. The number of significant digits is controlled
163 ;; by $fpprintprec.
214 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
215 ;; support printing of bfloats.
221 ;; Compute the (decimal exponent) of the bfloat number X.
225 ;; FIXME: do something better than printing and reading
226 ;; the result.
230 ;; Print the bfloat X with FDIGITS after the decimal
231 ;; point. This means, roughtly, FDIGITS+1 significant
232 ;; digits.
235 1
241 ;; Depending on the value of k, move the decimal
242 ;; point. DIGITS was printed assuming the decimal point
243 ;; is after the first digit. But if fdigits = 0, fpformat
244 ;; actually printed out one too many digits, so we need
245 ;; to remove that.
251 ;; Put the leading decimal and then some zeroes
256 ;; The number is scaled by 10^k. Do this by
257 ;; putting the decimal point in the right place,
258 ;; appending zeroes if needed.
262 digits
274 ;; Append some zeroes to get the desired number of digits
280 len
283 1
287 0
293 ;; Exponent overflow
297 ;; The hairy case
301 d)
302 nil))
307 nil)))
317 fstr flen lpoint tpoint)
321 ;; See CLHS 22.3.3.2. "If the parameter d is
322 ;; omitted, ... [and] if the fraction to be
323 ;; printed is zero then a single zero digit should
324 ;; appear after the decimal point." So we need to
325 ;; subtract one from here because we're going to
326 ;; add an extra 0 digit later.
338 w spaceleft overflowchar)
340 ;; Significand overflow; output the overflow char
354 ;; Add a zero if we need it. Which means
355 ;; we figured out we need one above, or
356 ;; another condition. Basically, append a
357 ;; zero if there are no width constraints
358 ;; and if the last char to print was a
359 ;; decimal (so the trailing fraction is
360 ;; zero.)
366 exponentchar
368 stream)
376 ;; NOTE: This is a modified version of FORMAT-FIXED-AUX from CMUCL to
377 ;; support printing of bfloats.
382 ;; Compute the (decimal exponent) of the bfloat number X.
386 ;; FIXME: do something better than printing and reading
387 ;; the result.
391 ;; Print the bfloat X with FDIGITS after the decimal
392 ;; point. To do this we need to know the exponent because
393 ;; fpformat always produces exponential output. If the
394 ;; exponent is E, and we want FDIGITS after the decimal
395 ;; point, we need FDIGITS + E digits printed.
397 #+nil
408 #+nil
413 ;; Depending on the value of scale, move the decimal
414 ;; point. DIGITS was printed assuming the decimal point
415 ;; is after the first digit. But if fdigits = 0, fpformat
416 ;; actually printed out one too many digits, so we need
417 ;; to remove that.
418 #+nil
420 #+nil
423 ;; Figure out where the decimal point should go. An
424 ;; exponent of 0 means the decimal is after the first
425 ;; digit.
431 #+nil
441 #+nil
444 ;; Append some zeroes to get the desired number of digits
450 len
453 1
460 ;;if caller specifically requested no fraction digits, suppress the
461 ;;optional trailing zero
465 ;;optional leading zero
470 ;;optional trailing zero
476 ;;field width overflow
478 t)
487 nil))))))
489 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
490 ;; support printing of bfloats.
496 ;; Compute the (decimal exponent) of the bfloat number X.
500 ;; FIXME: do something better than printing and reading
501 ;; the result.
505 ;; Print the bfloat X with FDIGITS after the decimal
506 ;; point. This means, roughtly, FDIGITS+1 significant
507 ;; digits.
510 1
516 ;; Depending on the value of k, move the decimal
517 ;; point. DIGITS was printed assuming the decimal point
518 ;; is after the first digit. But if fdigits = 0, fpformat
519 ;; actually printed out one too many digits, so we need
520 ;; to remove that.
526 ;; Put the leading decimal and then some zeroes
531 ;; The number is scaled by 10^k. Do this by
532 ;; putting the decimal point in the right place,
533 ;; appending zeroes if needed.
537 digits
549 ;; Append some zeroes to get the desired number of digits
555 len
558 1
562 ;; Default d if omitted. The procedure is taken directly from
563 ;; the definition given in the manual (CLHS 22.3.3.3), and is
564 ;; not very efficient, since we generate the digits twice.
565 ;; Future maintainers are encouraged to improve on this.
566 ;;
567 ;; It's also not very clear whether q in the spec is the
568 ;; number of significant digits or not. I (rtoy) think it
569 ;; makes more sense if q is the number of significant digits.
570 ;; That way 1d300 isn't printed as 1 followed by 300 zeroes.
571 ;; Exponential notation would be used instead.
586 ;; Use dd fraction digits, even if that would cause
587 ;; the width to be exceeded. We choose accuracy over
588 ;; width in this case.
590 ww
591 dd
592 0
593 ovf pad)
594 ovf
599 w
600 orig-d
602 ovf pad exponentchar)))))))
604 ;; Tells you if you have a bigfloat object. BUT, if it is a bigfloat,
605 ;; it will normalize it by making the precision of the bigfloat match
606 ;; the current precision setting in fpprec. And it will also convert
607 ;; bogus zeroes (mantissa is zero, but exponent is not) to a true
608 ;; zero.
610 ;; A bigfloat object looks like '((bigfloat simp <prec>) <mantissa> <exp>)
614 ;; Precision matches. (Should we fix up bogus bigfloat
615 ;; zeros?)
618 ;; Current precision is higher than bigfloat precision.
619 ;; Scale up mantissa and adjust exponent to get the
620 ;; correct precision.
624 ;; Current precision is LOWER than bigfloat precision.
625 ;; Round the number to the desired precision.
628 ;; Fix up any bogus zeros that we might have created.
655 exp))
660 $ratepsilon)))))
677 fprateps)))
694 ;; BLECHH, I HATE ENUMERATING CASES. IS THERE A BETTER WAY ??
700 ;; NOT SURE THE FOLLOWING REALLY WORKS
701 ;; #+(and cmu double-double)
702 ;; (kernel:double-double-float
703 ;; (values most-negative-double-double-float most-positive-double-double-float))
704 ))
706 ;; Convert a floating point number into a bigfloat.
715 ;; Need to check for zero because different lisps return different
716 ;; values for integer-decode-float of a 0. In particular CMUCL
717 ;; returns 0, -1075. A bigfloat zero needs to have an exponent and
718 ;; mantissa of zero.
723 ;; Scale frac to the desired number of bits, and adjust the
724 ;; exponent accordingly.
729 ;; Convert a bigfloat into a floating point number.
736 ;; Round the mantissa to the number of bits of precision of the
737 ;; machine, and then convert it to a floating point fraction. We
738 ;; have 0.5 <= mantissa < 1
740 ;; Multiply the mantissa by the exponent portion. I'm not sure
741 ;; why the exponent computation is so complicated.
742 ;;
743 ;; GCL doesn't signal overflow from scale-float if the number
744 ;; would overflow. We have to do it this way. 0.5 <= mantissa <
745 ;; 1. The largest double-float is .999999 * 2^1024. So if the
746 ;; exponent is 1025 or higher, we have an overflow.
752 ;; New machine-independent version of FIXFLOAT. This may be buggy. - CWH
753 ;; It is buggy! On the PDP10 it dies on (RATIONALIZE -1.16066076E-7)
754 ;; which calls FLOAT on some rather big numbers. ($RATEPSILON is approx.
755 ;; 7.45E-9) - JPG
761 ;; Takes a flonum arg and returns a rational number corresponding to the flonum
762 ;; in the form of a dotted pair of two integers. Since the denominator will
763 ;; always be a positive power of 2, this number will not always be in lowest
764 ;; terms.
780 x))
829 ;; R is a Maxima ratio, represented as a list of the numerator and
830 ;; denominator. FLOAT-RATIO doesn't like it if the numerator is 0,
831 ;; so handle that here.
836 ;; This is borrowed from CMUCL (float-ratio-float), and modified for
837 ;; converting ratios to Maxima's bfloat numbers.
846 ;;
847 ;; Strip any trailing zeros from the denominator and move it into the scale
848 ;; factor (to minimize the size of the operands.)
853 ;;
854 ;; Guess how much we need to scale by from the magnitudes of the numerator
855 ;; and denominator. We want one extra bit for a guard bit.
891 fraction-and-guard
915 ;; Bigfloat atan
920 ;; |x| > 1.
921 ;;
922 ;; Use A&S 4.4.5:
923 ;; atan(x) + acot(x) = +/- pi/2 (+ for x >= 0, - for x < 0)
924 ;;
925 ;; and A&S 4.4.8
926 ;; acot(z) = atan(1/z)
933 ;; |x| > 1/2
934 ;;
935 ;; Use A&S 4.4.42, third formula:
936 ;;
937 ;; atan(z) = z/(1+z^2)*[1 + 2/3*r + (2*4)/(3*5)*r^2 + ...]
938 ;;
939 ;; r = z^2/(1+z^2)
951 ;; |x| <= 1/2. Use Taylor series (A&S 4.4.42, first
952 ;; formula).
961 ;; atan(y/x) taking into account the quadrant. (Also equal to
962 ;; arg(x+%i*y).)
965 ;; atan(y/0) = atan(inf), but what sign?
969 ;; We're on the negative imaginary axis, so -pi/2.
972 ;; The positive imaginary axis, so +pi/2
975 ;; x > 0. atan(y/x) is the correct value.
978 ;; x < 0, and y > 0. We're in quadrant II, so the angle we
979 ;; want is pi+atan(y/x).
982 ;; x <= 0 and y <= 0. We're in quadrant III, so the angle we
983 ;; want is atan(y/x)-pi.
994 ;; Returns a list of a mantissa and an exponent.
1004 ;; It seems to me that this function gets called on an integer
1005 ;; and returns the mantissa portion of the mantissa/exponent pair.
1007 ;; "STICKY BIT" CALCULATION FIXED 10/14/75 --RJF
1008 ;; BASE must not get temporarily bound to NIL by being placed
1009 ;; in a PROG list as this will confuse stepping programs.
1015 ;;*M will be positive if the precision of the argument is greater than
1016 ;;the current precision being used.
1021 ;;FPSHIFT is essentially LSH.
1028 ;ONLY ZEROES SHIFTED OFF
1040 ;; Compute (* L (expt d n)) where D is 2 or 10 depending on
1041 ;; *decfp. Throw away an fractional part by truncating to zero.
1045 ;; Left shift of negative number requires some
1046 ;; care. (That is, (truncate l (expt 2 n)), but use
1047 ;; shifts instead.)
1057 ;; Bignum LSH -- N is assumed (and declared above) to be a fixnum.
1058 ;; This isn't really LSH, since the sign bit isn't propagated when
1059 ;; shifting to the right, i.e. (BIGLSH -100 -3) = -40, whereas
1060 ;; (LSH -100 -3) = 777777777770 (on a 36 bit machine).
1061 ;; This actually computes (* X (EXPT 2 N)). As of 12/21/80, this function
1062 ;; was only called by FPSHIFT. I would like to hear an argument as why this
1063 ;; is more efficient than simply writing (* X (EXPT 2 N)). Is the
1064 ;; intermediate result created by (EXPT 2 N) the problem? I assume that
1065 ;; EXPT tries to LSH when possible.
1071 ;; Either we are shifting a fixnum to the right, or shifting
1072 ;; a fixnum to the left, but not far enough left for it to become
1073 ;; a bignum.
1077 ;; The form which follows is nearly identical to (ASH X N), however
1078 ;; (ASH -100 -20) = -1, whereas (BIGLSH -100 -20) = 0.
1082 ;; If we get here, then either X is a bignum or our answer is
1083 ;; going to be a bignum.
1090 ;; Isn't this the kind of optimization that compilers are
1091 ;; supposed to make?
1096 ;; exp(x)
1097 ;;
1098 ;; For negative x, use exp(-x) = 1/exp(x)
1099 ;;
1100 ;; For x > 0, exp(x) = exp(r+y) = exp(r) * exp(y), where x = r + y and
1101 ;; r = floor(x).
1117 ;; exp(x) for small x, using Taylor series.
1128 ;; Does one higher precision to round correctly.
1129 ;; A and B are each a list of a mantissa and an exponent.
1149 min))
1154 max))
1156 ;; The following functions compute bigfloat values for %e, %pi,
1157 ;; %gamma, and log(2). For each precision, the computed value is
1158 ;; cached in a hash table so it doesn't need to be computed again.
1159 ;; There are functions to return the hash table or clear the hash
1160 ;; table, for debugging.
1161 ;;
1162 ;; Note that each of these return a bigfloat number, but without the
1163 ;; bigfloat tag.
1164 ;;
1165 ;; See
1166 ;; https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2910437&group_id=4933
1167 ;; for an explanation.
1172 value
1175 table)
1183 value
1186 table)
1194 value
1197 table)
1205 value
1208 table)
1212 ;; This doesn't need a hash table because there's never a problem with
1213 ;; using a high precision value and rounding to a lower precision
1214 ;; value because 1 is always an exact bfloat.
1220 ;;----------------------------------------------------------------------------;;
1221 ;;
1222 ;; The values of %e, %pi, %gamma and log(2) are computed by the technique of
1223 ;; binary splitting. See http://www.ginac.de/CLN/binsplit.pdf for details.
1224 ;;
1225 ;; Volker van Nek, Sept. 2014
1227 ;;
1228 ;; Euler's number E
1229 ;;
1233 ;;
1234 ;; Taylor: %e = sum(s[i] ,i,0,inf) where s[i] = 1/i!
1235 ;;
1240 ;;
1241 ;; binary splitting:
1242 ;;
1243 ;; 1
1244 ;; s[i] = ----------------------
1245 ;; q[0]*q[1]*q[2]*..*q[i]
1246 ;;
1247 ;; where q[0] = 1
1248 ;; q[i] = i
1249 ;;
1261 ;;
1262 ;; stop when i! > 2^fpprec
1263 ;;
1264 ;; log(i!) = sum(log(k), k,1,i) > fpprec * log(2)
1265 ;;
1272 ;;
1273 ;;----------------------------------------------------------------------------;;
1274 ;;
1275 ;; PI
1276 ;;
1280 ;;
1281 ;; Chudnovsky & Chudnovsky:
1282 ;;
1283 ;; C^(3/2)/(12*%pi) = sum(s[i], i,0,inf),
1284 ;;
1285 ;; where s[i] = (-1)^i*(6*i)!*(A*i+B) / (i!^3*(3*i)!*C^(3*i))
1286 ;;
1287 ;; and A = 545140134, B = 13591409, C = 640320
1288 ;;
1292 ;; STEP 1:
1293 ;; compute n/d = sqrt(10005) :
1294 ;;
1295 ;; n[0] n[i+1] = n[i]^2+a*d[i]^2 n[inf]
1296 ;; quadratic Heron: x[0] = ----, , sqrt(a) = ------
1297 ;; d[0] d[i+1] = 2*n[i]*d[i] d[inf]
1298 ;;
1304 ;; STEP 2:
1305 ;; divide C^(3/2)/12 = 3335*2^7*sqrt(10005)
1306 ;; by Chudnovsky-sum = tt/qq :
1307 ;;
1309 ;; fpprec*log(2)/log(C^3/(24*6*2*6))
1315 ;;
1316 ;; The returned n and d serve as start values for the iteration.
1317 ;; n/d = sqrt(10005) with a precision of p = ceiling(prec/2^nr) bits
1318 ;; where nr is the number of needed iterations.
1319 ;;
1334 ;;
1335 ;; binary splitting:
1336 ;;
1337 ;; a[i] * p[0]*p[1]*p[2]*..*p[i]
1338 ;; s[i] = -----------------------------
1339 ;; q[0]*q[1]*q[2]*..*q[i]
1340 ;;
1341 ;; where a[0] = B
1342 ;; p[0] = q[0] = 1
1343 ;; a[i] = A*i+B
1344 ;; p[i] = - (6*i-5)*(2*i-1)*(6*i-1)
1345 ;; q[i] = C^3/24*i^3
1346 ;;
1365 ;;
1366 ;;----------------------------------------------------------------------------;;
1367 ;;
1368 ;; Euler-Mascheroni constant GAMMA
1369 ;;
1373 ;;
1374 ;; Brent-McMillan algorithm
1375 ;;
1376 ;; Let
1377 ;; alpha = 4.970625759544
1378 ;;
1379 ;; n > 0 and N-1 >= alpha*n
1380 ;;
1381 ;; H(k) = sum(1/i, i,1,k)
1382 ;;
1383 ;; S = sum(H(k)*(n^k/k!)^2, k,0,N-1)
1384 ;;
1385 ;; I = sum((n^k/k!)^2, k,0,N-1)
1386 ;;
1387 ;; T = 1/(4*n)*sum((2*k)!^3/(k!^4*(16*n)^(2*k)), k,0,2*n-1)
1388 ;;
1389 ;; and
1390 ;; %gamma = S/I - T/I^2 - log(n)
1391 ;;
1392 ;; Then
1393 ;; |%gamma - gamma| < 24 * e^(-8*n)
1394 ;;
1395 ;; (Corollary 2, Remark 2, Brent/Johansson http://arxiv.org/pdf/1312.0039v1.pdf)
1396 ;;
1404 sumi sumi2 ;; I and I^2
1407 ;;
1408 ;; sums = vv/(qq*dd)
1409 ;; sumi = tt/qq
1410 ;; sums/sumi = vv/(qq*dd)*qq/tt = vv/(dd*tt)
1411 ;;
1415 ;;
1417 ;;
1418 ;; sumt = 1/(4*n)*ttt/qqq
1419 ;; sumt/sumi2 = ttt/(4*n*qqq*sumi2)
1420 ;;
1423 ;; %gamma :
1425 ;;
1426 ;; split S and I simultaneously:
1427 ;;
1428 ;; summands I[0] = 1, I[i]/I[i-1] = n^2/i^2
1429 ;;
1430 ;; S[0] = 0, S[i]/S[i-1] = n^2/i^2*H(i)/H(i-1)
1431 ;;
1432 ;; p[0]*p[1]*p[2]*..*p[i]
1433 ;; I[i] = ----------------------
1434 ;; q[0]*q[1]*q[2]*..*q[i]
1435 ;;
1436 ;; where p[0] = n^2
1437 ;; q[0] = 1
1438 ;; p[i] = n^2
1439 ;; q[i] = i^2
1440 ;; c[0] c[1] c[2] c[i]
1441 ;; S[i] = H[i] * I[i], where H[i] = ---- + ---- + ---- + .. + ----
1442 ;; d[0] d[1] d[2] d[i]
1443 ;; and c[0] = 0
1444 ;; d[0] = 1
1445 ;; c[i] = 1
1446 ;; d[i] = i
1447 ;;
1467 ;;
1468 ;; split 4*n*T:
1469 ;;
1470 ;; summands T[0] = 1, T[i]/T[i-1] = (2*i-1)^3/(32*i*n^2)
1471 ;;
1472 ;; p[0]*p[1]*p[2]*..*p[i]
1473 ;; T[i] = ----------------------
1474 ;; q[0]*q[1]*q[2]*..*q[i]
1475 ;;
1476 ;; where p[0] = q[0] = 1
1477 ;; p[i] = (2*i-1)^3
1478 ;; q[i] = 32*i*n^2
1479 ;;
1495 ;;
1496 ;;----------------------------------------------------------------------------;;
1497 ;;
1498 ;; log(2) = 18*L(26) - 2*L(4801) + 8*L(8749)
1499 ;;
1500 ;; where L(k) = atanh(1/k)
1501 ;;
1502 ;; see http://numbers.computation.free.fr/Constants/constants.html
1503 ;;
1504 ;;;(defun $log2 () (bcons (comp-log2))) ;; checked against reference table
1505 ;;
1513 ;;
1514 ;; Taylor: atanh(1/k) = sum(s[i], i,0,inf)
1515 ;;
1516 ;; where s[i] = 1/((2*i+1)*k^(2*i+1))
1517 ;;
1523 ;;
1524 ;; binary splitting:
1525 ;; 1
1526 ;; s[i] = -----------------------------
1527 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1528 ;;
1529 ;; where b[0] = 1
1530 ;; q[0] = k
1531 ;; b[i] = 2*i+1
1532 ;; q[i] = k^2
1533 ;;
1547 ;;
1548 ;;----------------------------------------------------------------------------;;
1549 ;;
1550 ;; log(n) = log(n/2^k) + k*log(2)
1551 ;;
1552 ;;;(defun $log10 () (bcons (log-n 10))) ;; checked against reference table
1553 ;;
1562 ;; choose k so that |n/2^k - 1| is as small as possible:
1564 ;; now |n/2^k - 1| <= 1/3
1568 ;;
1569 ;; log(1+u/v) = 2 * sum(s[i], i,0,inf)
1570 ;;
1571 ;; where s[i] = (u/(2*v+u))^(2*i+1)/(2*i+1)
1572 ;;
1587 ;;
1588 ;; binary splitting:
1589 ;;
1590 ;; p[0]*p[1]*p[2]*..*p[i]
1591 ;; s[i] = -----------------------------
1592 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1593 ;;
1594 ;; where b[0] = 1
1595 ;; p[0] = u
1596 ;; q[0] = w = 2*v+u
1597 ;; b[i] = 2*i+1
1598 ;; p[i] = u^2
1599 ;; q[i] = w^2
1600 ;;
1615 ;;
1616 ;;----------------------------------------------------------------------------;;
1617 \f
1624 x
1641 ; COMPUTE STICKY BIT
1643 ; MAKE ROOM FOR GUARD DIGIT & STICKY BIT
1663 ;; Don't use the symbol BASE since it is SPECIAL.
1671 bas)))
1673 ;; NN is positive or negative integer
1712 $%i)))))
1727 ;; The number of extra bits required
1733 ;; Special case for a = 0b0. General algorithm loops endlessly in that case.
1735 ;; Unlike many or maybe all of the other functions named FP-something,
1736 ;; FPROOT assumes it is called with an argument like
1737 ;; '((BIGFLOAT ...) FOO BAR) instead of '(FOO BAR).
1738 ;; However FPROOT does return something like '(FOO BAR).
1765 fans
1769 ;; If A is a bigfloat, be sure to round it to the current precision.
1770 ;; (See Bug 2543079 for one of the symptoms.)
1796 ;; FPPREC 16, to get rid of the miniscule imaginary
1797 ;; part of the a+bi answer.
1804 e))
1812 ;; THIS VERSION OF FPSIN COMPUTES SIN OR COS TO PRECISION FPPREC,
1813 ;; BUT CHECKS FOR THE POSSIBILITY OF CATASTROPHIC CANCELLATION DURING
1814 ;; ARGUMENT REDUCTION (E.G. SIN(N*%PI+EPSILON))
1815 ;; *FPSINCHECK* WILL CAUSE PRINTOUT OF ADDITIONAL INFO WHEN
1816 ;; EXTRA PRECISION IS NEEDED FOR SIN/COS CALCULATION. KNOWN
1817 ;; BAD FEATURES: IT IS NOT NECESSARY TO USE EXTRA PRECISION FOR, E.G.
1818 ;; SIN(PI/2), WHICH IS NOT NEAR ZERO, BUT EXTRA
1819 ;; PRECISION IS USED SINCE IT IS NEEDED FOR COS(PI/2).
1820 ;; PRECISION SEEMS TO BE 100% SATSIFACTORY FOR LARGE ARGUMENTS, E.G.
1821 ;; SIN(31415926.0B0), BUT LESS SO FOR SIN(3.1415926B0). EXPLANATION
1822 ;; NOT KNOWN. (9/12/75 RJF)
1869 ;; Compute SIN or COS in (0,PI/2). FL is T for SIN, NIL for COS.
1870 ;;
1871 ;; Use Taylor series
1889 x
1896 ;; Calculate the integer part of a floating point number that is represented as
1897 ;; a list
1898 ;;
1899 ;; (MANTISSA EXPONENT)
1900 ;;
1901 ;; The special variable fpprec should be bound to the precision (in bits) of the
1902 ;; number. This encodes how many bits are known of the result and also a right
1903 ;; shift. The pair denotes the number MANTISSA * 2^(EXPONENT - FPPREC), of which
1904 ;; FPPREC bits are known.
1905 ;;
1906 ;; If EXPONENT is large and positive then we might not have enough information
1907 ;; to calculate the integer part. Specifically, we only have enough information
1908 ;; if EXPONENT < FPPREC. If that isn't the case, we signal a Maxima error.
1923 ;;; Computes the log of a bigfloat number.
1924 ;;;
1925 ;;; Uses the series
1926 ;;;
1927 ;;; log(1+x) = sum((x/(x+2))^(2*n+1)/(2*n+1),n,0,inf);
1928 ;;;
1929 ;;;
1930 ;;; INF x 2 n + 1
1931 ;;; ==== (-----)
1932 ;;; \ x + 2
1933 ;;; = 2 > --------------
1934 ;;; / 2 n + 1
1935 ;;; ====
1936 ;;; n = 0
1937 ;;;
1938 ;;;
1939 ;;; which converges for x > 0.
1940 ;;;
1941 ;;; Note that FPLOG is given 1+X, not X.
1942 ;;;
1943 ;;; However, to aid convergence of the series, we scale 1+x until 1/e
1944 ;;; < 1+x <= e.
1945 ;;;
1953 ;; Scale X until 1/e < X <= E. ANS keeps track of how
1954 ;; many factors of E were used. Set X to NIL if X is E.
1967 ;; Prepare X for the series. The series is for 1 + x, so
1968 ;; get x from our X. TERM is (x/(x+2)). X becomes
1969 ;; (x/(x+2))^2.
1973 ;; Sum the series until the sum (in ANS) doesn't change
1974 ;; anymore.
1990 \f
1991 ;;;; Bigfloat implementations of special functions.
1992 ;;;;
1994 ;;; This is still a bit messy. Some functions here take bigfloat
1995 ;;; numbers, represented by ((bigfloat) <mant> <exp>), but others want
1996 ;;; just the FP number, represented by (<mant> <exp>). Likewise, some
1997 ;;; return a bigfloat, some return just the FP.
1998 ;;;
1999 ;;; This needs to be systemized somehow. It isn't helped by the fact
2000 ;;; that some of the routines above also do the samething.
2001 ;;;
2002 ;;; The implementation for the special functions for a complex
2003 ;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex
2004 ;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in
2005 ;;; Iserles and Powell (eds.) "The State of the Art in Numerical
2006 ;;; Analysis", pp 165-211, Clarendon Press, 1987
2008 ;; Compute exp(x) - 1, but do it carefully to preserve precision when
2009 ;; |x| is small. X is a FP number, and a FP number is returned. That
2010 ;; is, no bigfloat stuff.
2012 ;; What is the right breakpoint here? Is 1 ok? Perhaps 1/e is better?
2014 ;; exp(x) - 1
2017 ;; Use Taylor series for exp(x) - 1
2026 ans))))
2028 ;; log(1+x) for small x. X is FP number, and a FP number is returned.
2030 ;; Use the same series as given above for fplog. For small x we use
2031 ;; the series, otherwise fplog is accurate enough.
2046 ;; sinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2048 ;; X must be a maxima bigfloat
2050 ;; See, for example, Hart et al., Computer Approximations, 6.2.27:
2051 ;;
2052 ;; sinh(x) = 1/2*(D(x) + D(x)/(1+D(x)))
2053 ;;
2054 ;; where D(x) = exp(x) - 1.
2055 ;;
2056 ;; But for negative x, use sinh(x) = -sinh(-x) because D(x)
2057 ;; approaches -1 for large negative x.
2059 ;; Special case: x=0. Return immediately.
2062 ;; x is positive.
2067 ;; x is negative.
2072 ;; The rectform for sinh for complex args should be numerically
2073 ;; accurate, so return nil in that case.
2077 ;; asinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2079 ;; asinh(x) = sign(x) * log(|x| + sqrt(1+x*x))
2080 ;;
2081 ;; And
2082 ;;
2083 ;; asinh(x) = x, if 1+x*x = 1
2084 ;; = sign(x) * (log(2) + log(x)), large |x|
2085 ;; = sign(x) * log(2*|x| + 1/(|x|+sqrt(1+x*x))), if |x| > 2
2086 ;; = sign(x) * log1p(|x|+x^2/(1+sqrt(1+x*x))), otherwise.
2087 ;;
2088 ;; But I'm lazy right now and we only implement the last 2 cases.
2089 ;; We should implement all cases.
2095 result)
2096 ;; We only use two formulas here. |x| <= 2 and |x| > 2. Should
2097 ;; we add one for very big x and one for very small x, as given above.
2099 ;; |x| > 2
2100 ;;
2101 ;; log(2*|x| + 1/(|x|+sqrt(1+x^2)))
2109 ;; |x| <= 2
2110 ;;
2111 ;; log1p(|x|+x^2/(1+sqrt(1+x^2)))
2123 ;; asinh(z) = -%i * asin(%i*z)
2136 ;; asin(x) = atan(x/(sqrt(1-x^2))
2137 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2138 ;;
2139 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2140 ;;
2141 ;; If |x| > 1, we need to do something else.
2142 ;;
2143 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2144 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2145 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2146 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2150 ;; asin(-x) = -asin(x);
2153 ;; 0 <= x < 1/2
2154 ;; asin(x) = atan(x/sqrt(1-x^2))
2162 ;; x > 1
2163 ;; asin(x) = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2164 ;;
2165 ;; Should we try to do something a little fancier with the
2166 ;; argument to log and use log1p for better accuracy?
2175 ;; 1/2 <= x <= 1
2176 ;; asin(x) = %pi/2 - atan(sqrt(1-x^2)/x)
2184 fp-x)))))))))
2186 ;; asin(x) for real x. X is a bigfloat, and a maxima number (real or
2187 ;; complex) is returned.
2189 ;; asin(x) = atan(x/(sqrt(1-x^2))
2190 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2191 ;;
2192 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2193 ;;
2194 ;; If |x| > 1, we need to do something else.
2195 ;;
2196 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2197 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2198 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2199 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2203 ;; Square root of a complex number (xx, yy). Both are bigfloats. FP
2204 ;; (non-bigfloat) numbers are returned.
2222 rho))))
2225 ;; asin(z) for complex z = x + %i*y. X and Y are bigfloats. The real
2226 ;; and imaginary parts are returned as bigfloat numbers.
2236 ;; Realpart is atan(x/Re(sqrt(1-z)*sqrt(1+z)))
2237 ;; Imagpart is asinh(Im(conj(sqrt(1-z))*sqrt(1+z)))
2243 ;; Check for division by zero. If we would divide
2244 ;; by zero, return pi/2 or -pi/2 according to the
2245 ;; sign of X.
2265 ;; tanh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2267 ;; X is Maxima bigfloat
2268 ;; tanh(x) = D(2*x)/(2+D(2*x))
2274 ;; tanh(z), z = x + %i*y. X, Y are bigfloats, and a maxima number is
2275 ;; returned.
2292 ;; atanh(x) for real x, |x| <= 1. X is a bigfloat, and a bigfloat is
2293 ;; returned.
2295 ;; atanh(x) = -atanh(-x)
2296 ;; = 1/2*log1p(2*x/(1-x)), x >= 0.5
2297 ;; = 1/2*log1p(2*x+2*x*x/(1-x)), x <= 0.5
2301 ;; atanh(x) = -atanh(-x)
2304 ;; x > 1, so use complex version.
2309 ;; atanh(x) = 1/2*log1p(2*x/(1-x))
2315 ;; atanh(x) = 1/2*log1p(2*x + 2*x*x/(1-x))
2323 ;; Stuff which follows is derived from atanh z = (log(1 + z) - log(1 - z))/2
2324 ;; which apparently originates with Kahan's "Much ado" paper.
2326 ;; The formulas for eta and nu below can be easily derived from
2327 ;; rectform(atanh(x+%i*y)) =
2328 ;;
2329 ;; 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)))
2330 ;;
2331 ;; Expand the argument of log out and divide it out and we get
2332 ;;
2333 ;; log(((1+x)^2+y^2)/((1-x)^2+y^2)) = log(1+4*x/((1-x)^2+y^2))
2334 ;;
2335 ;; When y = 0, Im atanh z = 1/2 (arg(1 + x) - arg(1 - x))
2336 ;; = if x < -1 then %pi/2 else if x > 1 then -%pi/2 else <whatever>
2337 ;;
2338 ;; Otherwise, arg(1 - x + %i*(-y)) = - arg(1 - x + %i*y),
2339 ;; and Im atanh z = 1/2 (arg(1 + x + %i*y) + arg(1 - x + %i*y)).
2340 ;; Since arg(x)+arg(y) = arg(x*y) (almost), we can simplify the
2341 ;; imaginary part to
2342 ;;
2343 ;; arg((1+x+%i*y)*(1-x+%i*y)) = arg((1-x)*(1+x)-y^2+2*y*%i)
2344 ;; = atan2(2*y,((1-x)*(1+x)-y^2))
2345 ;;
2346 ;; These are the eta and nu forms below.
2358 ;; Kahan has rho = 4/most-positive-float. What should we do
2359 ;; here about that? Our big floats don't really have a
2360 ;; most-positive float value.
2365 ;; eta = log(1+4*x/((1-x)^2+y^2))/4
2371 ;; If y = 0, then Im atanh z = %pi/2 or -%pi/2.
2372 ;; Otherwise nu = 1/2*atan2(2*y,(1-x)*(1+x)-y^2)
2374 ;; EXTRA FPMINUS HERE TO COUNTERACT FPMINUS IN RETURN VALUE
2385 ;; WTF IS FPMINUS DOING HERE ??
2394 ;; acos(x) for real x. X is a bigfloat, and a maxima number is returned.
2396 ;; acos(x) = %pi/2 - asin(x)
2427 ;; No warning message, please.
2450 ;; x is (mantissa exp), where mantissa = frac*2^fpprec,
2451 ;; with 1/2 < frac <= 1 and x is frac*2^exp. To
2452 ;; compute log(x), use log(x) = log(frac)+ exp*log(2).
2468 ;; Special case for log(1). See Bug 3381301:
2469 ;; https://sourceforge.net/tracker/?func=detail&aid=3381301&group_id=4933&atid=104933
2472 ;; ??? Do we want to return an exact %i*%pi or a float
2473 ;; approximation?