;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*- ;; (c) 1986, 1991 Richard Fateman ;; Common Lisp Bigfloats, part II: ;; some elementary functions: sin, cos, exp,log,atan ;(proclaim '(optimize (speed 3))) ;(proclaim '(optimize (safety 0))) (eval-when (compile) (load "bf")) (proclaim '(special bigfloat-bin-prec bigfloatone bigfloatzero bigfloat%e bigfloat%pi)) (defmacro fpsin1 (x) `(fpsincos1 ,x t)) (defmacro fpcos1 (x) `(fpsincos1 ,x nil)) ;; initialization (defun bigfloat-init (q) ;; q is an integer: how many decimal digits needed. (setq bigfloat-bin-prec (+ 2 (integer-length (expt 10 q))) bigfloatone (intofp 1) bigfloatzero (bcons2 0 0) bigfloat%e (fpexp1 (intofp 1)) bigfloat%pi (fppi1 bigfloat-bin-prec)) q) (defun bigfloat-atan (x) (let*(term x2 tmp (one (intofp 1)) (two (intofp 2)) (ans one) (oans two)) (cond ((bigfloat-> (bigfloat-abs x) one) (setq tmp (bigfloat-/ (fppi) two)) (setq ans (bigfloat-- tmp (bigfloat-atan (bigfloat-/ one x)))) (cond ((bigfloat-> ans tmp) (setq ans (bigfloat-- ans (fppi)))))) ((bigfloat-> (bigfloat-abs x) (bigfloat-/ one two)) (setq tmp (bigfloat-/ x (bigfloat-+ (bigfloat-* x x) one))) (setq x2 (bigfloat-* x tmp) term (setq ans one)) (do ((n 0 (1+ n))) ((bigfloat-eql ans oans) ans) (setq term (bigfloat-* term (bigfloat-* x2 (bigfloat-/ (intofp (+ 2 (* 2 n))) (intofp (+ (* 2 n) 3)))))) (setq oans ans ans (bigfloat-+ term ans))) (setq ans (bigfloat-* tmp ans))) (t (setq ans x x2 (bigfloat-- (bigfloat-* x x)) term x) (do ((n 3 (+ n 2))) ((bigfloat-eql ans oans) ans) (setq term (bigfloat-* term x2)) (setq oans ans ans (bigfloat-+ ans (bigfloat-/ term (intofp n))))))) ans)) (defun bigfloat-atan2 (y x) ;atan(y/x) from -pi to pi (cond ((bigfloat-zerop x) ; atan(inf), but what sign? (cond ((bigfloat-zerop y) (error "atan(0/0)")) ((bigfloat-posp y) (bigfloat-/ (fppi) (intofp 2))) (t(bigfloat-/ (fppi) (intofp -2))))) ((bigfloat-posp x) (cond ((bigfloat-posp y) (bigfloat-atan (bigfloat-/ y x))) (t (bigfloat-- (bigfloat-atan (bigfloat-/ y x)))))) ((bigfloat-posp y) (bigfloat-+ (fppi) (bigfloat-atan (bigfloat-/ y x)))) (t (bigfloat-- (bigfloat-atan (bigfloat-/ y x)) (fppi))))) (defun bigfloat-tan (a) (bigfloat-convert (let* ((bigfloat-bin-prec (+ 8. bigfloat-bin-prec)) (y (bigfloat-convert a ))) (bigfloat-/ (bigfloat-sin1 y t) (bigfloat-sin1 y nil))))) (defun bigfloat-exp (x) ;; exponential for arbitrary arg (prog (r s) (cond ((bigfloat-posp x)) ;; for negative x, compute 1/e^(-x) (t (return (bigfloat-inv (bigfloat-exp (bigfloat-- x)))))) (setq r (bigfloat-intpart x)) ;; for large x, compute e^r * exp(x-r) (return (cond ((< r 2) (fpexp1 x)) (t (setq s (fpexp1 (bigfloat-- x (intofp r)))) (bigfloat-* s (bigfloat-convert (let ((bigfloat-bin-prec (+ bigfloat-bin-prec (integer-length r) -1))) (bigfloat-expt (fpe) r))))))))) ;; exponential for small arg (defun fpexp1 (x) (do (( n 1 (1+ n)) (ans (intofp 1)) (term (intofp 1)) (oans (intofp 0))) ((bigfloat-eql ans oans) ans) ;quit when two iterations same (setq term(bigfloat-/ (bigfloat-* x term) (intofp n)) oans ans ans(bigfloat-+ ans term)))) ;; (fpe) return big floating point %e. ;; it returns truncated bigfloat%e if possible, else recomputes. ;; in any case, bigfloat%e is set to last used value. (defun fpe nil (cond ((= bigfloat-bin-prec (bigfloat-precision bigfloat%e)) bigfloat%e) ((< bigfloat-bin-prec (bigfloat-precision bigfloat%e)) (setq bigfloat%e (bigfloat-convert bigfloat%e))) (t (setq bigfloat%e (fpexp1 (intofp 1)))))) (defun fppi nil (cond ((= bigfloat-bin-prec (bigfloat-precision bigfloat%pi)) bigfloat%pi) ((< bigfloat-bin-prec (bigfloat-precision bigfloat%pi)) (setq bigfloat%pi (bigfloat-convert bigfloat%pi))) (t (setq bigfloat%pi (fppi1 bigfloat-bin-prec))))) (defun fpone nil (cond ((= bigfloat-bin-prec (bigfloat-precision bigfloatone)) bigfloatone) (t (setq bigfloatone(intofp 1))))) ;; comppi computes pi to n bits. ;; that is, (comppi n)/(2^n) is an approximation to pi. (defun comppi (n) (let* ((a (expt 2 n)) (b (truncate a 8)) (c (+ (* 3 a) b))) (do ((i 4 (+ i 2))) ((zerop b) c) (setq b (truncate (* b (1- i) (1- i)) (* 4 i (1+ i)))) (setq c (+ c b))) )) (defun fppi1 (n) ;; compute pi to that many bits plus a few (bigfloat-/ (intofp (comppi (+ n 3))) (intofp (expt 2 (+ 3 n))))) ;; fl= t for sin, nil for cos (defun bigfloat-sin (x) (bigfloat-sin1 x t)) (defun bigfloat-cos (x) (bigfloat-sin1 x nil)) ;; Bigfloat-sin1 computes sin or cos to precision bigfloat-bin-prec, ;; It uses extra precision for the value of pi for large argument ;; values. It is still possible to get fooled somewhat by an argument ;; value that is surprisingly accurate in its approximation to an ;; integer multiple of pi/2. (More accurate than precision+exponent+4) ;; in which case the result will not be as accurate as the ;; precision would permit. (defun bigfloat-sin1 (x fl) (flet((insine (bigfloat-bin-prec xt) (let* ((x (bigfloat-convert xt)) (piby2 (bigfloat-/ (fppi) (intofp 2))) (r (bigfloat-intpart (bigfloat-/ x piby2)))) ; (format t "~a, ~a, ~a~%" x piby2 r) ;debug info (setq x (bigfloat-+ x (bigfloat-* (intofp (- r)) piby2))) (setq r (mod r 4)) (cond (fl (cond ((= r 0) (fpsin1 x)) ((= r 1) (fpcos1 x)) ((= r 2) (bigfloat-- (fpsin1 x))) ((= r 3) (bigfloat-- (fpcos1 x))))) (t (cond ((= r 0) (fpcos1 x)) ((= r 1) (bigfloat-- (fpsin1 x))) ((= r 2) (bigfloat-- (fpcos1 x))) ((= r 3) (fpsin1 x)))))))) (let ((sign (if fl (bigfloat-posp x) nil)) res) (setq x (bigfloat-abs x)) (setq res (cond ((bigfloat-zerop x) (if fl (intofp 0) ;sin(0) (intofp 1) ;cos(0) )) (t (bigfloat-convert (insine (max bigfloat-bin-prec (+ 4 bigfloat-bin-prec (bigfloat-exponent x))) x))))) (cond (sign res) (t (bigfloat-- res)))))) ;; compute sin or cos in (0,pi/2). fl is t for sin, nil for cos. (defun fpsincos1 (x fl) (prog (ans term oans x2) (setq ans (cond (fl x)(t(intofp 1))) x2 (bigfloat-- (bigfloat-* x x))) (setq term ans) (do ((n (cond (fl 3) (t 2)) (+ n 2))) ((and oans(bigfloat-eql ans oans)) ans) ;;oans is initial nil.. bug (setq term (bigfloat-* term (bigfloat-/ x2 (intofp (* n (1- n)))))) (setq oans ans ans (bigfloat-+ ans term))) (return ans))) (defun bigfloat-log (x) (prog (over two ans term e) (cond ((not (bigfloat-posp x)) (error "non-pos arg to bigfloat-log"))) (setq e (fpe) over (bigfloat-inv e) ans 0) (do () (nil) (cond ((bigfloat-eql x e) (setq x nil) (return nil)) ((and (bigfloat-< x e) (bigfloat-< over x)) (return nil)) ((bigfloat-< x over) (setq x (bigfloat-* x e)) (setq ans (1- ans))) (t (setq ans (1+ ans)) (setq x (bigfloat-/ x e))))) (cond ((null x) (return (intofp (1+ ans))))) (setq x (bigfloat-- x (fpone)) ans (intofp ans)) (setq x (bigfloat-expt (setq term (bigfloat-/ x (bigfloat-+ x (setq two (intofp 2))))) 2)) (do ((n 0 (1+ n)) (oldans ans) (ans ans (bigfloat-+ ans (bigfloat-* two (bigfloat-/ term (intofp (1+ (* 2 n))))))) (term (bigfloat-* term x))) ((bigfloat-eql ans oldans) ans))))