| blob | c236f4dc3b54ee9e80afbc21ea3e0836029123c3 |
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))))))
201 ;; If the rounded integer has more digits than
202 ;; expected, we need to adjust the exponent by
203 ;; this amount. This also means we need to remove
204 ;; these extra digits so that the result has the
205 ;; desired number of digits.
212 exploded-integer)))
218 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
219 ;; support printing of bfloats.
225 ;; Compute the (decimal exponent) of the bfloat number X.
229 ;; FIXME: do something better than printing and reading
230 ;; the result.
234 ;; Print the bfloat X with FDIGITS after the decimal
235 ;; point. This means, roughtly, FDIGITS+1 significant
236 ;; digits.
239 1
245 ;; Depending on the value of k, move the decimal
246 ;; point. DIGITS was printed assuming the decimal point
247 ;; is after the first digit. But if fdigits = 0, fpformat
248 ;; actually printed out one too many digits, so we need
249 ;; to remove that.
255 ;; Put the leading decimal and then some zeroes
260 ;; The number is scaled by 10^k. Do this by
261 ;; putting the decimal point in the right place,
262 ;; appending zeroes if needed.
266 digits
278 ;; Append some zeroes to get the desired number of digits
284 len
287 1
291 0
297 ;; Exponent overflow
301 ;; The hairy case
305 d)
306 nil))
311 nil)))
321 fstr flen lpoint tpoint)
325 ;; See CLHS 22.3.3.2. "If the parameter d is
326 ;; omitted, ... [and] if the fraction to be
327 ;; printed is zero then a single zero digit should
328 ;; appear after the decimal point." So we need to
329 ;; subtract one from here because we're going to
330 ;; add an extra 0 digit later.
342 w spaceleft overflowchar)
344 ;; Significand overflow; output the overflow char
358 ;; Add a zero if we need it. Which means
359 ;; we figured out we need one above, or
360 ;; another condition. Basically, append a
361 ;; zero if there are no width constraints
362 ;; and if the last char to print was a
363 ;; decimal (so the trailing fraction is
364 ;; zero.)
370 exponentchar
372 stream)
380 ;; NOTE: This is a modified version of FORMAT-FIXED-AUX from CMUCL to
381 ;; support printing of bfloats.
386 ;; Compute the (decimal exponent) of the bfloat number X.
390 ;; FIXME: do something better than printing and reading
391 ;; the result.
395 ;; Print the bfloat X with FDIGITS after the decimal
396 ;; point. To do this we need to know the exponent because
397 ;; fpformat always produces exponential output. If the
398 ;; exponent is E, and we want FDIGITS after the decimal
399 ;; point, we need FDIGITS + E digits printed.
401 #+nil
412 #+nil
417 ;; Depending on the value of scale, move the decimal
418 ;; point. DIGITS was printed assuming the decimal point
419 ;; is after the first digit. But if fdigits = 0, fpformat
420 ;; actually printed out one too many digits, so we need
421 ;; to remove that.
422 #+nil
424 #+nil
427 ;; Figure out where the decimal point should go. An
428 ;; exponent of 0 means the decimal is after the first
429 ;; digit.
435 #+nil
445 #+nil
448 ;; Append some zeroes to get the desired number of digits
454 len
457 1
464 ;;if caller specifically requested no fraction digits, suppress the
465 ;;optional trailing zero
469 ;;optional leading zero
474 ;;optional trailing zero
480 ;;field width overflow
482 t)
491 nil))))))
493 ;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to
494 ;; support printing of bfloats.
500 ;; Compute the (decimal exponent) of the bfloat number X.
504 ;; FIXME: do something better than printing and reading
505 ;; the result.
509 ;; Print the bfloat X with FDIGITS after the decimal
510 ;; point. This means, roughtly, FDIGITS+1 significant
511 ;; digits.
514 1
520 ;; Depending on the value of k, move the decimal
521 ;; point. DIGITS was printed assuming the decimal point
522 ;; is after the first digit. But if fdigits = 0, fpformat
523 ;; actually printed out one too many digits, so we need
524 ;; to remove that.
530 ;; Put the leading decimal and then some zeroes
535 ;; The number is scaled by 10^k. Do this by
536 ;; putting the decimal point in the right place,
537 ;; appending zeroes if needed.
541 digits
553 ;; Append some zeroes to get the desired number of digits
559 len
562 1
566 ;; Default d if omitted. The procedure is taken directly from
567 ;; the definition given in the manual (CLHS 22.3.3.3), and is
568 ;; not very efficient, since we generate the digits twice.
569 ;; Future maintainers are encouraged to improve on this.
570 ;;
571 ;; It's also not very clear whether q in the spec is the
572 ;; number of significant digits or not. I (rtoy) think it
573 ;; makes more sense if q is the number of significant digits.
574 ;; That way 1d300 isn't printed as 1 followed by 300 zeroes.
575 ;; Exponential notation would be used instead.
590 ;; Use dd fraction digits, even if that would cause
591 ;; the width to be exceeded. We choose accuracy over
592 ;; width in this case.
594 ww
595 dd
596 0
597 ovf pad)
598 ovf
603 w
604 orig-d
606 ovf pad exponentchar)))))))
608 ;; Tells you if you have a bigfloat object. BUT, if it is a bigfloat,
609 ;; it will normalize it by making the precision of the bigfloat match
610 ;; the current precision setting in fpprec. And it will also convert
611 ;; bogus zeroes (mantissa is zero, but exponent is not) to a true
612 ;; zero.
614 ;; A bigfloat object looks like '((bigfloat simp <prec>) <mantissa> <exp>)
615 ;; Note bene that the simp flag is optional -- don't count on its presence.
619 ;; Precision matches. (Should we fix up bogus bigfloat
620 ;; zeros?)
623 ;; Current precision is higher than bigfloat precision.
624 ;; Scale up mantissa and adjust exponent to get the
625 ;; correct precision.
629 ;; Current precision is LOWER than bigfloat precision.
630 ;; Round the number to the desired precision.
633 ;; Fix up any bogus zeros that we might have created.
661 exp))
666 $ratepsilon)))))
683 fprateps)))
700 ;; BLECHH, I HATE ENUMERATING CASES. IS THERE A BETTER WAY ??
706 ;; NOT SURE THE FOLLOWING REALLY WORKS
707 ;; #+(and cmu double-double)
708 ;; (kernel:double-double-float
709 ;; (values most-negative-double-double-float most-positive-double-double-float))
710 ))
712 ;; Convert a floating point number into a bigfloat.
723 ;; Need to check for zero because different lisps return different
724 ;; values for integer-decode-float of a 0. In particular CMUCL
725 ;; returns 0, -1075. A bigfloat zero needs to have an exponent and
726 ;; mantissa of zero.
731 ;; Scale frac to the desired number of bits, and adjust the
732 ;; exponent accordingly.
737 ;; Convert a bigfloat into a floating point number.
744 ;; Round the mantissa to the number of bits of precision of the
745 ;; machine, and then convert it to a floating point fraction. We
746 ;; have 0.5 <= mantissa < 1
748 ;; Multiply the mantissa by the exponent portion. I'm not sure
749 ;; why the exponent computation is so complicated.
750 ;;
751 ;; GCL doesn't signal overflow from scale-float if the number
752 ;; would overflow. We have to do it this way. 0.5 <= mantissa <
753 ;; 1. The largest double-float is .999999 * 2^1024. So if the
754 ;; exponent is 1025 or higher, we have an overflow.
760 ;; New machine-independent version of FIXFLOAT. This may be buggy. - CWH
761 ;; It is buggy! On the PDP10 it dies on (RATIONALIZE -1.16066076E-7)
762 ;; which calls FLOAT on some rather big numbers. ($RATEPSILON is approx.
763 ;; 7.45E-9) - JPG
769 ;; Takes a flonum arg and returns a rational number corresponding to the flonum
770 ;; in the form of a dotted pair of two integers. Since the denominator will
771 ;; always be a positive power of 2, this number will not always be in lowest
772 ;; terms.
788 x))
837 ;; R is a Maxima ratio, represented as a list of the numerator and
838 ;; denominator. FLOAT-RATIO doesn't like it if the numerator is 0,
839 ;; so handle that here.
844 ;; This is borrowed from CMUCL (float-ratio-float), and modified for
845 ;; converting ratios to Maxima's bfloat numbers.
854 ;;
855 ;; Strip any trailing zeros from the denominator and move it into the scale
856 ;; factor (to minimize the size of the operands.)
861 ;;
862 ;; Guess how much we need to scale by from the magnitudes of the numerator
863 ;; and denominator. We want one extra bit for a guard bit.
899 fraction-and-guard
923 ;; Bigfloat atan
928 ;; |x| > 1.
929 ;;
930 ;; Use A&S 4.4.5:
931 ;; atan(x) + acot(x) = +/- pi/2 (+ for x >= 0, - for x < 0)
932 ;;
933 ;; and A&S 4.4.8
934 ;; acot(z) = atan(1/z)
941 ;; |x| > 1/2
942 ;;
943 ;; Use A&S 4.4.42, third formula:
944 ;;
945 ;; atan(z) = z/(1+z^2)*[1 + 2/3*r + (2*4)/(3*5)*r^2 + ...]
946 ;;
947 ;; r = z^2/(1+z^2)
959 ;; |x| <= 1/2. Use Taylor series (A&S 4.4.42, first
960 ;; formula).
969 ;; atan(y/x) taking into account the quadrant. (Also equal to
970 ;; arg(x+%i*y).)
973 ;; atan(y/0) = atan(inf), but what sign?
977 ;; We're on the negative imaginary axis, so -pi/2.
980 ;; The positive imaginary axis, so +pi/2
983 ;; x > 0. atan(y/x) is the correct value.
986 ;; x < 0, and y > 0. We're in quadrant II, so the angle we
987 ;; want is pi+atan(y/x).
990 ;; x <= 0 and y <= 0. We're in quadrant III, so the angle we
991 ;; want is atan(y/x)-pi.
1002 ;; Returns a list of a mantissa and an exponent.
1012 ;; It seems to me that this function gets called on an integer
1013 ;; and returns the mantissa portion of the mantissa/exponent pair.
1015 ;; "STICKY BIT" CALCULATION FIXED 10/14/75 --RJF
1016 ;; BASE must not get temporarily bound to NIL by being placed
1017 ;; in a PROG list as this will confuse stepping programs.
1023 ;;*M will be positive if the precision of the argument is greater than
1024 ;;the current precision being used.
1029 ;;FPSHIFT is essentially LSH.
1036 ;ONLY ZEROES SHIFTED OFF
1048 ;; Compute (* L (expt d n)) where D is 2 or 10 depending on
1049 ;; *decfp. Throw away an fractional part by truncating to zero.
1053 ;; Left shift of negative number requires some
1054 ;; care. (That is, (truncate l (expt 2 n)), but use
1055 ;; shifts instead.)
1065 ;; Bignum LSH -- N is assumed (and declared above) to be a fixnum.
1066 ;; This isn't really LSH, since the sign bit isn't propagated when
1067 ;; shifting to the right, i.e. (BIGLSH -100 -3) = -40, whereas
1068 ;; (LSH -100 -3) = 777777777770 (on a 36 bit machine).
1069 ;; This actually computes (* X (EXPT 2 N)). As of 12/21/80, this function
1070 ;; was only called by FPSHIFT. I would like to hear an argument as why this
1071 ;; is more efficient than simply writing (* X (EXPT 2 N)). Is the
1072 ;; intermediate result created by (EXPT 2 N) the problem? I assume that
1073 ;; EXPT tries to LSH when possible.
1079 ;; Either we are shifting a fixnum to the right, or shifting
1080 ;; a fixnum to the left, but not far enough left for it to become
1081 ;; a bignum.
1085 ;; The form which follows is nearly identical to (ASH X N), however
1086 ;; (ASH -100 -20) = -1, whereas (BIGLSH -100 -20) = 0.
1090 ;; If we get here, then either X is a bignum or our answer is
1091 ;; going to be a bignum.
1098 ;; Isn't this the kind of optimization that compilers are
1099 ;; supposed to make?
1104 ;; exp(x)
1105 ;;
1106 ;; For negative x, use exp(-x) = 1/exp(x)
1107 ;;
1108 ;; For x > 0, exp(x) = exp(r+y) = exp(r) * exp(y), where x = r + y and
1109 ;; r = floor(x).
1125 ;; exp(x) for small x, using Taylor series.
1136 ;; Does one higher precision to round correctly.
1137 ;; A and B are each a list of a mantissa and an exponent.
1157 min))
1162 max))
1164 ;; The following functions compute bigfloat values for %e, %pi,
1165 ;; %gamma, and log(2). For each precision, the computed value is
1166 ;; cached in a hash table so it doesn't need to be computed again.
1167 ;; There are functions to return the hash table or clear the hash
1168 ;; table, for debugging.
1169 ;;
1170 ;; Note that each of these return a bigfloat number, but without the
1171 ;; bigfloat tag.
1172 ;;
1173 ;; See
1174 ;; https://sourceforge.net/p/maxima/bugs/1842/
1175 ;; for an explanation.
1180 value
1183 table)
1191 value
1194 table)
1202 value
1205 table)
1213 value
1216 table)
1220 ;; This doesn't need a hash table because there's never a problem with
1221 ;; using a high precision value and rounding to a lower precision
1222 ;; value because 1 is always an exact bfloat.
1228 ;;----------------------------------------------------------------------------;;
1229 ;;
1230 ;; The values of %e, %pi, %gamma and log(2) are computed by the technique of
1231 ;; binary splitting. See http://www.ginac.de/CLN/binsplit.pdf for details.
1232 ;;
1233 ;; Volker van Nek, Sept. 2014
1235 ;;
1236 ;; Euler's number E
1237 ;;
1241 ;;
1242 ;; Taylor: %e = sum(s[i] ,i,0,inf) where s[i] = 1/i!
1243 ;;
1248 ;;
1249 ;; binary splitting:
1250 ;;
1251 ;; 1
1252 ;; s[i] = ----------------------
1253 ;; q[0]*q[1]*q[2]*..*q[i]
1254 ;;
1255 ;; where q[0] = 1
1256 ;; q[i] = i
1257 ;;
1269 ;;
1270 ;; stop when i! > 2^fpprec
1271 ;;
1272 ;; log(i!) = sum(log(k), k,1,i) > fpprec * log(2)
1273 ;;
1280 ;;
1281 ;;----------------------------------------------------------------------------;;
1282 ;;
1283 ;; PI
1284 ;;
1288 ;;
1289 ;; Chudnovsky & Chudnovsky:
1290 ;;
1291 ;; C^(3/2)/(12*%pi) = sum(s[i], i,0,inf),
1292 ;;
1293 ;; where s[i] = (-1)^i*(6*i)!*(A*i+B) / (i!^3*(3*i)!*C^(3*i))
1294 ;;
1295 ;; and A = 545140134, B = 13591409, C = 640320
1296 ;;
1300 ;; STEP 1:
1301 ;; compute n/d = sqrt(10005) :
1302 ;;
1303 ;; n[0] n[i+1] = n[i]^2+a*d[i]^2 n[inf]
1304 ;; quadratic Heron: x[0] = ----, , sqrt(a) = ------
1305 ;; d[0] d[i+1] = 2*n[i]*d[i] d[inf]
1306 ;;
1312 ;; STEP 2:
1313 ;; divide C^(3/2)/12 = 3335*2^7*sqrt(10005)
1314 ;; by Chudnovsky-sum = tt/qq :
1315 ;;
1317 ;; fpprec*log(2)/log(C^3/(24*6*2*6))
1323 ;;
1324 ;; The returned n and d serve as start values for the iteration.
1325 ;; n/d = sqrt(10005) with a precision of p = ceiling(prec/2^nr) bits
1326 ;; where nr is the number of needed iterations.
1327 ;;
1342 ;;
1343 ;; binary splitting:
1344 ;;
1345 ;; a[i] * p[0]*p[1]*p[2]*..*p[i]
1346 ;; s[i] = -----------------------------
1347 ;; q[0]*q[1]*q[2]*..*q[i]
1348 ;;
1349 ;; where a[0] = B
1350 ;; p[0] = q[0] = 1
1351 ;; a[i] = A*i+B
1352 ;; p[i] = - (6*i-5)*(2*i-1)*(6*i-1)
1353 ;; q[i] = C^3/24*i^3
1354 ;;
1373 ;;
1374 ;;----------------------------------------------------------------------------;;
1375 ;;
1376 ;; Euler-Mascheroni constant GAMMA
1377 ;;
1381 ;;
1382 ;; Brent-McMillan algorithm
1383 ;;
1384 ;; Let
1385 ;; alpha = 4.970625759544
1386 ;;
1387 ;; n > 0 and N-1 >= alpha*n
1388 ;;
1389 ;; H(k) = sum(1/i, i,1,k)
1390 ;;
1391 ;; S = sum(H(k)*(n^k/k!)^2, k,0,N-1)
1392 ;;
1393 ;; I = sum((n^k/k!)^2, k,0,N-1)
1394 ;;
1395 ;; T = 1/(4*n)*sum((2*k)!^3/(k!^4*(16*n)^(2*k)), k,0,2*n-1)
1396 ;;
1397 ;; and
1398 ;; %gamma = S/I - T/I^2 - log(n)
1399 ;;
1400 ;; Then
1401 ;; |%gamma - gamma| < 24 * e^(-8*n)
1402 ;;
1403 ;; (Corollary 2, Remark 2, Brent/Johansson http://arxiv.org/pdf/1312.0039v1.pdf)
1404 ;;
1412 sumi sumi2 ;; I and I^2
1415 ;;
1416 ;; sums = vv/(qq*dd)
1417 ;; sumi = tt/qq
1418 ;; sums/sumi = vv/(qq*dd)*qq/tt = vv/(dd*tt)
1419 ;;
1423 ;;
1425 ;;
1426 ;; sumt = 1/(4*n)*ttt/qqq
1427 ;; sumt/sumi2 = ttt/(4*n*qqq*sumi2)
1428 ;;
1431 ;; %gamma :
1433 ;;
1434 ;; split S and I simultaneously:
1435 ;;
1436 ;; summands I[0] = 1, I[i]/I[i-1] = n^2/i^2
1437 ;;
1438 ;; S[0] = 0, S[i]/S[i-1] = n^2/i^2*H(i)/H(i-1)
1439 ;;
1440 ;; p[0]*p[1]*p[2]*..*p[i]
1441 ;; I[i] = ----------------------
1442 ;; q[0]*q[1]*q[2]*..*q[i]
1443 ;;
1444 ;; where p[0] = n^2
1445 ;; q[0] = 1
1446 ;; p[i] = n^2
1447 ;; q[i] = i^2
1448 ;; c[0] c[1] c[2] c[i]
1449 ;; S[i] = H[i] * I[i], where H[i] = ---- + ---- + ---- + .. + ----
1450 ;; d[0] d[1] d[2] d[i]
1451 ;; and c[0] = 0
1452 ;; d[0] = 1
1453 ;; c[i] = 1
1454 ;; d[i] = i
1455 ;;
1475 ;;
1476 ;; split 4*n*T:
1477 ;;
1478 ;; summands T[0] = 1, T[i]/T[i-1] = (2*i-1)^3/(32*i*n^2)
1479 ;;
1480 ;; p[0]*p[1]*p[2]*..*p[i]
1481 ;; T[i] = ----------------------
1482 ;; q[0]*q[1]*q[2]*..*q[i]
1483 ;;
1484 ;; where p[0] = q[0] = 1
1485 ;; p[i] = (2*i-1)^3
1486 ;; q[i] = 32*i*n^2
1487 ;;
1503 ;;
1504 ;;----------------------------------------------------------------------------;;
1505 ;;
1506 ;; log(2) = 18*L(26) - 2*L(4801) + 8*L(8749)
1507 ;;
1508 ;; where L(k) = atanh(1/k)
1509 ;;
1510 ;; see http://numbers.computation.free.fr/Constants/constants.html
1511 ;;
1512 ;;;(defun $log2 () (bcons (comp-log2))) ;; checked against reference table
1513 ;;
1521 ;;
1522 ;; Taylor: atanh(1/k) = sum(s[i], i,0,inf)
1523 ;;
1524 ;; where s[i] = 1/((2*i+1)*k^(2*i+1))
1525 ;;
1531 ;;
1532 ;; binary splitting:
1533 ;; 1
1534 ;; s[i] = -----------------------------
1535 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1536 ;;
1537 ;; where b[0] = 1
1538 ;; q[0] = k
1539 ;; b[i] = 2*i+1
1540 ;; q[i] = k^2
1541 ;;
1555 ;;
1556 ;;----------------------------------------------------------------------------;;
1557 ;;
1558 ;; log(n) = log(n/2^k) + k*log(2)
1559 ;;
1560 ;;;(defun $log10 () (bcons (log-n 10))) ;; checked against reference table
1561 ;;
1570 ;; choose k so that |n/2^k - 1| is as small as possible:
1572 ;; now |n/2^k - 1| <= 1/3
1576 ;;
1577 ;; log(1+u/v) = 2 * sum(s[i], i,0,inf)
1578 ;;
1579 ;; where s[i] = (u/(2*v+u))^(2*i+1)/(2*i+1)
1580 ;;
1595 ;;
1596 ;; binary splitting:
1597 ;;
1598 ;; p[0]*p[1]*p[2]*..*p[i]
1599 ;; s[i] = -----------------------------
1600 ;; b[i] * q[0]*q[1]*q[2]*..*q[i]
1601 ;;
1602 ;; where b[0] = 1
1603 ;; p[0] = u
1604 ;; q[0] = w = 2*v+u
1605 ;; b[i] = 2*i+1
1606 ;; p[i] = u^2
1607 ;; q[i] = w^2
1608 ;;
1623 ;;
1624 ;;----------------------------------------------------------------------------;;
1625 \f
1632 x
1649 ; COMPUTE STICKY BIT
1651 ; MAKE ROOM FOR GUARD DIGIT & STICKY BIT
1671 ;; Don't use the symbol BASE since it is SPECIAL.
1679 bas)))
1681 ;; NN is positive or negative integer
1720 $%i)))))
1735 ;; The number of extra bits required
1741 ;; Special case for a = 0b0. General algorithm loops endlessly in that case.
1743 ;; Unlike many or maybe all of the other functions named FP-something,
1744 ;; FPROOT assumes it is called with an argument like
1745 ;; '((BIGFLOAT ...) FOO BAR) instead of '(FOO BAR).
1746 ;; However FPROOT does return something like '(FOO BAR).
1773 fans
1777 ;; If A is a bigfloat, be sure to round it to the current precision.
1778 ;; (See Bug 2543079 for one of the symptoms.)
1804 ;; FPPREC 16, to get rid of the miniscule imaginary
1805 ;; part of the a+bi answer.
1812 e))
1820 ;; THIS VERSION OF FPSIN COMPUTES SIN OR COS TO PRECISION FPPREC,
1821 ;; BUT CHECKS FOR THE POSSIBILITY OF CATASTROPHIC CANCELLATION DURING
1822 ;; ARGUMENT REDUCTION (E.G. SIN(N*%PI+EPSILON))
1823 ;; *FPSINCHECK* WILL CAUSE PRINTOUT OF ADDITIONAL INFO WHEN
1824 ;; EXTRA PRECISION IS NEEDED FOR SIN/COS CALCULATION. KNOWN
1825 ;; BAD FEATURES: IT IS NOT NECESSARY TO USE EXTRA PRECISION FOR, E.G.
1826 ;; SIN(PI/2), WHICH IS NOT NEAR ZERO, BUT EXTRA
1827 ;; PRECISION IS USED SINCE IT IS NEEDED FOR COS(PI/2).
1828 ;; PRECISION SEEMS TO BE 100% SATSIFACTORY FOR LARGE ARGUMENTS, E.G.
1829 ;; SIN(31415926.0B0), BUT LESS SO FOR SIN(3.1415926B0). EXPLANATION
1830 ;; NOT KNOWN. (9/12/75 RJF)
1834 ;; FL is a T for sin and NIL for cos.
1878 ;; Compute SIN or COS in (0,PI/2). FL is T for SIN, NIL for COS.
1879 ;;
1880 ;; Use Taylor series
1898 x
1905 ;; Calculate the integer part of a floating point number that is represented as
1906 ;; a list
1907 ;;
1908 ;; (MANTISSA EXPONENT)
1909 ;;
1910 ;; The special variable fpprec should be bound to the precision (in bits) of the
1911 ;; number. This encodes how many bits are known of the result and also a right
1912 ;; shift. The pair denotes the number MANTISSA * 2^(EXPONENT - FPPREC), of which
1913 ;; FPPREC bits are known.
1914 ;;
1915 ;; If EXPONENT is large and positive then we might not have enough
1916 ;; information to calculate the integer part. Specifically, we only
1917 ;; have enough information if EXPONENT < FPPREC. If that isn't the
1918 ;; case, we signal a Maxima error. However, if SKIP-EXPONENT-CHECK-P
1919 ;; is non-NIL, this check is skipped, and we compute the integer part
1920 ;; as requested.
1921 ;;
1922 ;; For the bigfloat code here, skip-exponent-check-p should be true.
1923 ;; For other uses (see commit 576c7508 and bug #2784), this should be
1924 ;; nil, which is the default.
1927 f
1943 ;;; Computes the log of a bigfloat number.
1944 ;;;
1945 ;;; Uses the series
1946 ;;;
1947 ;;; log(1+x) = sum((x/(x+2))^(2*n+1)/(2*n+1),n,0,inf);
1948 ;;;
1949 ;;;
1950 ;;; INF x 2 n + 1
1951 ;;; ==== (-----)
1952 ;;; \ x + 2
1953 ;;; = 2 > --------------
1954 ;;; / 2 n + 1
1955 ;;; ====
1956 ;;; n = 0
1957 ;;;
1958 ;;;
1959 ;;; which converges for x > 0.
1960 ;;;
1961 ;;; Note that FPLOG is given 1+X, not X.
1962 ;;;
1963 ;;; However, to aid convergence of the series, we scale 1+x until 1/e
1964 ;;; < 1+x <= e.
1965 ;;;
1973 ;; Scale X until 1/e < X <= E. ANS keeps track of how
1974 ;; many factors of E were used. Set X to NIL if X is E.
1987 ;; Prepare X for the series. The series is for 1 + x, so
1988 ;; get x from our X. TERM is (x/(x+2)). X becomes
1989 ;; (x/(x+2))^2.
1993 ;; Sum the series until the sum (in ANS) doesn't change
1994 ;; anymore.
2010 \f
2011 ;;;; Bigfloat implementations of special functions.
2012 ;;;;
2014 ;;; This is still a bit messy. Some functions here take bigfloat
2015 ;;; numbers, represented by ((bigfloat) <mant> <exp>), but others want
2016 ;;; just the FP number, represented by (<mant> <exp>). Likewise, some
2017 ;;; return a bigfloat, some return just the FP.
2018 ;;;
2019 ;;; This needs to be systemized somehow. It isn't helped by the fact
2020 ;;; that some of the routines above also do the samething.
2021 ;;;
2022 ;;; The implementation for the special functions for a complex
2023 ;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex
2024 ;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in
2025 ;;; Iserles and Powell (eds.) "The State of the Art in Numerical
2026 ;;; Analysis", pp 165-211, Clarendon Press, 1987
2028 ;; Compute exp(x) - 1, but do it carefully to preserve precision when
2029 ;; |x| is small. X is a FP number, and a FP number is returned. That
2030 ;; is, no bigfloat stuff.
2032 ;; What is the right breakpoint here? Is 1 ok? Perhaps 1/e is better?
2034 ;; exp(x) - 1
2037 ;; Use Taylor series for exp(x) - 1
2046 ans))))
2048 ;; log(1+x) for small x. X is FP number, and a FP number is returned.
2050 ;; Use the same series as given above for fplog. For small x we use
2051 ;; the series, otherwise fplog is accurate enough.
2066 ;; sinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2068 ;; X must be a maxima bigfloat
2070 ;; See, for example, Hart et al., Computer Approximations, 6.2.27:
2071 ;;
2072 ;; sinh(x) = 1/2*(D(x) + D(x)/(1+D(x)))
2073 ;;
2074 ;; where D(x) = exp(x) - 1.
2075 ;;
2076 ;; But for negative x, use sinh(x) = -sinh(-x) because D(x)
2077 ;; approaches -1 for large negative x.
2079 ;; Special case: x=0. Return immediately.
2082 ;; x is positive.
2087 ;; x is negative.
2092 ;; The rectform for sinh for complex args should be numerically
2093 ;; accurate, so return nil in that case.
2097 ;; asinh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2099 ;; asinh(x) = sign(x) * log(|x| + sqrt(1+x*x))
2100 ;;
2101 ;; And
2102 ;;
2103 ;; asinh(x) = x, if 1+x*x = 1
2104 ;; = sign(x) * (log(2) + log(x)), large |x|
2105 ;; = sign(x) * log(2*|x| + 1/(|x|+sqrt(1+x*x))), if |x| > 2
2106 ;; = sign(x) * log1p(|x|+x^2/(1+sqrt(1+x*x))), otherwise.
2107 ;;
2108 ;; But I'm lazy right now and we only implement the last 2 cases.
2109 ;; We should implement all cases.
2115 result)
2116 ;; We only use two formulas here. |x| <= 2 and |x| > 2. Should
2117 ;; we add one for very big x and one for very small x, as given above.
2119 ;; |x| > 2
2120 ;;
2121 ;; log(2*|x| + 1/(|x|+sqrt(1+x^2)))
2129 ;; |x| <= 2
2130 ;;
2131 ;; log1p(|x|+x^2/(1+sqrt(1+x^2)))
2143 ;; asinh(z) = -%i * asin(%i*z)
2156 ;; asin(x) = atan(x/(sqrt(1-x^2))
2157 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2158 ;;
2159 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2160 ;;
2161 ;; If |x| > 1, we need to do something else.
2162 ;;
2163 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2164 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2165 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2166 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2170 ;; asin(-x) = -asin(x);
2173 ;; 0 <= x < 1/2
2174 ;; asin(x) = atan(x/sqrt(1-x^2))
2182 ;; x > 1
2183 ;; asin(x) = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2184 ;;
2185 ;; Should we try to do something a little fancier with the
2186 ;; argument to log and use log1p for better accuracy?
2195 ;; 1/2 <= x <= 1
2196 ;; asin(x) = %pi/2 - atan(sqrt(1-x^2)/x)
2204 fp-x)))))))))
2206 ;; asin(x) for real x. X is a bigfloat, and a maxima number (real or
2207 ;; complex) is returned.
2209 ;; asin(x) = atan(x/(sqrt(1-x^2))
2210 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))]
2211 ;;
2212 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1.
2213 ;;
2214 ;; If |x| > 1, we need to do something else.
2215 ;;
2216 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x)
2217 ;; = -%i*log(%i*x + %i*sqrt(x^2-1))
2218 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2]
2219 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|)
2223 ;; Square root of a complex number (xx, yy). Both are bigfloats. FP
2224 ;; (non-bigfloat) numbers are returned.
2242 rho))))
2245 ;; asin(z) for complex z = x + %i*y. X and Y are bigfloats. The real
2246 ;; and imaginary parts are returned as bigfloat numbers.
2256 ;; Realpart is atan(x/Re(sqrt(1-z)*sqrt(1+z)))
2257 ;; Imagpart is asinh(Im(conj(sqrt(1-z))*sqrt(1+z)))
2263 ;; Check for division by zero. If we would divide
2264 ;; by zero, return pi/2 or -pi/2 according to the
2265 ;; sign of X.
2285 ;; tanh(x) for real x. X is a bigfloat, and a bigfloat is returned.
2287 ;; X is Maxima bigfloat
2288 ;; tanh(x) = D(2*x)/(2+D(2*x))
2294 ;; tanh(z), z = x + %i*y. X, Y are bigfloats, and a maxima number is
2295 ;; returned.
2312 ;; atanh(x) for real x, |x| <= 1. X is a bigfloat, and a bigfloat is
2313 ;; returned.
2315 ;; atanh(x) = -atanh(-x)
2316 ;; = 1/2*log1p(2*x/(1-x)), x >= 0.5
2317 ;; = 1/2*log1p(2*x+2*x*x/(1-x)), x <= 0.5
2321 ;; atanh(x) = -atanh(-x)
2324 ;; x > 1, so use complex version.
2329 ;; atanh(x) = 1/2*log1p(2*x/(1-x))
2335 ;; atanh(x) = 1/2*log1p(2*x + 2*x*x/(1-x))
2343 ;; Stuff which follows is derived from atanh z = (log(1 + z) - log(1 - z))/2
2344 ;; which apparently originates with Kahan's "Much ado" paper.
2346 ;; The formulas for eta and nu below can be easily derived from
2347 ;; rectform(atanh(x+%i*y)) =
2348 ;;
2349 ;; 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)))
2350 ;;
2351 ;; Expand the argument of log out and divide it out and we get
2352 ;;
2353 ;; log(((1+x)^2+y^2)/((1-x)^2+y^2)) = log(1+4*x/((1-x)^2+y^2))
2354 ;;
2355 ;; When y = 0, Im atanh z = 1/2 (arg(1 + x) - arg(1 - x))
2356 ;; = if x < -1 then %pi/2 else if x > 1 then -%pi/2 else <whatever>
2357 ;;
2358 ;; Otherwise, arg(1 - x + %i*(-y)) = - arg(1 - x + %i*y),
2359 ;; and Im atanh z = 1/2 (arg(1 + x + %i*y) + arg(1 - x + %i*y)).
2360 ;; Since arg(x)+arg(y) = arg(x*y) (almost), we can simplify the
2361 ;; imaginary part to
2362 ;;
2363 ;; arg((1+x+%i*y)*(1-x+%i*y)) = arg((1-x)*(1+x)-y^2+2*y*%i)
2364 ;; = atan2(2*y,((1-x)*(1+x)-y^2))
2365 ;;
2366 ;; These are the eta and nu forms below.
2378 ;; Kahan has rho = 4/most-positive-float. What should we do
2379 ;; here about that? Our big floats don't really have a
2380 ;; most-positive float value.
2385 ;; eta = log(1+4*x/((1-x)^2+y^2))/4
2391 ;; If y = 0, then Im atanh z = %pi/2 or -%pi/2.
2392 ;; Otherwise nu = 1/2*atan2(2*y,(1-x)*(1+x)-y^2)
2394 ;; EXTRA FPMINUS HERE TO COUNTERACT FPMINUS IN RETURN VALUE
2405 ;; WTF IS FPMINUS DOING HERE ??
2414 ;; acos(x) for real x. X is a bigfloat, and a maxima number is returned.
2416 ;; acos(x) = %pi/2 - asin(x)
2447 ;; No warning message, please.
2470 ;; x is (mantissa exp), where mantissa = frac*2^fpprec,
2471 ;; with 1/2 < frac <= 1 and x is frac*2^exp. To
2472 ;; compute log(x), use log(x) = log(frac)+ exp*log(2).
2488 ;; Special case for log(1). See:
2489 ;; https://sourceforge.net/p/maxima/bugs/2240/
2492 ;; ??? Do we want to return an exact %i*%pi or a float
2493 ;; approximation?