1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10-*- ;;;; 2 3;;; This package contains a numeric class for use with Maxima. The 4;;; purpose is to allow users to write numerical algorithms that 5;;; support double-float, (complex double-float) and Maxima bfloat and 6;;; complex bfloat arithmetic, without having to write separate 7;;; versions for each. Of course, specially written versions for 8;;; double-float and (complex double-float) will probably be much 9;;; faster, but this allows users to write just one routine that can 10;;; handle all of the data types in a more "natural" Lisp style. 11 12#+cmu 13(eval-when (:compile-toplevel :load-toplevel :execute) 14 (setf lisp::*enable-package-locked-errors* nil)) 15 16(in-package #-gcl #:bigfloat #+gcl "BIGFLOAT") 17 18(defun intofp (re) 19 ;; Kind of like Maxima's INTOFP, but we only handle numeric types. 20 ;; We should return a Maxima bigfloat object (list of bigfloat 21 ;; marker, mantissa, and exponent). 22 (cond ((floatp re) 23 (maxima::bcons (maxima::floattofp re))) 24 ((eql re 0) 25 (maxima::bcons '(0 0))) 26 ((integerp re) 27 (maxima::bcons (list (maxima::fpround re) (cl:+ maxima::*m maxima::fpprec)))) 28 ((typep re 'ratio) 29 ;; Should we do something better than converting the 30 ;; numerator and denominator to floats and dividing? 31 (maxima::bcons (maxima::fpquotient (cdr (intofp (numerator re))) 32 (cdr (intofp (denominator re)))))) 33 ((maxima::$bfloatp re) 34 ;; Call bigfloatp to make sure we round/scale the bigfloat to 35 ;; the correct precision! 36 (maxima::bigfloatp re)) 37 (t 38 (error "Don't know how to convert ~S to a BIGFLOAT" re)))) 39 40(defclass numeric () 41 () 42 (:documentation "Basic class, like CL's NUMBER type")) 43 44(defclass bigfloat (numeric) 45 ;; We store the Maxima internal bigfloat format because we also need 46 ;; the precision in case we have mixed size bigfloat operations. 47 ;; (We could recompute it from the size of the mantissa part, but 48 ;; why bother? 49 ((real :initform (intofp 0) 50 :initarg :real 51 :documentation "A Maxima bigfloat. This contains the full 52 Maxima bigfloat object including the mantissa, the exponent 53 and the bigfloat marker and precision." )) 54 (:documentation "Big float, equivalent to a Maxima bfloat object")) 55 56;; Extract the internal representation of a bigfloat, and adjust the 57;; precision to the current value of fpprec. 58(defmethod real-value ((b bigfloat)) 59 (maxima::bigfloatp (slot-value b 'real))) 60 61(defclass complex-bigfloat (numeric) 62 ;; Currently, the real and imaginary parts contain a Maxima bigfloat 63 ;; including the bigfloat marker and the mantissa and exponent. 64 ;; Should they be BIGFLOAT objects instead? 65 ((real :initform (intofp 0) 66 :initarg :real 67 :documentation "Real part of a complex bigfloat") 68 (imag :initform (intofp 0) 69 :initarg :imag 70 :documentation "Imaginary part of a complex bigfloat")) 71 (:documentation "Complex bigfloat")) 72 73;; Extract the internal representation of the real part of a 74;; complex-bigfloat, and adjust the precision to the current value of 75;; fpprec. 76(defmethod real-value ((b complex-bigfloat)) 77 (maxima::bigfloatp (slot-value b 'real))) 78 79;; Extract the internal representation of the imaginary part of a 80;; complex-bigfloat, and adjust the precision to the current value of 81;; fpprec. 82(defmethod imag-value ((b complex-bigfloat)) 83 (maxima::bigfloatp (slot-value b 'imag))) 84 85(defmethod make-load-form ((x bigfloat) &optional environment) 86 (declare (ignore environment)) 87 `(make-instance ',(class-of x) 88 :real ',(real-value x))) 89 90;;; BIGFLOAT - External 91;;; 92;;; BIGFLOAT converts a number to a BIGFLOAT or COMPLEX-BIGFLOAT. 93;;; This is intended to convert CL numbers or Maxima (internal) 94;;; numbers to a bigfloat object. 95(defun bigfloat (re &optional im) 96 "Convert RE to a BIGFLOAT. If IM is given, return a COMPLEX-BIGFLOAT" 97 (cond (im 98 (make-instance 'complex-bigfloat 99 :real (intofp re) 100 :imag (intofp im))) 101 ((cl:realp re) 102 (make-instance 'bigfloat :real (intofp re))) 103 ((cl:complexp re) 104 (make-instance 'complex-bigfloat 105 :real (intofp (cl:realpart re)) 106 :imag (intofp (cl:imagpart re)))) 107 ((maxima::$bfloatp re) 108 (make-instance 'bigfloat :real (intofp re))) 109 ((maxima::complex-number-p re 'maxima::bigfloat-or-number-p) 110 (make-instance 'complex-bigfloat 111 :real (intofp (maxima::$realpart re)) 112 :imag (intofp (maxima::$imagpart re)))) 113 ((typep re 'bigfloat) 114 ;; Done this way so the new bigfloat updates the precision of 115 ;; the given bigfloat, if necessary. 116 (make-instance 'bigfloat :real (real-value re))) 117 ((typep re 'complex-bigfloat) 118 ;; Done this way so the new bigfloat updates the precision of 119 ;; the given bigfloat, if necessary. 120 (make-instance 'complex-bigfloat 121 :real (real-value re) 122 :imag (imag-value re))) 123 (t 124 (make-instance 'bigfloat :real (intofp re))))) 125 126 127;;; MAXIMA::TO - External 128;;; 129;;; Convert a CL number, a BIGFLOAT, or a COMPLEX-BIGFLOAT to 130;;; Maxima's internal representation of the number. 131(defmethod maxima::to ((z cl:float)) 132 z) 133 134(defmethod maxima::to ((z cl:rational)) 135 (if (typep z 'ratio) 136 (list '(maxima::rat maxima::simp) (numerator z) (denominator z)) 137 z)) 138 139(defmethod maxima::to ((z cl:complex)) 140 (maxima::add (maxima::to (cl:realpart z)) 141 (maxima::mul 'maxima::$%i 142 (maxima::to (cl:imagpart z))))) 143 144(defmethod maxima::to ((z bigfloat)) 145 "Convert BIGFLOAT object to a maxima number" 146 (real-value z)) 147 148(defmethod maxima::to ((z complex-bigfloat)) 149 "Convert COMPLEX-BIGFLOAT object to a maxima number" 150 (maxima::add (real-value z) 151 (maxima::mul 'maxima::$%i 152 (imag-value z)))) 153 154(defmethod maxima::to ((z t)) 155 z) 156 157;; MAX-EXPONENT roughly computes the log2(|x|). If x is real and x = 158;; 2^n*f, with |f| < 1, MAX-EXPONENT returns |n|. For complex 159;; numbers, we return one more than the max of the exponent of the 160;; real and imaginary parts. 161(defmethod max-exponent ((x bigfloat)) 162 ;; The third element is the exponent of a bigfloat. 163 (cl:abs (third (slot-value x 'real)))) 164 165(defmethod max-exponent ((x complex-bigfloat)) 166 (cl:1+ (cl:max (cl:abs (third (slot-value x 'real))) 167 (cl:abs (third (slot-value x 'imag)))))) 168 169(defmethod max-exponent ((x cl:float)) 170 (if (zerop x) 171 0 172 (cl:abs (nth-value 1 (cl:decode-float x))))) 173 174(defmethod max-exponent ((x cl:rational)) 175 (if (zerop x) 176 0 177 (cl:ceiling (cl:log (cl:abs x) 2)))) 178 179(defmethod max-exponent ((x cl:complex)) 180 (cl:1+ (cl:max (max-exponent (cl:realpart x)) 181 (max-exponent (cl:imagpart x))))) 182 183;; When computing x^a using exp(a*log(x)), we need extra bits because 184;; the integer part of a*log(x) doesn't contribute to the accuracy of 185;; the result. The number of extra bits needed is basically the 186;; "size" of a plus the number of bits for ceiling(log(x)). We need 187;; ceiling(log(x)) extra bits because that's how many bits are taken 188;; up by the log(x). The "size" of a is, basically, the exponent of 189;; a. If a = 2^n*f where |f| < 1, then the size is abs(n) because 190;; that's how many extra bits are added to the integer part of 191;; a*log(x). If either |x| or |a| < 1, the size is 0, since no 192;; additional bits are taken up. 193(defun expt-extra-bits (x a) 194 (max 1 (+ (integer-length (max-exponent x)) 195 (max-exponent a)))) 196 197;;; WITH-EXTRA-PRECISION - Internal 198;;; 199;;; Executes the body BODY with extra precision. The precision is 200;;; increased by EXTRA, and the list of variables given in VARLIST have 201;;; the precision increased. The precision of the first value of the 202;;; body is then reduced back to the normal precision. 203(defmacro with-extra-precision ((extra (&rest varlist)) &body body) 204 (let ((result (gensym)) 205 (old-fpprec (gensym))) 206 `(let ((,result 207 (let ((,old-fpprec maxima::fpprec)) 208 (unwind-protect 209 (let ((maxima::fpprec (cl:+ maxima::fpprec ,extra))) 210 (let ,(mapcar #'(lambda (v) 211 ;; Could probably do this in a faster 212 ;; way, but conversion to a maxima 213 ;; form automatically increases the 214 ;; precision of the bigfloat to the 215 ;; new precision. Conversion of that 216 ;; to a bigfloat object preserves the 217 ;; precision. 218 `(,v (bigfloat:to (maxima::to ,v)))) 219 varlist) 220 ,@body)) 221 (setf maxima::fpprec ,old-fpprec))))) 222 ;; Conversion of the result to a maxima number adjusts the 223 ;; precision appropriately. 224 (bigfloat:to (maxima::to ,result))))) 225 226;;; REALP 227 228;; GCL doesn't have the REAL class! But since a real is a rational or 229;; float, we can fake it by defining methods on rationals and floats 230;; for gcl. 231#-gcl 232(defmethod realp ((x cl:real)) 233 t) 234 235#+gcl 236(progn 237 (defmethod realp ((x cl:rational)) 238 t) 239 (defmethod realp ((x cl:float)) 240 t) 241 ) 242 243 244(defmethod realp ((x bigfloat)) 245 t) 246 247(defmethod realp ((x t)) 248 nil) 249 250;;; COMPLEXP 251(defmethod complexp ((x cl:complex)) 252 t) 253 254(defmethod complexp ((x complex-bigfloat)) 255 t) 256 257(defmethod complexp ((x t)) 258 nil) 259 260;;; NUMBERP 261(defmethod numberp ((x cl:number)) 262 t) 263 264(defmethod numberp ((x bigfloat)) 265 t) 266 267(defmethod numberp ((x complex-bigfloat)) 268 t) 269 270(defmethod numberp ((x t)) 271 nil) 272 273(defmethod make-load-form ((x complex-bigfloat) &optional environment) 274 (declare (ignore environment)) 275 `(make-instance ',(class-of x) 276 :real ',(real-value x) 277 :imag ',(imag-value x))) 278 279;; The print-object and describe-object methods are mostly for 280;; debugging purposes. Maxima itself shouldn't ever see such objects. 281(defmethod print-object ((x bigfloat) stream) 282 (let ((r (cdr (real-value x)))) 283 (multiple-value-bind (sign output-list) 284 (if (cl:minusp (first r)) 285 (values "-" (maxima::fpformat (maxima::bcons (list (cl:- (first r)) (second r))))) 286 (values "+" (maxima::fpformat (maxima::bcons r)))) 287 (format stream "~A~{~D~}" sign output-list)))) 288 289(defmethod print-object ((x complex-bigfloat) stream) 290 (format stream "~A~A*%i" (realpart x) (imagpart x))) 291 292(defmethod describe-object ((x bigfloat) stream) 293 (let ((r (slot-value x 'real))) 294 (format stream "~&~S is a ~D-bit BIGFLOAT with mantissa ~D and exponent ~D~%" 295 x (third (first r)) (second r) (third r)))) 296 297(defmethod describe-object ((x complex-bigfloat) stream) 298 (format stream "~S is a COMPLEX-BIGFLOAT~%" x) 299 (describe-object (make-instance 'bigfloat :real (slot-value x 'real)) stream) 300 (describe-object (make-instance 'bigfloat :real (slot-value x 'imag)) stream)) 301 302 303(defgeneric add1 (a) 304 (:documentation "Add 1")) 305 306(defgeneric sub1 (a) 307 (:documentation "Subtract 1")) 308 309 310(defgeneric two-arg-+ (a b) 311 (:documentation "A + B")) 312 313(defgeneric two-arg-- (a b) 314 (:documentation "A - B")) 315 316(defgeneric two-arg-* (a b) 317 (:documentation "A * B")) 318 319(defgeneric two-arg-/ (a b) 320 (:documentation "A / B")) 321 322(defgeneric two-arg-< (a b) 323 (:documentation "A < B")) 324 325(defgeneric two-arg-> (a b) 326 (:documentation "A > B")) 327 328(defgeneric two-arg-<= (a b) 329 (:documentation "A <= B")) 330 331(defgeneric two-arg->= (a b) 332 (:documentation "A >= B")) 333 334(defgeneric two-arg-= (a b) 335 (:documentation "A = B?")) 336 337 338(defgeneric unary-minus (a) 339 (:documentation "-A")) 340 341(defgeneric unary-divide (a) 342 (:documentation "1 / A")) 343 344 345;;; Basic arithmetic operations 346 347;;; 1+ and 1- 348 349(defmethod add1 ((a number)) 350 (cl::1+ a)) 351 352(defmethod add1 ((a bigfloat)) 353 (make-instance 'bigfloat 354 :real (maxima::bcons 355 (maxima::fpplus (cdr (real-value a)) 356 (maxima::fpone))))) 357 358(defmethod add1 ((a complex-bigfloat)) 359 (make-instance 'complex-bigfloat 360 :real (maxima::bcons 361 (maxima::fpplus (cdr (real-value a)) 362 (maxima::fpone))) 363 :imag (imag-value a))) 364 365(defmethod sub1 ((a number)) 366 (cl::1- a)) 367 368(defmethod sub1 ((a bigfloat)) 369 (make-instance 'bigfloat 370 :real (maxima::bcons 371 (maxima::fpdifference (cdr (real-value a)) 372 (maxima::fpone))))) 373 374(defmethod sub1 ((a complex-bigfloat)) 375 (make-instance 'complex-bigfloat 376 :real (maxima::bcons 377 (maxima::fpdifference (cdr (real-value a)) 378 (maxima::fpone))) 379 :imag (imag-value a))) 380 381(declaim (inline 1+ 1-)) 382 383(defun 1+ (x) 384 (add1 x)) 385 386(defun 1- (x) 387 (sub1 x)) 388 389;; Add two numbers 390(defmethod two-arg-+ ((a cl:number) (b cl:number)) 391 (cl:+ a b)) 392 393(defmethod two-arg-+ ((a bigfloat) (b bigfloat)) 394 (make-instance 'bigfloat 395 :real (maxima::bcons 396 (maxima::fpplus (cdr (real-value a)) 397 (cdr (real-value b)))))) 398 399(defmethod two-arg-+ ((a complex-bigfloat) (b complex-bigfloat)) 400 (make-instance 'complex-bigfloat 401 :real (maxima::bcons 402 (maxima::fpplus (cdr (real-value a)) 403 (cdr (real-value b)))) 404 :imag (maxima::bcons 405 (maxima::fpplus (cdr (imag-value a)) 406 (cdr (imag-value b)))))) 407 408;; Handle contagion for two-arg-+ 409(defmethod two-arg-+ ((a bigfloat) (b cl:float)) 410 (make-instance 'bigfloat 411 :real (maxima::bcons 412 (maxima::fpplus (cdr (real-value a)) 413 (cdr (intofp b)))))) 414 415(defmethod two-arg-+ ((a bigfloat) (b cl:rational)) 416 (make-instance 'bigfloat 417 :real (maxima::bcons 418 (maxima::fpplus (cdr (real-value a)) 419 (cdr (intofp b)))))) 420 421(defmethod two-arg-+ ((a bigfloat) (b cl:complex)) 422 (make-instance 'complex-bigfloat 423 :real (maxima::bcons 424 (maxima::fpplus (cdr (real-value a)) 425 (cdr (intofp (realpart b))))) 426 :imag (intofp (imagpart b)))) 427 428(defmethod two-arg-+ ((a cl:number) (b bigfloat)) 429 (two-arg-+ b a)) 430 431(defmethod two-arg-+ ((a complex-bigfloat) (b bigfloat)) 432 (make-instance 'complex-bigfloat 433 :real (maxima::bcons 434 (maxima::fpplus (cdr (real-value a)) 435 (cdr (real-value b)))) 436 :imag (imag-value a))) 437 438(defmethod two-arg-+ ((a complex-bigfloat) (b number)) 439 (make-instance 'complex-bigfloat 440 :real (maxima::bcons 441 (maxima::fpplus (cdr (real-value a)) 442 (cdr (intofp (cl:realpart b))))) 443 :imag (maxima::bcons 444 (maxima::fpplus (cdr (imag-value a)) 445 (cdr (intofp (cl:imagpart b))))))) 446 447(defmethod two-arg-+ ((a bigfloat) (b complex-bigfloat)) 448 (two-arg-+ b a)) 449 450(defmethod two-arg-+ ((a number) (b complex-bigfloat)) 451 (two-arg-+ b a)) 452 453(defun + (&rest args) 454 (if (null args) 455 0 456 (do ((args (cdr args) (cdr args)) 457 (res (car args) 458 (two-arg-+ res (car args)))) 459 ((null args) res)))) 460 461;; Negate a number 462(defmethod unary-minus ((a number)) 463 (cl:- a)) 464 465(defmethod unary-minus ((a bigfloat)) 466 (make-instance 'bigfloat 467 :real (maxima::bcons 468 (maxima::fpminus (cdr (real-value a)))))) 469 470(defmethod unary-minus ((a complex-bigfloat)) 471 (make-instance 'complex-bigfloat 472 :real (maxima::bcons 473 (maxima::fpminus (cdr (real-value a)))) 474 :imag (maxima::bcons 475 (maxima::fpminus (cdr (imag-value a)))))) 476 477;;; Subtract two numbers 478(defmethod two-arg-- ((a number) (b number)) 479 (cl:- a b)) 480 481(defmethod two-arg-- ((a bigfloat) (b bigfloat)) 482 (make-instance 'bigfloat 483 :real (maxima::bcons 484 (maxima::fpdifference (cdr (real-value a)) 485 (cdr (real-value b)))))) 486 487(defmethod two-arg-- ((a complex-bigfloat) (b complex-bigfloat)) 488 (make-instance 'complex-bigfloat 489 :real (maxima::bcons 490 (maxima::fpdifference (cdr (real-value a)) 491 (cdr (real-value b)))) 492 :imag (maxima::bcons 493 (maxima::fpdifference (cdr (imag-value a)) 494 (cdr (imag-value b)))))) 495 496;; Handle contagion for two-arg-- 497(defmethod two-arg-- ((a bigfloat) (b cl:float)) 498 (make-instance 'bigfloat 499 :real (maxima::bcons 500 (maxima::fpdifference (cdr (real-value a)) 501 (cdr (intofp b)))))) 502 503(defmethod two-arg-- ((a bigfloat) (b cl:rational)) 504 (make-instance 'bigfloat 505 :real (maxima::bcons 506 (maxima::fpdifference (cdr (real-value a)) 507 (cdr (intofp b)))))) 508 509(defmethod two-arg-- ((a bigfloat) (b cl:complex)) 510 (make-instance 'complex-bigfloat 511 :real (maxima::bcons 512 (maxima::fpdifference (cdr (real-value a)) 513 (cdr (intofp (realpart b))))) 514 :imag (maxima::bcons (maxima::fpminus (cdr (intofp (imagpart b))))))) 515 516(defmethod two-arg-- ((a cl:float) (b bigfloat)) 517 (make-instance 'bigfloat 518 :real (maxima::bcons 519 (maxima::fpdifference (cdr (intofp a)) 520 (cdr (real-value b)))))) 521 522(defmethod two-arg-- ((a cl:rational) (b bigfloat)) 523 (make-instance 'bigfloat 524 :real (maxima::bcons 525 (maxima::fpdifference (cdr (intofp a)) 526 (cdr (real-value b)))))) 527 528(defmethod two-arg-- ((a cl:complex) (b bigfloat)) 529 (two-arg-- (bigfloat (cl:realpart a) (cl:imagpart a)) b)) 530 531(defmethod two-arg-- ((a complex-bigfloat) (b bigfloat)) 532 (make-instance 'complex-bigfloat 533 :real (maxima::bcons 534 (maxima::fpdifference (cdr (real-value a)) 535 (cdr (real-value b)))) 536 :imag (imag-value a))) 537 538(defmethod two-arg-- ((a complex-bigfloat) (b number)) 539 (if (cl:complexp b) 540 (two-arg-- a (bigfloat (cl:realpart b) (cl:imagpart b))) 541 (two-arg-- a (bigfloat b)))) 542 543(defmethod two-arg-- ((a bigfloat) (b complex-bigfloat)) 544 (make-instance 'complex-bigfloat 545 :real (maxima::bcons 546 (maxima::fpdifference (cdr (real-value a)) 547 (cdr (real-value b)))) 548 :imag (maxima::bcons (maxima::fpminus (cdr (imag-value b)))))) 549 550(defmethod two-arg-- ((a number) (b complex-bigfloat)) 551 (if (cl:complexp a) 552 (two-arg-- (bigfloat (cl:realpart a) (cl:imagpart a)) 553 b) 554 (two-arg-- (bigfloat a) b))) 555 556(defun - (number &rest more-numbers) 557 (if more-numbers 558 (do ((nlist more-numbers (cdr nlist)) 559 (result number)) 560 ((atom nlist) result) 561 (declare (list nlist)) 562 (setq result (two-arg-- result (car nlist)))) 563 (unary-minus number))) 564 565;;; Multiply two numbers 566(defmethod two-arg-* ((a number) (b number)) 567 (cl:* a b)) 568 569(defmethod two-arg-* ((a bigfloat) (b bigfloat)) 570 (make-instance 'bigfloat 571 :real (maxima::bcons 572 (maxima::fptimes* (cdr (real-value a)) 573 (cdr (real-value b)))))) 574 575(defmethod two-arg-* ((a complex-bigfloat) (b complex-bigfloat)) 576 (let ((a-re (cdr (real-value a))) 577 (a-im (cdr (imag-value a))) 578 (b-re (cdr (real-value b))) 579 (b-im (cdr (imag-value b)))) 580 (make-instance 'complex-bigfloat 581 :real (maxima::bcons 582 (maxima::fpdifference (maxima::fptimes* a-re b-re) 583 (maxima::fptimes* a-im b-im))) 584 :imag (maxima::bcons 585 (maxima::fpplus (maxima::fptimes* a-re b-im) 586 (maxima::fptimes* a-im b-re)))))) 587 588;; Handle contagion for two-arg-* 589(defmethod two-arg-* ((a bigfloat) (b cl:float)) 590 (make-instance 'bigfloat 591 :real (maxima::bcons 592 (maxima::fptimes* (cdr (real-value a)) 593 (cdr (intofp b)))))) 594 595(defmethod two-arg-* ((a bigfloat) (b cl:rational)) 596 (make-instance 'bigfloat 597 :real (maxima::bcons 598 (maxima::fptimes* (cdr (real-value a)) 599 (cdr (intofp b)))))) 600 601(defmethod two-arg-* ((a bigfloat) (b cl:complex)) 602 (make-instance 'complex-bigfloat 603 :real (maxima::bcons 604 (maxima::fptimes* (cdr (real-value a)) 605 (cdr (intofp (realpart b))))) 606 :imag (maxima::bcons 607 (maxima::fptimes* (cdr (real-value a)) 608 (cdr (intofp (imagpart b))))))) 609 610(defmethod two-arg-* ((a cl:number) (b bigfloat)) 611 (two-arg-* b a)) 612 613(defmethod two-arg-* ((a complex-bigfloat) (b bigfloat)) 614 (make-instance 'complex-bigfloat 615 :real (maxima::bcons 616 (maxima::fptimes* (cdr (real-value a)) 617 (cdr (real-value b)))) 618 :imag (maxima::bcons 619 (maxima::fptimes* (cdr (imag-value a)) 620 (cdr (real-value b)))))) 621 622(defmethod two-arg-* ((a complex-bigfloat) (b number)) 623 (if (cl:complexp b) 624 (two-arg-* a (bigfloat (cl:realpart b) (cl:imagpart b))) 625 (two-arg-* a (bigfloat b)))) 626 627(defmethod two-arg-* ((a bigfloat) (b complex-bigfloat)) 628 (two-arg-* b a)) 629 630(defmethod two-arg-* ((a number) (b complex-bigfloat)) 631 (two-arg-* b a)) 632 633(defun * (&rest args) 634 (if (null args) 635 1 636 (do ((args (cdr args) (cdr args)) 637 (res (car args) 638 (two-arg-* res (car args)))) 639 ((null args) res)))) 640 641;;; Reciprocal of a number 642(defmethod unary-divide ((a number)) 643 (cl:/ a)) 644 645(defmethod unary-divide ((a bigfloat)) 646 (make-instance 'bigfloat 647 :real (maxima::bcons 648 (maxima::fpquotient (maxima::fpone) 649 (cdr (real-value a)))))) 650 651(defmethod unary-divide ((b complex-bigfloat)) 652 ;; Could just call two-arg-/, but let's optimize this a little 653 (let ((a-re (maxima::fpone)) 654 (b-re (cdr (real-value b))) 655 (b-im (cdr (imag-value b)))) 656 (if (maxima::fpgreaterp (maxima::fpabs b-re) (maxima::fpabs b-im)) 657 (let* ((r (maxima::fpquotient b-im b-re)) 658 (dn (maxima::fpplus b-re (maxima::fptimes* r b-im)))) 659 (make-instance 'complex-bigfloat 660 :real (maxima::bcons (maxima::fpquotient a-re dn)) 661 :imag (maxima::bcons 662 (maxima::fpquotient (maxima::fpminus r) 663 dn)))) 664 (let* ((r (maxima::fpquotient b-re b-im)) 665 (dn (maxima::fpplus b-im (maxima::fptimes* r b-re)))) 666 (make-instance 'complex-bigfloat 667 :real (maxima::bcons (maxima::fpquotient r dn)) 668 :imag (maxima::bcons 669 (maxima::fpquotient (maxima::fpminus a-re) 670 dn))))))) 671;;; Divide two numbers 672(defmethod two-arg-/ ((a number) (b number)) 673 (cl:/ a b)) 674 675(defmethod two-arg-/ ((a bigfloat) (b bigfloat)) 676 (make-instance 'bigfloat 677 :real (maxima::bcons 678 (maxima::fpquotient (cdr (real-value a)) 679 (cdr (real-value b)))))) 680 681(defmethod two-arg-/ ((a complex-bigfloat) (b complex-bigfloat)) 682 (let ((a-re (cdr (real-value a))) 683 (a-im (cdr (imag-value a))) 684 (b-re (cdr (real-value b))) 685 (b-im (cdr (imag-value b)))) 686 (if (maxima::fpgreaterp (maxima::fpabs b-re) (maxima::fpabs b-im)) 687 (let* ((r (maxima::fpquotient b-im b-re)) 688 (dn (maxima::fpplus b-re (maxima::fptimes* r b-im)))) 689 (make-instance 'complex-bigfloat 690 :real (maxima::bcons 691 (maxima::fpquotient 692 (maxima::fpplus a-re 693 (maxima::fptimes* a-im r)) 694 dn)) 695 :imag (maxima::bcons 696 (maxima::fpquotient 697 (maxima::fpdifference a-im 698 (maxima::fptimes* a-re r)) 699 dn)))) 700 (let* ((r (maxima::fpquotient b-re b-im)) 701 (dn (maxima::fpplus b-im (maxima::fptimes* r b-re)))) 702 (make-instance 'complex-bigfloat 703 :real (maxima::bcons 704 (maxima::fpquotient 705 (maxima::fpplus (maxima::fptimes* a-re r) 706 a-im) 707 dn)) 708 :imag (maxima::bcons 709 (maxima::fpquotient (maxima::fpdifference 710 (maxima::fptimes* a-im r) 711 a-re) 712 dn))))))) 713;; Handle contagion for two-arg-/ 714(defmethod two-arg-/ ((a bigfloat) (b cl:float)) 715 (make-instance 'bigfloat 716 :real (maxima::bcons 717 (maxima::fpquotient (cdr (real-value a)) 718 (cdr (intofp b)))))) 719 720(defmethod two-arg-/ ((a bigfloat) (b cl:rational)) 721 (make-instance 'bigfloat 722 :real (maxima::bcons 723 (maxima::fpquotient (cdr (real-value a)) 724 (cdr (intofp b)))))) 725 726(defmethod two-arg-/ ((a bigfloat) (b cl:complex)) 727 (two-arg-/ (complex a) b)) 728 729(defmethod two-arg-/ ((a cl:float) (b bigfloat)) 730 (make-instance 'bigfloat 731 :real (maxima::bcons 732 (maxima::fpquotient (cdr (intofp a)) 733 (cdr (real-value b)))))) 734 735(defmethod two-arg-/ ((a cl:rational) (b bigfloat)) 736 (make-instance 'bigfloat 737 :real (maxima::bcons 738 (maxima::fpquotient (cdr (intofp a)) 739 (cdr (real-value b)))))) 740 741(defmethod two-arg-/ ((a cl:complex) (b bigfloat)) 742 (two-arg-/ (bigfloat a) b)) 743 744 745(defmethod two-arg-/ ((a complex-bigfloat) (b bigfloat)) 746 (make-instance 'complex-bigfloat 747 :real (maxima::bcons 748 (maxima::fpquotient (cdr (real-value a)) 749 (cdr (real-value b)))) 750 :imag (maxima::bcons 751 (maxima::fpquotient (cdr (imag-value a)) 752 (cdr (real-value b)))))) 753 754(defmethod two-arg-/ ((a complex-bigfloat) (b number)) 755 (if (cl:complexp b) 756 (two-arg-/ a (bigfloat (cl:realpart b) (cl:imagpart b))) 757 (two-arg-/ a (bigfloat b)))) 758 759(defmethod two-arg-/ ((a bigfloat) (b complex-bigfloat)) 760 (two-arg-/ (make-instance 'complex-bigfloat :real (real-value a)) 761 b)) 762 763(defmethod two-arg-/ ((a number) (b complex-bigfloat)) 764 (if (cl:complexp a) 765 (two-arg-/ (bigfloat (cl:realpart a) (cl:imagpart a)) 766 b) 767 (two-arg-/ (bigfloat a) b))) 768 769 770(defun / (number &rest more-numbers) 771 (if more-numbers 772 (do ((nlist more-numbers (cdr nlist)) 773 (result number)) 774 ((atom nlist) result) 775 (declare (list nlist)) 776 (setq result (two-arg-/ result (car nlist)))) 777 (unary-divide number))) 778 779;;; Compare against zero (zerop, minusp, plusp) 780(macrolet 781 ((frob (name) 782 (let ((cl-name (intern (symbol-name name) :cl))) 783 `(progn 784 (defmethod ,name ((x cl:float)) 785 (,cl-name x)) 786 (defmethod ,name ((x cl:rational)) 787 (,cl-name x)))))) 788 (frob plusp) 789 (frob minusp)) 790 791(defmethod zerop ((x number)) 792 (cl:zerop x)) 793 794(defmethod zerop ((x bigfloat)) 795 (let ((r (cdr (real-value x)))) 796 (and (zerop (first r)) 797 (zerop (second r))))) 798 799(defmethod zerop ((a complex-bigfloat)) 800 (and (equal (cdr (real-value a)) '(0 0)) 801 (equal (cdr (imag-value a)) '(0 0)))) 802 803(defmethod plusp ((x bigfloat)) 804 (cl:plusp (first (cdr (real-value x))))) 805 806(defmethod minusp ((x bigfloat)) 807 (cl:minusp (first (cdr (real-value x))))) 808 809 810 811;;; Equality 812(defmethod two-arg-= ((a number) (b number)) 813 (cl:= a b)) 814 815(defmethod two-arg-= ((a bigfloat) (b bigfloat)) 816 (zerop (two-arg-- a b))) 817 818(defmethod two-arg-= ((a complex-bigfloat) (b complex-bigfloat)) 819 (zerop (two-arg-- a b))) 820 821;; Handle contagion for two-arg-=. This needs some work. CL says 822;; floats and rationals are compared by converting the float to a 823;; rational before converting. 824(defmethod two-arg-= ((a bigfloat) (b number)) 825 (zerop (two-arg-- a b))) 826 827(defmethod two-arg-= ((a number) (b bigfloat)) 828 (two-arg-= b a)) 829 830(defmethod two-arg-= ((a complex-bigfloat) (b number)) 831 (zerop (two-arg-- a b))) 832 833(defmethod two-arg-= ((a number) (b complex-bigfloat)) 834 (two-arg-= b a)) 835 836(defun = (number &rest more-numbers) 837 "Returns T if all of its arguments are numerically equal, NIL otherwise." 838 (declare (optimize (safety 2)) 839 #-gcl (dynamic-extent more-numbers)) 840 (do ((nlist more-numbers (cdr nlist))) 841 ((atom nlist) t) 842 (declare (list nlist)) 843 (if (not (two-arg-= (car nlist) number)) 844 (return nil)))) 845 846(defun /= (number &rest more-numbers) 847 "Returns T if no two of its arguments are numerically equal, NIL otherwise." 848 (declare (optimize (safety 2)) 849 #-gcl (dynamic-extent more-numbers)) 850 (do* ((head number (car nlist)) 851 (nlist more-numbers (cdr nlist))) 852 ((atom nlist) t) 853 (declare (list nlist)) 854 (unless (do* ((nl nlist (cdr nl))) 855 ((atom nl) t) 856 (declare (list nl)) 857 (if (two-arg-= head (car nl)) 858 (return nil))) 859 (return nil)))) 860 861;;; Comparison operations 862(macrolet 863 ((frob (op) 864 (let ((method (intern (concatenate 'string 865 (string '#:two-arg-) 866 (symbol-name op)))) 867 (cl-fun (find-symbol (symbol-name op) :cl))) 868 `(progn 869 (defmethod ,method ((a cl:float) (b cl:float)) 870 (,cl-fun a b)) 871 (defmethod ,method ((a cl:float) (b cl:rational)) 872 (,cl-fun a b)) 873 (defmethod ,method ((a cl:rational) (b cl:float)) 874 (,cl-fun a b)) 875 (defmethod ,method ((a cl:rational) (b cl:rational)) 876 (,cl-fun a b)) 877 (defun ,op (number &rest more-numbers) 878 "Returns T if its arguments are in strictly increasing order, NIL otherwise." 879 (declare (optimize (safety 2)) 880 #-gcl (dynamic-extent more-numbers)) 881 (do* ((n number (car nlist)) 882 (nlist more-numbers (cdr nlist))) 883 ((atom nlist) t) 884 (declare (list nlist)) 885 (if (not (,method n (car nlist))) (return nil)))))))) 886 (frob <) 887 (frob >) 888 (frob <=) 889 (frob >=)) 890 891(defmethod two-arg-< ((x bigfloat) (y bigfloat)) 892 (maxima::fplessp (cdr (real-value x)) (cdr (real-value y)))) 893 894(defmethod two-arg-< ((x bigfloat) (y cl:float)) 895 (maxima::fplessp (cdr (real-value x)) (cdr (intofp y)))) 896 897(defmethod two-arg-< ((x bigfloat) (y cl:rational)) 898 (maxima::fplessp (cdr (real-value x)) (cdr (intofp y)))) 899 900(defmethod two-arg-< ((x cl:float) (y bigfloat)) 901 (maxima::fplessp (cdr (intofp x)) (cdr (real-value y)))) 902 903(defmethod two-arg-< ((x cl:rational) (y bigfloat)) 904 (maxima::fplessp (cdr (intofp x)) (cdr (real-value y)))) 905 906(defmethod two-arg-> ((x bigfloat) (y bigfloat)) 907 (maxima::fpgreaterp (cdr (real-value x)) (cdr (real-value y)))) 908 909(defmethod two-arg-> ((x bigfloat) (y cl:float)) 910 (maxima::fpgreaterp (cdr (real-value x)) (cdr (intofp y)))) 911 912(defmethod two-arg-> ((x bigfloat) (y cl:rational)) 913 (maxima::fpgreaterp (cdr (real-value x)) (cdr (intofp y)))) 914 915(defmethod two-arg-> ((x cl:float) (y bigfloat)) 916 (maxima::fpgreaterp (cdr (intofp x)) (cdr (real-value y)))) 917 918(defmethod two-arg-> ((x cl:rational) (y bigfloat)) 919 (maxima::fpgreaterp (cdr (intofp x)) (cdr (real-value y)))) 920 921(defmethod two-arg-<= ((x bigfloat) (y bigfloat)) 922 (or (zerop (two-arg-- x y)) 923 (maxima::fplessp (cdr (real-value x)) (cdr (real-value y))))) 924 925(defmethod two-arg-<= ((x bigfloat) (y cl:float)) 926 (or (zerop (two-arg-- x y)) 927 (maxima::fplessp (cdr (real-value x)) (cdr (intofp y))))) 928 929(defmethod two-arg-<= ((x bigfloat) (y cl:rational)) 930 (or (zerop (two-arg-- x y)) 931 (maxima::fplessp (cdr (real-value x)) (cdr (intofp y))))) 932 933(defmethod two-arg-<= ((x cl:float) (y bigfloat)) 934 (or (zerop (two-arg-- x y)) 935 (maxima::fplessp (cdr (intofp x)) (cdr (real-value y))))) 936 937(defmethod two-arg-<= ((x cl:rational) (y bigfloat)) 938 (or (zerop (two-arg-- x y)) 939 (maxima::fplessp (cdr (intofp x)) (cdr (real-value y))))) 940 941(defmethod two-arg->= ((x bigfloat) (y bigfloat)) 942 (or (zerop (two-arg-- x y)) 943 (maxima::fpgreaterp (cdr (real-value x)) (cdr (real-value y))))) 944 945(defmethod two-arg->= ((x bigfloat) (y cl:float)) 946 (or (zerop (two-arg-- x y)) 947 (maxima::fpgreaterp (cdr (real-value x)) (cdr (intofp y))))) 948 949(defmethod two-arg->= ((x bigfloat) (y cl:rational)) 950 (or (zerop (two-arg-- x y)) 951 (maxima::fpgreaterp (cdr (real-value x)) (cdr (intofp y))))) 952 953(defmethod two-arg->= ((x cl:float) (y bigfloat)) 954 (or (zerop (two-arg-- x y)) 955 (maxima::fpgreaterp (cdr (intofp x)) (cdr (real-value y))))) 956 957(defmethod two-arg->= ((x cl:rational) (y bigfloat)) 958 (or (zerop (two-arg-- x y)) 959 (maxima::fpgreaterp (cdr (intofp x)) (cdr (real-value y))))) 960 961;; Need to define incf and decf to call our generic +/- methods. 962(defmacro incf (place &optional (delta 1) &environment env) 963 "The first argument is some location holding a number. This number is 964 incremented by the second argument, DELTA, which defaults to 1." 965 (multiple-value-bind (dummies vals newval setter getter) 966 (get-setf-expansion place env) 967 (let ((d (gensym))) 968 `(let* (,@(mapcar #'list dummies vals) 969 (,d ,delta) 970 (,(car newval) (+ ,getter ,d))) 971 ,setter)))) 972 973(defmacro decf (place &optional (delta 1) &environment env) 974 "The first argument is some location holding a number. This number is 975 decremented by the second argument, DELTA, which defaults to 1." 976 (multiple-value-bind (dummies vals newval setter getter) 977 (get-setf-expansion place env) 978 (let ((d (gensym))) 979 `(let* (,@(mapcar #'list dummies vals) 980 (,d ,delta) 981 (,(car newval) (- ,getter ,d))) 982 ,setter)))) 983 984 985 986;;; Special functions for real-valued arguments 987(macrolet 988 ((frob (name) 989 (let ((cl-name (intern (symbol-name name) :cl))) 990 `(progn 991 (defmethod ,name ((x number)) 992 (,cl-name x)))))) 993 (frob sqrt) 994 (frob abs) 995 (frob exp) 996 (frob sin) 997 (frob cos) 998 (frob tan) 999 (frob asin) 1000 (frob acos) 1001 (frob sinh) 1002 (frob cosh) 1003 (frob tanh) 1004 (frob asinh) 1005 (frob acosh) 1006 (frob atanh) 1007 ) 1008 1009(defmethod abs ((x bigfloat)) 1010 (make-instance 'bigfloat 1011 :real (maxima::bcons (maxima::fpabs (cdr (real-value x)))))) 1012 1013(defmethod abs ((z complex-bigfloat)) 1014 (let ((x (make-instance 'bigfloat :real (real-value z))) 1015 (y (make-instance 'bigfloat :real (imag-value z)))) 1016 ;; Bigfloats don't overflow, so we don't need anything special. 1017 (sqrt (+ (* x x) (* y y))))) 1018 1019(defmethod exp ((x bigfloat)) 1020 (make-instance 'bigfloat 1021 :real (maxima::bcons (maxima::fpexp (cdr (real-value x)))))) 1022 1023(defmethod sin ((x bigfloat)) 1024 (make-instance 'bigfloat 1025 :real (maxima::bcons (maxima::fpsin (cdr (real-value x)) t)))) 1026 1027(defmethod cos ((x bigfloat)) 1028 (make-instance 'bigfloat 1029 :real (maxima::bcons (maxima::fpsin (cdr (real-value x)) nil)))) 1030 1031(defmethod tan ((x bigfloat)) 1032 (let ((r (cdr (real-value x)))) 1033 (make-instance 'bigfloat 1034 :real (maxima::bcons 1035 (maxima::fpquotient (maxima::fpsin r t) 1036 (maxima::fpsin r nil)))))) 1037 1038(defmethod asin ((x bigfloat)) 1039 (bigfloat (maxima::fpasin (real-value x)))) 1040 1041(defmethod acos ((x bigfloat)) 1042 (bigfloat (maxima::fpacos (real-value x)))) 1043 1044 1045(defmethod sqrt ((x bigfloat)) 1046 (if (minusp x) 1047 (make-instance 'complex-bigfloat 1048 :real (intofp 0) 1049 :imag (maxima::bcons 1050 (maxima::fproot (maxima::bcons (maxima::fpabs (cdr (real-value x)))) 2))) 1051 (make-instance 'bigfloat 1052 :real (maxima::bcons 1053 (maxima::fproot (real-value x) 2))))) 1054 1055(defmethod sqrt ((z complex-bigfloat)) 1056 (multiple-value-bind (rx ry) 1057 (maxima::complex-sqrt (real-value z) 1058 (imag-value z)) 1059 (make-instance 'complex-bigfloat 1060 :real (maxima::bcons rx) 1061 :imag (maxima::bcons ry)))) 1062 1063(defmethod one-arg-log ((a number)) 1064 (cl:log a)) 1065 1066(defmethod one-arg-log ((a bigfloat)) 1067 (if (minusp a) 1068 (make-instance 'complex-bigfloat 1069 :real (maxima::bcons 1070 (maxima::fplog (maxima::fpabs (cdr (real-value a))))) 1071 :imag (maxima::bcons (maxima::fppi))) 1072 (make-instance 'bigfloat 1073 :real (maxima::bcons (maxima::fplog (cdr (real-value a))))))) 1074 1075(defmethod one-arg-log ((a complex-bigfloat)) 1076 (let* ((res (maxima::big-float-log (real-value a) 1077 (imag-value a)))) 1078 (bigfloat res))) 1079 1080(defmethod two-arg-log ((a number) (b number)) 1081 (cl:log a b)) 1082 1083(defmethod two-arg-log ((a numeric) (b numeric)) 1084 (two-arg-/ (one-arg-log a) (one-arg-log b))) 1085 1086(defmethod two-arg-log ((a numeric) (b cl:number)) 1087 (two-arg-/ (one-arg-log a) (one-arg-log (bigfloat b)))) 1088 1089(defmethod two-arg-log ((a cl:number) (b numeric)) 1090 (two-arg-/ (one-arg-log (bigfloat a)) (one-arg-log b))) 1091 1092(defun log (a &optional b) 1093 (if b 1094 (two-arg-log a b) 1095 (one-arg-log a))) 1096 1097(defmethod sinh ((x bigfloat)) 1098 (let ((r (real-value x))) 1099 (make-instance 'bigfloat :real (maxima::fpsinh r)))) 1100 1101(defmethod cosh ((x bigfloat)) 1102 (let ((r (real-value x))) 1103 (make-instance 'bigfloat :real (maxima::$bfloat `((maxima::%cosh) ,r))))) 1104 1105(defmethod tanh ((x bigfloat)) 1106 (let ((r (real-value x))) 1107 (make-instance 'bigfloat :real (maxima::fptanh r)))) 1108 1109(defmethod asinh ((x bigfloat)) 1110 (let ((r (real-value x))) 1111 (make-instance 'bigfloat :real (maxima::fpasinh r)))) 1112 1113(defmethod atanh ((x bigfloat)) 1114 (let ((r (maxima::big-float-atanh (real-value x)))) 1115 (if (maxima::bigfloatp r) 1116 (make-instance 'bigfloat :real r) 1117 (make-instance 'complex-bigfloat 1118 :real (maxima::$realpart r) 1119 :imag (maxima::$imagpart r))))) 1120 1121(defmethod acosh ((x bigfloat)) 1122 (let* ((r (real-value x)) 1123 (value (maxima::mevalp `((maxima::%acosh maxima::simp) 1124 ,r)))) 1125 (if (maxima::bigfloatp value) 1126 (make-instance 'bigfloat :real value) 1127 (make-instance 'complex-bigfloat 1128 :real (maxima::$realpart value) 1129 :imag (maxima::$imagpart value))))) 1130 1131;;; Complex arguments 1132 1133;;; Special functions for complex args 1134(macrolet 1135 ((frob (name &optional big-float-op-p) 1136 (if big-float-op-p 1137 (let ((big-op (intern (concatenate 'string 1138 (string '#:big-float-) 1139 (string name)) 1140 '#:maxima))) 1141 `(defmethod ,name ((a complex-bigfloat)) 1142 (let ((res (,big-op (real-value a) 1143 (imag-value a)))) 1144 (to res)))) 1145 (let ((max-op (intern (concatenate 'string "$" (string name)) '#:maxima))) 1146 `(defmethod ,name ((a complex-bigfloat)) 1147 ;; We should do something better than calling mevalp 1148 (let* ((arg (maxima::add (real-value a) 1149 (maxima::mul 'maxima::$%i (imag-value a)))) 1150 (result (maxima::mevalp `((,',max-op maxima::simp) ,arg)))) 1151 (to result))))))) 1152 (frob exp) 1153 (frob sin) 1154 (frob cos) 1155 (frob tan) 1156 (frob asin t) 1157 (frob acos t) 1158 (frob sinh) 1159 (frob cosh) 1160 (frob tanh t) 1161 (frob asinh t) 1162 (frob acosh) 1163 (frob atanh t)) 1164 1165(defmethod one-arg-atan ((a number)) 1166 (cl:atan a)) 1167 1168(defmethod one-arg-atan ((a bigfloat)) 1169 (make-instance 'bigfloat 1170 :real (maxima::bcons (maxima::fpatan (cdr (real-value a)))))) 1171 1172(defmethod one-arg-atan ((a complex-bigfloat)) 1173 ;; We should do something better than calling mevalp 1174 (let* ((arg (maxima::add (real-value a) 1175 (maxima::mul 'maxima::$%i (imag-value a)))) 1176 (result (maxima::mevalp `((maxima::%atan maxima::simp) ,arg)))) 1177 (make-instance 'complex-bigfloat 1178 :real (maxima::$realpart result) 1179 :imag (maxima::$imagpart result)))) 1180 1181;; Really want type real, but gcl doesn't like that. Define methods for rational and float 1182#-gcl 1183(defmethod two-arg-atan ((a real) (b real)) 1184 (cl:atan a b)) 1185 1186#+gcl 1187(progn 1188 (defmethod two-arg-atan ((a cl:float) (b cl:float)) 1189 (cl:atan a b)) 1190 (defmethod two-arg-atan ((a cl:rational) (b cl:rational)) 1191 (cl:atan a b)) 1192 (defmethod two-arg-atan ((a cl:float) (b cl:rational)) 1193 (cl:atan a b)) 1194 (defmethod two-arg-atan ((a cl:rational) (b cl:float)) 1195 (cl:atan a b)) 1196 ) 1197 1198(defmethod two-arg-atan ((a bigfloat) (b bigfloat)) 1199 (make-instance 'bigfloat 1200 :real (maxima::bcons 1201 (maxima::fpatan2 (cdr (real-value a)) 1202 (cdr (real-value b)))))) 1203 1204(defmethod two-arg-atan ((a bigfloat) (b cl:float)) 1205 (make-instance 'bigfloat 1206 :real (maxima::bcons (maxima::fpatan2 (cdr (real-value a)) 1207 (cdr (intofp b)))))) 1208 1209(defmethod two-arg-atan ((a bigfloat) (b cl:rational)) 1210 (make-instance 'bigfloat 1211 :real (maxima::bcons (maxima::fpatan2 (cdr (real-value a)) 1212 (cdr (intofp b)))))) 1213 1214(defmethod two-arg-atan ((a cl:float) (b bigfloat)) 1215 (make-instance 'bigfloat 1216 :real (maxima::bcons (maxima::fpatan2 (cdr (intofp a)) 1217 (cdr (real-value b)))))) 1218 1219(defmethod two-arg-atan ((a cl:rational) (b bigfloat)) 1220 (make-instance 'bigfloat 1221 :real (maxima::bcons (maxima::fpatan2 (cdr (intofp a)) 1222 (cdr (real-value b)))))) 1223 1224(defun atan (a &optional b) 1225 (if b 1226 (two-arg-atan a b) 1227 (one-arg-atan a))) 1228 1229(defmethod scale-float ((a cl:float) (n integer)) 1230 (cl:scale-float a n)) 1231 1232(defmethod scale-float ((a bigfloat) (n integer)) 1233 (if (cl:zerop (second (real-value a))) 1234 (make-instance 'bigfloat :real (maxima::bcons (list 0 0))) 1235 (destructuring-bind (marker mantissa exp) 1236 (real-value a) 1237 (declare (ignore marker)) 1238 (make-instance 'bigfloat :real (maxima::bcons (list mantissa (+ exp n))))))) 1239 1240(macrolet 1241 ((frob (name) 1242 (let ((cl-name (intern (string name) '#:cl))) 1243 `(defmethod ,name ((a number)) 1244 (,cl-name a))))) 1245 (frob realpart) 1246 (frob imagpart) 1247 (frob conjugate) 1248 (frob phase)) 1249 1250(macrolet 1251 ((frob (name) 1252 (let ((cl-name (intern (string name) '#:cl))) 1253 `(defmethod ,name ((a number) &optional (divisor 1)) 1254 (,cl-name a divisor))))) 1255 (frob floor) 1256 (frob ffloor) 1257 (frob ceiling) 1258 (frob fceiling) 1259 (frob truncate) 1260 (frob ftruncate) 1261 (frob round) 1262 (frob fround)) 1263 1264 1265(defmethod realpart ((a bigfloat)) 1266 (make-instance 'bigfloat :real (real-value a))) 1267 1268(defmethod realpart ((a complex-bigfloat)) 1269 (make-instance 'bigfloat :real (real-value a))) 1270 1271(defmethod imagpart ((a bigfloat)) 1272 (make-instance 'bigfloat :real (intofp 0))) 1273 1274(defmethod imagpart ((a complex-bigfloat)) 1275 (make-instance 'bigfloat :real (imag-value a))) 1276 1277(defmethod conjugate ((a bigfloat)) 1278 (make-instance 'bigfloat :real (real-value a))) 1279 1280(defmethod conjugate ((a complex-bigfloat)) 1281 (make-instance 'complex-bigfloat 1282 :real (real-value a) 1283 :imag (maxima::bcons (maxima::fpminus (cdr (imag-value a)))))) 1284 1285(defmethod cis ((a cl:float)) 1286 (cl:cis a)) 1287 1288(defmethod cis ((a cl:rational)) 1289 (cl:cis a)) 1290 1291(defmethod cis ((a bigfloat)) 1292 (make-instance 'complex-bigfloat 1293 :real (maxima::bcons (maxima::fpsin (cdr (real-value a)) nil)) 1294 :imag (maxima::bcons (maxima::fpsin (cdr (real-value a)) t)))) 1295 1296(defmethod phase ((a bigfloat)) 1297 (let ((r (cdr (real-value a)))) 1298 (if (cl:>= (car r) 0) 1299 (make-instance 'bigfloat :real (maxima::bcons (list 0 0))) 1300 (make-instance 'bigfloat :real (maxima::bcons (maxima::fppi)))))) 1301 1302(defmethod phase ((a complex-bigfloat)) 1303 (make-instance 'bigfloat 1304 :real (maxima::bcons (maxima::fpatan2 (cdr (imag-value a)) 1305 (cdr (real-value a)))))) 1306 1307(defun max (number &rest more-numbers) 1308 "Returns the greatest of its arguments." 1309 (declare (optimize (safety 2)) (type (or real bigfloat) number) 1310 #-gcl (dynamic-extent more-numbers)) 1311 (dolist (real more-numbers) 1312 (when (> real number) 1313 (setq number real))) 1314 number) 1315 1316(defun min (number &rest more-numbers) 1317 "Returns the least of its arguments." 1318 (declare (optimize (safety 2)) (type (or real bigfloat) number) 1319 #-gcl (dynamic-extent more-numbers)) 1320 (do ((nlist more-numbers (cdr nlist)) 1321 (result (the (or real bigfloat) number))) 1322 ((null nlist) (return result)) 1323 (declare (list nlist)) 1324 (if (< (car nlist) result) 1325 (setq result (car nlist))))) 1326 1327;; We really want a real type, but gcl doesn't like it, so use number 1328;; instead. 1329#-gcl 1330(defmethod one-arg-complex ((a real)) 1331 (cl:complex a)) 1332 1333#+gcl 1334(progn 1335(defmethod one-arg-complex ((a cl:float)) 1336 (cl:complex a)) 1337(defmethod one-arg-complex ((a cl:rational)) 1338 (cl:complex a)) 1339) 1340 1341(defmethod one-arg-complex ((a bigfloat)) 1342 (make-instance 'complex-bigfloat 1343 :real (real-value a) 1344 :imag (intofp 0))) 1345 1346#-gcl 1347(defmethod two-arg-complex ((a real) (b real)) 1348 (cl:complex a b)) 1349 1350#+gcl 1351(progn 1352(defmethod two-arg-complex ((a cl:float) (b cl:float)) 1353 (cl:complex a b)) 1354(defmethod two-arg-complex ((a cl:rational) (b cl:rational)) 1355 (cl:complex a b)) 1356(defmethod two-arg-complex ((a cl:float) (b cl:rational)) 1357 (cl:complex a b)) 1358(defmethod two-arg-complex ((a cl:rational) (b cl:float)) 1359 (cl:complex a b)) 1360) 1361 1362(defmethod two-arg-complex ((a bigfloat) (b bigfloat)) 1363 (make-instance 'complex-bigfloat 1364 :real (real-value a) 1365 :imag (real-value b))) 1366 1367(defmethod two-arg-complex ((a cl:float) (b bigfloat)) 1368 (make-instance 'complex-bigfloat 1369 :real (intofp a) 1370 :imag (real-value b))) 1371 1372(defmethod two-arg-complex ((a cl:rational) (b bigfloat)) 1373 (make-instance 'complex-bigfloat 1374 :real (intofp a) 1375 :imag (real-value b))) 1376 1377(defmethod two-arg-complex ((a bigfloat) (b cl:float)) 1378 (make-instance 'complex-bigfloat 1379 :real (real-value a) 1380 :imag (intofp b))) 1381 1382(defmethod two-arg-complex ((a bigfloat) (b cl:rational)) 1383 (make-instance 'complex-bigfloat 1384 :real (real-value a) 1385 :imag (intofp b))) 1386 1387(defun complex (a &optional b) 1388 (if b 1389 (two-arg-complex a b) 1390 (one-arg-complex a))) 1391 1392(defmethod unary-floor ((a bigfloat)) 1393 ;; fpentier truncates to zero, so adjust for negative numbers 1394 (let ((trunc (maxima::fpentier (real-value a)))) 1395 (cond ((minusp a) 1396 ;; If the truncated value is the same as the original, 1397 ;; there's nothing to do because A was an integer. 1398 ;; Otherwise, we need to subtract 1 to make it the floor. 1399 (if (= trunc a) 1400 trunc 1401 (1- trunc))) 1402 (t 1403 trunc)))) 1404 1405(defmethod unary-ffloor ((a bigfloat)) 1406 ;; We can probably do better than converting to an integer and 1407 ;; converting back to a float. 1408 (make-instance 'bigfloat :real (intofp (unary-floor a)))) 1409 1410(defmethod floor ((a bigfloat) &optional (divisor 1)) 1411 (if (= divisor 1) 1412 (let ((int (unary-floor a))) 1413 (values int (- a int))) 1414 (let ((q (unary-floor (/ a divisor)))) 1415 (values q (- a (* q divisor)))))) 1416 1417(defmethod ffloor ((a bigfloat) &optional (divisor 1)) 1418 (if (= divisor 1) 1419 (let ((int (unary-ffloor a))) 1420 (values int (- a int))) 1421 (let ((q (unary-ffloor (/ a divisor)))) 1422 (values q (- a (* q divisor)))))) 1423 1424(defmethod unary-truncate ((a bigfloat)) 1425 (maxima::fpentier (real-value a))) 1426 1427(defmethod unary-ftruncate ((a bigfloat)) 1428 ;; We can probably do better than converting to an integer and 1429 ;; converting back to a float. 1430 (make-instance 'bigfloat :real (intofp (unary-truncate a)))) 1431 1432(defmethod truncate ((a bigfloat) &optional (divisor 1)) 1433 (if (eql divisor 1) 1434 (let ((int (unary-truncate a))) 1435 (values int (- a int))) 1436 (let ((q (unary-truncate (/ a divisor)))) 1437 (values q (- a (* q divisor)))))) 1438 1439(defmethod ftruncate ((a bigfloat) &optional (divisor 1)) 1440 (if (eql divisor 1) 1441 (let ((int (unary-ftruncate a))) 1442 (values int (- a int))) 1443 (let ((q (unary-ftruncate (/ a divisor)))) 1444 (values q (- a (* q divisor)))))) 1445 1446(defmethod unary-ceiling ((a bigfloat)) 1447 ;; fpentier truncates to zero, so adjust for positive numbers. 1448 (if (minusp a) 1449 (maxima::fpentier (real-value a)) 1450 (maxima::fpentier (real-value (+ a 1))))) 1451 1452(defmethod unary-fceiling ((a bigfloat)) 1453 ;; We can probably do better than converting to an integer and 1454 ;; converting back to a float. 1455 (make-instance 'bigfloat :real (intofp (unary-ceiling a)))) 1456 1457(defmethod ceiling ((a bigfloat) &optional (divisor 1)) 1458 (if (eql divisor 1) 1459 (let ((int (unary-ceiling a))) 1460 (values int (- a int))) 1461 (let ((q (unary-ceiling (/ a divisor)))) 1462 (values q (- a (* q divisor)))))) 1463 1464(defmethod fceiling ((a bigfloat) &optional (divisor 1)) 1465 (if (eql divisor 1) 1466 (let ((int (unary-fceiling a))) 1467 (values int (- a int))) 1468 (let ((q (unary-fceiling (/ a divisor)))) 1469 (values q (- a (* q divisor)))))) 1470 1471;; Stolen from CMUCL. 1472(defmethod round ((a bigfloat) &optional (divisor 1)) 1473 (multiple-value-bind (tru rem) 1474 (truncate a divisor) 1475 (if (zerop rem) 1476 (values tru rem) 1477 (let ((thresh (/ (abs divisor) 2))) 1478 (cond ((or (> rem thresh) 1479 (and (= rem thresh) (oddp tru))) 1480 (if (minusp divisor) 1481 (values (- tru 1) (+ rem divisor)) 1482 (values (+ tru 1) (- rem divisor)))) 1483 ((let ((-thresh (- thresh))) 1484 (or (< rem -thresh) 1485 (and (= rem -thresh) (oddp tru)))) 1486 (if (minusp divisor) 1487 (values (+ tru 1) (- rem divisor)) 1488 (values (- tru 1) (+ rem divisor)))) 1489 (t (values tru rem))))))) 1490 1491(defmethod fround ((number bigfloat) &optional (divisor 1)) 1492 "Same as ROUND, but returns first value as a float." 1493 (multiple-value-bind (res rem) 1494 (round number divisor) 1495 (values (bigfloat res) rem))) 1496 1497(defmethod expt ((a number) (b number)) 1498 (cl:expt a b)) 1499 1500;; This needs more work 1501(defmethod expt ((a numeric) (b numeric)) 1502 (if (zerop b) 1503 ;; CLHS says if the power is 0, the answer is 1 of the appropriate type. 1504 (if (or (typep a 'complex-bigfloat) 1505 (typep b 'complex-bigfloat)) 1506 (complex (bigfloat 1)) 1507 (bigfloat 1)) 1508 (cond ((and (zerop a) (plusp (realpart b))) 1509 (* a b)) 1510 ((and (typep b 'bigfloat) (= b (truncate b))) 1511 ;; Use the numeric^number method because it can be much 1512 ;; more accurate when b is an integer. 1513 (expt a (truncate b))) 1514 (t 1515 (with-extra-precision ((expt-extra-bits a b) 1516 (a b)) 1517 (exp (* b (log a)))))))) 1518 1519(defmethod expt ((a cl:number) (b numeric)) 1520 (if (zerop b) 1521 ;; CLHS says if the power is 0, the answer is 1 of the appropriate type. 1522 (if (or (typep a 'cl:complex) 1523 (typep b 'complex-bigfloat)) 1524 (complex (bigfloat 1)) 1525 (bigfloat 1)) 1526 (cond ((and (zerop a) (plusp (realpart b))) 1527 (* a b)) 1528 ((= b (truncate b)) 1529 (with-extra-precision ((expt-extra-bits a b) 1530 (a b)) 1531 (intofp (expt a (truncate b))))) 1532 (t 1533 (with-extra-precision ((expt-extra-bits a b) 1534 (a b)) 1535 (exp (* b (log (bigfloat a))))))))) 1536 1537(defmethod expt ((a numeric) (b cl:number)) 1538 (if (zerop b) 1539 ;; CLHS says if the power is 0, the answer is 1 of the appropriate type. 1540 (if (or (typep a 'complex-bigfloat) 1541 (typep b 'cl:complex)) 1542 (complex (bigfloat 1)) 1543 (bigfloat 1)) 1544 (if (and (zerop a) (plusp (realpart b))) 1545 (* a b) 1546 ;; Handle a few special cases using multiplication. 1547 (cond ((= b 1) 1548 a) 1549 ((= b -1) 1550 (/ a)) 1551 ((= b 2) 1552 (* a a)) 1553 ((= b -2) 1554 (/ (* a a))) 1555 ((= b 3) (* a a a)) 1556 ((= b -3) (/ (* a a a))) 1557 ((= b 4) 1558 (let ((a2 (* a a))) 1559 (* a2 a2))) 1560 ((= b -4) 1561 (let ((a2 (* a a))) 1562 (/ (* a2 a2)))) 1563 (t 1564 (with-extra-precision ((expt-extra-bits a b) 1565 (a b)) 1566 (exp (* (bigfloat b) (log a))))))))) 1567 1568;; Handle a^b a little more carefully because the result is known to 1569;; be real when a is real and b is an integer. 1570(defmethod expt ((a bigfloat) (b integer)) 1571 (cond ((zerop b) 1572 (bigfloat 1)) 1573 ((and (zerop a) (plusp b)) 1574 ;; 0^b, for positive b 1575 (* a b)) 1576 ;; Handle a few special cases using multiplication. 1577 ((eql b 1) a) 1578 ((eql b -1) (/ a)) 1579 ((eql b 2) (* a a)) 1580 ((eql b -2) (/ (* a a))) 1581 ((eql b 3) (* a a a)) 1582 ((eql b -3) (/ (* a a a))) 1583 ((eql b 4) 1584 (let ((a2 (* a a))) 1585 (* a2 a2))) 1586 ((eql b -4) 1587 (let ((a2 (* a a))) 1588 (/ (* a2 a2)))) 1589 ((minusp a) 1590 ;; a^b = exp(b*log(|a|) + %i*%pi*b) 1591 ;; = exp(b*log(|a|))*exp(%i*%pi*b) 1592 ;; = (-1)^b*exp(b*log(|a|)) 1593 (with-extra-precision ((expt-extra-bits a b) 1594 (a b)) 1595 (* (exp (* b (log (abs a)))) 1596 (if (oddp b) -1 1)))) 1597 (t 1598 (with-extra-precision ((expt-extra-bits a b) 1599 (a b)) 1600 (exp (* b (log a))))))) 1601 1602;;; TO - External 1603;;; 1604;;; TO takes a maxima number and converts it. Floats remain 1605;;; floats, maxima cl:rationals are converted to CL cl:rationals. Maxima 1606;;; bigfloats are convert to BIGFLOATS. Maxima complex numbers are 1607;;; converted to CL complex numbers or to COMPLEX-BIGFLOAT's. 1608(defun to (maxima-num &optional imag) 1609 (let ((result (ignore-errors (%to maxima-num imag)))) 1610 (or result 1611 (maxima::merror (intl:gettext "BIGFLOAT: unable to convert ~M to a CL or BIGFLOAT number.") 1612 (if imag 1613 (maxima::add maxima-num (maxima::mul imag 'maxima::$%i)) 1614 maxima-num))))) 1615 1616;;; MAYBE-TO - External 1617;;; 1618;;; Like TO, but if the maxima number can't be converted to a CL 1619;;; number or BIGFLOAT, just return the maxima number. 1620(defun maybe-to (maxima-num &optional imag) 1621 (let ((result (ignore-errors (%to maxima-num imag)))) 1622 (or result 1623 (if imag 1624 (maxima::add maxima-num imag) 1625 maxima-num)))) 1626 1627(defun %to (maxima-num &optional imag) 1628 (cond (imag 1629 ;; Clisp has a "feature" that (complex rat float) does not 1630 ;; make the both components of the complex number a float. 1631 ;; Sometimes this is nice, but other times it's annoying 1632 ;; because it is non-ANSI behavior. For our code, we really 1633 ;; want both components to be a float. 1634 #-clisp 1635 (complex (to maxima-num) (to imag)) 1636 #+clisp 1637 (let ((re (to maxima-num)) 1638 (im (to imag))) 1639 (cond ((and (rationalp re) (floatp im)) 1640 (setf re (float re im))) 1641 ((and (rational im) (floatp re)) 1642 (setf im (float im re)))) 1643 (complex re im))) 1644 (t 1645 (cond ((cl:numberp maxima-num) 1646 maxima-num) 1647 ((eq maxima-num 'maxima::$%i) 1648 ;; Convert %i to an exact complex cl:rational. 1649 #c(0 1)) 1650 ((consp maxima-num) 1651 ;; Some kind of maxima number 1652 (cond ((maxima::ratnump maxima-num) 1653 ;; Maxima cl:rational ((mrat ...) num den) 1654 (/ (second maxima-num) (third maxima-num))) 1655 ((maxima::$bfloatp maxima-num) 1656 (bigfloat maxima-num)) 1657 ((maxima::complex-number-p maxima-num #'(lambda (x) 1658 (or (cl:realp x) 1659 (maxima::$bfloatp x) 1660 (and (consp x) 1661 (eq (caar x) 'maxima::rat))))) 1662 ;; We have some kind of complex number whose 1663 ;; parts are a cl:real, a bfloat, or a Maxima 1664 ;; cl:rational. 1665 (let ((re (maxima::$realpart maxima-num)) 1666 (im (maxima::$imagpart maxima-num))) 1667 (to re im))))) 1668 ((or (typep maxima-num 'bigfloat) 1669 (typep maxima-num 'complex-bigfloat)) 1670 maxima-num) 1671 (t 1672 (error "BIGFLOAT: unable to convert to a CL or BIGFLOAT number.")))))) 1673 1674;;; EPSILON - External 1675;;; 1676;;; Return the float epsilon value for the given float type. 1677(defmethod epsilon ((x cl:float)) 1678 (etypecase x 1679 (short-float short-float-epsilon) 1680 (single-float single-float-epsilon) 1681 (double-float double-float-epsilon) 1682 (long-float long-float-epsilon))) 1683 1684(defmethod epsilon ((x cl:complex)) 1685 (epsilon (cl:realpart x))) 1686 1687(defmethod epsilon ((x bigfloat)) 1688 ;; epsilon is just above 2^(-fpprec). 1689 (make-instance 'bigfloat 1690 :real (maxima::bcons (list (1+ (ash 1 (1- maxima::fpprec))) 1691 (- (1- maxima::fpprec)))))) 1692 1693(defmethod epsilon ((x complex-bigfloat)) 1694 (epsilon (realpart x))) 1695 1696 1697 1698;; Compiler macros to convert + to multiple calls to two-arg-+. Same 1699;; for -, *, and /. 1700(define-compiler-macro + (&whole form &rest args) 1701 (declare (ignore form)) 1702 (if (null args) 1703 0 1704 (do ((args (cdr args) (cdr args)) 1705 (res (car args) 1706 `(two-arg-+ ,res ,(car args)))) 1707 ((null args) res)))) 1708 1709(define-compiler-macro - (&whole form number &rest more-numbers) 1710 (declare (ignore form)) 1711 (if more-numbers 1712 (do ((nlist more-numbers (cdr nlist)) 1713 (result number)) 1714 ((atom nlist) result) 1715 (declare (list nlist)) 1716 (setq result `(two-arg-- ,result ,(car nlist)))) 1717 `(unary-minus ,number))) 1718 1719(define-compiler-macro * (&whole form &rest args) 1720 (declare (ignore form)) 1721 (if (null args) 1722 1 1723 (do ((args (cdr args) (cdr args)) 1724 (res (car args) 1725 `(two-arg-* ,res ,(car args)))) 1726 ((null args) res)))) 1727 1728(define-compiler-macro / (number &rest more-numbers) 1729 (if more-numbers 1730 (do ((nlist more-numbers (cdr nlist)) 1731 (result number)) 1732 ((atom nlist) result) 1733 (declare (list nlist)) 1734 (setq result `(two-arg-/ ,result ,(car nlist)))) 1735 `(unary-divide ,number))) 1736 1737(define-compiler-macro /= (&whole form number &rest more-numbers) 1738 ;; Convert (/= x y) to (not (two-arg-= x y)). Should we try to 1739 ;; handle a few more cases? 1740 (if (cdr more-numbers) 1741 form 1742 `(not (two-arg-= ,number ,(car more-numbers))))) 1743 1744;; Compiler macros to convert <, >, <=, and >= into multiple calls of 1745;; the corresponding two-arg-<foo> function. 1746(macrolet 1747 ((frob (op) 1748 (let ((method (intern (concatenate 'string 1749 (string '#:two-arg-) 1750 (symbol-name op))))) 1751 `(define-compiler-macro ,op (number &rest more-numbers) 1752 (do* ((n number (car nlist)) 1753 (nlist more-numbers (cdr nlist)) 1754 (res nil)) 1755 ((atom nlist) 1756 `(and ,@(nreverse res))) 1757 (push `(,',method ,n ,(car nlist)) res)))))) 1758 (frob <) 1759 (frob >) 1760 (frob <=) 1761 (frob >=)) 1762 1763(defmethod integer-decode-float ((x cl:float)) 1764 (cl:integer-decode-float x)) 1765 1766(defmethod integer-decode-float ((x bigfloat)) 1767 (let ((r (real-value x))) 1768 (values (abs (second r)) 1769 (- (third r) (third (first r))) 1770 (signum (second r))))) 1771 1772(defmethod decode-float ((x cl:float)) 1773 (cl:decode-float x)) 1774 1775(defmethod decode-float ((x bigfloat)) 1776 (let ((r (real-value x))) 1777 (values (make-instance 'bigfloat 1778 :real (maxima::bcons (list (abs (second r)) 0))) 1779 (third r) 1780 (bigfloat (signum (second r)))))) 1781 1782;; GCL doesn't have a REAL class! 1783#+gcl 1784(progn 1785(defmethod float ((x cl:float) (y cl:float)) 1786 (cl:float x y)) 1787 1788(defmethod float ((x cl:rational) (y cl:float)) 1789 (cl:float x y)) 1790 1791(defmethod float ((x cl:float) (y bigfloat)) 1792 (bigfloat x)) 1793 1794(defmethod float ((x cl:rational) (y bigfloat)) 1795 (bigfloat x)) 1796) 1797 1798#-gcl 1799(progn 1800(defmethod float ((x real) (y cl:float)) 1801 (cl:float x y)) 1802 1803(defmethod float ((x real) (y bigfloat)) 1804 (bigfloat x)) 1805) 1806 1807;; Like Maxima's fp2flo, but for single-float numbers. 1808(defun fp2single (l) 1809 (let ((precision (caddar l)) 1810 (mantissa (cadr l)) 1811 (exponent (caddr l)) 1812 (fpprec (float-digits 1f0)) 1813 (maxima::*m 0)) 1814 ;; Round the mantissa to the number of bits of precision of the 1815 ;; machine, and then convert it to a floating point fraction. We 1816 ;; have 0.5 <= mantissa < 1 1817 (setq mantissa (maxima::quotient (maxima::fpround mantissa) 1818 (expt 2f0 fpprec))) 1819 ;; Multiply the mantissa by the exponent portion. I'm not sure 1820 ;; why the exponent computation is so complicated. 1821 ;; 1822 ;; GCL doesn't signal overflow from scale-float if the number 1823 ;; would overflow. We have to do it this way. 0.5 <= mantissa < 1824 ;; 1. The largest double-float is .999999 * 2^128. So if the 1825 ;; exponent is 128 or higher, we have an overflow. 1826 (let ((e (+ exponent (- precision) maxima::*m fpprec))) 1827 (if (>= (abs e) 129) 1828 (maxima::merror (intl:gettext "FP2SINGLE: floating point overflow converting ~:M to float.") l) 1829 (cl:scale-float mantissa e))))) 1830 1831 1832(defmethod float ((x bigfloat) (y cl:float)) 1833 (if (typep y 'maxima::flonum) 1834 (maxima::fp2flo (real-value x)) 1835 (fp2single (real-value x)))) 1836 1837(defmethod random ((x cl:float) &optional (state cl:*random-state*)) 1838 (cl:random x state)) 1839(defmethod random ((x integer) &optional (state cl:*random-state*)) 1840 (cl:random x state)) 1841 1842(defmethod random ((x bigfloat) &optional (state cl:*random-state*)) 1843 ;; Generate an integer with fpprec bits, and convert to a bigfloat 1844 ;; by making the exponent 0. Then multiply by the arg to get the 1845 ;; correct range. 1846 (if (plusp x) 1847 (let ((int (cl:random (ash 1 maxima::fpprec) state))) 1848 (* x (bigfloat (maxima::bcons (list int 0))))) 1849 (error "Argument is not a positive bigfloat: ~A~%" x))) 1850 1851(defmethod signum ((x number)) 1852 (cl:signum x)) 1853 1854(defmethod signum ((x bigfloat)) 1855 (cond ((minusp x) 1856 (bigfloat -1)) 1857 ((plusp x) 1858 (bigfloat 1)) 1859 (t 1860 x))) 1861 1862(defmethod signum ((x complex-bigfloat)) 1863 (/ x (abs x))) 1864 1865(defmethod float-sign ((x cl:float)) 1866 (cl:float-sign x)) 1867 1868(defmethod float-sign ((x bigfloat)) 1869 (if (minusp x) 1870 (bigfloat -1) 1871 (bigfloat 1))) 1872 1873(defmethod float-digits ((x cl:float)) 1874 (cl:float-digits x)) 1875 1876(defmethod float-digits ((x bigfloat)) 1877 ;; Should we just return fpprec or should we get the actual number 1878 ;; of bits in the bigfloat number? We choose the latter in case the 1879 ;; number and fpprec don't match. 1880 (let ((r (slot-value x 'real))) 1881 (third (first r)))) 1882 1883#-gcl 1884(defmethod rational ((x real)) 1885 (cl:rational x)) 1886 1887#+gcl 1888(progn 1889(defmethod rational ((x cl:float)) 1890 (cl:rational x)) 1891(defmethod rational ((x cl:rational)) 1892 (cl:rational x)) 1893) 1894 1895(defmethod rational ((x bigfloat)) 1896 (destructuring-bind ((marker simp prec) mantissa exp) 1897 (real-value x) 1898 (declare (ignore marker simp)) 1899 (* mantissa (expt 2 (- exp prec))))) 1900 1901#-gcl 1902(defmethod rationalize ((x real)) 1903 (cl:rationalize x)) 1904 1905#+gcl 1906(progn 1907(defmethod rationalize ((x cl:float)) 1908 (cl:rationalize x)) 1909(defmethod rationalize ((x cl:rational)) 1910 (cl:rationalize x)) 1911) 1912 1913 1914;;; This routine taken from CMUCL, which, in turn is a routine from 1915;;; CLISP, which is GPL. 1916;;; 1917;;; I (rtoy) have modified it from CMUCL so that it only handles bigfloats. 1918;;; 1919;;; RATIONALIZE -- Public 1920;;; 1921;;; The algorithm here is the method described in CLISP. Bruno Haible has 1922;;; graciously given permission to use this algorithm. He says, "You can use 1923;;; it, if you present the following explanation of the algorithm." 1924;;; 1925;;; Algorithm (recursively presented): 1926;;; If x is a rational number, return x. 1927;;; If x = 0.0, return 0. 1928;;; If x < 0.0, return (- (rationalize (- x))). 1929;;; If x > 0.0: 1930;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa, 1931;;; exponent, sign). 1932;;; If m = 0 or e >= 0: return x = m*2^e. 1933;;; Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e 1934;;; with smallest possible numerator and denominator. 1935;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e. 1936;;; But in this case the result will be x itself anyway, regardless of 1937;;; the choice of a. Therefore we can simply ignore this case. 1938;;; Note 2: At first, we need to consider the closed interval [a,b]. 1939;;; but since a and b have the denominator 2^(|e|+1) whereas x itself 1940;;; has a denominator <= 2^|e|, we can restrict the seach to the open 1941;;; interval (a,b). 1942;;; So, for given a and b (0 < a < b) we are searching a rational number 1943;;; y with a <= y <= b. 1944;;; Recursive algorithm fraction_between(a,b): 1945;;; c := (ceiling a) 1946;;; if c < b 1947;;; then return c ; because a <= c < b, c integer 1948;;; else 1949;;; ; a is not integer (otherwise we would have had c = a < b) 1950;;; k := c-1 ; k = floor(a), k < a < b <= k+1 1951;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k)) 1952;;; ; note 1 <= 1/(b-k) < 1/(a-k) 1953;;; 1954;;; You can see that we are actually computing a continued fraction expansion. 1955;;; 1956;;; Algorithm (iterative): 1957;;; If x is rational, return x. 1958;;; Call (integer-decode-float x). It returns a m,e,s (mantissa, 1959;;; exponent, sign). 1960;;; If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.) 1961;;; Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1) 1962;;; (positive and already in lowest terms because the denominator is a 1963;;; power of two and the numerator is odd). 1964;;; Start a continued fraction expansion 1965;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0. 1966;;; Loop 1967;;; c := (ceiling a) 1968;;; if c >= b 1969;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)), 1970;;; goto Loop 1971;;; finally partial_quotient(c). 1972;;; Here partial_quotient(c) denotes the iteration 1973;;; i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2]. 1974;;; At the end, return s * (p[i]/q[i]). 1975;;; This rational number is already in lowest terms because 1976;;; p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i. 1977;;; 1978(defmethod rationalize ((x bigfloat)) 1979 (multiple-value-bind (frac expo sign) 1980 (integer-decode-float x) 1981 (cond ((or (zerop frac) (>= expo 0)) 1982 (if (minusp sign) 1983 (- (ash frac expo)) 1984 (ash frac expo))) 1985 (t 1986 ;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e), 1987 ;; so build the fraction up immediately, without having to do 1988 ;; a gcd. 1989 (let ((a (/ (- (* 2 frac) 1) (ash 1 (- 1 expo)))) 1990 (b (/ (+ (* 2 frac) 1) (ash 1 (- 1 expo)))) 1991 (p0 0) 1992 (q0 1) 1993 (p1 1) 1994 (q1 0)) 1995 (do ((c (ceiling a) (ceiling a))) 1996 ((< c b) 1997 (let ((top (+ (* c p1) p0)) 1998 (bot (+ (* c q1) q0))) 1999 (/ (if (minusp sign) 2000 (- top) 2001 top) 2002 bot))) 2003 (let* ((k (- c 1)) 2004 (p2 (+ (* k p1) p0)) 2005 (q2 (+ (* k q1) q0))) 2006 (psetf a (/ (- b k)) 2007 b (/ (- a k))) 2008 (setf p0 p1 2009 q0 q1 2010 p1 p2 2011 q1 q2)))))))) 2012 2013(defun coerce (obj type) 2014 (flet ((coerce-error () 2015 (error "Cannot coerce ~A to type ~S" obj type))) 2016 (cond ((typep obj type) 2017 obj) 2018 ((subtypep type 'bigfloat) 2019 ;; (coerce foo 'bigfloat). Foo has to be a real 2020 (cond ((typep obj 'real) 2021 (bigfloat obj)) 2022 (t 2023 (coerce-error)))) 2024 ((subtypep type 'complex-bigfloat) 2025 ;; (coerce foo 'complex-bigfloat). Foo has to be a real or complex 2026 (cond ((typep obj 'real) 2027 (bigfloat obj 0)) 2028 ((typep obj 'cl:complex) 2029 (bigfloat obj)) 2030 ((typep obj 'bigfloat) 2031 (bigfloat obj 0)) 2032 (t 2033 (coerce-error)))) 2034 ((typep obj 'bigfloat) 2035 ;; (coerce bigfloat foo) 2036 (cond ((subtypep type 'cl:float) 2037 (float obj (cl:coerce 0 type))) 2038 ((subtypep type '(cl:complex long-float)) 2039 (cl:complex (float (realpart obj) 1l0) 2040 (float (imagpart obj) 1l0))) 2041 ((subtypep type '(cl:complex double-float)) 2042 (cl:complex (float (realpart obj) 1d0) 2043 (float (imagpart obj) 1d0))) 2044 ((subtypep type '(cl:complex single-float)) 2045 (cl:complex (float (realpart obj) 1f0) 2046 (float (imagpart obj) 1f0))) 2047 ((subtypep type 'cl:complex) 2048 ;; What should we do here? Return a 2049 ;; complex-bigfloat? A complex double-float? 2050 ;; complex single-float? I arbitrarily select 2051 ;; complex maxima:flonum for now. 2052 (cl:complex (float (realpart obj) 1.0) 2053 (float (imagpart obj) 1.0))) 2054 (t 2055 (coerce-error)))) 2056 ((typep obj 'complex-bigfloat) 2057 ;; (coerce complex-bigfloat foo) 2058 (cond ((subtypep type 'complex-bigfloat) 2059 obj) 2060 ((subtypep type '(cl:complex long-float)) 2061 (cl:complex (float (realpart obj) 1l0) 2062 (float (imagpart obj) 1l0))) 2063 ((subtypep type '(cl:complex double-float)) 2064 (cl:complex (float (realpart obj) 1d0) 2065 (float (imagpart obj) 1d0))) 2066 ((subtypep type '(cl:complex single-float)) 2067 (cl:complex (float (realpart obj) 1f0) 2068 (float (imagpart obj) 1f0))) 2069 (t 2070 (coerce-error)))) 2071 (t 2072 (cl:coerce obj type))))) 2073 2074;;; %PI - External 2075;;; 2076;;; Return a value of pi with the same precision as the argument. 2077;;; For rationals, we return a single-float approximation. 2078(defmethod %pi ((x cl:rational)) 2079 (cl:coerce cl:pi 'single-float)) 2080 2081(defmethod %pi ((x cl:float)) 2082 (cl:float cl:pi x)) 2083 2084(defmethod %pi ((x bigfloat)) 2085 (to (maxima::bcons (maxima::fppi)))) 2086 2087(defmethod %pi ((x cl:complex)) 2088 (cl:float cl:pi (realpart x))) 2089 2090(defmethod %pi ((x complex-bigfloat)) 2091 (to (maxima::bcons (maxima::fppi)))) 2092 2093;;; %e - External 2094;;; 2095;;; Return a value of e with the same precision as the argument. 2096;;; For rationals, we return a single-float approximation. 2097(defmethod %e ((x cl:rational)) 2098 (cl:coerce maxima::%e-val 'single-float)) 2099 2100(defmethod %e ((x cl:float)) 2101 (cl:float maxima::%e-val x)) 2102 2103(defmethod %e ((x bigfloat)) 2104 (to (maxima::bcons (maxima::fpe)))) 2105 2106(defmethod %e ((x cl:complex)) 2107 (cl:float maxima::%e-val (realpart x))) 2108 2109(defmethod %e ((x complex-bigfloat)) 2110 (to (maxima::bcons (maxima::fpe)))) 2111 2112;;;; Useful routines 2113 2114;;; Evaluation of continued fractions 2115 2116(defvar *debug-cf-eval* 2117 nil 2118 "When true, enable some debugging prints when evaluating a 2119 continued fraction.") 2120 2121;; Max number of iterations allowed when evaluating the continued 2122;; fraction. When this is reached, we assume that the continued 2123;; fraction did not converge. 2124(defvar *max-cf-iterations* 2125 10000 2126 "Max number of iterations allowed when evaluating the continued 2127 fraction. When this is reached, we assume that the continued 2128 fraction did not converge.") 2129 2130;;; LENTZ - External 2131;;; 2132;;; Lentz's algorithm for evaluating continued fractions. 2133;;; 2134;;; Let the continued fraction be: 2135;;; 2136;;; a1 a2 a3 2137;;; b0 + ---- ---- ---- 2138;;; b1 + b2 + b3 + 2139;;; 2140;;; 2141;;; Then LENTZ expects two functions, each taking a single fixnum 2142;;; index. The first returns the b term and the second returns the a 2143;;; terms as above for a give n. 2144(defun lentz (bf af) 2145 (let ((tiny-value-count 0)) 2146 (flet ((value-or-tiny (v) 2147 ;; If v is zero, return a "tiny" number. 2148 (if (zerop v) 2149 (progn 2150 (incf tiny-value-count) 2151 (etypecase v 2152 ((or double-float cl:complex) 2153 (sqrt least-positive-normalized-double-float)) 2154 ((or bigfloat complex-bigfloat) 2155 ;; What is a "tiny" bigfloat? Bigfloats have 2156 ;; unbounded exponents, so we need something 2157 ;; small, but not zero. Arbitrarily choose an 2158 ;; exponent of 50 times the precision. 2159 (expt 10 (- (* 50 maxima::$fpprec)))))) 2160 v))) 2161 (let* ((f (value-or-tiny (funcall bf 0))) 2162 (c f) 2163 (d 0) 2164 (eps (epsilon f))) 2165 (loop 2166 for j from 1 upto *max-cf-iterations* 2167 for an = (funcall af j) 2168 for bn = (funcall bf j) 2169 do (progn 2170 (setf d (value-or-tiny (+ bn (* an d)))) 2171 (setf c (value-or-tiny (+ bn (/ an c)))) 2172 (when *debug-cf-eval* 2173 (format t "~&j = ~d~%" j) 2174 (format t " an = ~s~%" an) 2175 (format t " bn = ~s~%" bn) 2176 (format t " c = ~s~%" c) 2177 (format t " d = ~s~%" d)) 2178 (let ((delta (/ c d))) 2179 (setf d (/ d)) 2180 (setf f (* f delta)) 2181 (when *debug-cf-eval* 2182 (format t " dl= ~S (|dl - 1| = ~S)~%" delta (abs (1- delta))) 2183 (format t " f = ~S~%" f)) 2184 (when (<= (abs (- delta 1)) eps) 2185 (return-from lentz (values f j tiny-value-count))))) 2186 finally 2187 (error 'simple-error 2188 :format-control "~<Continued fraction failed to converge after ~D iterations.~% Delta = ~S~>" 2189 :format-arguments (list *max-cf-iterations* (/ c d)))))))) 2190 2191;;; SUM-POWER-SERIES - External 2192;;; 2193;;; SUM-POWER-SERIES sums the given power series, adding terms until 2194;;; the next term would not change the sum. 2195;;; 2196;;; The series to be summed is 2197;;; 2198;;; S = 1 + sum(c[k]*x^k, k, 1, inf) 2199;;; = 1 + sum(prod(f[n]*x, n, 1, k), k, 1, inf) 2200;;; 2201;;; where f[n] = c[n]/c[n-1]. 2202;;; 2203(defun sum-power-series (x f) 2204 (let ((eps (epsilon x))) 2205 (do* ((k 1 (+ 1 k)) 2206 (sum 1 (+ sum term)) 2207 (term (* x (funcall f 1)) 2208 (* term x (funcall f k)))) 2209 ((< (abs term) (* eps (abs sum))) 2210 sum) 2211 #+nil 2212 (format t "~4d: ~S ~S ~S~%" k sum term (funcall f k))))) 2213 2214;; Format bigfloats using ~E format. This is suitable as a ~// format. 2215;; 2216;; NOTE: This is a modified version of FORMAT-EXPONENTIAL from CMUCL to 2217;; support printing of bfloats. 2218 2219(defun format-e (stream number colonp atp 2220 &optional w d e k 2221 overflowchar padchar exponentchar) 2222 (typecase number 2223 (bigfloat 2224 (maxima::bfloat-format-e stream (real-value number) colonp atp 2225 w d e (or k 1) 2226 overflowchar 2227 (or padchar #\space) 2228 (or exponentchar #\b))) 2229 (complex-bigfloat 2230 ;; FIXME: Do something better than this since this doesn't honor 2231 ;; any of the parameters. 2232 (princ number stream)) 2233 (otherwise 2234 ;; We were given some other kind of object. Just use CL's normal 2235 ;; ~E printer to print it. 2236 (let ((f 2237 (with-output-to-string (s) 2238 ;; Construct a suitable ~E format string from the given 2239 ;; parameters. First, handle w,d,e,k. 2240 (write-string "~V,V,V,V," s) 2241 (if overflowchar 2242 (format s "'~C," overflowchar) 2243 (write-string "," s)) 2244 (if padchar 2245 (format s "'~C," padchar) 2246 (write-string "," s)) 2247 (when exponentchar 2248 (format s "'~C" exponentchar)) 2249 (when colonp 2250 (write-char #\: s)) 2251 (when atp 2252 (write-char #\@ s)) 2253 (write-char #\E s)))) 2254 (format stream f w d e k number))))) 2255 2256#| 2257(defmacro assert-equal (expected form) 2258 (let ((result (gensym)) 2259 (e (gensym))) 2260 `(let ((,e ,expected) 2261 (,result ,form)) 2262 (unless (equal ,e ,result) 2263 (format *debug-io* "Assertion failed: Expected ~S but got ~S~%" ,e ,result))))) 2264 2265(assert-equal " 0.990E+00" (format nil 2266 "~11,3,2,0,'*,,'E/bigfloat::format-e/" 2267 (bigfloat:bigfloat 99/100))) 2268(assert-equal " 0.999E+00" (format nil 2269 "~11,3,2,0,'*,,'E/bigfloat::format-e/" 2270 (bigfloat:bigfloat 999/1000))) 2271;; Actually " 0.100E+01", but format-e doesn't round the output. 2272(assert-equal " 0.999E+00" (format nil 2273 "~11,3,2,0,'*,,'E/bigfloat::format-e/" 2274 (bigfloat:bigfloat 9999/10000))) 2275(assert-equal " 0.999E-04" (format nil 2276 "~11,3,2,0,'*,,'E/bigfloat::format-e/" 2277 (bigfloat:bigfloat 0000999/10000000))) 2278;; Actually " 0.100E-03", but format-e doesn't round the output. 2279(assert-equal " 0.999E-0e" (format nil 2280 "~11,3,2,0,'*,,'E/bigfloat::format-e/" 2281 (bigfloat:bigfloat 00009999/100000000))) 2282(assert-equal " 9.999E-05" (format nil 2283 "~11,3,2,,'*,,'E/bigfloat::format-e/" 2284 (bigfloat:bigfloat 00009999/100000000))) 2285;; Actually " 1.000E-04", but format-e doesn't round the output. 2286(assert-equal " 9.999E-05" (format nil 2287 "~11,3,2,,'*,,'E/bigfloat::format-e/" 2288 (bigfloat:bigfloat 000099999/1000000000))) 2289;; All of these currently fail. 2290(assert-equal ".00123d+6" (format nil 2291 "~9,,,-2/bigfloat::format-e/" 2292 (bigfloat:bigfloat 1.2345689d3))) 2293(assert-equal "-.0012d+6" (format nil 2294 "~9,,,-2/bigfloat::format-e/" 2295 (bigfloat:bigfloat -1.2345689d3))) 2296(assert-equal ".00123d+0" (format nil 2297 "~9,,,-2/bigfloat::format-e/" 2298 (bigfloat:bigfloat 1.2345689d-3))) 2299(assert-equal "-.0012d+0" (format nil 2300 "~9,,,-2/bigfloat::format-e/" 2301 (bigfloat:bigfloat -1.2345689d-3))) 2302 2303;; These fail because too many digits are printed and because the 2304;; scale factor isn't properly applied. 2305(assert-equal ".00000003d+8" (format nil 2306 "~9,4,,-7E" 2307 (bigfloat:bigfloat pi))) 2308(assert-equal ".000003d+6" (format nil 2309 "~9,4,,-5E" 2310 (bigfloat:bigfloat pi))) 2311(assert-equal "3141600.d-6" (format nil 2312 "~5,4,,7E" 2313 (bigfloat:bigfloat pi))) 2314(assert-equal " 314.16d-2" (format nil 2315 "~11,4,,3E" 2316 (bigfloat:bigfloat pi))) 2317(assert-equal " 31416.d-4" (format nil 2318 "~11,4,,5E" 2319 (bigfloat:bigfloat pi))) 2320(assert-equal " 0.3142d+1" (format nil 2321 "~11,4,,0E" 2322 (bigfloat:bigfloat pi))) 2323(assert-equal ".03142d+2" (format nil 2324 "~9,,,-1E" 2325 (bigfloat:bigfloat pi))) 2326(assert-equal "0.003141592653589793d+3" (format nil 2327 "~,,,-2E" 2328 (bigfloat:bigfloat pi))) 2329(assert-equal "31.41592653589793d-1" (format nil 2330 "~,,,2E" 2331 (bigfloat:bigfloat pi))) 2332;; Fails because exponent is printed as "b0" instead of "b+0" 2333(assert-equal "3.141592653589793b+0" (format nil "~E" (bigfloat:bigfloat pi))) 2334 2335 2336;; These fail because too many digits are printed and because the 2337;; scale factor isn't properly applied. 2338(assert-equal ".03142d+2" (format nil "~9,5,,-1E" (bigfloat:bigfloat pi))) 2339(assert-equal " 0.03142d+2" (format nil "~11,5,,-1E" (bigfloat:bigfloat pi))) 2340(assert-equal "| 3141593.d-06|" (format nil "|~13,6,2,7E|" (bigfloat:bigfloat pi))) 2341(assert-equal "0.314d+01" (format nil "~9,3,2,0,'%E" (bigfloat:bigfloat pi))) 2342(assert-equal "+.003d+03" (format nil "~9,3,2,-2,'%@E" (bigfloat:bigfloat pi))) 2343(assert-equal "+0.003d+03" (format nil "~10,3,2,-2,'%@E" (bigfloat:bigfloat pi))) 2344(assert-equal "=====+0.003d+03" (format nil "~15,3,2,-2,'%,'=@E" (bigfloat:bigfloat pi))) 2345(assert-equal "0.003d+03" (format nil "~9,3,2,-2,'%E" (bigfloat:bigfloat pi))) 2346(assert-equal "%%%%%%%%" (format nil "~8,3,2,-2,'%@E" (bigfloat:bigfloat pi))) 2347 2348;; Works 2349(assert-equal "0.0f+0" (format nil "~e" 0)) 2350 2351;; Fails because exponent is printed as "b0" instead of "b+0' 2352(assert-equal "0.0b+0" (format nil "~e" (bigfloat:bigfloat 0d0))) 2353;; Fails because exponent is printed as "b0 " instead of "b+0000" 2354(assert-equal "0.0b+0000" (format nil "~9,,4e" (bigfloat:bigfloat 0d0))) 2355;; Fails because exponent is printed as "b4" isntead of "b+4" 2356(assert-equal "1.2345678901234567b+4" (format nil "~E" 2357 (bigfloat:bigfloat 1.234567890123456789d4))) 2358 2359;; Fails because exponent is printed as "b36" instead of "b+36" 2360(assert-equal "1.32922799578492b+36" (format nil "~20E" 2361 (bigfloat:bigfloat (expt 2d0 120)))) 2362;; Fails because too many digits are printed and the exponent doesn't include "+". 2363(assert-equal " 1.32922800b+36" (format nil "~21,8E" 2364 (bigfloat:bigfloat (expt 2d0 120)))) 2365|# 2366 2367 2368;; Format bigfloats using ~F format. This is suitable as a ~// format. 2369;; 2370;; NOTE: This is a modified version of FORMAT-FIXED from CMUCL to 2371;; support printing of bfloats. 2372 2373(defun format-f (stream number colonp atp 2374 &optional w d k overflowchar padchar) 2375 (typecase number 2376 (bigfloat 2377 (maxima::bfloat-format-f stream (real-value number) colonp atp 2378 w d (or k 0) 2379 overflowchar 2380 (or padchar #\space))) 2381 (complex-bigfloat 2382 ;; FIXME: Do something better than this since this doesn't honor 2383 ;; any of the parameters. 2384 (princ number stream)) 2385 (otherwise 2386 ;; We were given some other kind of object. Just use CL's normal 2387 ;; ~F printer to print it. 2388 (let ((f 2389 (with-output-to-string (s) 2390 ;; Construct a suitable ~F format string from the given 2391 ;; parameters. First handle w,d,k. 2392 (write-string "~V,V,V," s) 2393 (if overflowchar 2394 (format s "'~C," overflowchar) 2395 (write-string "," s)) 2396 (if (char= padchar #\space) 2397 (write-string "," s) 2398 (format s "'~C," padchar)) 2399 (when colonp 2400 (write-char #\: s)) 2401 (when atp 2402 (write-char #\@ s)) 2403 (write-char #\F s)))) 2404 (format stream f w d k number))))) 2405 2406;; Format bigfloats using ~G format. This is suitable as a ~// format. 2407;; 2408;; NOTE: This is a modified version of FORMAT-GENERAL from CMUCL to 2409;; support printing of bfloats. 2410 2411(defun format-g (stream number colonp atp 2412 &optional w d e k overflowchar padchar marker) 2413 (typecase number 2414 (bigfloat 2415 (maxima::bfloat-format-g stream (real-value number) colonp atp 2416 w d e (or k 1) 2417 overflowchar 2418 (or padchar #\space) 2419 (or marker #\b))) 2420 (complex-bigfloat 2421 ;; FIXME: Do something better than this since this doesn't honor 2422 ;; any of the parameters. 2423 (princ number stream)) 2424 (otherwise 2425 ;; We were given some other kind of object. Just use CL's normal 2426 ;; ~G printer to print it. 2427 (let ((f 2428 (with-output-to-string (s) 2429 ;; Construct a suitable ~E format string from the given 2430 ;; parameters. First, handle w,d,e,k. 2431 (write-string "~V,V,V,V," s) 2432 (if overflowchar 2433 (format s "'~C," overflowchar) 2434 (write-string "," s)) 2435 (if padchar 2436 (format s "'~C," padchar) 2437 (write-string "," s)) 2438 (when marker 2439 (format s "'~C" marker)) 2440 (when colonp 2441 (write-char #\: s)) 2442 (when atp 2443 (write-char #\@ s)) 2444 (write-char #\G s)))) 2445 (format stream f w d e k number))))) 2446 2447 2448