;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*- (in-package :mma) ;;; Common Lisp Bigfloat Package, Part I: Basic arithmetic only. ;; copyright (c) 1990 Richard Fateman, UC Berkeley ;;; Bigfloats are stored as structures: precision, fraction and exponent are ;;; integers, with an implied binary point to the left of the fraction. ;;; The fraction is normalized (except for zero). ;;; We use the common lisp data structure type discrimination mechanism ;;; for tagging, printing, etc. bigfloats. ;;; Another file, bfelem provides sin/cos/log/exp of bigfloats. ;;; (provide 'bf) (defvar bigfloat-bin-prec 10 "The number of bits in the fraction of a new bigfloat") (defvar bigfloat-chop nil "If T makes bigfloat arithmetic faster and less accurate by disabling rounding") (defvar bigfloat-print-trunc t "If T means drop trailing 0's in the fraction when printing bigfloats") (defvar bigfloat-printprec 0 "If >0, then restrict printing to that many digits when printing bigfloats") (proclaim '(special bigfloatone bigfloatzero)) ;; declare the structure of bigfloats (defstruct (bigfloat (:constructor bcons2 (fraction exponent &optional (precision bigfloat-bin-prec))) (:print-function bigfloatprintfunction)) (fraction 0 :type integer :read-only t) (exponent 0 :type integer :read-only t) (precision bigfloat-bin-prec :type integer :read-only t) ) ;; test for equality of two bigfloat (defun bigfloat-eql(x y) (and (eql (bigfloat-fraction x)(bigfloat-fraction y)) (eql (bigfloat-exponent x)(bigfloat-exponent y)))) ;; how many bits does it take to represent x decimal digits? (defun bits-in-decimal(x) (ceiling (* x #.(log 10 2)))) ;; initialization routines ;; Any time bigfloat precision is changed, one of these routine should ;; be called: either bigfloat-init-dec or bigfloat-init-bin. ;; bigfloat-init-dec is ;; given: a certain number of decimal digits to be carried (approximately) (defun bigfloat-init-dec (q) (bigfloat-init-bin (bits-in-decimal q))) ;; bigfloat-init-bin is ;; given: a certain number of bits to be carried in the fraction (defun bigfloat-init-bin (q) (setq bigfloat-bin-prec q bigfloatone (intofp 1) ;;initialize one bigfloatzero (bcons2 0 0)) ;; initialize zero q) (defun intlen(x)(integer-length (abs x))) ;; into-bigfloat is called by the reader to make a bigfloat. ;; The two arguments are integers representing frac * 10^exp. ;; The input routines should convert 1.23*10^5 to 123 * 10^(-2), ;; and then call into-bigfloat on 123 and -2. ;; the return value is an accurately converted number. (defun into-bigfloat (frac exp) ;; log[2](10) is 3.321928 ... ;; jack up precision if it is needed for good conversion (let ((bigfloat-bin-prec (max bigfloat-bin-prec (+ (ceiling 3.321193 (intlen exp)) (intlen frac))))) ;; compute frac*10^exp (bigfloat-* (intofp frac) (bigfloat-expt(intofp 10)exp)))) (defun bigfloat-convert(x) ;; takes whatever x is, and returns a bigfloat. (cond ((typep x 'bigfloat) (cond((= bigfloat-bin-prec (bigfloat-precision x)) x) (t (bfnormal (bigfloat-fraction x) (bigfloat-exponent x))))) ((typep x 'integer)(intofp x)) ;bignum, fixnum, ;; this is the common-lisp rational number type ((typep x 'ratio) (bigfloat-/ ;; can't lose on under/overflow (intofp (numerator x)) (intofp (denominator x)))) ;; for a flooating point number, convert to rational, first. ((typep x 'float) ;single, double or other float (bigfloat-convert (rational x))) (t x))) ;; if it can't be converted, return x? ;; integer to bigfloat conversion (defun intofp(l)(if (equal 0 l) bigfloatzero (bfnormal l (intlen l)))) ;; bfnormal "normalizes" a floating point fraction. ;; the proposed number has fraction frac (with binary bit at the left) ;; and exponent exp. ;; But frac may have the wrong number of bits. To be ;; normalized, (intlen frac) must be = to bigfloat-bin-prec ;; unless frac =0. (defun bfnormal(frac exp &optional (il (intlen frac))) (cond ((= il bigfloat-bin-prec) (bcons2 frac exp)) (bigfloat-chop (bfnormal-chop frac exp il)) (t(bfnormal-round frac exp il)))) (defun bfnormal-chop (frac exp il); (if (= frac 0) bigfloatzero (bcons2 (ash frac (- bigfloat-bin-prec il)) exp))) (defun bfnormal-round (frac exp il); (if (= frac 0) bigfloatzero (let ((s (- bigfloat-bin-prec il))) (cond ((> s 0) ;; not enough bits, so shift left what you have (bcons2 (ash frac s) exp)) ;; s= 0 taken care of in bfnormal. (t ;; that is, (< s 0) ;; too many bits. Must round. (let ((r (round frac (ash 1 (- s))))) (if (= (intlen r) bigfloat-bin-prec) (bcons2 r exp) ;; usual case ;; Take care of rare case when rounding up bumps the length. ;; Analogy in decimal: round 999 x 10^0 to 2 digits. ;; You get 100 x 10^3. (bcons2 (ash r -1) (1+ exp))))))))) ;; add two bigfloats of any precision; converts each to global precision. ;; c = a+b. (defun bigfloat-+ (a b) (setq a (bigfloat-convert a)) (setq b (bigfloat-convert b)) (let* ((fa (bigfloat-fraction a)) (fb (bigfloat-fraction b)) (ea (bigfloat-exponent a)) (eb (bigfloat-exponent b)) fc (dif (- ea eb))) (cond ((> dif bigfloat-bin-prec) a) ;b is insignificant ((> (- dif) bigfloat-bin-prec) b) ;a is insignificant ((>= dif 0) (setq fc (+ (ash fa dif) fb)) (bfnormal fc (+ eb (- (intlen fc) bigfloat-bin-prec)))) (t (setq fc (+ (ash fb (- dif)) fa)) (bfnormal fc (+ ea (- (intlen fc) bigfloat-bin-prec))))))) ;; bigfloat-* multiplies two bigfloats of arbitrary (perhaps different) ;; precisions, and return a bigfloat of global (bigfloat-bin-prec) ;; precision. (defun bigfloat-* (a b) (let* ((f (* (bigfloat-fraction a)(bigfloat-fraction b))) (il (intlen f)) (exp (+ (bigfloat-exponent a) (bigfloat-exponent b) (- (bigfloat-precision a)) (- (bigfloat-precision b)) il))) (bfnormal f exp il))) ;; bigfloat-/ divides two bigfloats of arbitrary (perhaps different) ;; precisions, and return a bigfloat of global (bigfloat-bin-prec) ;; precision. (defun bigfloat-/ (a b) ;; division by zero will be noticed in the "round" below (let* ((f (round (ash (bigfloat-fraction a) (+ bigfloat-bin-prec (bigfloat-precision b))) (bigfloat-fraction b))) (il (intlen f)) (exp (+ (bigfloat-exponent a) (- (bigfloat-exponent b)) (- (bigfloat-precision a)) (- bigfloat-bin-prec) il))) (bfnormal f exp il))) ;; compute 1/x for bigfloat x (defun bigfloat-inv(x) (bigfloat-/ (bfone) x)) ;; coerce-bigfloat is like the CL function coerce, but allows ;; the first argument to be of type bigfloat. It converts the ;; first argument to ratio, double-float, single-float (= float), ;; integer (= bignum or fixnum) (defun coerce-bigfloat(x typ) (cond ;; convert bigfloat to bigfloat of global precision ((eq typ 'bigfloat) (bigfloat-convert x)) ((eq typ 'ratio) ;; convert a bigfloat to a rational fraction. This is exact in ;; the sense that for any represented bigfloat x in bigfloat-bin-prec, ;; (bigfloat-convert( bigfloattorat x)) = x (* (bigfloat-fraction x) (expt 2 (- (bigfloat-exponent x) (bigfloat-precision x))))) ((eq typ 'double-float) ;; convert a bigfloat into a floating point double when it is ;; simple to do. Won't work for very high precision, since ;; the fraction part will overflow. (* (coerce (bigfloat-fraction x) 'double-float) (expt 2.0d0 (- (bigfloat-exponent x)(bigfloat-precision x))))) ((member typ '(float single-float) :test #'eq) (coerce (coerce-bigfloat x 'double-float) 'single-float)) ((member typ '(integer bignum fixnum)) (round (coerce-bigfloat x 'ratio))) (t (error "Can't coerce ~s to type ~s~%" x typ)))) ;; how many decimal digits in an integer x? ;; slow, but accurate. Used in i/o conversions. (defun decimalsin(x) (length (format nil "~s" x))) ;; various useful predicates (defun bigfloat-> (a b) (bigfloat-posp (bigfloat-diff a b))) (defun bigfloat-< (a b) (bigfloat-posp (bigfloat-diff b a))) (defun bigfloat-posp (x)(> (bigfloat-fraction x) 0)) (defun bigfloat-zerop(x)(= (bigfloat-fraction x) 0)) (defun bigfloat-min(&rest args) (let (min) (cond ((null args) (error "bigfloat-min with no arguments"))) (setq min (car args)) (do ((a (cdr args) (cdr a))) ((null a) min) (cond ((bigfloat-< (car a) min) (setq min (car a))))))) (defun bigfloat-max(&rest args) (let (max) (cond ((null args) (error "bigfloat-max with no arguments"))) (setq max (car args)) (do ((a (cdr args) (cdr a))) ((null a) max) (cond ((bigfloat-> (car a) max) (setq max (car a))))))) ;; return an existing 1 or make a new one. (defun bfone nil (cond ((= bigfloat-bin-prec (bigfloat-precision bigfloatone)) bigfloatone) (t (setq bigfloatone (intofp 1))))) ;; compute x-y or -x (if y is missing) (defun bigfloat-- (x &optional (y nil)) (cond (y (bigfloat-+ x (bigfloat-- y)));; compute x-y ((bigfloat-zerop x) x) (t (bcons2 (- (bigfloat-fraction x)) (bigfloat-exponent x))))) ;; compute p^n for bigfloat p. nn is positive or negative integer. ;; Note that we do NOT allow bigfloat exponent here. (use log/exp for that) (defun bigfloat-expt (p nn) (cond ((eql nn 0) (bfone)) ((eql nn 1) p) ((< nn 0) (bigfloat-inv (bigfloat-expt p (- nn)))) (t (do ((n (floor nn 2)(floor n 2)) (s (cond ((oddp nn) p) (t (bfone))))) ((zerop n) s) (setq p (bigfloat-* p p)) (and (oddp n) (setq s (bigfloat-* s p))) ) ))) ;; compute square-root of bigfloat x (defun bigfloat-sqrt(x)(bigfloat-root x 2)) ;; compute a^(1/n) see Fitch, SIGSAM Bull Nov 74 (defun bigfloat-root (a n) (let* ((ofprec bigfloat-bin-prec) (bigfloat-bin-prec (+ bigfloat-bin-prec 2)) ;assumes a>0 n>=2 (bk (bigfloat-expt (intofp 2) (1+ (truncate (bigfloat-exponent (setq a (bigfloat-convert a))) n))))) (do ((x bk (bigfloat-diff x (setq bk (bigfloat-/ (bigfloat-diff x (bigfloat-/ a (bigfloat-expt x n1))) n)))) (n1 (1- n)) (n (intofp n))) ((or (bigfloat-zerop bk) (> (- (bigfloat-exponent x) (bigfloat-exponent bk)) ofprec)) (setq a x)))) (bigfloat-convert a)) ;return to previous precision (defun bigfloat-abs (x) (cond ((>= (bigfloat-fraction x) 0) x) (t (bcons2 (- (bigfloat-fraction x)) (bigfloat-exponent x))))) ;; integer part of a bigfloat (defun bigfloat-intpart (f) (let ((m (- (bigfloat-precision f) (bigfloat-exponent f)))) (if (> m 0) (truncate (bigfloat-fraction f) (ash 1 m)) (* (bigfloat-fraction f) (ash 1 (- m)))))) ;;set up user's precision ;;(bigfloat-init-dec 3);; 3 decimal digits for debugging ;; some tests, some utilities (defun bc (x &optional (bigfloat-bin-prec bigfloat-bin-prec )) (bigfloat-convert x)) ;; this definition helps for tracing functions so you can read everything.. ;; at least while you are using bignums that are exactly representable (defvar bigfloat-print-trace nil "t if you want to see internals of bigfloats") (defun bigfloatprintfunction (x s pl) ;pl, print-level, is not used. (cond (bigfloat-print-trace (format s "[~s*2^~s]=~a" (bigfloat-fraction x) (bigfloat-exponent x) ;; (coerce-bigfloat x 'double-float) ; test (bigfloat-print x nil) )) (t (bigfloat-print x s)))) ;;(bigfloat-print x) normally prints the bigfloat to the standard output. ;; it looks like d.dd*10^ee where the number of digits, d, is specified ;; by the width w, or if w is not provided, by the number of decimal ;; digits assumed plausible by the binary precision. ;; ee, the exponent, is computed as appropriate. ;;(bigfloat-print x stream) prints to a stream, and if stream=nil, ;; it returns a string. ;;(bigfloat-print x stream w) prints to a stream (if stream=t, standard ;; output; and also uses a width for the fraction part of w digits. ;; The total field width will be at most w + 4 + width of exponent. ;; optionally there will be a negative sign taking one space, also. ;; this version of bigfloat-print prints d.ddd..d *10^n ;; except if n=0, when it prints d.ddd..d (no *10^n). Also ;; trailing zeros in the fraction are dropped off if bigfloat-print-trunc ;; is t, and (defun bigfloat-print (x &optional (stream t) (w bigfloat-printprec)&aux sign) (flet ;; utility functions ((t0(s); change "12300" to "123" (if bigfloat-print-trunc (string-right-trim "0" s) s)) (putindot(s) ; change "123" to "1.23" (concatenate 'string (subseq s 0 1) "." (subseq s 1)))) (cond((bigfloat-zerop x) (format stream "0.0")) (t (multiple-value-bind (frac exp) (bf-format x w) (setq sign (< frac 0)) (setq frac (format nil "~s" (abs frac))) (format stream "~:[~;-~]~a~:[*10^~s~;~]" sign (t0(putindot frac)) (zerop (setq exp(+ (length frac) -1 exp))) exp)))))) ;; bf-format returns a possible pair of decimal fraction and exponent ;; values for printing x. ;; In particular, two integer values f and e are returned, and the ;; number x is closely approximated by f x 10^e. (defun bf-format( x w) (cond ((bigfloat-zerop x) (values 0 0)) (t (let*((f (abs(bigfloat-fraction x))) (e (bigfloat-exponent x)) (p (bigfloat-precision x)) (bigfloat-bin-prec p) ;use the precision of the number (width (if (> w 0) w (truncate (* bigfloat-bin-prec 0.30103)))) ;; to convert .f x 2^e to decimal ;; compute the log[10] of the number, approximately (d (floor (+ ;; log of the exponent (* (+ e (- p) (intlen f)) #.(/ 1.0d0(log 10.0d0 2))) ;; log of the fraction ;; we can do this in double-float since it is ;; unreasonable to have an exponent exceeding +/- 10^307. ;; Numbers of the size (2 ^(10^307)) won't ;; print. (log (coerce (/ f (expt 2 (intlen f))) 'double-float) 10)))) (s (if (>= d 0) (bigfloat-/ x ;;note: x still has its sign ;;this is cheap ;;(bigfloat-expt (intofp 10) d) (intofp (expt 10 d));; this is accurate ) (bigfloat-* x (intofp (expt 10 (- d)))))) ;; now 1 < |s| <= 10 ? , x = s * 10^d (left(round (/(* (expt 10 (1- width))(bigfloat-fraction s)) (expt 2 (- (bigfloat-precision s) (bigfloat-exponent s)))))) ) (cond ((> (decimalsin left) width) (setq d (1+ d) left (round left 10)))) ;; this could be set up as a string to return.. ;; The two items of interest that could be returned are ;; the fraction and exponent (values left (- d (1- width)))))))