| blob | e3e8875a2e49bfa1c4ea70ae913553ae4c4d1bec |
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 ;; Assume that X has the form ((BIGFLOAT ... <prec>) ...).
162 ;; Return <prec>.
166 ;; Converts the bigfloat L to list of digits including |.| and the
167 ;; exponent marker |b|. The number of significant digits is controlled
168 ;; by $fpprintprec.
178 (let* ((extradigs (floor (1+ (quotient (integer-length (caddr l)) #.(/ (log 10.0) (log 2.0))))))
203 exploded-integer)))
209 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
210 ;; support printing of bfloats.
216 ;; Compute the (decimal exponent) of the bfloat number X.
220 ;; FIXME: do something better than printing and reading
221 ;; the result.
225 ;; Print the bfloat X with FDIGITS after the decimal
226 ;; point. This means, roughtly, FDIGITS+1 significant
227 ;; digits.
230 1
236 ;; Depending on the value of k, move the decimal
237 ;; point. DIGITS was printed assuming the decimal point
238 ;; is after the first digit. But if fdigits = 0, fpformat
239 ;; actually printed out one too many digits, so we need
240 ;; to remove that.
246 ;; Put the leading decimal and then some zeroes
251 ;; The number is scaled by 10^k. Do this by
252 ;; putting the decimal point in the right place,
253 ;; appending zeroes if needed.
257 digits
269 ;; Append some zeroes to get the desired number of digits
275 len
278 1
282 0
288 ;; Exponent overflow
292 ;; The hairy case
296 d)
297 nil))
302 nil)))
312 fstr flen lpoint tpoint)
316 ;; See CLHS 22.3.3.2. "If the parameter d is
317 ;; omitted, ... [and] if the fraction to be
318 ;; printed is zero then a single zero digit should
319 ;; appear after the decimal point." So we need to
320 ;; subtract one from here because we're going to
321 ;; add an extra 0 digit later.
333 w spaceleft overflowchar)
335 ;; Significand overflow; output the overflow char
349 ;; Add a zero if we need it. Which means
350 ;; we figured out we need one above, or
351 ;; another condition. Basically, append a
352 ;; zero if there are no width constraints
353 ;; and if the last char to print was a
354 ;; decimal (so the trailing fraction is
355 ;; zero.)
361 exponentchar
363 stream)
371 ;; NOTE: This is a modified version of FORMAT-FIXED-AUX from CMUCL to
372 ;; support printing of bfloats.
377 ;; Compute the (decimal exponent) of the bfloat number X.
381 ;; FIXME: do something better than printing and reading
382 ;; the result.
386 ;; Print the bfloat X with FDIGITS after the decimal
387 ;; point. To do this we need to know the exponent because
388 ;; fpformat always produces exponential output. If the
389 ;; exponent is E, and we want FDIGITS after the decimal
390 ;; point, we need FDIGITS + E digits printed.
392 #+nil
403 #+nil
408 ;; Depending on the value of scale, move the decimal
409 ;; point. DIGITS was printed assuming the decimal point
410 ;; is after the first digit. But if fdigits = 0, fpformat
411 ;; actually printed out one too many digits, so we need
412 ;; to remove that.
413 #+nil
415 #+nil
418 ;; Figure out where the decimal point should go. An
419 ;; exponent of 0 means the decimal is after the first
420 ;; digit.
426 #+nil
436 #+nil
439 ;; Append some zeroes to get the desired number of digits
445 len
448 1
455 ;;if caller specifically requested no fraction digits, suppress the
456 ;;optional trailing zero
460 ;;optional leading zero
465 ;;optional trailing zero
471 ;;field width overflow
473 t)
482 nil))))))
484 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
485 ;; support printing of bfloats.
491 ;; Compute the (decimal exponent) of the bfloat number X.
495 ;; FIXME: do something better than printing and reading
496 ;; the result.
500 ;; Print the bfloat X with FDIGITS after the decimal
501 ;; point. This means, roughtly, FDIGITS+1 significant
502 ;; digits.
505 1
511 ;; Depending on the value of k, move the decimal
512 ;; point. DIGITS was printed assuming the decimal point
513 ;; is after the first digit. But if fdigits = 0, fpformat
514 ;; actually printed out one too many digits, so we need
515 ;; to remove that.
521 ;; Put the leading decimal and then some zeroes
526 ;; The number is scaled by 10^k. Do this by
527 ;; putting the decimal point in the right place,
528 ;; appending zeroes if needed.
532 digits
544 ;; Append some zeroes to get the desired number of digits
550 len
553 1
557 ;; Default d if omitted. The procedure is taken directly from
558 ;; the definition given in the manual (CLHS 22.3.3.3), and is
559 ;; not very efficient, since we generate the digits twice.
560 ;; Future maintainers are encouraged to improve on this.
561 ;;
562 ;; It's also not very clear whether q in the spec is the
563 ;; number of significant digits or not. I (rtoy) think it
564 ;; makes more sense if q is the number of significant digits.
565 ;; That way 1d300 isn't printed as 1 followed by 300 zeroes.
566 ;; Exponential notation would be used instead.
581 ;; Use dd fraction digits, even if that would cause
582 ;; the width to be exceeded. We choose accuracy over
583 ;; width in this case.
585 ww
586 dd
587 0
588 ovf pad)
589 ovf
594 w
595 orig-d
597 ovf pad exponentchar)))))))
599 ;; Tells you if you have a bigfloat object. BUT, if it is a bigfloat,
600 ;; it will normalize it by making the precision of the bigfloat match
601 ;; the current precision setting in fpprec. And it will also convert
602 ;; bogus zeroes (mantissa is zero, but exponent is not) to a true
603 ;; zero.
605 ;; A bigfloat object looks like '((bigfloat simp <prec>) <mantissa> <exp>)
606 ;; Note bene that the simp flag is optional -- don't count on its presence.
610 ;; Precision matches. (Should we fix up bogus bigfloat
611 ;; zeros?)
614 ;; Current precision is higher than bigfloat precision.
615 ;; Scale up mantissa and adjust exponent to get the
616 ;; correct precision.
620 ;; Current precision is LOWER than bigfloat precision.
621 ;; Round the number to the desired precision.
624 ;; Fix up any bogus zeros that we might have created.
652 exp))
657 $ratepsilon)))))
674 fprateps)))
691 ;; BLECHH, I HATE ENUMERATING CASES. IS THERE A BETTER WAY ??
697 ;; NOT SURE THE FOLLOWING REALLY WORKS
698 ;; #+(and cmu double-double)
699 ;; (kernel:double-double-float
700 ;; (values most-negative-double-double-float most-positive-double-double-float))
701 ))
703 ;; Convert a floating point number into a bigfloat.
714 ;; Need to check for zero because different lisps return different
715 ;; values for integer-decode-float of a 0. In particular CMUCL
716 ;; returns 0, -1075. A bigfloat zero needs to have an exponent and
717 ;; mantissa of zero.
722 ;; Scale frac to the desired number of bits, and adjust the
723 ;; exponent accordingly.
728 ;; Convert a bigfloat into a floating point number.
735 ;; Round the mantissa to the number of bits of precision of the
736 ;; machine, and then convert it to a floating point fraction. We
737 ;; have 0.5 <= mantissa < 1
739 ;; Multiply the mantissa by the exponent portion. I'm not sure
740 ;; why the exponent computation is so complicated.
741 ;;
742 ;; GCL doesn't signal overflow from scale-float if the number
743 ;; would overflow. We have to do it this way. 0.5 <= mantissa <
744 ;; 1. The largest double-float is .999999 * 2^1024. So if the
745 ;; exponent is 1025 or higher, we have an overflow.
751 ;; New machine-independent version of FIXFLOAT. This may be buggy. - CWH
752 ;; It is buggy! On the PDP10 it dies on (RATIONALIZE -1.16066076E-7)
753 ;; which calls FLOAT on some rather big numbers. ($RATEPSILON is approx.
754 ;; 7.45E-9) - JPG
760 ;; Takes a flonum arg and returns a rational number corresponding to the flonum
761 ;; in the form of a dotted pair of two integers. Since the denominator will
762 ;; always be a positive power of 2, this number will not always be in lowest
763 ;; terms.
779 x))
828 ;; R is a Maxima ratio, represented as a list of the numerator and
829 ;; denominator. FLOAT-RATIO doesn't like it if the numerator is 0,
830 ;; so handle that here.
835 ;; This is borrowed from CMUCL (float-ratio-float), and modified for
836 ;; converting ratios to Maxima's bfloat numbers.
845 ;;
846 ;; Strip any trailing zeros from the denominator and move it into the scale
847 ;; factor (to minimize the size of the operands.)
852 ;;
853 ;; Guess how much we need to scale by from the magnitudes of the numerator
854 ;; and denominator. We want one extra bit for a guard bit.
890 fraction-and-guard
914 ;; Bigfloat atan
919 ;; |x| > 1.
920 ;;
921 ;; Use A&S 4.4.5:
922 ;; atan(x) + acot(x) = +/- pi/2 (+ for x >= 0, - for x < 0)
923 ;;
924 ;; and A&S 4.4.8
925 ;; acot(z) = atan(1/z)
932 ;; |x| > 1/2
933 ;;
934 ;; Use A&S 4.4.42, third formula:
935 ;;
936 ;; atan(z) = z/(1+z^2)*[1 + 2/3*r + (2*4)/(3*5)*r^2 + ...]
937 ;;
938 ;; r = z^2/(1+z^2)
950 ;; |x| <= 1/2. Use Taylor series (A&S 4.4.42, first
951 ;; formula).
960 ;; atan(y/x) taking into account the quadrant. (Also equal to
961 ;; arg(x+%i*y).)
964 ;; atan(y/0) = atan(inf), but what sign?
968 ;; We're on the negative imaginary axis, so -pi/2.
971 ;; The positive imaginary axis, so +pi/2
974 ;; x > 0. atan(y/x) is the correct value.
977 ;; x < 0, and y > 0. We're in quadrant II, so the angle we
978 ;; want is pi+atan(y/x).
981 ;; x <= 0 and y <= 0. We're in quadrant III, so the angle we
982 ;; want is atan(y/x)-pi.
993 ;; Returns a list of a mantissa and an exponent.
1003 ;; It seems to me that this function gets called on an integer
1004 ;; and returns the mantissa portion of the mantissa/exponent pair.
1006 ;; "STICKY BIT" CALCULATION FIXED 10/14/75 --RJF
1007 ;; BASE must not get temporarily bound to NIL by being placed
1008 ;; in a PROG list as this will confuse stepping programs.
1014 ;;*M will be positive if the precision of the argument is greater than
1015 ;;the current precision being used.
1020 ;;FPSHIFT is essentially LSH.
1027 ;ONLY ZEROES SHIFTED OFF
1039 ;; Compute (* L (expt d n)) where D is 2 or 10 depending on
1040 ;; *decfp. Throw away an fractional part by truncating to zero.
1044 ;; Left shift of negative number requires some
1045 ;; care. (That is, (truncate l (expt 2 n)), but use
1046 ;; shifts instead.)
1056 ;; Bignum LSH -- N is assumed (and declared above) to be a fixnum.
1057 ;; This isn't really LSH, since the sign bit isn't propagated when
1058 ;; shifting to the right, i.e. (BIGLSH -100 -3) = -40, whereas
1059 ;; (LSH -100 -3) = 777777777770 (on a 36 bit machine).
1060 ;; This actually computes (* X (EXPT 2 N)). As of 12/21/80, this function
1061 ;; was only called by FPSHIFT. I would like to hear an argument as why this
1062 ;; is more efficient than simply writing (* X (EXPT 2 N)). Is the
1063 ;; intermediate result created by (EXPT 2 N) the problem? I assume that
1064 ;; EXPT tries to LSH when possible.
1070 ;; Either we are shifting a fixnum to the right, or shifting
1071 ;; a fixnum to the left, but not far enough left for it to become
1072 ;; a bignum.
1076 ;; The form which follows is nearly identical to (ASH X N), however
1077 ;; (ASH -100 -20) = -1, whereas (BIGLSH -100 -20) = 0.
1081 ;; If we get here, then either X is a bignum or our answer is
1082 ;; going to be a bignum.
1089 ;; Isn't this the kind of optimization that compilers are
1090 ;; supposed to make?
1095 ;; exp(x)
1096 ;;
1097 ;; For negative x, use exp(-x) = 1/exp(x)
1098 ;;
1099 ;; For x > 0, exp(x) = exp(r+y) = exp(r) * exp(y), where x = r + y and
1100 ;; r = floor(x).
1116 ;; exp(x) for small x, using Taylor series.
1127 ;; Does one higher precision to round correctly.
1128 ;; A and B are each a list of a mantissa and an exponent.
1148 min))
1153 max))
1155 ;; The following functions compute bigfloat values for %e, %pi,
1156 ;; %gamma, and log(2). For each precision, the computed value is
1157 ;; cached in a hash table so it doesn't need to be computed again.
1158 ;; There are functions to return the hash table or clear the hash
1159 ;; table, for debugging.
1160 ;;
1161 ;; Note that each of these return a bigfloat number, but without the
1162 ;; bigfloat tag.
1163 ;;
1164 ;; See
1165 ;; https://sourceforge.net/p/maxima/bugs/1842/
1166 ;; for an explanation.
1171 value
1174 table)
1182 value
1185 table)
1193 value
1196 table)
1204 value
1207 table)
1211 ;; This doesn't need a hash table because there's never a problem with
1212 ;; using a high precision value and rounding to a lower precision
1213 ;; value because 1 is always an exact bfloat.
1219 ;;----------------------------------------------------------------------------;;
1220 ;;
1221 ;; The values of %e, %pi, %gamma and log(2) are computed by the technique of
1222 ;; binary splitting. See http://www.ginac.de/CLN/binsplit.pdf for details.
1223 ;;
1224 ;; Volker van Nek, Sept. 2014
1226 ;;
1227 ;; Euler's number E
1228 ;;
1232 ;;
1233 ;; Taylor: %e = sum(s[i] ,i,0,inf) where s[i] = 1/i!
1234 ;;
1239 ;;
1240 ;; binary splitting:
1241 ;;
1242 ;; 1
1243 ;; s[i] = ----------------------
1244 ;; q[0]*q[1]*q[2]*..*q[i]
1245 ;;
1246 ;; where q[0] = 1
1247 ;; q[i] = i
1248 ;;
1260 ;;
1261 ;; stop when i! > 2^fpprec
1262 ;;
1263 ;; log(i!) = sum(log(k), k,1,i) > fpprec * log(2)
1264 ;;
1271 ;;
1272 ;;----------------------------------------------------------------------------;;
1273 ;;
1274 ;; PI
1275 ;;
1279 ;;
1280 ;; Chudnovsky & Chudnovsky:
1281 ;;
1282 ;; C^(3/2)/(12*%pi) = sum(s[i], i,0,inf),
1283 ;;
1284 ;; where s[i] = (-1)^i*(6*i)!*(A*i+B) / (i!^3*(3*i)!*C^(3*i))
1285 ;;
1286 ;; and A = 545140134, B = 13591409, C = 640320
1287 ;;
1291 ;; STEP 1:
1292 ;; compute n/d = sqrt(10005) :
1293 ;;
1294 ;; n[0] n[i+1] = n[i]^2+a*d[i]^2 n[inf]
1295 ;; quadratic Heron: x[0] = ----, , sqrt(a) = ------
1296 ;; d[0] d[i+1] = 2*n[i]*d[i] d[inf]
1297 ;;
1303 ;; STEP 2:
1304 ;; divide C^(3/2)/12 = 3335*2^7*sqrt(10005)
1305 ;; by Chudnovsky-sum = tt/qq :
1306 ;;
1308 ;; fpprec*log(2)/log(C^3/(24*6*2*6))
1314 ;;
1315 ;; The returned n and d serve as start values for the iteration.
1316 ;; n/d = sqrt(10005) with a precision of p = ceiling(prec/2^nr) bits
1317 ;; where nr is the number of needed iterations.
1318 ;;
1333 ;;
1334 ;; binary splitting:
1335 ;;
1336 ;; a[i] * p[0]*p[1]*p[2]*..*p[i]
1337 ;; s[i] = -----------------------------
1338 ;; q[0]*q[1]*q[2]*..*q[i]
1339 ;;
1340 ;; where a[0] = B
1341 ;; p[0] = q[0] = 1
1342 ;; a[i] = A*i+B
1343 ;; p[i] = - (6*i-5)*(2*i-1)*(6*i-1)
1344 ;; q[i] = C^3/24*i^3
1345 ;;
1364 ;;
1365 ;;----------------------------------------------------------------------------;;
1366 ;;
1367 ;; Euler-Mascheroni constant GAMMA
1368 ;;
1372 ;;
1373 ;; Brent-McMillan algorithm
1374 ;;
1375 ;; Let
1376 ;; alpha = 4.970625759544
1377 ;;
1378 ;; n > 0 and N-1 >= alpha*n
1379 ;;
1380 ;; H(k) = sum(1/i, i,1,k)
1381 ;;
1382 ;; S = sum(H(k)*(n^k/k!)^2, k,0,N-1)
1383 ;;
1384 ;; I = sum((n^k/k!)^2, k,0,N-1)
1385 ;;
1386 ;; T = 1/(4*n)*sum((2*k)!^3/(k!^4*(16*n)^(2*k)), k,0,2*n-1)
1387 ;;
1388 ;; and
1389 ;; %gamma = S/I - T/I^2 - log(n)
1390 ;;
1391 ;; Then
1392 ;; |%gamma - gamma| < 24 * e^(-8*n)
1393 ;;
1394 ;; (Corollary 2, Remark 2, Brent/Johansson http://arxiv.org/pdf/1312.0039v1.pdf)
1395 ;;
1403 sumi sumi2 ;; I and I^2
1406 ;;
1407 ;; sums = vv/(qq*dd)
1408 ;; sumi = tt/qq
1409 ;; sums/sumi = vv/(qq*dd)*qq/tt = vv/(dd*tt)
1410 ;;
1414 ;;
1416 ;;
1417 ;; sumt = 1/(4*n)*ttt/qqq
1418 ;; sumt/sumi2 = ttt/(4*n*qqq*sumi2)
1419 ;;
1422 ;; %gamma :
1424 ;;
1425 ;; split S and I simultaneously:
1426 ;;
1427 ;; summands I[0] = 1, I[i]/I[i-1] = n^2/i^2
1428 ;;
1429 ;; S[0] = 0, S[i]/S[i-1] = n^2/i^2*H(i)/H(i-1)
1430 ;;
1431 ;; p[0]*p[1]*p[2]*..*p[i]
1432 ;; I[i] = ----------------------
1433 ;; q[0]*q[1]*q[2]*..*q[i]
1434 ;;
1435 ;; where p[0] = n^2
1436 ;; q[0] = 1
1437 ;; p[i] = n^2
1438 ;; q[i] = i^2
1439 ;; c[0] c[1] c[2] c[i]
1440 ;; S[i] = H[i] * I[i], where H[i] = ---- + ---- + ---- + .. + ----
1441 ;; d[0] d[1] d[2] d[i]
1442 ;; and c[0] = 0
1443 ;; d[0] = 1
1444 ;; c[i] = 1
1445 ;; d[i] = i
1446 ;;
1466 ;;
1467 ;; split 4*n*T:
1468 ;;
1469 ;; summands T[0] = 1, T[i]/T[i-1] = (2*i-1)^3/(32*i*n^2)
1470 ;;
1471 ;; p[0]*p[1]*p[2]*..*p[i]
1472 ;; T[i] = ----------------------
1473 ;; q[0]*q[1]*q[2]*..*q[i]
1474 ;;
1475 ;; where p[0] = q[0] = 1
1476 ;; p[i] = (2*i-1)^3
1477 ;; q[i] = 32*i*n^2
1478 ;;
1494 ;;
1495 ;;----------------------------------------------------------------------------;;
1496 ;;
1497 ;; log(2) = 18*L(26) - 2*L(4801) + 8*L(8749)
1498 ;;
1499 ;; where L(k) = atanh(1/k)
1500 ;;
1501 ;; see http://numbers.computation.free.fr/Constants/constants.html
1502 ;;
1503 ;;;(defun $log2 () (bcons (comp-log2))) ;; checked against reference table
1504 ;;
1512 ;;
1513 ;; Taylor: atanh(1/k) = sum(s[i], i,0,inf)
1514 ;;
1515 ;; where s[i] = 1/((2*i+1)*k^(2*i+1))
1516 ;;
1522 ;;
1523 ;; binary splitting:
1524 ;; 1
1525 ;; s[i] = -----------------------------
1526 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1527 ;;
1528 ;; where b[0] = 1
1529 ;; q[0] = k
1530 ;; b[i] = 2*i+1
1531 ;; q[i] = k^2
1532 ;;
1546 ;;
1547 ;;----------------------------------------------------------------------------;;
1548 ;;
1549 ;; log(n) = log(n/2^k) + k*log(2)
1550 ;;
1551 ;;;(defun $log10 () (bcons (log-n 10))) ;; checked against reference table
1552 ;;
1561 ;; choose k so that |n/2^k - 1| is as small as possible:
1563 ;; now |n/2^k - 1| <= 1/3
1567 ;;
1568 ;; log(1+u/v) = 2 * sum(s[i], i,0,inf)
1569 ;;
1570 ;; where s[i] = (u/(2*v+u))^(2*i+1)/(2*i+1)
1571 ;;
1586 ;;
1587 ;; binary splitting:
1588 ;;
1589 ;; p[0]*p[1]*p[2]*..*p[i]
1590 ;; s[i] = -----------------------------
1591 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1592 ;;
1593 ;; where b[0] = 1
1594 ;; p[0] = u
1595 ;; q[0] = w = 2*v+u
1596 ;; b[i] = 2*i+1
1597 ;; p[i] = u^2
1598 ;; q[i] = w^2
1599 ;;
1614 ;;
1615 ;;----------------------------------------------------------------------------;;
1616 \f
1623 x
1640 ; COMPUTE STICKY BIT
1642 ; MAKE ROOM FOR GUARD DIGIT & STICKY BIT
1662 ;; Don't use the symbol BASE since it is SPECIAL.
1670 bas)))
1672 ;; NN is positive or negative integer
1711 $%i)))))
1726 ;; The number of extra bits required
1732 ;; Special case for a = 0b0. General algorithm loops endlessly in that case.
1734 ;; Unlike many or maybe all of the other functions named FP-something,
1735 ;; FPROOT assumes it is called with an argument like
1736 ;; '((BIGFLOAT ...) FOO BAR) instead of '(FOO BAR).
1737 ;; However FPROOT does return something like '(FOO BAR).
1764 fans
1768 ;; If A is a bigfloat, be sure to round it to the current precision.
1769 ;; (See Bug 2543079 for one of the symptoms.)
1795 ;; FPPREC 16, to get rid of the miniscule imaginary
1796 ;; part of the a+bi answer.
1803 e))
1811 ;; THIS VERSION OF FPSIN COMPUTES SIN OR COS TO PRECISION FPPREC,
1812 ;; BUT CHECKS FOR THE POSSIBILITY OF CATASTROPHIC CANCELLATION DURING
1813 ;; ARGUMENT REDUCTION (E.G. SIN(N*%PI+EPSILON))
1814 ;; *FPSINCHECK* WILL CAUSE PRINTOUT OF ADDITIONAL INFO WHEN
1815 ;; EXTRA PRECISION IS NEEDED FOR SIN/COS CALCULATION. KNOWN
1816 ;; BAD FEATURES: IT IS NOT NECESSARY TO USE EXTRA PRECISION FOR, E.G.
1817 ;; SIN(PI/2), WHICH IS NOT NEAR ZERO, BUT EXTRA
1818 ;; PRECISION IS USED SINCE IT IS NEEDED FOR COS(PI/2).
1819 ;; PRECISION SEEMS TO BE 100% SATSIFACTORY FOR LARGE ARGUMENTS, E.G.
1820 ;; SIN(31415926.0B0), BUT LESS SO FOR SIN(3.1415926B0). EXPLANATION
1821 ;; NOT KNOWN. (9/12/75 RJF)
1825 ;; FL is a T for sin and NIL for cos.
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
1907 ;; information to calculate the integer part. Specifically, we only
1908 ;; have enough information if EXPONENT < FPPREC. If that isn't the
1909 ;; case, we signal a Maxima error. However, if SKIP-EXPONENT-CHECK-P
1910 ;; is non-NIL, this check is skipped, and we compute the integer part
1911 ;; as requested.
1912 ;;
1913 ;; For the bigfloat code here, skip-exponent-check-p should be true.
1914 ;; For other uses (see commit 576c7508 and bug #2784), this should be
1915 ;; nil, which is the default.
1918 f
1934 ;;; Computes the log of a bigfloat number.
1935 ;;;
1936 ;;; Uses the series
1937 ;;;
1938 ;;; log(1+x) = sum((x/(x+2))^(2*n+1)/(2*n+1),n,0,inf);
1939 ;;;
1940 ;;;
1941 ;;; INF x 2 n + 1
1942 ;;; ==== (-----)
1943 ;;; \ x + 2
1944 ;;; = 2 > --------------
1945 ;;; / 2 n + 1
1946 ;;; ====
1947 ;;; n = 0
1948 ;;;
1949 ;;;
1950 ;;; which converges for x > 0.
1951 ;;;
1952 ;;; Note that FPLOG is given 1+X, not X.
1953 ;;;
1954 ;;; However, to aid convergence of the series, we scale 1+x until 1/e
1955 ;;; < 1+x <= e.
1956 ;;;
1964 ;; Scale X until 1/e < X <= E. ANS keeps track of how
1965 ;; many factors of E were used. Set X to NIL if X is E.
1978 ;; Prepare X for the series. The series is for 1 + x, so
1979 ;; get x from our X. TERM is (x/(x+2)). X becomes
1980 ;; (x/(x+2))^2.
1984 ;; Sum the series until the sum (in ANS) doesn't change
1985 ;; anymore.
2001 \f
2002 ;;;; Bigfloat implementations of special functions.
2003 ;;;;
2005 ;;; This is still a bit messy. Some functions here take bigfloat
2006 ;;; numbers, represented by ((bigfloat) <mant> <exp>), but others want
2007 ;;; just the FP number, represented by (<mant> <exp>). Likewise, some
2008 ;;; return a bigfloat, some return just the FP.
2009 ;;;
2010 ;;; This needs to be systemized somehow. It isn't helped by the fact
2011 ;;; that some of the routines above also do the samething.
2012 ;;;
2013 ;;; The implementation for the special functions for a complex
2014 ;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex
2015 ;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in
2016 ;;; Iserles and Powell (eds.) "The State of the Art in Numerical
2017 ;;; Analysis", pp 165-211, Clarendon Press, 1987
2019 ;; Compute exp(x) - 1, but do it carefully to preserve precision when
2020 ;; |x| is small. X is a FP number, and a FP number is returned. That
2021 ;; is, no bigfloat stuff.
2023 ;; What is the right breakpoint here? Is 1 ok? Perhaps 1/e is better?
2025 ;; exp(x) - 1
2028 ;; Use Taylor series for exp(x) - 1
2037 ans))))
2039 ;; log(1+x) for small x. X is FP number, and a FP number is returned.
2041 ;; Use the same series as given above for fplog. For small x we use
2042 ;; the series, otherwise fplog is accurate enough.
2057 ;; sinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2059 ;; X must be a maxima bigfloat
2061 ;; See, for example, Hart et al., Computer Approximations, 6.2.27:
2062 ;;
2063 ;; sinh(x) = 1/2*(D(x) + D(x)/(1+D(x)))
2064 ;;
2065 ;; where D(x) = exp(x) - 1.
2066 ;;
2067 ;; But for negative x, use sinh(x) = -sinh(-x) because D(x)
2068 ;; approaches -1 for large negative x.
2070 ;; Special case: x=0. Return immediately.
2073 ;; x is positive.
2078 ;; x is negative.
2083 ;; The rectform for sinh for complex args should be numerically
2084 ;; accurate, so return nil in that case.
2088 ;; asinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2090 ;; asinh(x) = sign(x) * log(|x| + sqrt(1+x*x))
2091 ;;
2092 ;; And
2093 ;;
2094 ;; asinh(x) = x, if 1+x*x = 1
2095 ;; = sign(x) * (log(2) + log(x)), large |x|
2096 ;; = sign(x) * log(2*|x| + 1/(|x|+sqrt(1+x*x))), if |x| > 2
2097 ;; = sign(x) * log1p(|x|+x^2/(1+sqrt(1+x*x))), otherwise.
2098 ;;
2099 ;; But I'm lazy right now and we only implement the last 2 cases.
2100 ;; We should implement all cases.
2106 result)
2107 ;; We only use two formulas here. |x| <= 2 and |x| > 2. Should
2108 ;; we add one for very big x and one for very small x, as given above.
2110 ;; |x| > 2
2111 ;;
2112 ;; log(2*|x| + 1/(|x|+sqrt(1+x^2)))
2120 ;; |x| <= 2
2121 ;;
2122 ;; log1p(|x|+x^2/(1+sqrt(1+x^2)))
2134 ;; asinh(z) = -%i * asin(%i*z)
2147 ;; asin(x) = atan(x/(sqrt(1-x^2))
2148 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2149 ;;
2150 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2151 ;;
2152 ;; If |x| > 1, we need to do something else.
2153 ;;
2154 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2155 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2156 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2157 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2161 ;; asin(-x) = -asin(x);
2164 ;; 0 <= x < 1/2
2165 ;; asin(x) = atan(x/sqrt(1-x^2))
2173 ;; x > 1
2174 ;; asin(x) = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2175 ;;
2176 ;; Should we try to do something a little fancier with the
2177 ;; argument to log and use log1p for better accuracy?
2186 ;; 1/2 <= x <= 1
2187 ;; asin(x) = %pi/2 - atan(sqrt(1-x^2)/x)
2195 fp-x)))))))))
2197 ;; asin(x) for real x. X is a bigfloat, and a maxima number (real or
2198 ;; complex) is returned.
2200 ;; asin(x) = atan(x/(sqrt(1-x^2))
2201 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2202 ;;
2203 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2204 ;;
2205 ;; If |x| > 1, we need to do something else.
2206 ;;
2207 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2208 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2209 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2210 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2214 ;; Square root of a complex number (xx, yy). Both are bigfloats. FP
2215 ;; (non-bigfloat) numbers are returned.
2233 rho))))
2236 ;; asin(z) for complex z = x + %i*y. X and Y are bigfloats. The real
2237 ;; and imaginary parts are returned as bigfloat numbers.
2247 ;; Realpart is atan(x/Re(sqrt(1-z)*sqrt(1+z)))
2248 ;; Imagpart is asinh(Im(conj(sqrt(1-z))*sqrt(1+z)))
2254 ;; Check for division by zero. If we would divide
2255 ;; by zero, return pi/2 or -pi/2 according to the
2256 ;; sign of X.
2276 ;; tanh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2278 ;; X is Maxima bigfloat
2279 ;; tanh(x) = D(2*x)/(2+D(2*x))
2285 ;; tanh(z), z = x + %i*y. X, Y are bigfloats, and a maxima number is
2286 ;; returned.
2303 ;; atanh(x) for real x, |x| <= 1. X is a bigfloat, and a bigfloat is
2304 ;; returned.
2306 ;; atanh(x) = -atanh(-x)
2307 ;; = 1/2*log1p(2*x/(1-x)), x >= 0.5
2308 ;; = 1/2*log1p(2*x+2*x*x/(1-x)), x <= 0.5
2312 ;; atanh(x) = -atanh(-x)
2315 ;; x > 1, so use complex version.
2320 ;; atanh(x) = 1/2*log1p(2*x/(1-x))
2326 ;; atanh(x) = 1/2*log1p(2*x + 2*x*x/(1-x))
2334 ;; Stuff which follows is derived from atanh z = (log(1 + z) - log(1 - z))/2
2335 ;; which apparently originates with Kahan's "Much ado" paper.
2337 ;; The formulas for eta and nu below can be easily derived from
2338 ;; rectform(atanh(x+%i*y)) =
2339 ;;
2340 ;; 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)))
2341 ;;
2342 ;; Expand the argument of log out and divide it out and we get
2343 ;;
2344 ;; log(((1+x)^2+y^2)/((1-x)^2+y^2)) = log(1+4*x/((1-x)^2+y^2))
2345 ;;
2346 ;; When y = 0, Im atanh z = 1/2 (arg(1 + x) - arg(1 - x))
2347 ;; = if x < -1 then %pi/2 else if x > 1 then -%pi/2 else <whatever>
2348 ;;
2349 ;; Otherwise, arg(1 - x + %i*(-y)) = - arg(1 - x + %i*y),
2350 ;; and Im atanh z = 1/2 (arg(1 + x + %i*y) + arg(1 - x + %i*y)).
2351 ;; Since arg(x)+arg(y) = arg(x*y) (almost), we can simplify the
2352 ;; imaginary part to
2353 ;;
2354 ;; arg((1+x+%i*y)*(1-x+%i*y)) = arg((1-x)*(1+x)-y^2+2*y*%i)
2355 ;; = atan2(2*y,((1-x)*(1+x)-y^2))
2356 ;;
2357 ;; These are the eta and nu forms below.
2369 ;; Kahan has rho = 4/most-positive-float. What should we do
2370 ;; here about that? Our big floats don't really have a
2371 ;; most-positive float value.
2376 ;; eta = log(1+4*x/((1-x)^2+y^2))/4
2382 ;; If y = 0, then Im atanh z = %pi/2 or -%pi/2.
2383 ;; Otherwise nu = 1/2*atan2(2*y,(1-x)*(1+x)-y^2)
2385 ;; EXTRA FPMINUS HERE TO COUNTERACT FPMINUS IN RETURN VALUE
2396 ;; WTF IS FPMINUS DOING HERE ??
2405 ;; acos(x) for real x. X is a bigfloat, and a maxima number is returned.
2407 ;; acos(x) = %pi/2 - asin(x)
2438 ;; No warning message, please.
2461 ;; x is (mantissa exp), where mantissa = frac*2^fpprec,
2462 ;; with 1/2 < frac <= 1 and x is frac*2^exp. To
2463 ;; compute log(x), use log(x) = log(frac)+ exp*log(2).
2479 ;; Special case for log(1). See:
2480 ;; https://sourceforge.net/p/maxima/bugs/2240/
2483 ;; ??? Do we want to return an exact %i*%pi or a float
2484 ;; approximation?