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;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module float) 14 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. 24 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. 29 30(eval-when 31 #+gcl (compile load eval) 32 #-gcl (:compile-toplevel :load-toplevel :execute) 33 (defconstant +machine-fixnum-precision+ (integer-length most-positive-fixnum))) 34 35;; External variables 36 37(defmvar $float2bf t 38 "If TRUE, no MAXIMA-ERROR message is printed when a floating point number is 39converted to a bigfloat number.") 40 41(defmvar $bftorat nil 42 "Controls the conversion of bigfloat numbers to rational numbers. If 43FALSE, RATEPSILON will be used to control the conversion (this results in 44relatively small rational numbers). If TRUE, the rational number generated 45will accurately represent the bigfloat.") 46 47(defmvar $bftrunc t 48 "If TRUE, printing of bigfloat numbers will truncate trailing zeroes. 49 Otherwise, all trailing zeroes are printed.") 50 51(defmvar $fpprintprec 0 52 "Controls the number of significant digits printed for floats. If 53 0, then full precision is used." 54 fixnum) 55 56(defmvar $maxfpprintprec (ceiling (log (expt 2 (float-digits 1.0)) 10.0)) 57 "The maximum number of significant digits printed for floats.") 58 59(defmvar $fpprec $maxfpprintprec 60 "Number of decimal digits of precision to use when creating new bigfloats. 61One extra decimal digit in actual representation for rounding purposes.") 62 63(defmvar bigfloatzero '((bigfloat simp 56.) 0 0) 64 "Bigfloat representation of 0" in-core) 65 66(defmvar bigfloatone '((bigfloat simp 56.) #.(expt 2 55.) 1) 67 "Bigfloat representation of 1" in-core) 68 69(defmvar bfhalf '((bigfloat simp 56.) #.(expt 2 55.) 0) 70 "Bigfloat representation of 1/2") 71 72(defmvar bfmhalf '((bigfloat simp 56.) #.(- (expt 2 55.)) 0) 73 "Bigfloat representation of -1/2") 74 75(defmvar bigfloat%e '((bigfloat simp 56.) 48968212118944587. 2) 76 "Bigfloat representation of %E") 77 78(defmvar bigfloat%pi '((bigfloat simp 56.) 56593902016227522. 2) 79 "Bigfloat representation of %pi") 80 81(defmvar bigfloat%gamma '((bigfloat simp 56.) 41592772053807304. 0) 82 "Bigfloat representation of %gamma") 83 84(defmvar bigfloat_log2 '((bigfloat simp 56.) 49946518145322874. 0) 85 "Bigfloat representation of log(2)") 86 87;; Internal specials 88 89;; Number of bits of precision in the mantissa of newly created bigfloats. 90;; FPPREC = ($FPPREC+1)*(Log base 2 of 10) 91 92(defvar fpprec) 93 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. 99 100(defvar *m) 101 102;; *DECFP = T if the computation is being done in decimal radix. NIL implies 103;; base 2. Decimal radix is used only during output. 104 105(defvar *decfp nil) 106 107(defvar max-bfloat-%pi bigfloat%pi) 108(defvar max-bfloat-%e bigfloat%e) 109(defvar max-bfloat-%gamma bigfloat%gamma) 110(defvar max-bfloat-log2 bigfloat_log2) 111 112 113(declare-top (special *cancelled $float $bfloat $ratprint $ratepsilon $domain $m1pbranch)) 114 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)). 122 123(defun hipart (x nn) 124 (if (bignump nn) 125 (abs x) 126 (haipart x nn))) 127 128(defun fpprec1 (assign-var q) 129 (declare (ignore assign-var)) 130 (if (or (not (fixnump q)) (< q 1)) 131 (merror (intl:gettext "fpprec: value must be a positive integer; found: ~M") q)) 132 (setq fpprec (+ 2 (integer-length (expt 10. q))) 133 bigfloatone ($bfloat 1) 134 bigfloatzero ($bfloat 0) 135 bfhalf (list (car bigfloatone) (cadr bigfloatone) 0) 136 bfmhalf (list (car bigfloatone) (- (cadr bigfloatone)) 0)) 137 q) 138 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. 144 145(defun fpscan (lft rt exp &aux (*read-base* 10.) (*m 1) (*cancelled 0)) 146 (setq exp (readlist exp)) 147 (bigfloatp 148 (let ((fpprec (+ 4 fpprec (integer-length exp) 149 (floor (1+ (* #.(/ (log 10.0) (log 2.0)) (length lft)))))) 150 $float temp) 151 (setq temp (add (readlist lft) 152 (div (readlist rt) (expt 10. (length rt))))) 153 ($bfloat (cond ((> (abs exp) 1000.) 154 (cons '(mtimes) (list temp (list '(mexpt) 10. exp)))) 155 (t (mul2 temp (power 10. exp)))))))) 156 157(defun dim-bigfloat (form result) 158 (let (($lispdisp nil)) 159 (dimension-atom (maknam (fpformat form)) result))) 160 161;; Assume that X has the form ((BIGFLOAT ... <prec>) ...). 162;; Return <prec>. 163(defun bigfloat-prec (x) 164 (car (last (car x)))) 165 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. 169(defun fpformat (l) 170 (if (not (member 'simp (cdar l) :test #'eq)) 171 (setq l (cons (cons (caar l) (cons 'simp (cdar l))) (cdr l)))) 172 (cond ((equal (cadr l) 0) 173 (if (not (equal (caddr l) 0)) 174 (mtell "FPFORMAT: warning: detected an incorrect form of 0.0b0: ~M, ~M~%" 175 (cadr l) (caddr l))) 176 (list '|0| '|.| '|0| '|b| '|0|)) 177 (t ;; L IS ALWAYS POSITIVE FP NUMBER 178 (let* ((extradigs (floor (1+ (quotient (integer-length (caddr l)) #.(/ (log 10.0) (log 2.0)))))) 179 (fpprec (+ extradigs (decimalsin (- (bigfloat-prec l) 2)))) 180 (*m 1) 181 (*cancelled 0)) 182 (setq l 183 (let ((*decfp t) 184 (of (bigfloat-prec l)) 185 (l (cdr l)) 186 (expon nil)) 187 (setq expon (- (cadr l) of)) 188 (setq l (if (minusp expon) 189 (fpquotient (intofp (car l)) (fpintexpt 2 (- expon) of)) 190 (fptimes* (intofp (car l)) (fpintexpt 2 expon of)))) 191 (incf fpprec (- extradigs)) 192 (list (fpround (car l)) (+ (- extradigs) *m (cadr l))))) 193 (let ((*print-base* 10.) 194 *print-radix* 195 (l1 nil)) 196 (setq l1 (let* 197 ((effective-printprec (if (or (= $fpprintprec 0) (> $fpprintprec fpprec)) fpprec $fpprintprec)) 198 (integer-to-explode (round (car l) (expt 10 (- fpprec effective-printprec)))) 199 (exploded-integer (explodec integer-to-explode))) 200 (if $bftrunc 201 (do ((l (nreverse exploded-integer) (cdr l))) 202 ((not (eq '|0| (car l))) (nreverse l))) 203 exploded-integer))) 204 (nconc (ncons (car l1)) (ncons '|.|) 205 (or (cdr l1) (ncons '|0|)) 206 (ncons '|b|) 207 (explodec (1- (cadr l))))))))) 208 209;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to 210;; support printing of bfloats. 211(defun bfloat-format-e (stream arg colonp atp 212 &optional w d e (k 1) 213 overflowchar (padchar #\space) exponentchar) 214 (declare (ignore colonp)) 215 (flet ((exponent-value (x) 216 ;; Compute the (decimal exponent) of the bfloat number X. 217 (let* (($fpprintprec 1) 218 (f (fpformat x)) 219 (marker (position '|b| f))) 220 ;; FIXME: do something better than printing and reading 221 ;; the result. 222 (read-from-string 223 (format nil "~{~A~}" (nthcdr (1+ marker) f))))) 224 (bfloat-to-string (x fdigits scale) 225 ;; Print the bfloat X with FDIGITS after the decimal 226 ;; point. This means, roughtly, FDIGITS+1 significant 227 ;; digits. 228 (let* (($fpprintprec (if fdigits 229 (if (zerop fdigits) 230 1 231 (+ fdigits scale)) 232 0)) 233 (f (fpformat (bcons (fpabs (cdr x))))) 234 (marker (position '|b| f)) 235 (digits (remove '|.| (subseq f 0 marker)))) 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. 241 (when (and fdigits (zerop fdigits)) 242 (setf digits (butlast digits))) 243 (cond ((zerop k) 244 (push '|.| digits)) 245 ((minusp k) 246 ;; Put the leading decimal and then some zeroes 247 (dotimes (i (abs k)) 248 (push #\0 digits)) 249 (push '|.| digits)) 250 (t 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. 254 (setf digits 255 (cond ((> k (length digits)) 256 (concatenate 'list 257 digits 258 (make-list (- k (length digits)) 259 :initial-element #\0) 260 (list '|.|))) 261 (t 262 (concatenate 'list 263 (subseq digits 0 k) 264 (list '|.|) 265 (subseq digits k))))))) 266 (let* ((str (format nil "~{~A~}" digits)) 267 (len (length str))) 268 (when (and fdigits (>= fdigits len)) 269 ;; Append some zeroes to get the desired number of digits 270 (setf str (concatenate 'string str 271 (make-string (+ 1 k (- fdigits len)) 272 :initial-element #\0))) 273 (setf len (length str))) 274 (values str 275 len 276 (char= (aref str 0) #\.) 277 (char= (aref str (1- (length str))) #\.) 278 1 279 0))))) 280 (let* ((num-expt (exponent-value arg)) 281 (expt (if (zerop (second arg)) 282 0 283 (1+ (- num-expt k)))) 284 (estr (format nil "~D" (abs expt))) 285 (elen (if e (max (length estr) e) (length estr))) 286 (add-zero-p nil)) 287 (cond ((and w overflowchar e (> elen e)) 288 ;; Exponent overflow 289 (dotimes (i w) 290 (write-char overflowchar stream))) 291 (t 292 ;; The hairy case 293 (let* ((fdig (if d 294 (if (plusp k) 295 (1+ (- d k)) 296 d) 297 nil)) 298 (spaceleft (if w 299 (- w 2 elen 300 (if (or atp (minusp (second arg))) 301 1 0)) 302 nil))) 303 #+(or) 304 (progn 305 (format t "d, k = ~D ~D~%" d k) 306 (format t "fdig = ~D, spaceleft = ~D~%" fdig spaceleft)) 307 308 (multiple-value-bind (fstr flen lpoint tpoint) 309 (bfloat-to-string arg fdig (or k 1)) 310 #+(or) 311 (format t "fstr flen lpoint tpoint = ~S ~S ~S ~S~%" 312 fstr flen lpoint tpoint) 313 (when (and d (zerop d)) (setq tpoint nil)) 314 (when w 315 (decf spaceleft flen) 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. 322 (when (and (null d) (char= (aref fstr (1- flen)) #\.)) 323 (setf add-zero-p t) 324 (decf spaceleft)) 325 (when lpoint 326 (if (or (> spaceleft 0) tpoint) 327 (decf spaceleft) 328 (setq lpoint nil))) 329 (when (and tpoint (<= spaceleft 0)) 330 (setq tpoint nil))) 331 #+(or) 332 (format t "w, spaceleft overflowchar = ~S ~S ~S~%" 333 w spaceleft overflowchar) 334 (cond ((and w (< spaceleft 0) overflowchar) 335 ;; Significand overflow; output the overflow char 336 (dotimes (i w) 337 (write-char overflowchar stream))) 338 (t 339 (when w 340 (dotimes (i spaceleft) 341 (write-char padchar stream))) 342 (if (minusp (second arg)) 343 (write-char #\- stream) 344 (when atp (write-char #\+ stream))) 345 (when lpoint 346 (write-char #\0 stream)) 347 348 (write-string fstr stream) 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.) 356 (when (or add-zero-p 357 (and (null w) 358 (char= (aref fstr (1- flen)) #\.))) 359 (write-char #\0 stream)) 360 (write-char (if exponentchar 361 exponentchar 362 #\b) 363 stream) 364 (write-char (if (minusp expt) #\- #\+) stream) 365 (when e 366 (dotimes (i (- e (length estr))) 367 (write-char #\0 stream))) 368 (write-string estr stream))))))))) 369 (values)) 370 371;; NOTE: This is a modified version of FORMAT-FIXED-AUX from CMUCL to 372;; support printing of bfloats. 373(defun bfloat-format-f (stream number colonp atsign &optional w d (k 0) ovf (pad #\space)) 374 (declare (ignore colonp)) 375 (labels 376 ((exponent-value (x) 377 ;; Compute the (decimal exponent) of the bfloat number X. 378 (let* (($fpprintprec 1) 379 (f (fpformat x)) 380 (marker (position '|b| f))) 381 ;; FIXME: do something better than printing and reading 382 ;; the result. 383 (read-from-string 384 (format nil "~{~A~}" (nthcdr (1+ marker) f))))) 385 (bfloat-to-string (x fdigits scale spaceleft) 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. 391 (flet ((compute-prec (exp spaceleft) 392 #+nil 393 (format t "compute-prec ~D ~D~%" exp spaceleft) 394 (cond (fdigits 395 (+ fdigits exp 1)) 396 (spaceleft 397 (max (1- spaceleft) (1+ exp))) 398 (t 399 (max (1+ exp) 0))))) 400 (let* ((exp (+ k (exponent-value x))) 401 ($fpprintprec (compute-prec exp spaceleft)) 402 (f (let ((maxima::$bftrunc nil)) 403 #+nil 404 (format t "printprec = ~D~%" $fpprintprec) 405 (fpformat (bcons (fpabs (cdr x)))))) 406 (marker (position '|b| f)) 407 (digits (remove '|.| (subseq f 0 marker)))) 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 414 (format t "exp, fdigits = ~D ~D, digits = ~S~%" exp fdigits digits) 415 #+nil 416 (when (and fdigits (zerop fdigits)) 417 (setf digits (butlast digits))) 418 ;; Figure out where the decimal point should go. An 419 ;; exponent of 0 means the decimal is after the first 420 ;; digit. 421 (cond ((minusp exp) 422 (dotimes (k (1- (abs exp))) 423 (push '|0| digits)) 424 (push '|.| digits)) 425 ((< exp (length digits)) 426 #+nil 427 (format t "exp, len = ~D ~D~%" exp (length digits)) 428 (setf digits (concatenate 'list 429 (subseq digits 0 (1+ exp)) 430 (list '|.|) 431 (subseq digits (1+ exp))))) 432 (t 433 (setf digits (append digits (list '|.|))))) 434 (let* ((str (format nil "~{~A~}" digits)) 435 (len (length str))) 436 #+nil 437 (format t "str = ~S~%" str) 438 (when (and fdigits (>= fdigits len)) 439 ;; Append some zeroes to get the desired number of digits 440 (setf str (concatenate 'string str 441 (make-string (+ 1 scale (- fdigits len)) 442 :initial-element #\0))) 443 (setf len (length str))) 444 (values str 445 len 446 (char= (aref str 0) #\.) 447 (char= (aref str (1- (length str))) #\.) 448 1 449 0)))))) 450 (let ((spaceleft w)) 451 (when (and w (or atsign (minusp (second number)))) 452 (decf spaceleft)) 453 (multiple-value-bind (str len lpoint tpoint) 454 (bfloat-to-string number d k spaceleft) 455 ;;if caller specifically requested no fraction digits, suppress the 456 ;;optional trailing zero 457 (when (and d (zerop d)) (setq tpoint nil)) 458 (when w 459 (decf spaceleft len) 460 ;;optional leading zero 461 (when lpoint 462 (if (or (> spaceleft 0) tpoint) ;force at least one digit 463 (decf spaceleft) 464 (setq lpoint nil))) 465 ;;optional trailing zero 466 (when tpoint 467 (if (> spaceleft 0) 468 (decf spaceleft) 469 (setq tpoint nil)))) 470 (cond ((and w (< spaceleft 0) ovf) 471 ;;field width overflow 472 (dotimes (i w) (write-char ovf stream)) 473 t) 474 (t 475 (when w (dotimes (i spaceleft) (write-char pad stream))) 476 (if (minusp (second number)) 477 (write-char #\- stream) 478 (if atsign (write-char #\+ stream))) 479 (when lpoint (write-char #\0 stream)) 480 (write-string str stream) 481 (when tpoint (write-char #\0 stream)) 482 nil)))))) 483 484;; NOTE: This is a modified version of FORMAT-EXP-AUX from CMUCL to 485;; support printing of bfloats. 486(defun bfloat-format-g (stream arg colonp atsign 487 &optional w d e (k 1) 488 ovf (pad #\space) exponentchar) 489 (declare (ignore colonp)) 490 (flet ((exponent-value (x) 491 ;; Compute the (decimal exponent) of the bfloat number X. 492 (let* (($fpprintprec 1) 493 (f (fpformat x)) 494 (marker (position '|b| f))) 495 ;; FIXME: do something better than printing and reading 496 ;; the result. 497 (read-from-string 498 (format nil "~{~A~}" (nthcdr (1+ marker) f))))) 499 (bfloat-to-string (x fdigits) 500 ;; Print the bfloat X with FDIGITS after the decimal 501 ;; point. This means, roughtly, FDIGITS+1 significant 502 ;; digits. 503 (let* (($fpprintprec (if fdigits 504 (if (zerop fdigits) 505 1 506 (1+ fdigits)) 507 0)) 508 (f (fpformat (bcons (fpabs (cdr x))))) 509 (marker (position '|b| f)) 510 (digits (remove '|.| (subseq f 0 marker)))) 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. 516 (when (and fdigits (zerop fdigits)) 517 (setf digits (butlast digits))) 518 (cond ((zerop k) 519 (push '|.| digits)) 520 ((minusp k) 521 ;; Put the leading decimal and then some zeroes 522 (dotimes (i (abs k)) 523 (push #\0 digits)) 524 (push '|.| digits)) 525 (t 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. 529 (setf digits 530 (cond ((> k (length digits)) 531 (concatenate 'list 532 digits 533 (make-list (- k (length digits)) 534 :initial-element #\0) 535 (list '|.|))) 536 (t 537 (concatenate 'list 538 (subseq digits 0 k) 539 (list '|.|) 540 (subseq digits k))))))) 541 (let* ((str (format nil "~{~A~}" digits)) 542 (len (length str))) 543 (when (and fdigits (>= fdigits len)) 544 ;; Append some zeroes to get the desired number of digits 545 (setf str (concatenate 'string str 546 (make-string (+ 1 k (- fdigits len)) 547 :initial-element #\0))) 548 (setf len (length str))) 549 (values str 550 len 551 (char= (aref str 0) #\.) 552 (char= (aref str (1- (length str))) #\.) 553 1 554 0))))) 555 (let* ((n (1+ (exponent-value arg))) 556 (orig-d d)) 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. 567 (unless d 568 (let* ((q (1- (nth-value 1 (bfloat-to-string arg nil))))) 569 (setq d (max q (min n 7))))) 570 (let* ((ee (if e (+ e 2) 4)) 571 (ww (if w (- w ee) nil)) 572 (dd (- d n))) 573 #+(or) 574 (progn 575 (format t "d = ~A~%" d) 576 (format t "ee = ~A~%" ee) 577 (format t "ww = ~A~%" ww) 578 (format t "dd = ~A~%" dd) 579 (format t "n = ~A~%" n)) 580 (cond ((<= 0 dd d) 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. 584 (let* ((fill-char (if (bfloat-format-f stream arg nil atsign 585 ww 586 dd 587 0 588 ovf pad) 589 ovf 590 #\space))) 591 (dotimes (i ee) (write-char fill-char stream)))) 592 (t 593 (bfloat-format-e stream arg nil atsign 594 w 595 orig-d 596 e (or k 1) 597 ovf pad exponentchar))))))) 598 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. 604(defun bigfloatp (x) 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. 607 (prog (x-prec) 608 (cond ((not ($bfloatp x)) (return nil)) 609 ((= fpprec (setq x-prec (bigfloat-prec x))) 610 ;; Precision matches. (Should we fix up bogus bigfloat 611 ;; zeros?) 612 (return x)) 613 ((> fpprec x-prec) 614 ;; Current precision is higher than bigfloat precision. 615 ;; Scale up mantissa and adjust exponent to get the 616 ;; correct precision. 617 (setq x (bcons (list (fpshift (cadr x) (- fpprec x-prec)) 618 (caddr x))))) 619 (t 620 ;; Current precision is LOWER than bigfloat precision. 621 ;; Round the number to the desired precision. 622 (setq x (bcons (list (fpround (cadr x)) 623 (+ (caddr x) *m fpprec (- x-prec))))))) 624 ;; Fix up any bogus zeros that we might have created. 625 (return (if (equal (cadr x) 0) (bcons (list 0 0)) x)))) 626 627(defun bigfloat2rat (x) 628 (setq x (bigfloatp x)) 629 (let (($float2bf t) 630 (exp nil) 631 (y nil) 632 (sign nil)) 633 (setq exp (cond ((minusp (cadr x)) 634 (setq sign t 635 y (fpration1 (cons (car x) (fpabs (cdr x))))) 636 (rplaca y (* -1 (car y)))) 637 (t (fpration1 x)))) 638 (when $ratprint 639 (princ "`rat' replaced ") 640 (when sign (princ "-")) 641 (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) 642 (princ " by ") 643 (princ (car exp)) 644 (write-char #\/) 645 (princ (cdr exp)) 646 (princ " = ") 647 (setq x ($bfloat (list '(rat simp) (car exp) (cdr exp)))) 648 (when sign (princ "-")) 649 (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) 650 (terpri) 651 (finish-output)) 652 exp)) 653 654(defun fpration1 (x) 655 (let ((fprateps (cdr ($bfloat (if $bftorat 656 (list '(rat simp) 1 (exptrl 2 (1- fpprec))) 657 $ratepsilon))))) 658 (or (and (equal x bigfloatzero) (cons 0 1)) 659 (prog (y a) 660 (return (do ((xx x (setq y (invertbigfloat 661 (bcons (fpdifference (cdr xx) (cdr ($bfloat a))))))) 662 (num (setq a (fpentier x)) 663 (+ (* (setq a (fpentier y)) num) onum)) 664 (den 1 (+ (* a den) oden)) 665 (onum 1 num) 666 (oden 0 den)) 667 ((and (not (zerop den)) 668 (not (fpgreaterp 669 (fpabs (fpquotient 670 (fpdifference (cdr x) 671 (fpquotient (cdr ($bfloat num)) 672 (cdr ($bfloat den)))) 673 (cdr x))) 674 fprateps))) 675 (cons num den)))))))) 676 677(defun float-nan-p (x) 678 (and (floatp x) (not (= x x)))) 679 680(defun float-inf-p (x) 681 (and (floatp x) (not (float-nan-p x)) (beyond-extreme-values x))) 682 683(defun beyond-extreme-values (x) 684 (multiple-value-bind (most-negative most-positive) (extreme-float-values x) 685 (cond 686 ((< x 0) (< x most-negative)) 687 ((> x 0) (> x most-positive)) 688 (t nil)))) 689 690(defun extreme-float-values (x) 691 ;; BLECHH, I HATE ENUMERATING CASES. IS THERE A BETTER WAY ?? 692 (case (type-of x) 693 (short-float (values most-negative-short-float most-positive-short-float)) 694 (single-float (values most-negative-single-float most-positive-single-float)) 695 (double-float (values most-negative-double-float most-positive-double-float)) 696 (long-float (values most-negative-long-float most-positive-long-float)) 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 )) 702 703;; Convert a floating point number into a bigfloat. 704(defun floattofp (x) 705 (if (float-nan-p x) 706 (merror (intl:gettext "bfloat: attempted conversion of floating point NaN (not-a-number).~%"))) 707 (if (float-inf-p x) 708 (merror (intl:gettext "bfloat: attempted conversion of floating-point infinity.~%"))) 709 (unless $float2bf 710 (mtell (intl:gettext "bfloat: converting float ~S to bigfloat.~%") x)) 711 712 ;; Need to check for zero because different lisps return different 713 ;; values for integer-decode-float of a 0. In particular CMUCL 714 ;; returns 0, -1075. A bigfloat zero needs to have an exponent and 715 ;; mantissa of zero. 716 (if (zerop x) 717 (list 0 0) 718 (multiple-value-bind (frac exp sign) 719 (integer-decode-float x) 720 ;; Scale frac to the desired number of bits, and adjust the 721 ;; exponent accordingly. 722 (let ((scale (- fpprec (integer-length frac)))) 723 (list (ash (* sign frac) scale) 724 (+ fpprec (- exp scale))))))) 725 726;; Convert a bigfloat into a floating point number. 727(defun fp2flo (l) 728 (let ((precision (bigfloat-prec l)) 729 (mantissa (cadr l)) 730 (exponent (caddr l)) 731 (fpprec machine-mantissa-precision) 732 (*m 0)) 733 ;; Round the mantissa to the number of bits of precision of the 734 ;; machine, and then convert it to a floating point fraction. We 735 ;; have 0.5 <= mantissa < 1 736 (setq mantissa (quotient (fpround mantissa) (expt 2.0 machine-mantissa-precision))) 737 ;; Multiply the mantissa by the exponent portion. I'm not sure 738 ;; why the exponent computation is so complicated. 739 ;; 740 ;; GCL doesn't signal overflow from scale-float if the number 741 ;; would overflow. We have to do it this way. 0.5 <= mantissa < 742 ;; 1. The largest double-float is .999999 * 2^1024. So if the 743 ;; exponent is 1025 or higher, we have an overflow. 744 (let ((e (+ exponent (- precision) *m machine-mantissa-precision))) 745 (if (>= e 1025) 746 (merror (intl:gettext "float: floating point overflow converting ~:M") l) 747 (scale-float mantissa e))))) 748 749;; New machine-independent version of FIXFLOAT. This may be buggy. - CWH 750;; It is buggy! On the PDP10 it dies on (RATIONALIZE -1.16066076E-7) 751;; which calls FLOAT on some rather big numbers. ($RATEPSILON is approx. 752;; 7.45E-9) - JPG 753 754(defun fixfloat (x) 755 (let (($ratepsilon (expt 2.0 (- machine-mantissa-precision)))) 756 (maxima-rationalize x))) 757 758;; Takes a flonum arg and returns a rational number corresponding to the flonum 759;; in the form of a dotted pair of two integers. Since the denominator will 760;; always be a positive power of 2, this number will not always be in lowest 761;; terms. 762 763(defun bcons (s) 764 `((bigfloat simp ,fpprec) . ,s)) 765 766(defmfun $bfloat (x) 767 (let (y) 768 (cond ((bigfloatp x)) 769 ((or (numberp x) 770 (member x '($%e $%pi $%gamma) :test #'eq)) 771 (bcons (intofp x))) 772 ((or (atom x) (member 'array (cdar x) :test #'eq)) 773 (if (eq x '$%phi) 774 ($bfloat '((mtimes simp) 775 ((rat simp) 1 2) 776 ((mplus simp) 1 ((mexpt simp) 5 ((rat simp) 1 2))))) 777 x)) 778 ((eq (caar x) 'mexpt) 779 (if (equal (cadr x) '$%e) 780 (*fpexp ($bfloat (caddr x))) 781 (exptbigfloat ($bfloat (cadr x)) (caddr x)))) 782 ((eq (caar x) 'mncexpt) 783 (list '(mncexpt) ($bfloat (cadr x)) (caddr x))) 784 ((eq (caar x) 'rat) 785 (ratbigfloat (cdr x))) 786 ((setq y (safe-get (caar x) 'floatprog)) 787 (funcall y (mapcar #'$bfloat (cdr x)))) 788 ((or (trigp (caar x)) (arcp (caar x)) (eq (caar x) '$entier)) 789 (setq y ($bfloat (cadr x))) 790 (if ($bfloatp y) 791 (cond ((eq (caar x) '$entier) ($entier y)) 792 ((arcp (caar x)) 793 (setq y ($bfloat (logarc (caar x) y))) 794 (if (free y '$%i) 795 y (let ($ratprint) (fparcsimp ($rectform y))))) 796 ((member (caar x) '(%cot %sec %csc) :test #'eq) 797 (invertbigfloat 798 ($bfloat (list (ncons (safe-get (caar x) 'recip)) y)))) 799 (t ($bfloat (exponentialize (caar x) y)))) 800 (subst0 (list (ncons (caar x)) y) x))) 801 (t (recur-apply #'$bfloat x))))) 802 803(defprop mplus addbigfloat floatprog) 804(defprop mtimes timesbigfloat floatprog) 805(defprop %sin sinbigfloat floatprog) 806(defprop %cos cosbigfloat floatprog) 807(defprop rat ratbigfloat floatprog) 808(defprop %atan atanbigfloat floatprog) 809(defprop %tan tanbigfloat floatprog) 810(defprop %log logbigfloat floatprog) 811(defprop mabs mabsbigfloat floatprog) 812 813(defun addbigfloat (h) 814 (prog (fans tst r nfans) 815 (setq fans (setq tst bigfloatzero) nfans 0) 816 (do ((l h (cdr l))) 817 ((null l)) 818 (cond ((setq r (bigfloatp (car l))) 819 (setq fans (bcons (fpplus (cdr r) (cdr fans))))) 820 (t (setq nfans (list '(mplus) (car l) nfans))))) 821 (return (cond ((equal nfans 0) fans) 822 ((equal fans tst) nfans) 823 (t (simplify (list '(mplus) fans nfans))))))) 824 825(defun ratbigfloat (r) 826 ;; R is a Maxima ratio, represented as a list of the numerator and 827 ;; denominator. FLOAT-RATIO doesn't like it if the numerator is 0, 828 ;; so handle that here. 829 (if (zerop (car r)) 830 (bcons (list 0 0)) 831 (bcons (float-ratio r)))) 832 833;; This is borrowed from CMUCL (float-ratio-float), and modified for 834;; converting ratios to Maxima's bfloat numbers. 835(defun float-ratio (x) 836 (let* ((signed-num (first x)) 837 (plusp (plusp signed-num)) 838 (num (if plusp signed-num (- signed-num))) 839 (den (second x)) 840 (digits fpprec) 841 (scale 0)) 842 (declare (fixnum digits scale)) 843 ;; 844 ;; Strip any trailing zeros from the denominator and move it into the scale 845 ;; factor (to minimize the size of the operands.) 846 (let ((den-twos (1- (integer-length (logxor den (1- den)))))) 847 (declare (fixnum den-twos)) 848 (decf scale den-twos) 849 (setq den (ash den (- den-twos)))) 850 ;; 851 ;; Guess how much we need to scale by from the magnitudes of the numerator 852 ;; and denominator. We want one extra bit for a guard bit. 853 (let* ((num-len (integer-length num)) 854 (den-len (integer-length den)) 855 (delta (- den-len num-len)) 856 (shift (1+ (the fixnum (+ delta digits)))) 857 (shifted-num (ash num shift))) 858 (declare (fixnum delta shift)) 859 (decf scale delta) 860 (labels ((float-and-scale (bits) 861 (let* ((bits (ash bits -1)) 862 (len (integer-length bits))) 863 (cond ((> len digits) 864 (assert (= len (the fixnum (1+ digits)))) 865 (multiple-value-bind (f0) 866 (floatit (ash bits -1)) 867 (list (first f0) (+ (second f0) 868 (1+ scale))))) 869 (t 870 (multiple-value-bind (f0) 871 (floatit bits) 872 (list (first f0) (+ (second f0) scale))))))) 873 (floatit (bits) 874 (let ((sign (if plusp 1 -1))) 875 (list (* sign bits) 0)))) 876 (loop 877 (multiple-value-bind (fraction-and-guard rem) 878 (truncate shifted-num den) 879 (let ((extra (- (integer-length fraction-and-guard) digits))) 880 (declare (fixnum extra)) 881 (cond ((/= extra 1) 882 (assert (> extra 1))) 883 ((oddp fraction-and-guard) 884 (return 885 (if (zerop rem) 886 (float-and-scale 887 (if (zerop (logand fraction-and-guard 2)) 888 fraction-and-guard 889 (1+ fraction-and-guard))) 890 (float-and-scale (1+ fraction-and-guard))))) 891 (t 892 (return (float-and-scale fraction-and-guard))))) 893 (setq shifted-num (ash shifted-num -1)) 894 (incf scale))))))) 895 896(defun decimalsin (x) 897 (do ((i (quotient (* 59. x) 196.) (1+ i))) ;log[10](2)=.301029 898 (nil) 899 (when (> (integer-length (expt 10. i)) x) 900 (return (1- i))))) 901 902(defun atanbigfloat (x) 903 (*fpatan (car x) (cdr x))) 904 905(defun *fpatan (a y) 906 (fpend (let ((fpprec (+ 8. fpprec))) 907 (if (null y) 908 (if ($bfloatp a) (fpatan (cdr ($bfloat a))) 909 (list '(%atan) a)) 910 (fpatan2 (cdr ($bfloat a)) (cdr ($bfloat (car y)))))))) 911 912;; Bigfloat atan 913(defun fpatan (x) 914 (prog (term x2 ans oans one two tmp) 915 (setq one (intofp 1) two (intofp 2)) 916 (cond ((fpgreaterp (fpabs x) one) 917 ;; |x| > 1. 918 ;; 919 ;; Use A&S 4.4.5: 920 ;; atan(x) + acot(x) = +/- pi/2 (+ for x >= 0, - for x < 0) 921 ;; 922 ;; and A&S 4.4.8 923 ;; acot(z) = atan(1/z) 924 (setq tmp (fpquotient (fppi) two)) 925 (setq ans (fpdifference tmp (fpatan (fpquotient one x)))) 926 (return (cond ((fplessp x (intofp 0)) 927 (fpdifference ans (fppi))) 928 (t ans)))) 929 ((fpgreaterp (fpabs x) (fpquotient one two)) 930 ;; |x| > 1/2 931 ;; 932 ;; Use A&S 4.4.42, third formula: 933 ;; 934 ;; atan(z) = z/(1+z^2)*[1 + 2/3*r + (2*4)/(3*5)*r^2 + ...] 935 ;; 936 ;; r = z^2/(1+z^2) 937 (setq tmp (fpquotient x (fpplus (fptimes* x x) one))) 938 (setq x2 (fptimes* x tmp) term (setq ans one)) 939 (do ((n 0 (1+ n))) 940 ((equal ans oans)) 941 (setq term 942 (fptimes* term (fptimes* x2 (fpquotient 943 (intofp (+ 2 (* 2 n))) 944 (intofp (+ (* 2 n) 3)))))) 945 (setq oans ans ans (fpplus term ans))) 946 (setq ans (fptimes* tmp ans))) 947 (t 948 ;; |x| <= 1/2. Use Taylor series (A&S 4.4.42, first 949 ;; formula). 950 (setq ans x x2 (fpminus (fptimes* x x)) term x) 951 (do ((n 3 (+ n 2))) 952 ((equal ans oans)) 953 (setq term (fptimes* term x2)) 954 (setq oans ans 955 ans (fpplus ans (fpquotient term (intofp n))))))) 956 (return ans))) 957 958;; atan(y/x) taking into account the quadrant. (Also equal to 959;; arg(x+%i*y).) 960(defun fpatan2 (y x) 961 (cond ((equal (car x) 0) 962 ;; atan(y/0) = atan(inf), but what sign? 963 (cond ((equal (car y) 0) 964 (merror (intl:gettext "atan2: atan2(0, 0) is undefined."))) 965 ((minusp (car y)) 966 ;; We're on the negative imaginary axis, so -pi/2. 967 (fpquotient (fppi) (intofp -2))) 968 (t 969 ;; The positive imaginary axis, so +pi/2 970 (fpquotient (fppi) (intofp 2))))) 971 ((signp g (car x)) 972 ;; x > 0. atan(y/x) is the correct value. 973 (fpatan (fpquotient y x))) 974 ((signp g (car y)) 975 ;; x < 0, and y > 0. We're in quadrant II, so the angle we 976 ;; want is pi+atan(y/x). 977 (fpplus (fppi) (fpatan (fpquotient y x)))) 978 (t 979 ;; x <= 0 and y <= 0. We're in quadrant III, so the angle we 980 ;; want is atan(y/x)-pi. 981 (fpdifference (fpatan (fpquotient y x)) (fppi))))) 982 983(defun tanbigfloat (a) 984 (setq a (car a)) 985 (fpend (let ((fpprec (+ 8. fpprec))) 986 (cond (($bfloatp a) 987 (setq a (cdr ($bfloat a))) 988 (fpquotient (fpsin a t) (fpsin a nil))) 989 (t (list '(%tan) a)))))) 990 991;; Returns a list of a mantissa and an exponent. 992(defun intofp (l) 993 (cond ((not (atom l)) ($bfloat l)) 994 ((floatp l) (floattofp l)) 995 ((equal 0 l) '(0 0)) 996 ((eq l '$%pi) (fppi)) 997 ((eq l '$%e) (fpe)) 998 ((eq l '$%gamma) (fpgamma)) 999 (t (list (fpround l) (+ *m fpprec))))) 1000 1001;; It seems to me that this function gets called on an integer 1002;; and returns the mantissa portion of the mantissa/exponent pair. 1003 1004;; "STICKY BIT" CALCULATION FIXED 10/14/75 --RJF 1005;; BASE must not get temporarily bound to NIL by being placed 1006;; in a PROG list as this will confuse stepping programs. 1007 1008(defun fpround (l &aux (*print-base* 10.) *print-radix*) 1009 (prog (adjust) 1010 (cond 1011 ((null *decfp) 1012 ;;*M will be positive if the precision of the argument is greater than 1013 ;;the current precision being used. 1014 (setq *m (- (integer-length l) fpprec)) 1015 (when (= *m 0) 1016 (setq *cancelled 0) 1017 (return l)) 1018 ;;FPSHIFT is essentially LSH. 1019 (setq adjust (fpshift 1 (1- *m))) 1020 (when (minusp l) (setq adjust (- adjust))) 1021 (incf l adjust) 1022 (setq *m (- (integer-length l) fpprec)) 1023 (setq *cancelled (abs *m)) 1024 (cond ((zerop (hipart l (- *m))) 1025 ;ONLY ZEROES SHIFTED OFF 1026 (return (fpshift (fpshift l (- -1 *m)) 1027 1))) ; ROUND TO MAKE EVEN 1028 (t (return (fpshift l (- *m)))))) 1029 (t 1030 (setq *m (- (flatsize (abs l)) fpprec)) 1031 (setq adjust (fpshift 1 (1- *m))) 1032 (when (minusp l) (setq adjust (- adjust))) 1033 (setq adjust (* 5 adjust)) 1034 (setq *m (- (flatsize (abs (setq l (+ l adjust)))) fpprec)) 1035 (return (fpshift l (- *m))))))) 1036 1037;; Compute (* L (expt d n)) where D is 2 or 10 depending on 1038;; *decfp. Throw away an fractional part by truncating to zero. 1039(defun fpshift (l n) 1040 (cond ((null *decfp) 1041 (cond ((and (minusp n) (minusp l)) 1042 ;; Left shift of negative number requires some 1043 ;; care. (That is, (truncate l (expt 2 n)), but use 1044 ;; shifts instead.) 1045 (- (ash (- l) n))) 1046 (t 1047 (ash l n)))) 1048 ((> n 0) 1049 (* l (expt 10. n))) 1050 ((< n 0.) 1051 (quotient l (expt 10. (- n)))) 1052 (t l))) 1053 1054;; Bignum LSH -- N is assumed (and declared above) to be a fixnum. 1055;; This isn't really LSH, since the sign bit isn't propagated when 1056;; shifting to the right, i.e. (BIGLSH -100 -3) = -40, whereas 1057;; (LSH -100 -3) = 777777777770 (on a 36 bit machine). 1058;; This actually computes (* X (EXPT 2 N)). As of 12/21/80, this function 1059;; was only called by FPSHIFT. I would like to hear an argument as why this 1060;; is more efficient than simply writing (* X (EXPT 2 N)). Is the 1061;; intermediate result created by (EXPT 2 N) the problem? I assume that 1062;; EXPT tries to LSH when possible. 1063 1064(defun biglsh (x n) 1065 (cond ((and (not (bignump x)) 1066 (< n #.(- +machine-fixnum-precision+))) 1067 0) 1068 ;; Either we are shifting a fixnum to the right, or shifting 1069 ;; a fixnum to the left, but not far enough left for it to become 1070 ;; a bignum. 1071 ((and (not (bignump x)) 1072 (or (<= n 0) 1073 (< (+ (integer-length x) n) #.+machine-fixnum-precision+))) 1074 ;; The form which follows is nearly identical to (ASH X N), however 1075 ;; (ASH -100 -20) = -1, whereas (BIGLSH -100 -20) = 0. 1076 (if (>= x 0) 1077 (ash x n) 1078 (- (biglsh (- x) n)))) ;(- x) may be a bignum even is x is a fixnum. 1079 ;; If we get here, then either X is a bignum or our answer is 1080 ;; going to be a bignum. 1081 ((< n 0) 1082 (cond ((> (abs n) (integer-length x)) 0) 1083 ((> x 0) 1084 (hipart x (+ (integer-length x) n))) 1085 (t (- (hipart x (+ (integer-length x) n)))))) 1086 ((= n 0) x) 1087 ;; Isn't this the kind of optimization that compilers are 1088 ;; supposed to make? 1089 ((< n #.(1- +machine-fixnum-precision+)) (* x (ash 1 n))) 1090 (t (* x (expt 2 n))))) 1091 1092 1093;; exp(x) 1094;; 1095;; For negative x, use exp(-x) = 1/exp(x) 1096;; 1097;; For x > 0, exp(x) = exp(r+y) = exp(r) * exp(y), where x = r + y and 1098;; r = floor(x). 1099(defun fpexp (x) 1100 (prog (r s) 1101 (unless (signp ge (car x)) 1102 (return (fpquotient (fpone) (fpexp (fpabs x))))) 1103 (setq r (fpintpart x :skip-exponent-check-p t)) 1104 (return (cond ((< r 2) 1105 (fpexp1 x)) 1106 (t 1107 (setq s (fpexp1 (fpdifference x (intofp r)))) 1108 (fptimes* s 1109 (cdr (bigfloatp 1110 (let ((fpprec (+ fpprec (integer-length r) -1)) 1111 (r r)) 1112 (bcons (fpexpt (fpe) r))))))))))) ; patch for full precision %E 1113 1114;; exp(x) for small x, using Taylor series. 1115(defun fpexp1 (x) 1116 (prog (term ans oans) 1117 (setq ans (setq term (fpone))) 1118 (do ((n 1 (1+ n))) 1119 ((equal ans oans)) 1120 (setq term (fpquotient (fptimes* x term) (intofp n))) 1121 (setq oans ans) 1122 (setq ans (fpplus ans term))) 1123 (return ans))) 1124 1125;; Does one higher precision to round correctly. 1126;; A and B are each a list of a mantissa and an exponent. 1127(defun fpquotient (a b) 1128 (cond ((equal (car b) 0) 1129 (merror (intl:gettext "pquotient: attempted quotient by zero."))) 1130 ((equal (car a) 0) '(0 0)) 1131 (t (list (fpround (quotient (fpshift (car a) (+ 3 fpprec)) (car b))) 1132 (+ -3 (- (cadr a) (cadr b)) *m))))) 1133 1134(defun fpgreaterp (a b) 1135 (fpposp (fpdifference a b))) 1136 1137(defun fplessp (a b) 1138 (fpposp (fpdifference b a))) 1139 1140(defun fpposp (x) 1141 (> (car x) 0)) 1142 1143(defun fpmin (arg1 &rest args) 1144 (let ((min arg1)) 1145 (mapc #'(lambda (u) (if (fplessp u min) (setq min u))) args) 1146 min)) 1147 1148(defun fpmax (arg1 &rest args) 1149 (let ((max arg1)) 1150 (mapc #'(lambda (u) (if (fpgreaterp u max) (setq max u))) args) 1151 max)) 1152 1153;; The following functions compute bigfloat values for %e, %pi, 1154;; %gamma, and log(2). For each precision, the computed value is 1155;; cached in a hash table so it doesn't need to be computed again. 1156;; There are functions to return the hash table or clear the hash 1157;; table, for debugging. 1158;; 1159;; Note that each of these return a bigfloat number, but without the 1160;; bigfloat tag. 1161;; 1162;; See 1163;; https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2910437&group_id=4933 1164;; for an explanation. 1165(let ((table (make-hash-table))) 1166 (defun fpe () 1167 (let ((value (gethash fpprec table))) 1168 (if value 1169 value 1170 (setf (gethash fpprec table) (cdr (fpe1)))))) 1171 (defun fpe-table () 1172 table) 1173 (defun clear_fpe_table () 1174 (clrhash table))) 1175 1176(let ((table (make-hash-table))) 1177 (defun fppi () 1178 (let ((value (gethash fpprec table))) 1179 (if value 1180 value 1181 (setf (gethash fpprec table) (cdr (fppi1)))))) 1182 (defun fppi-table () 1183 table) 1184 (defun clear_fppi_table () 1185 (clrhash table))) 1186 1187(let ((table (make-hash-table))) 1188 (defun fpgamma () 1189 (let ((value (gethash fpprec table))) 1190 (if value 1191 value 1192 (setf (gethash fpprec table) (cdr (fpgamma1)))))) 1193 (defun fpgamma-table () 1194 table) 1195 (defun clear_fpgamma_table () 1196 (clrhash table))) 1197 1198(let ((table (make-hash-table))) 1199 (defun fplog2 () 1200 (let ((value (gethash fpprec table))) 1201 (if value 1202 value 1203 (setf (gethash fpprec table) (comp-log2))))) 1204 (defun fplog2-table () 1205 table) 1206 (defun clear_fplog2_table () 1207 (clrhash table))) 1208 1209;; This doesn't need a hash table because there's never a problem with 1210;; using a high precision value and rounding to a lower precision 1211;; value because 1 is always an exact bfloat. 1212(defun fpone () 1213 (cond (*decfp (intofp 1)) 1214 ((= fpprec (bigfloat-prec bigfloatone)) (cdr bigfloatone)) 1215 (t (intofp 1)))) 1216 1217;;----------------------------------------------------------------------------;; 1218;; 1219;; The values of %e, %pi, %gamma and log(2) are computed by the technique of 1220;; binary splitting. See http://www.ginac.de/CLN/binsplit.pdf for details. 1221;; 1222;; Volker van Nek, Sept. 2014 1223 1224;; 1225;; Euler's number E 1226;; 1227(defun fpe1 () 1228 (let ((e (compe (+ fpprec 12)))) ;; compute additional bits 1229 (bcons (list (fpround (car e)) (cadr e))) )) ;; round to fpprec 1230;; 1231;; Taylor: %e = sum(s[i] ,i,0,inf) where s[i] = 1/i! 1232;; 1233(defun compe (prec) 1234 (let ((fpprec prec)) 1235 (multiple-value-bind (tt qq) (split-taylor-e 0 (taylor-e-size prec)) 1236 (fpquotient (intofp tt) (intofp qq)) ))) 1237;; 1238;; binary splitting: 1239;; 1240;; 1 1241;; s[i] = ---------------------- 1242;; q[0]*q[1]*q[2]*..*q[i] 1243;; 1244;; where q[0] = 1 1245;; q[i] = i 1246;; 1247(defun split-taylor-e (i j) 1248 (let (qq tt) 1249 (if (= (- j i) 1) 1250 (setq qq (if (= i 0) 1 i) 1251 tt 1 ) 1252 (let ((m (ash (+ i j) -1))) 1253 (multiple-value-bind (tl ql) (split-taylor-e i m) 1254 (multiple-value-bind (tr qr) (split-taylor-e m j) 1255 (setq qq (* ql qr) 1256 tt (+ (* qr tl) tr) ))))) 1257 (values tt qq) )) 1258;; 1259;; stop when i! > 2^fpprec 1260;; 1261;; log(i!) = sum(log(k), k,1,i) > fpprec * log(2) 1262;; 1263(defun taylor-e-size (prec) 1264 (let ((acc 0) 1265 (lim (* prec (log 2))) ) 1266 (do ((i 1 (1+ i))) 1267 ((> acc lim) i) 1268 (incf acc (log i)) ))) 1269;; 1270;;----------------------------------------------------------------------------;; 1271;; 1272;; PI 1273;; 1274(defun fppi1 () 1275 (let ((pi1 (comppi (+ fpprec 10)))) 1276 (bcons (list (fpround (car pi1)) (cadr pi1))) )) 1277;; 1278;; Chudnovsky & Chudnovsky: 1279;; 1280;; C^(3/2)/(12*%pi) = sum(s[i], i,0,inf), 1281;; 1282;; where s[i] = (-1)^i*(6*i)!*(A*i+B) / (i!^3*(3*i)!*C^(3*i)) 1283;; 1284;; and A = 545140134, B = 13591409, C = 640320 1285;; 1286(defun comppi (prec) 1287 (let ((fpprec prec) 1288 nr n d oldn tt qq n*qq ) 1289 ;; STEP 1: 1290 ;; compute n/d = sqrt(10005) : 1291 ;; 1292 ;; n[0] n[i+1] = n[i]^2+a*d[i]^2 n[inf] 1293 ;; quadratic Heron: x[0] = ----, , sqrt(a) = ------ 1294 ;; d[0] d[i+1] = 2*n[i]*d[i] d[inf] 1295 ;; 1296 (multiple-value-setq (nr n d) (sqrt-10005-constants fpprec)) 1297 (dotimes (i nr) 1298 (setq oldn n 1299 n (+ (* n n) (* 10005 d d)) 1300 d (* 2 oldn d) )) 1301 ;; STEP 2: 1302 ;; divide C^(3/2)/12 = 3335*2^7*sqrt(10005) 1303 ;; by Chudnovsky-sum = tt/qq : 1304 ;; 1305 (setq nr (ceiling (* fpprec 0.021226729578153))) ;; nr of summands 1306 ;; fpprec*log(2)/log(C^3/(24*6*2*6)) 1307 (multiple-value-setq (tt qq) (split-chudnovsky 0 (1+ nr))) 1308 (setq n (* 3335 n) 1309 n*qq (intofp (* n qq)) ) 1310 (fpquotient (list (car n*qq) (+ (cadr n*qq) 7)) 1311 (intofp (* d tt)) ))) 1312;; 1313;; The returned n and d serve as start values for the iteration. 1314;; n/d = sqrt(10005) with a precision of p = ceiling(prec/2^nr) bits 1315;; where nr is the number of needed iterations. 1316;; 1317(defun sqrt-10005-constants (prec) 1318 (let (ilen p nr n d) 1319 (if (< prec 128) 1320 (setq nr 0 p prec) 1321 (setq ilen (integer-length prec) 1322 nr (- ilen 7) 1323 p (ceiling (* prec (expt 2.0 (- nr)))) )) 1324 (cond 1325 ((<= p 76) (setq n 256192036001 d 2561280120)) 1326 ((<= p 89) (setq n 51244811200700 d 512320048001)) 1327 ((<= p 102) (setq n 2050048640064001 d 20495363200160)) 1328 ((<= p 115) (setq n 410060972824000900 d 4099584960080001)) 1329 (t (setq n 16404488961600100001 d 164003893766400200)) ) 1330 (values nr n d) )) 1331;; 1332;; binary splitting: 1333;; 1334;; a[i] * p[0]*p[1]*p[2]*..*p[i] 1335;; s[i] = ----------------------------- 1336;; q[0]*q[1]*q[2]*..*q[i] 1337;; 1338;; where a[0] = B 1339;; p[0] = q[0] = 1 1340;; a[i] = A*i+B 1341;; p[i] = - (6*i-5)*(2*i-1)*(6*i-1) 1342;; q[i] = C^3/24*i^3 1343;; 1344(defun split-chudnovsky (i j) 1345 (let (aa pp/qq pp qq tt) 1346 (if (= (- j i) 1) 1347 (if (= i 0) 1348 (setq aa 13591409 pp 1 qq 1 tt aa) 1349 (setq aa (+ (* i 545140134) 13591409) 1350 pp/qq (/ (* (- 5 (* 6 i)) (- (* 2 i) 1) (- (* 6 i) 1)) 1351 10939058860032000 ) ; C^3/24 1352 pp (numerator pp/qq) 1353 qq (* (denominator pp/qq) (expt i 3)) 1354 tt (* aa pp) )) 1355 (let ((m (ash (+ i j) -1))) 1356 (multiple-value-bind (tl ql pl) (split-chudnovsky i m) 1357 (multiple-value-bind (tr qr pr) (split-chudnovsky m j) 1358 (setq pp (* pl pr) 1359 qq (* ql qr) 1360 tt (+ (* qr tl) (* pl tr)) ))))) 1361 (values tt qq pp) )) 1362;; 1363;;----------------------------------------------------------------------------;; 1364;; 1365;; Euler-Mascheroni constant GAMMA 1366;; 1367(defun fpgamma1 () 1368 (let ((res (comp-bf%gamma (+ fpprec 14)))) 1369 (bcons (list (fpround (car res)) (cadr res))) )) 1370;; 1371;; Brent-McMillan algorithm 1372;; 1373;; Let 1374;; alpha = 4.970625759544 1375;; 1376;; n > 0 and N-1 >= alpha*n 1377;; 1378;; H(k) = sum(1/i, i,1,k) 1379;; 1380;; S = sum(H(k)*(n^k/k!)^2, k,0,N-1) 1381;; 1382;; I = sum((n^k/k!)^2, k,0,N-1) 1383;; 1384;; T = 1/(4*n)*sum((2*k)!^3/(k!^4*(16*n)^(2*k)), k,0,2*n-1) 1385;; 1386;; and 1387;; %gamma = S/I - T/I^2 - log(n) 1388;; 1389;; Then 1390;; |%gamma - gamma| < 24 * e^(-8*n) 1391;; 1392;; (Corollary 2, Remark 2, Brent/Johansson http://arxiv.org/pdf/1312.0039v1.pdf) 1393;; 1394(defun comp-bf%gamma (prec) 1395 (let* ((fpprec prec) 1396 (n (ceiling (* 1/8 (+ (* prec (log 2.0)) (log 24.0))))) 1397 (n2 (* n n)) 1398 (alpha 4.970625759544) 1399 (lim (ceiling (* alpha n))) 1400 sums/sumi ;; S/I 1401 sumi sumi2 ;; I and I^2 1402 sumt/sumi2 ) ;; T/I^2 1403 (multiple-value-bind (vv tt qq dd) (split-gamma-1 1 (1+ lim) n2) 1404 ;; 1405 ;; sums = vv/(qq*dd) 1406 ;; sumi = tt/qq 1407 ;; sums/sumi = vv/(qq*dd)*qq/tt = vv/(dd*tt) 1408 ;; 1409 (setq sums/sumi (fpquotient (intofp vv) (intofp (* dd tt))) 1410 sumi (fpquotient (intofp tt) (intofp qq)) 1411 sumi2 (fptimes* sumi sumi) ) 1412 ;; 1413 (multiple-value-bind (ttt qqq) (split-gamma-2 0 (* 2 n) (* 32 n2)) 1414 ;; 1415 ;; sumt = 1/(4*n)*ttt/qqq 1416 ;; sumt/sumi2 = ttt/(4*n*qqq*sumi2) 1417 ;; 1418 (setq sumt/sumi2 (fpquotient (intofp ttt) 1419 (fptimes* (intofp (* 4 n qqq)) sumi2) )) 1420 ;; %gamma : 1421 (fpdifference sums/sumi (fpplus sumt/sumi2 (log-n n)) ))))) 1422;; 1423;; split S and I simultaneously: 1424;; 1425;; summands I[0] = 1, I[i]/I[i-1] = n^2/i^2 1426;; 1427;; S[0] = 0, S[i]/S[i-1] = n^2/i^2*H(i)/H(i-1) 1428;; 1429;; p[0]*p[1]*p[2]*..*p[i] 1430;; I[i] = ---------------------- 1431;; q[0]*q[1]*q[2]*..*q[i] 1432;; 1433;; where p[0] = n^2 1434;; q[0] = 1 1435;; p[i] = n^2 1436;; q[i] = i^2 1437;; c[0] c[1] c[2] c[i] 1438;; S[i] = H[i] * I[i], where H[i] = ---- + ---- + ---- + .. + ---- 1439;; d[0] d[1] d[2] d[i] 1440;; and c[0] = 0 1441;; d[0] = 1 1442;; c[i] = 1 1443;; d[i] = i 1444;; 1445(defun split-gamma-1 (i j n2) 1446 (let (pp cc dd qq tt vv) 1447 (cond 1448 ((= (- j i) 1) 1449 (if (= i 1) ;; S[0] is 0 -> start with i=1 and add I[0]=1 to tt : 1450 (setq pp n2 cc 1 dd 1 qq 1 tt (1+ n2) vv n2) 1451 (setq pp n2 cc 1 dd i qq (* i i) tt pp vv tt) )) 1452 (t 1453 (let* ((m (ash (+ i j) -1)) tmp) 1454 (multiple-value-bind (vl tl ql dl cl pl) (split-gamma-1 i m n2) 1455 (multiple-value-bind (vr tr qr dr cr pr) (split-gamma-1 m j n2) 1456 (setq pp (* pl pr) 1457 cc (+ (* cl dr) (* dl cr)) 1458 dd (* dl dr) 1459 qq (* ql qr) 1460 tmp (* pl tr) 1461 tt (+ (* tl qr) tmp) 1462 vv (+ (* dr (+ (* vl qr) (* cl tmp))) (* dl pl vr)) )))))) 1463 (values vv tt qq dd cc pp) )) 1464;; 1465;; split 4*n*T: 1466;; 1467;; summands T[0] = 1, T[i]/T[i-1] = (2*i-1)^3/(32*i*n^2) 1468;; 1469;; p[0]*p[1]*p[2]*..*p[i] 1470;; T[i] = ---------------------- 1471;; q[0]*q[1]*q[2]*..*q[i] 1472;; 1473;; where p[0] = q[0] = 1 1474;; p[i] = (2*i-1)^3 1475;; q[i] = 32*i*n^2 1476;; 1477(defun split-gamma-2 (i j n2*32) 1478 (let (pp qq tt) 1479 (cond 1480 ((= (- j i) 1) 1481 (if (= i 0) 1482 (setq pp 1 qq 1 tt 1) 1483 (setq pp (expt (1- (* 2 i)) 3) qq (* i n2*32) tt pp) )) 1484 (t 1485 (let* ((m (ash (+ i j) -1))) 1486 (multiple-value-bind (tl ql pl) (split-gamma-2 i m n2*32) 1487 (multiple-value-bind (tr qr pr) (split-gamma-2 m j n2*32) 1488 (setq pp (* pl pr) 1489 qq (* ql qr) 1490 tt (+ (* tl qr) (* pl tr)) )))))) 1491 (values tt qq pp) )) 1492;; 1493;;----------------------------------------------------------------------------;; 1494;; 1495;; log(2) = 18*L(26) - 2*L(4801) + 8*L(8749) 1496;; 1497;; where L(k) = atanh(1/k) 1498;; 1499;; see http://numbers.computation.free.fr/Constants/constants.html 1500;; 1501;;;(defun $log2 () (bcons (comp-log2))) ;; checked against reference table 1502;; 1503(defun comp-log2 () 1504 (let ((res 1505 (let ((fpprec (+ fpprec 12))) 1506 (fpplus 1507 (fpdifference (n*atanh-1/k 18 26) (n*atanh-1/k 2 4801)) 1508 (n*atanh-1/k 8 8749) )))) 1509 (list (fpround (car res)) (cadr res)) )) 1510;; 1511;; Taylor: atanh(1/k) = sum(s[i], i,0,inf) 1512;; 1513;; where s[i] = 1/((2*i+1)*k^(2*i+1)) 1514;; 1515(defun n*atanh-1/k (n k) ;; integer n,k 1516 (let* ((k2 (* k k)) 1517 (nr (ceiling (* fpprec (/ (log 2) (log k2))))) ) 1518 (multiple-value-bind (tt qq bb) (split-atanh-1/k 0 (1+ nr) k k2) 1519 (fpquotient (intofp (* n tt)) (intofp (* bb qq))) ))) 1520;; 1521;; binary splitting: 1522;; 1 1523;; s[i] = ----------------------------- 1524;; b[i] * q[0]*q[1]*q[2]*..*q[i] 1525;; 1526;; where b[0] = 1 1527;; q[0] = k 1528;; b[i] = 2*i+1 1529;; q[i] = k^2 1530;; 1531(defun split-atanh-1/k (i j k k2) 1532 (let (bb qq tt) 1533 (if (= (- j i) 1) 1534 (if (= i 0) 1535 (setq bb 1 qq k tt 1) 1536 (setq bb (1+ (* 2 i)) qq k2 tt 1) ) 1537 (let ((m (ash (+ i j) -1))) 1538 (multiple-value-bind (tl ql bl) (split-atanh-1/k i m k k2) 1539 (multiple-value-bind (tr qr br) (split-atanh-1/k m j k k2) 1540 (setq bb (* bl br) 1541 qq (* ql qr) 1542 tt (+ (* br qr tl) (* bl tr)) ))))) 1543 (values tt qq bb) )) 1544;; 1545;;----------------------------------------------------------------------------;; 1546;; 1547;; log(n) = log(n/2^k) + k*log(2) 1548;; 1549;;;(defun $log10 () (bcons (log-n 10))) ;; checked against reference table 1550;; 1551(defun log-n (n) ;; integer n > 0 1552 (cond 1553 ((= 1 n) (list 0 0)) 1554 ((= 2 n) (comp-log2)) 1555 (t 1556 (let ((res 1557 (let ((fpprec (+ fpprec 10)) 1558 (k (integer-length n)) ) 1559 ;; choose k so that |n/2^k - 1| is as small as possible: 1560 (when (< n (* (coerce 2/3 'flonum) (ash 1 k))) (decf k)) 1561 ;; now |n/2^k - 1| <= 1/3 1562 (fpplus (log-u/2^k n k fpprec) 1563 (fptimes* (intofp k) (comp-log2)) )))) 1564 (list (fpround (car res)) (cadr res)) )))) 1565;; 1566;; log(1+u/v) = 2 * sum(s[i], i,0,inf) 1567;; 1568;; where s[i] = (u/(2*v+u))^(2*i+1)/(2*i+1) 1569;; 1570(defun log-u/2^k (u k prec) ;; integer u k; x = u/2^k; |x - 1| < 1 1571 (setq u (- u (ash 1 k))) ;; x <-- x - 1 1572 (cond 1573 ((= 0 u) (list 0 0)) 1574 (t 1575 (while (evenp u) (setq u (ash u -1)) (decf k)) 1576 (let* ((u2 (* u u)) 1577 (w (+ u (ash 2 k))) 1578 (w2 (* w w)) 1579 (nr (ceiling (* prec (/ (log 2) 2 (log (abs (/ w u))))))) 1580 lg/2 ) 1581 (multiple-value-bind (tt qq bb) (split-log-1+u/v 0 (1+ nr) u u2 w w2) 1582 (setq lg/2 (fpquotient (intofp tt) (intofp (* bb qq)))) ;; sum 1583 (list (car lg/2) (1+ (cadr lg/2))) ))))) ;; 2*sum 1584;; 1585;; binary splitting: 1586;; 1587;; p[0]*p[1]*p[2]*..*p[i] 1588;; s[i] = ----------------------------- 1589;; b[i] * q[0]*q[1]*q[2]*..*q[i] 1590;; 1591;; where b[0] = 1 1592;; p[0] = u 1593;; q[0] = w = 2*v+u 1594;; b[i] = 2*i+1 1595;; p[i] = u^2 1596;; q[i] = w^2 1597;; 1598(defun split-log-1+u/v (i j u u2 w w2) 1599 (let (pp bb qq tt) 1600 (if (= (- j i) 1) 1601 (if (= i 0) 1602 (setq pp u bb 1 qq w tt u) 1603 (setq pp u2 bb (1+ (* 2 i)) qq w2 tt pp) ) 1604 (let ((m (ash (+ i j) -1))) 1605 (multiple-value-bind (tl ql bl pl) (split-log-1+u/v i m u u2 w w2) 1606 (multiple-value-bind (tr qr br pr) (split-log-1+u/v m j u u2 w w2) 1607 (setq bb (* bl br) 1608 pp (* pl pr) 1609 qq (* ql qr) 1610 tt (+ (* br qr tl) (* bl pl tr)) ))))) 1611 (values tt qq bb pp) )) 1612;; 1613;;----------------------------------------------------------------------------;; 1614 1615 1616(defun fpdifference (a b) 1617 (fpplus a (fpminus b))) 1618 1619(defun fpminus (x) 1620 (if (equal (car x) 0) 1621 x 1622 (list (- (car x)) (cadr x)))) 1623 1624(defun fpplus (a b) 1625 (prog (*m exp man sticky) 1626 (setq *cancelled 0) 1627 (cond ((equal (car a) 0) (return b)) 1628 ((equal (car b) 0) (return a))) 1629 (setq exp (- (cadr a) (cadr b))) 1630 (setq man (cond ((equal exp 0) 1631 (setq sticky 0) 1632 (fpshift (+ (car a) (car b)) 2)) 1633 ((> exp 0) 1634 (setq sticky (hipart (car b) (- 1 exp))) 1635 (setq sticky (cond ((signp e sticky) 0) 1636 ((signp l (car b)) -1) 1637 (t 1))) 1638 ; COMPUTE STICKY BIT 1639 (+ (fpshift (car a) 2) 1640 ; MAKE ROOM FOR GUARD DIGIT & STICKY BIT 1641 (fpshift (car b) (- 2 exp)))) 1642 (t (setq sticky (hipart (car a) (1+ exp))) 1643 (setq sticky (cond ((signp e sticky) 0) 1644 ((signp l (car a)) -1) 1645 (t 1))) 1646 (+ (fpshift (car b) 2) 1647 (fpshift (car a) (+ 2 exp)))))) 1648 (setq man (+ man sticky)) 1649 (return (cond ((equal man 0) '(0 0)) 1650 (t (setq man (fpround man)) 1651 (setq exp (+ -2 *m (max (cadr a) (cadr b)))) 1652 (list man exp)))))) 1653 1654(defun fptimes* (a b) 1655 (if (or (zerop (car a)) (zerop (car b))) 1656 '(0 0) 1657 (list (fpround (* (car a) (car b))) 1658 (+ *m (cadr a) (cadr b) (- fpprec))))) 1659 1660;; Don't use the symbol BASE since it is SPECIAL. 1661 1662(defun fpintexpt (int nn fixprec) ;INT is integer 1663 (setq fixprec (truncate fixprec (1- (integer-length int)))) ;NN is pos 1664 (let ((bas (intofp (expt int (min nn fixprec))))) 1665 (if (> nn fixprec) 1666 (fptimes* (intofp (expt int (rem nn fixprec))) 1667 (fpexpt bas (quotient nn fixprec))) 1668 bas))) 1669 1670;; NN is positive or negative integer 1671 1672(defun fpexpt (p nn) 1673 (cond ((zerop nn) (fpone)) 1674 ((eql nn 1) p) 1675 ((< nn 0) (fpquotient (fpone) (fpexpt p (- nn)))) 1676 (t (prog (u) 1677 (if (oddp nn) 1678 (setq u p) 1679 (setq u (fpone))) 1680 (do ((ii (quotient nn 2) (quotient ii 2))) 1681 ((zerop ii)) 1682 (setq p (fptimes* p p)) 1683 (when (oddp ii) 1684 (setq u (fptimes* u p)))) 1685 (return u))))) 1686 1687(defun exptbigfloat (p n) 1688 (cond ((equal n 1) p) 1689 ((equal n 0) ($bfloat 1)) 1690 ((not ($bfloatp p)) (list '(mexpt) p n)) 1691 ((equal (cadr p) 0) ($bfloat 0)) 1692 ((and (< (cadr p) 0) (ratnump n)) 1693 (mul2 (let ($numer $float $keepfloat $ratprint) 1694 (power -1 n)) 1695 (exptbigfloat (bcons (fpminus (cdr p))) n))) 1696 ((and (< (cadr p) 0) (not (integerp n))) 1697 (cond ((or (equal n 0.5) (equal n bfhalf)) 1698 (exptbigfloat p '((rat simp) 1 2))) 1699 ((or (equal n -0.5) (equal n bfmhalf)) 1700 (exptbigfloat p '((rat simp) -1 2))) 1701 (($bfloatp (setq n ($bfloat n))) 1702 (cond ((equal n ($bfloat (fpentier n))) 1703 (exptbigfloat p (fpentier n))) 1704 (t ;; for P<0: P^N = (-P)^N*cos(pi*N) + i*(-P)^N*sin(pi*N) 1705 (setq p (exptbigfloat (bcons (fpminus (cdr p))) n) 1706 n ($bfloat `((mtimes) $%pi ,n))) 1707 (add2 ($bfloat `((mtimes) ,p ,(*fpsin n nil))) 1708 `((mtimes simp) ,($bfloat `((mtimes) ,p ,(*fpsin n t))) 1709 $%i))))) 1710 (t (list '(mexpt) p n)))) 1711 ((and (ratnump n) (< (caddr n) 10.)) 1712 (bcons (fpexpt (fproot p (caddr n)) (cadr n)))) 1713 ((not (integerp n)) 1714 (setq n ($bfloat n)) 1715 (cond 1716 ((not ($bfloatp n)) (list '(mexpt) p n)) 1717 (t 1718 (let ((extrabits (max 1 (+ (caddr n) (integer-length (caddr p)))))) 1719 (setq p 1720 (let ((fpprec (+ extrabits fpprec))) 1721 (fpexp (fptimes* (cdr (bigfloatp n)) (fplog (cdr (bigfloatp p))))))) 1722 (setq p (list (fpround (car p)) (+ (- extrabits) *m (cadr p)))) 1723 (bcons p))))) 1724 ;; The number of extra bits required 1725 ((< n 0) (invertbigfloat (exptbigfloat p (- n)))) 1726 (t (bcons (fpexpt (cdr p) n))))) 1727 1728(defun fproot (a n) ; computes a^(1/n) see Fitch, SIGSAM Bull Nov 74 1729 1730 ;; Special case for a = 0b0. General algorithm loops endlessly in that case. 1731 1732 ;; Unlike many or maybe all of the other functions named FP-something, 1733 ;; FPROOT assumes it is called with an argument like 1734 ;; '((BIGFLOAT ...) FOO BAR) instead of '(FOO BAR). 1735 ;; However FPROOT does return something like '(FOO BAR). 1736 1737 (if (eql (cadr a) 0) 1738 '(0 0) 1739 (progn 1740 (let* ((ofprec fpprec) 1741 (fpprec (+ fpprec 2)) ;assumes a>0 n>=2 1742 (bk (fpexpt (intofp 2) (1+ (quotient (cadr (setq a (cdr (bigfloatp a)))) n))))) 1743 (do ((x bk (fpdifference x 1744 (setq bk (fpquotient (fpdifference 1745 x (fpquotient a (fpexpt x n1))) n)))) 1746 (n1 (1- n)) 1747 (n (intofp n))) 1748 ((or (equal bk '(0 0)) 1749 (> (- (cadr x) (cadr bk)) ofprec)) 1750 (setq a x)))) 1751 (list (fpround (car a)) (+ -2 *m (cadr a)))))) 1752 1753(defun timesbigfloat (h) 1754 (prog (fans r nfans) 1755 (setq fans (bcons (fpone)) nfans 1) 1756 (do ((l h (cdr l))) 1757 ((null l)) 1758 (if (setq r (bigfloatp (car l))) 1759 (setq fans (bcons (fptimes* (cdr r) (cdr fans)))) 1760 (setq nfans (list '(mtimes) (car l) nfans)))) 1761 (return (if (equal nfans 1) 1762 fans 1763 (simplify (list '(mtimes) fans nfans)))))) 1764 1765(defun invertbigfloat (a) 1766 ;; If A is a bigfloat, be sure to round it to the current precision. 1767 ;; (See Bug 2543079 for one of the symptoms.) 1768 (let ((b (bigfloatp a))) 1769 (if b 1770 (bcons (fpquotient (fpone) (cdr b))) 1771 (simplify (list '(mexpt) a -1))))) 1772 1773(defun *fpexp (a) 1774 (fpend (let ((fpprec (+ 8. fpprec))) 1775 (if ($bfloatp a) 1776 (fpexp (cdr (bigfloatp a))) 1777 (list '(mexpt) '$%e a))))) 1778 1779(defun *fpsin (a fl) 1780 (fpend (let ((fpprec (+ 8. fpprec))) 1781 (cond (($bfloatp a) (fpsin (cdr ($bfloat a)) fl)) 1782 (fl (list '(%sin) a)) 1783 (t (list '(%cos) a)))))) 1784 1785(defun fpend (a) 1786 (cond ((equal (car a) 0) (bcons a)) 1787 ((numberp (car a)) 1788 (setq a (list (fpround (car a)) (+ -8. *m (cadr a)))) 1789 (bcons a)) 1790 (t a))) 1791 1792(defun fparcsimp (e) ; needed for e.g. ASIN(.123567812345678B0) with 1793 ;; FPPREC 16, to get rid of the miniscule imaginary 1794 ;; part of the a+bi answer. 1795 (if (and (mplusp e) (null (cdddr e)) 1796 (mtimesp (caddr e)) (null (cdddr (caddr e))) 1797 ($bfloatp (cadr (caddr e))) 1798 (eq (caddr (caddr e)) '$%i) 1799 (< (caddr (cadr (caddr e))) (+ (- fpprec) 2))) 1800 (cadr e) 1801 e)) 1802 1803(defun sinbigfloat (x) 1804 (*fpsin (car x) t)) 1805 1806(defun cosbigfloat (x) 1807 (*fpsin (car x) nil)) 1808 1809;; THIS VERSION OF FPSIN COMPUTES SIN OR COS TO PRECISION FPPREC, 1810;; BUT CHECKS FOR THE POSSIBILITY OF CATASTROPHIC CANCELLATION DURING 1811;; ARGUMENT REDUCTION (E.G. SIN(N*%PI+EPSILON)) 1812;; *FPSINCHECK* WILL CAUSE PRINTOUT OF ADDITIONAL INFO WHEN 1813;; EXTRA PRECISION IS NEEDED FOR SIN/COS CALCULATION. KNOWN 1814;; BAD FEATURES: IT IS NOT NECESSARY TO USE EXTRA PRECISION FOR, E.G. 1815;; SIN(PI/2), WHICH IS NOT NEAR ZERO, BUT EXTRA 1816;; PRECISION IS USED SINCE IT IS NEEDED FOR COS(PI/2). 1817;; PRECISION SEEMS TO BE 100% SATSIFACTORY FOR LARGE ARGUMENTS, E.G. 1818;; SIN(31415926.0B0), BUT LESS SO FOR SIN(3.1415926B0). EXPLANATION 1819;; NOT KNOWN. (9/12/75 RJF) 1820 1821(defvar *fpsincheck* nil) 1822 1823;; FL is a T for sin and NIL for cos. 1824(defun fpsin (x fl) 1825 (prog (piby2 r sign res k *cancelled) 1826 (setq sign (cond (fl (signp g (car x))) 1827 (t)) 1828 x (fpabs x)) 1829 (when (equal (car x) 0) 1830 (return (if fl (intofp 0) (intofp 1)))) 1831 (return 1832 (cdr 1833 (bigfloatp 1834 (let ((fpprec (max fpprec (+ fpprec (cadr x)))) 1835 (xt (bcons x)) 1836 (*cancelled 0) 1837 (oldprec fpprec)) 1838 (prog (x) 1839 loop (setq x (cdr (bigfloatp xt))) 1840 (setq piby2 (fpquotient (fppi) (intofp 2))) 1841 (setq r (fpintpart (fpquotient x piby2) :skip-exponent-check-p t)) 1842 (setq x (fpplus x (fptimes* (intofp (- r)) piby2))) 1843 (setq k *cancelled) 1844 (fpplus x (fpminus piby2)) 1845 (setq *cancelled (max k *cancelled)) 1846 (when *fpsincheck* 1847 (print `(*canc= ,*cancelled fpprec= ,fpprec oldprec= ,oldprec))) 1848 (cond ((not (> oldprec (- fpprec *cancelled))) 1849 (setq r (rem r 4)) 1850 (setq res 1851 (cond (fl (cond ((= r 0) (fpsin1 x)) 1852 ((= r 1) (fpcos1 x)) 1853 ((= r 2) (fpminus (fpsin1 x))) 1854 ((= r 3) (fpminus (fpcos1 x))))) 1855 (t (cond ((= r 0) (fpcos1 x)) 1856 ((= r 1) (fpminus (fpsin1 x))) 1857 ((= r 2) (fpminus (fpcos1 x))) 1858 ((= r 3) (fpsin1 x)))))) 1859 (return (bcons (if sign res (fpminus res))))) 1860 (t 1861 (incf fpprec *cancelled) 1862 (go loop)))))))))) 1863 1864(defun fpcos1 (x) 1865 (fpsincos1 x nil)) 1866 1867;; Compute SIN or COS in (0,PI/2). FL is T for SIN, NIL for COS. 1868;; 1869;; Use Taylor series 1870(defun fpsincos1 (x fl) 1871 (prog (ans term oans x2) 1872 (setq ans (if fl x (intofp 1)) 1873 x2 (fpminus(fptimes* x x))) 1874 (setq term ans) 1875 (do ((n (if fl 3 2) (+ n 2))) 1876 ((equal ans oans)) 1877 (setq term (fptimes* term (fpquotient x2 (intofp (* n (1- n)))))) 1878 (setq oans ans 1879 ans (fpplus ans term))) 1880 (return ans))) 1881 1882(defun fpsin1(x) 1883 (fpsincos1 x t)) 1884 1885(defun fpabs (x) 1886 (if (signp ge (car x)) 1887 x 1888 (cons (- (car x)) (cdr x)))) 1889 1890(defun fpentier (f) 1891 (let ((fpprec (bigfloat-prec f))) 1892 (fpintpart (cdr f)))) 1893 1894;; Calculate the integer part of a floating point number that is represented as 1895;; a list 1896;; 1897;; (MANTISSA EXPONENT) 1898;; 1899;; The special variable fpprec should be bound to the precision (in bits) of the 1900;; number. This encodes how many bits are known of the result and also a right 1901;; shift. The pair denotes the number MANTISSA * 2^(EXPONENT - FPPREC), of which 1902;; FPPREC bits are known. 1903;; 1904;; If EXPONENT is large and positive then we might not have enough 1905;; information to calculate the integer part. Specifically, we only 1906;; have enough information if EXPONENT < FPPREC. If that isn't the 1907;; case, we signal a Maxima error. However, if SKIP-EXPONENT-CHECK-P 1908;; is non-NIL, this check is skipped, and we compute the integer part 1909;; as requested. 1910;; 1911;; For the bigfloat code here, skip-exponent-check-p should be true. 1912;; For other uses (see commit 576c7508 and bug #2784), this should be 1913;; nil, which is the default. 1914(defun fpintpart (f &key skip-exponent-check-p) 1915 (destructuring-bind (mantissa exponent) 1916 f 1917 (let ((m (- fpprec exponent))) 1918 (if (plusp m) 1919 (quotient mantissa (expt 2 (- fpprec exponent))) 1920 (if (and (not skip-exponent-check-p) (< exponent fpprec)) 1921 (merror "~M doesn't have enough precision to compute its integer part" 1922 `((bigfloat ,fpprec) ,mantissa ,exponent)) 1923 (* mantissa (expt 2 (- m)))))))) 1924 1925(defun logbigfloat (a) 1926 (cond (($bfloatp (car a)) 1927 (big-float-log ($bfloat (car a)))) 1928 (t 1929 (list '(%log) (car a))))) 1930 1931 1932;;; Computes the log of a bigfloat number. 1933;;; 1934;;; Uses the series 1935;;; 1936;;; log(1+x) = sum((x/(x+2))^(2*n+1)/(2*n+1),n,0,inf); 1937;;; 1938;;; 1939;;; INF x 2 n + 1 1940;;; ==== (-----) 1941;;; \ x + 2 1942;;; = 2 > -------------- 1943;;; / 2 n + 1 1944;;; ==== 1945;;; n = 0 1946;;; 1947;;; 1948;;; which converges for x > 0. 1949;;; 1950;;; Note that FPLOG is given 1+X, not X. 1951;;; 1952;;; However, to aid convergence of the series, we scale 1+x until 1/e 1953;;; < 1+x <= e. 1954;;; 1955(defun fplog (x) 1956 (prog (over two ans oldans term e sum) 1957 (unless (> (car x) 0) 1958 (merror (intl:gettext "fplog: argument must be positive; found: ~M") (car x))) 1959 (setq e (fpe) 1960 over (fpquotient (fpone) e) 1961 ans 0) 1962 ;; Scale X until 1/e < X <= E. ANS keeps track of how 1963 ;; many factors of E were used. Set X to NIL if X is E. 1964 (do () 1965 (nil) 1966 (cond ((equal x e) (setq x nil) (return nil)) 1967 ((and (fplessp x e) (fplessp over x)) 1968 (return nil)) 1969 ((fplessp x over) 1970 (setq x (fptimes* x e)) 1971 (decf ans)) 1972 (t 1973 (incf ans) 1974 (setq x (fpquotient x e))))) 1975 (when (null x) (return (intofp (1+ ans)))) 1976 ;; Prepare X for the series. The series is for 1 + x, so 1977 ;; get x from our X. TERM is (x/(x+2)). X becomes 1978 ;; (x/(x+2))^2. 1979 (setq x (fpdifference x (fpone)) 1980 ans (intofp ans)) 1981 (setq x (fpexpt (setq term (fpquotient x (fpplus x (setq two (intofp 2))))) 2)) 1982 ;; Sum the series until the sum (in ANS) doesn't change 1983 ;; anymore. 1984 (setq sum (intofp 0)) 1985 (do ((n 1 (+ n 2))) 1986 ((equal sum oldans)) 1987 (setq oldans sum) 1988 (setq sum (fpplus sum (fpquotient term (intofp n)))) 1989 (setq term (fptimes* term x))) 1990 (return (fpplus ans (fptimes* two sum))))) 1991 1992(defun mabsbigfloat (l) 1993 (prog (r) 1994 (setq r (bigfloatp (car l))) 1995 (return (if (null r) 1996 (list '(mabs) (car l)) 1997 (bcons (fpabs (cdr r))))))) 1998 1999 2000;;;; Bigfloat implementations of special functions. 2001;;;; 2002 2003;;; This is still a bit messy. Some functions here take bigfloat 2004;;; numbers, represented by ((bigfloat) <mant> <exp>), but others want 2005;;; just the FP number, represented by (<mant> <exp>). Likewise, some 2006;;; return a bigfloat, some return just the FP. 2007;;; 2008;;; This needs to be systemized somehow. It isn't helped by the fact 2009;;; that some of the routines above also do the samething. 2010;;; 2011;;; The implementation for the special functions for a complex 2012;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex 2013;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in 2014;;; Iserles and Powell (eds.) "The State of the Art in Numerical 2015;;; Analysis", pp 165-211, Clarendon Press, 1987 2016 2017;; Compute exp(x) - 1, but do it carefully to preserve precision when 2018;; |x| is small. X is a FP number, and a FP number is returned. That 2019;; is, no bigfloat stuff. 2020(defun fpexpm1 (x) 2021 ;; What is the right breakpoint here? Is 1 ok? Perhaps 1/e is better? 2022 (cond ((fpgreaterp (fpabs x) (fpone)) 2023 ;; exp(x) - 1 2024 (fpdifference (fpexp x) (fpone))) 2025 (t 2026 ;; Use Taylor series for exp(x) - 1 2027 (let ((ans x) 2028 (oans nil) 2029 (term x)) 2030 (do ((n 2 (1+ n))) 2031 ((equal ans oans)) 2032 (setf term (fpquotient (fptimes* x term) (intofp n))) 2033 (setf oans ans) 2034 (setf ans (fpplus ans term))) 2035 ans)))) 2036 2037;; log(1+x) for small x. X is FP number, and a FP number is returned. 2038(defun fplog1p (x) 2039 ;; Use the same series as given above for fplog. For small x we use 2040 ;; the series, otherwise fplog is accurate enough. 2041 (cond ((fpgreaterp (fpabs x) (fpone)) 2042 (fplog (fpplus x (fpone)))) 2043 (t 2044 (let* ((sum (intofp 0)) 2045 (term (fpquotient x (fpplus x (intofp 2)))) 2046 (f (fptimes* term term)) 2047 (oldans nil)) 2048 (do ((n 1 (+ n 2))) 2049 ((equal sum oldans)) 2050 (setq oldans sum) 2051 (setq sum (fpplus sum (fpquotient term (intofp n)))) 2052 (setq term (fptimes* term f))) 2053 (fptimes* sum (intofp 2)))))) 2054 2055;; sinh(x) for real x. X is a bigfloat, and a bigfloat is returned. 2056(defun fpsinh (x) 2057 ;; X must be a maxima bigfloat 2058 2059 ;; See, for example, Hart et al., Computer Approximations, 6.2.27: 2060 ;; 2061 ;; sinh(x) = 1/2*(D(x) + D(x)/(1+D(x))) 2062 ;; 2063 ;; where D(x) = exp(x) - 1. 2064 ;; 2065 ;; But for negative x, use sinh(x) = -sinh(-x) because D(x) 2066 ;; approaches -1 for large negative x. 2067 (cond ((equal 0 (cadr x)) 2068 ;; Special case: x=0. Return immediately. 2069 (bigfloatp x)) 2070 ((fpposp (cdr x)) 2071 ;; x is positive. 2072 (let ((d (fpexpm1 (cdr (bigfloatp x))))) 2073 (bcons (fpquotient (fpplus d (fpquotient d (fpplus d (fpone)))) 2074 (intofp 2))))) 2075 (t 2076 ;; x is negative. 2077 (bcons 2078 (fpminus (cdr (fpsinh (bcons (fpminus (cdr (bigfloatp x))))))))))) 2079 2080(defun big-float-sinh (x &optional y) 2081 ;; The rectform for sinh for complex args should be numerically 2082 ;; accurate, so return nil in that case. 2083 (unless y 2084 (fpsinh x))) 2085 2086;; asinh(x) for real x. X is a bigfloat, and a bigfloat is returned. 2087(defun fpasinh (x) 2088 ;; asinh(x) = sign(x) * log(|x| + sqrt(1+x*x)) 2089 ;; 2090 ;; And 2091 ;; 2092 ;; asinh(x) = x, if 1+x*x = 1 2093 ;; = sign(x) * (log(2) + log(x)), large |x| 2094 ;; = sign(x) * log(2*|x| + 1/(|x|+sqrt(1+x*x))), if |x| > 2 2095 ;; = sign(x) * log1p(|x|+x^2/(1+sqrt(1+x*x))), otherwise. 2096 ;; 2097 ;; But I'm lazy right now and we only implement the last 2 cases. 2098 ;; We should implement all cases. 2099 (let* ((fp-x (cdr (bigfloatp x))) 2100 (absx (fpabs fp-x)) 2101 (one (fpone)) 2102 (two (intofp 2)) 2103 (minus (minusp (car fp-x))) 2104 result) 2105 ;; We only use two formulas here. |x| <= 2 and |x| > 2. Should 2106 ;; we add one for very big x and one for very small x, as given above. 2107 (cond ((fpgreaterp absx two) 2108 ;; |x| > 2 2109 ;; 2110 ;; log(2*|x| + 1/(|x|+sqrt(1+x^2))) 2111 (setf result (fplog (fpplus (fptimes* absx two) 2112 (fpquotient one 2113 (fpplus absx 2114 (fproot (bcons (fpplus one 2115 (fptimes* absx absx))) 2116 2))))))) 2117 (t 2118 ;; |x| <= 2 2119 ;; 2120 ;; log1p(|x|+x^2/(1+sqrt(1+x^2))) 2121 (let ((x*x (fptimes* absx absx))) 2122 (setq result (fplog1p (fpplus absx 2123 (fpquotient x*x 2124 (fpplus one 2125 (fproot (bcons (fpplus one x*x)) 2126 2))))))))) 2127 (if minus 2128 (bcons (fpminus result)) 2129 (bcons result)))) 2130 2131(defun complex-asinh (x y) 2132 ;; asinh(z) = -%i * asin(%i*z) 2133 (multiple-value-bind (u v) 2134 (complex-asin (mul -1 y) x) 2135 (values v (bcons (fpminus (cdr u)))))) 2136 2137(defun big-float-asinh (x &optional y) 2138 (if y 2139 (multiple-value-bind (u v) 2140 (complex-asinh x y) 2141 (add u (mul '$%i v))) 2142 (fpasinh x))) 2143 2144(defun fpasin-core (x) 2145 ;; asin(x) = atan(x/(sqrt(1-x^2)) 2146 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))] 2147 ;; 2148 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1. 2149 ;; 2150 ;; If |x| > 1, we need to do something else. 2151 ;; 2152 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x) 2153 ;; = -%i*log(%i*x + %i*sqrt(x^2-1)) 2154 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2] 2155 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|) 2156 2157 (let ((fp-x (cdr (bigfloatp x)))) 2158 (cond ((minusp (car fp-x)) 2159 ;; asin(-x) = -asin(x); 2160 (mul -1 (fpasin (bcons (fpminus fp-x))))) 2161 ((fplessp fp-x (cdr bfhalf)) 2162 ;; 0 <= x < 1/2 2163 ;; asin(x) = atan(x/sqrt(1-x^2)) 2164 (bcons 2165 (fpatan (fpquotient fp-x 2166 (fproot (bcons 2167 (fptimes* (fpdifference (fpone) fp-x) 2168 (fpplus (fpone) fp-x))) 2169 2))))) 2170 ((fpgreaterp fp-x (fpone)) 2171 ;; x > 1 2172 ;; asin(x) = %pi/2 - %i*log(|x+sqrt(x^2-1)|) 2173 ;; 2174 ;; Should we try to do something a little fancier with the 2175 ;; argument to log and use log1p for better accuracy? 2176 (let ((arg (fpplus fp-x 2177 (fproot (bcons (fptimes* (fpdifference fp-x (fpone)) 2178 (fpplus fp-x (fpone)))) 2179 2)))) 2180 (add (div '$%pi 2) 2181 (mul -1 '$%i (bcons (fplog arg)))))) 2182 2183 (t 2184 ;; 1/2 <= x <= 1 2185 ;; asin(x) = %pi/2 - atan(sqrt(1-x^2)/x) 2186 (add (div '$%pi 2) 2187 (mul -1 2188 (bcons 2189 (fpatan 2190 (fpquotient (fproot (bcons (fptimes* (fpdifference (fpone) fp-x) 2191 (fpplus (fpone) fp-x))) 2192 2) 2193 fp-x))))))))) 2194 2195;; asin(x) for real x. X is a bigfloat, and a maxima number (real or 2196;; complex) is returned. 2197(defun fpasin (x) 2198 ;; asin(x) = atan(x/(sqrt(1-x^2)) 2199 ;; = sgn(x)*[%pi/2 - atan(sqrt(1-x^2)/abs(x))] 2200 ;; 2201 ;; Use the first for 0 <= x < 1/2 and the latter for 1/2 < x <= 1. 2202 ;; 2203 ;; If |x| > 1, we need to do something else. 2204 ;; 2205 ;; asin(x) = -%i*log(sqrt(1-x^2)+%i*x) 2206 ;; = -%i*log(%i*x + %i*sqrt(x^2-1)) 2207 ;; = -%i*[log(|x + sqrt(x^2-1)|) + %i*%pi/2] 2208 ;; = %pi/2 - %i*log(|x+sqrt(x^2-1)|) 2209 2210 ($bfloat (fpasin-core x))) 2211 2212;; Square root of a complex number (xx, yy). Both are bigfloats. FP 2213;; (non-bigfloat) numbers are returned. 2214(defun complex-sqrt (xx yy) 2215 (let* ((x (cdr (bigfloatp xx))) 2216 (y (cdr (bigfloatp yy))) 2217 (rho (fpplus (fptimes* x x) 2218 (fptimes* y y)))) 2219 (setf rho (fpplus (fpabs x) (fproot (bcons rho) 2))) 2220 (setf rho (fpplus rho rho)) 2221 (setf rho (fpquotient (fproot (bcons rho) 2) (intofp 2))) 2222 2223 (let ((eta rho) 2224 (nu y)) 2225 (when (fpgreaterp rho (intofp 0)) 2226 (setf nu (fpquotient (fpquotient nu rho) (intofp 2))) 2227 (when (fplessp x (intofp 0)) 2228 (setf eta (fpabs nu)) 2229 (setf nu (if (minusp (car y)) 2230 (fpminus rho) 2231 rho)))) 2232 (values eta nu)))) 2233 2234;; asin(z) for complex z = x + %i*y. X and Y are bigfloats. The real 2235;; and imaginary parts are returned as bigfloat numbers. 2236(defun complex-asin (x y) 2237 (let ((x (cdr (bigfloatp x))) 2238 (y (cdr (bigfloatp y)))) 2239 (multiple-value-bind (re-sqrt-1-z im-sqrt-1-z) 2240 (complex-sqrt (bcons (fpdifference (intofp 1) x)) 2241 (bcons (fpminus y))) 2242 (multiple-value-bind (re-sqrt-1+z im-sqrt-1+z) 2243 (complex-sqrt (bcons (fpplus (intofp 1) x)) 2244 (bcons y)) 2245 ;; Realpart is atan(x/Re(sqrt(1-z)*sqrt(1+z))) 2246 ;; Imagpart is asinh(Im(conj(sqrt(1-z))*sqrt(1+z))) 2247 (values (bcons 2248 (let ((d (fpdifference (fptimes* re-sqrt-1-z 2249 re-sqrt-1+z) 2250 (fptimes* im-sqrt-1-z 2251 im-sqrt-1+z)))) 2252 ;; Check for division by zero. If we would divide 2253 ;; by zero, return pi/2 or -pi/2 according to the 2254 ;; sign of X. 2255 (cond ((equal d '(0 0)) 2256 (if (fplessp x '(0 0)) 2257 (fpminus (fpquotient (fppi) (intofp 2))) 2258 (fpquotient (fppi) (intofp 2)))) 2259 (t 2260 (fpatan (fpquotient x d)))))) 2261 (fpasinh (bcons 2262 (fpdifference (fptimes* re-sqrt-1-z 2263 im-sqrt-1+z) 2264 (fptimes* im-sqrt-1-z 2265 re-sqrt-1+z))))))))) 2266 2267(defun big-float-asin (x &optional y) 2268 (if y 2269 (multiple-value-bind (u v) (complex-asin x y) 2270 (add u (mul '$%i v))) 2271 (fpasin x))) 2272 2273 2274;; tanh(x) for real x. X is a bigfloat, and a bigfloat is returned. 2275(defun fptanh (x) 2276 ;; X is Maxima bigfloat 2277 ;; tanh(x) = D(2*x)/(2+D(2*x)) 2278 (let* ((two (intofp 2)) 2279 (fp (cdr (bigfloatp x))) 2280 (d (fpexpm1 (fptimes* fp two)))) 2281 (bcons (fpquotient d (fpplus d two))))) 2282 2283;; tanh(z), z = x + %i*y. X, Y are bigfloats, and a maxima number is 2284;; returned. 2285(defun complex-tanh (x y) 2286 (let* ((tv (cdr (tanbigfloat (list y)))) 2287 (beta (fpplus (fpone) (fptimes* tv tv))) 2288 (s (cdr (fpsinh x))) 2289 (s^2 (fptimes* s s)) 2290 (rho (fproot (bcons (fpplus (fpone) s^2)) 2)) 2291 (den (fpplus (fpone) (fptimes* beta s^2)))) 2292 (values (bcons (fpquotient (fptimes* beta (fptimes* rho s)) den)) 2293 (bcons (fpquotient tv den))))) 2294 2295(defun big-float-tanh (x &optional y) 2296 (if y 2297 (multiple-value-bind (u v) (complex-tanh x y) 2298 (add u (mul '$%i v))) 2299 (fptanh x))) 2300 2301;; atanh(x) for real x, |x| <= 1. X is a bigfloat, and a bigfloat is 2302;; returned. 2303(defun fpatanh (x) 2304 ;; atanh(x) = -atanh(-x) 2305 ;; = 1/2*log1p(2*x/(1-x)), x >= 0.5 2306 ;; = 1/2*log1p(2*x+2*x*x/(1-x)), x <= 0.5 2307 2308 (let* ((fp-x (cdr (bigfloatp x)))) 2309 (cond ((fplessp fp-x (intofp 0)) 2310 ;; atanh(x) = -atanh(-x) 2311 (mul -1 (fpatanh (bcons (fpminus fp-x))))) 2312 ((fpgreaterp fp-x (fpone)) 2313 ;; x > 1, so use complex version. 2314 (multiple-value-bind (u v) 2315 (complex-atanh x (bcons (intofp 0))) 2316 (add u (mul '$%i v)))) 2317 ((fpgreaterp fp-x (cdr bfhalf)) 2318 ;; atanh(x) = 1/2*log1p(2*x/(1-x)) 2319 (bcons 2320 (fptimes* (cdr bfhalf) 2321 (fplog1p (fpquotient (fptimes* (intofp 2) fp-x) 2322 (fpdifference (fpone) fp-x)))))) 2323 (t 2324 ;; atanh(x) = 1/2*log1p(2*x + 2*x*x/(1-x)) 2325 (let ((2x (fptimes* (intofp 2) fp-x))) 2326 (bcons 2327 (fptimes* (cdr bfhalf) 2328 (fplog1p (fpplus 2x 2329 (fpquotient (fptimes* 2x fp-x) 2330 (fpdifference (fpone) fp-x))))))))))) 2331 2332;; Stuff which follows is derived from atanh z = (log(1 + z) - log(1 - z))/2 2333;; which apparently originates with Kahan's "Much ado" paper. 2334 2335;; The formulas for eta and nu below can be easily derived from 2336;; rectform(atanh(x+%i*y)) = 2337;; 2338;; 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))) 2339;; 2340;; Expand the argument of log out and divide it out and we get 2341;; 2342;; log(((1+x)^2+y^2)/((1-x)^2+y^2)) = log(1+4*x/((1-x)^2+y^2)) 2343;; 2344;; When y = 0, Im atanh z = 1/2 (arg(1 + x) - arg(1 - x)) 2345;; = if x < -1 then %pi/2 else if x > 1 then -%pi/2 else <whatever> 2346;; 2347;; Otherwise, arg(1 - x + %i*(-y)) = - arg(1 - x + %i*y), 2348;; and Im atanh z = 1/2 (arg(1 + x + %i*y) + arg(1 - x + %i*y)). 2349;; Since arg(x)+arg(y) = arg(x*y) (almost), we can simplify the 2350;; imaginary part to 2351;; 2352;; arg((1+x+%i*y)*(1-x+%i*y)) = arg((1-x)*(1+x)-y^2+2*y*%i) 2353;; = atan2(2*y,((1-x)*(1+x)-y^2)) 2354;; 2355;; These are the eta and nu forms below. 2356(defun complex-atanh (x y) 2357 (let* ((fpx (cdr (bigfloatp x))) 2358 (fpy (cdr (bigfloatp y))) 2359 (beta (if (minusp (car fpx)) 2360 (fpminus (fpone)) 2361 (fpone))) 2362 (x-lt-minus-1 (mevalp `((mlessp) ,x -1))) 2363 (x-gt-plus-1 (mevalp `((mgreaterp) ,x 1))) 2364 (y-equals-0 (like y '((bigfloat) 0 0))) 2365 (x (fptimes* beta fpx)) 2366 (y (fptimes* beta (fpminus fpy))) 2367 ;; Kahan has rho = 4/most-positive-float. What should we do 2368 ;; here about that? Our big floats don't really have a 2369 ;; most-positive float value. 2370 (rho (intofp 0)) 2371 (t1 (fpplus (fpabs y) rho)) 2372 (t1^2 (fptimes* t1 t1)) 2373 (1-x (fpdifference (fpone) x)) 2374 ;; eta = log(1+4*x/((1-x)^2+y^2))/4 2375 (eta (fpquotient 2376 (fplog1p (fpquotient (fptimes* (intofp 4) x) 2377 (fpplus (fptimes* 1-x 1-x) 2378 t1^2))) 2379 (intofp 4))) 2380 ;; If y = 0, then Im atanh z = %pi/2 or -%pi/2. 2381 ;; Otherwise nu = 1/2*atan2(2*y,(1-x)*(1+x)-y^2) 2382 (nu (if y-equals-0 2383 ;; EXTRA FPMINUS HERE TO COUNTERACT FPMINUS IN RETURN VALUE 2384 (fpminus (if x-lt-minus-1 2385 (cdr ($bfloat '((mquotient) $%pi 2))) 2386 (if x-gt-plus-1 2387 (cdr ($bfloat '((mminus) ((mquotient) $%pi 2)))) 2388 (merror "COMPLEX-ATANH: HOW DID I GET HERE?")))) 2389 (fptimes* (cdr bfhalf) 2390 (fpatan2 (fptimes* (intofp 2) y) 2391 (fpdifference (fptimes* 1-x (fpplus (fpone) x)) 2392 t1^2)))))) 2393 (values (bcons (fptimes* beta eta)) 2394 ;; WTF IS FPMINUS DOING HERE ?? 2395 (bcons (fpminus (fptimes* beta nu)))))) 2396 2397(defun big-float-atanh (x &optional y) 2398 (if y 2399 (multiple-value-bind (u v) (complex-atanh x y) 2400 (add u (mul '$%i v))) 2401 (fpatanh x))) 2402 2403;; acos(x) for real x. X is a bigfloat, and a maxima number is returned. 2404(defun fpacos (x) 2405 ;; acos(x) = %pi/2 - asin(x) 2406 ($bfloat (add (div '$%pi 2) (mul -1 (fpasin-core x))))) 2407 2408(defun complex-acos (x y) 2409 (let ((x (cdr (bigfloatp x))) 2410 (y (cdr (bigfloatp y)))) 2411 (multiple-value-bind (re-sqrt-1-z im-sqrt-1-z) 2412 (complex-sqrt (bcons (fpdifference (intofp 1) x)) 2413 (bcons (fpminus y))) 2414 (multiple-value-bind (re-sqrt-1+z im-sqrt-1+z) 2415 (complex-sqrt (bcons (fpplus (intofp 1) x)) 2416 (bcons y)) 2417 (values (bcons 2418 (fptimes* (intofp 2) 2419 (fpatan (fpquotient re-sqrt-1-z re-sqrt-1+z)))) 2420 (fpasinh (bcons 2421 (fpdifference 2422 (fptimes* re-sqrt-1+z im-sqrt-1-z) 2423 (fptimes* im-sqrt-1+z re-sqrt-1-z))))))))) 2424 2425 2426(defun big-float-acos (x &optional y) 2427 (if y 2428 (multiple-value-bind (u v) (complex-acos x y) 2429 (add u (mul '$%i v))) 2430 (fpacos x))) 2431 2432(defun complex-log (x y) 2433 (let* ((x (cdr (bigfloatp x))) 2434 (y (cdr (bigfloatp y))) 2435 (t1 (let (($float2bf t)) 2436 ;; No warning message, please. 2437 (floattofp 1.2))) 2438 (t2 (intofp 3)) 2439 (rho (fpplus (fptimes* x x) 2440 (fptimes* y y))) 2441 (abs-x (fpabs x)) 2442 (abs-y (fpabs y)) 2443 (beta (fpmax abs-x abs-y)) 2444 (theta (fpmin abs-x abs-y))) 2445 (values (if (or (fpgreaterp t1 beta) 2446 (fplessp rho t2)) 2447 (fpquotient (fplog1p (fpplus (fptimes* (fpdifference beta (fpone)) 2448 (fpplus beta (fpone))) 2449 (fptimes* theta theta))) 2450 (intofp 2)) 2451 (fpquotient (fplog rho) (intofp 2))) 2452 (fpatan2 y x)))) 2453 2454(defun big-float-log (x &optional y) 2455 (if y 2456 (multiple-value-bind (u v) (complex-log x y) 2457 (add (bcons u) (mul '$%i (bcons v)))) 2458 (flet ((%log (x) 2459 ;; x is (mantissa exp), where mantissa = frac*2^fpprec, 2460 ;; with 1/2 < frac <= 1 and x is frac*2^exp. To 2461 ;; compute log(x), use log(x) = log(frac)+ exp*log(2). 2462 (cdr 2463 (let* ((extra 8) 2464 (fpprec (+ fpprec extra)) 2465 (log-frac 2466 (fplog #+nil 2467 (cdr ($bfloat 2468 (cl-rat-to-maxima (/ (car x) 2469 (ash 1 (- fpprec 8)))))) 2470 (list (ash (car x) extra) 0))) 2471 (log-exp (fptimes* (intofp (second x)) (fplog2))) 2472 (result (bcons (fpplus log-frac log-exp)))) 2473 (let ((fpprec (- fpprec extra))) 2474 (bigfloatp result)))))) 2475 (let ((fp-x (cdr (bigfloatp x)))) 2476 (cond ((onep1 x) 2477 ;; Special case for log(1). See Bug 3381301: 2478 ;; https://sourceforge.net/tracker/?func=detail&aid=3381301&group_id=4933&atid=104933 2479 (bcons (intofp 0))) 2480 ((fplessp fp-x (intofp 0)) 2481 ;; ??? Do we want to return an exact %i*%pi or a float 2482 ;; approximation? 2483 (add (big-float-log (bcons (fpminus fp-x))) 2484 (mul '$%i (bcons (fppi))))) 2485 (t 2486 (bcons (%log fp-x)))))))) 2487 2488(defun big-float-sqrt (x &optional y) 2489 (if y 2490 (multiple-value-bind (u v) (complex-sqrt x y) 2491 (add (bcons u) (mul '$%i (bcons v)))) 2492 (let ((fp-x (cdr (bigfloatp x)))) 2493 (if (fplessp fp-x (intofp 0)) 2494 (mul '$%i (bcons (fproot (bcons (fpminus fp-x)) 2))) 2495 (bcons (fproot x 2)))))) 2496 2497(eval-when 2498 #+gcl (load eval) 2499 #-gcl (:load-toplevel :execute) 2500 (fpprec1 nil $fpprec)) ; Set up user's precision 2501