1;; Copyright 2009 by Barton Willis 2 3;; This is free software; you can redistribute it and/or 4;; modify it under the terms of the GNU General Public License, 5;; http://www.gnu.org/copyleft/gpl.html. 6 7;; This software has NO WARRANTY, not even the implied warranty of 8;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 9 10;; The last time I tried, the file "hypergeometric.lisp" must be loaded before compiling nfloat; so 11 12;(eval-when (compile) 13; ($load "hypergeometric.lisp")) 14 15(in-package :maxima) 16 17 18(in-package #-gcl #:bigfloat #+gcl "BIGFLOAT") 19 20;; Each member of the CL list l is a two member list (running error form). 21(defun running-error-plus (l) 22 (let ((acc 0) (err 0)) 23 (dolist (lk l) 24 (setq acc (+ acc (first lk))) 25 (setq err (+ err (second lk) (abs acc)))) 26 (list acc err))) 27 28;;(%i20) ah * lh * (1 + eps); 29;;(%o20) ah*(eps+1)*lh 30;;(%i21) %-a*l; 31;;(%o21) ah*(eps+1)*lh-a*l 32;;(%i22) subst([a = ah-e,l = lh-w],%); 33;;(%o22) ah*(eps+1)*lh-(ah-e)*(lh-w) 34;;(%i23) expand(%); 35;;(%o23) -e*w+ah*w+ah*eps*lh+e*lh 36 37(defun running-error-prod (l) 38 (let ((acc 1) (err 0) (z)) 39 (dolist (lk l) 40 (setq z acc) 41 (setq acc (* acc (first lk))) 42 (setq err (+ (* err (abs (first lk))) (* (abs z) (abs (second lk)))))) 43 (list acc err))) 44 45;; (%i1) (a+ae*eps)*(1+eps)/(b+be*eps)-a/b$ 46;; (%i2) expand(limit(%/eps,eps,0))$ 47;; (%i3) expand(ratsubst(Q,a/b,%))$ 48;; (%i4) map('abs,%)$ 49;; (%i5) facsum(%,abs(Q)); 50;; (%o5) ((abs(be)+abs(b))*abs(Q)+abs(ae))/abs(b) 51 52(defun running-error-quotient (l) 53 (let* ((a (first l)) (b (second l)) (s)) 54 (setq s (/ (first a) (first b))) 55 (list s (+ (* (abs s) (+ 1 (abs (/ (second b) (first b))))) (abs (/ (second a) (first b))))))) 56 57;; unary negation. 58(defun running-error-negate (x) 59 (setq x (first x)) 60 (list (- (first x)) (second x))) 61 62;;(%i46) (x*(1+ex))^(n *(1+en)); 63;;(%o46) ((ex+1)*x)^((en+1)*n) 64;;(%i47) taylor(%,[ex,en],0,1); 65;;(%o47) x^n+(x^n*n*ex+x^n*n*log(x)*en)+... 66;;(%i48) factor(%); 67;;(%o48) x^n*(en*n*log(x)+ex*n+1) 68 69(defun running-error-expt (l) 70 (let* ((s) (x (first l)) (n (second l)) (ex) (en)) 71 (setq ex (second x)) 72 (setq en (second n)) 73 (setq x (first x)) 74 (setq n (first n)) 75 (setq s (bigfloat::expt x n)) 76 (list s (+ (abs (* s en n (log x))) (abs (* s ex n)))))) 77 78;; sqrt(x + ex) = sqrt(x)+(sqrt(x)*ex)/(2*x)+... = sqrt(x) + ex / (2 * sqrt(x)) + ... 79(defun running-error-sqrt (x) 80 (setq x (first x)) 81 (let ((s (sqrt (first x)))) 82 (list s (abs (/ (second x) (* 2 s)))))) 83 84;; log(x + ex) = log(x) + ex / x + ... 85(defun running-error-log (x) 86 (setq x (first x)) 87 (list (log (first x)) (abs (/ (second x) (first x))))) 88 89(defun psi0 (x) 90 (bigfloat (maxima::$rectform (maxima::mfuncall 'maxima::$bfpsi0 (maxima::$bfloat (maxima::to x)) 91 maxima::$fpprec)))) 92 93(defun gamma (x) 94 (bigfloat (maxima::$bfloat (maxima::take '(maxima::%gamma) (maxima::to x))))) 95 96;; gamma(x + ex) = gamma(x) + ex * gamma(x) * psi[0](x) + .. 97(defun running-error-gamma (x) 98 (setq x (first x)) 99 (let ((s (gamma (first x)))) 100 (list s (abs (* s (second x) (psi0 (first x))))))) 101 102(defun running-error-hypergeometric (a b x subs bits) 103 (let ((dig) (d) (f)) 104 105 ;; To do this correctly, we'd need the partial derivatives of the hypergeometric functions 106 ;; with respect the the parameters. Ouch! 107 108 (setq a (mapcar #'(lambda (s) (car (running-error-eval s subs bits))) (maxima::margs a))) 109 (setq b (mapcar #'(lambda (s) (car (running-error-eval s subs bits))) (maxima::margs b))) 110 (setq x (car (running-error-eval x subs bits))) 111 (cond ((< (abs x) 0.99) 112 (multiple-value-setq (f d) (hypergeometric-by-series a b x)) 113 (list f (* (abs f) (expt 10 (- d))))) 114 (t 115 (setq dig (ceiling (* bits #.(/ (log 2.0) (log 10.0))))) 116 (setq a (mapcar 'maxima::to a)) 117 (setq b (mapcar 'maxima::to b)) 118 (setq x (maxima::to x)) 119 (multiple-value-setq (f d) (hypergeometric-float-eval a b x dig t)) 120 (list f (* (abs f) (expt 10 (- d)))))))) 121 122(defun running-error-sum (l subs bits) 123 (let ((sumand (first l)) 124 (v (second l)) 125 (lo (third l)) 126 (hi (fourth l)) 127 (acc 0) (err 0) (x) (q)) 128 129 (cond ((and (integerp lo) (integerp hi)) 130 (maxima::while (<= lo hi) 131 (setq q (maxima::$sublis `((maxima::mlist) ((maxima::mequal) ,v ,lo)) sumand)) 132 (setq q (maxima::simplify q)) 133 (setq x (running-error-eval q subs bits)) 134 (incf lo) 135 (setq acc (+ acc (first x))) 136 (setq err (+ err (second x) (abs acc)))) 137 (list acc err)) 138 (t (throw 'maxima::nfloat-nounform-return 'return-nounform))))) 139 140(defun running-error-product (l subs bits) 141 (let ((prodand (first l)) ;; a sum has a summand, so a product has a ... 142 (v (second l)) 143 (lo (third l)) 144 (hi (fourth l)) 145 (acc 1) (err 0) (x)) 146 147 (cond ((and (integerp lo) (integerp hi)) 148 (maxima::while (<= lo hi) 149 (setq x (maxima::$sublis `((maxima::mlist) ((maxima::mequal) ,v ,lo)) prodand)) 150 (setq x (maxima::simplify x)) 151 (setq x (running-error-eval x subs bits)) 152 (incf lo) 153 (setq acc (* acc (first x))) 154 (setq err (+ err (second x) (abs acc)))) 155 (list acc err)) 156 (t (throw 'maxima::nfloat-nounform-return 'return-nounform))))) 157 158(defun running-error-abs (l) 159 (setq l (first l)) 160 (list (abs (first l)) (second l))) 161 162(defun running-error-conjugate (l) 163 (setq l (first l)) 164 (list (conjugate (first l)) (second l))) 165 166;; untested!!!!! 167(defun running-error-factorial (l) 168 (setq l (first l)) 169 (if (integerp l) 170 (list (maxima::take (list 'maxima::mfactorial) l) 0) 171 (running-error-gamma (list (list (+ 1 (first l)) (second l)))))) 172 173(defun running-error-atan2 (l) 174 (let* ((y (first l)) 175 (x (second l)) 176 (d (/ 1 (+ (* (first x) (first x)) (* (first y) (first y)))))) 177 (list (atan (first y) (first x)) 178 (* (+ (* (abs (second y)) (first x)) (* (abs (second x)) (first y))) d)))) 179 180(defun running-error-realpart (l) 181 (setq l (first l)) 182 (list (realpart (first l)) (second l))) 183 184(defun running-error-imagpart (l) 185 (setq l (first l)) 186 (list (imagpart (first l)) (second l))) 187 188 189;; For a similar hashtable mechanism, see trig.lisp. 190(defvar *running-error-op* (make-hash-table :size 16) 191 "Hash table mapping a maxima function to a corresponding Lisp 192 function to evaluate the maxima function numerically using a running error.") 193 194(setf (gethash 'maxima::mplus *running-error-op*) #'running-error-plus) 195(setf (gethash 'maxima::mtimes *running-error-op*) #'running-error-prod) 196(setf (gethash 'maxima::mquotient *running-error-op*) #'running-error-quotient) 197(setf (gethash 'maxima::mminus *running-error-op*) #'running-error-negate) 198(setf (gethash 'maxima::mexpt *running-error-op*) #'running-error-expt) 199(setf (gethash 'maxima::%sqrt *running-error-op*) #'running-error-sqrt) 200(setf (gethash 'maxima::%log *running-error-op*) #'running-error-log) 201(setf (gethash 'maxima::%gamma *running-error-op*) #'running-error-gamma) 202(setf (gethash 'maxima::mabs *running-error-op*) #'running-error-abs) 203(setf (gethash 'maxima::$cabs *running-error-op*) #'running-error-abs) 204(setf (gethash 'maxima::$conjugate *running-error-op*) #'running-error-conjugate) 205(setf (gethash 'maxima::mfactorial *running-error-op*) #'running-error-factorial) 206(setf (gethash 'maxima::$atan2 *running-error-op*) #'running-error-atan2) 207(setf (gethash 'maxima::$realpart *running-error-op*) #'running-error-realpart) 208(setf (gethash 'maxima::$imagpart *running-error-op*) #'running-error-imagpart) 209 210(defun running-error-eval (e subs bits) 211 (let ((f)) 212 213 (cond ((eq e 'maxima::$%i) 214 (setq e (bigfloat::to (if (> bits #.(float-digits 1.0e0)) (maxima::$bfloat 1) (maxima::$float 1)))) 215 (list (bigfloat::to 0 e) (abs e))) 216 217 ((maxima::complex-number-p e #'(lambda (s) (or (maxima::$ratnump s) (maxima::$numberp s)))) 218 (setq e (bigfloat::to (if (> bits #.(float-digits 1.0e0)) (maxima::$bfloat e) (maxima::$float e)))) 219 (list e (abs e))) 220 221 ((and (atom e) (maxima::mget e '$numer)) 222 (running-error-eval (maxima::mget e 'maxima::$numer) '((mlist)) bits)) 223 224 ((and (atom e) (get e 'maxima::sysconst)) 225 (running-error-eval (maxima::$bfloat e) '((mlist)) bits)) 226 227 ((atom e) 228 (setq e (maxima::$sublis subs e)) 229 (if (maxima::complex-number-p e 'maxima::bigfloat-or-number-p) 230 (running-error-eval e nil bits) 231 (throw 'maxima::nfloat-nounform-return 'return-nounform))) 232 233 ;; Return a nounform for expressions (arrays, CRE expressions) that don't 234 ;; appear to be Maxima expressions of the form ((op) a1 a2 ...). 235 ((not (and (consp e) (consp (car e)))) 236 (throw 'maxima::nfloat-nounform-return 'return-nounform)) 237 238 ;; Special case exp(x) (more efficient & accurate than sending this through mexpt). 239 ((and (eq 'maxima::mexpt (caar e)) (eq (second e) 'maxima::$%e)) 240 (setq e (running-error-eval (third e) subs bits)) 241 (let ((z (exp (first e)))) 242 (list z (abs (* (second e) z))))) 243 244 ;; Special case x^n, where n is an integer. For this case, we do not want to 245 ;; convert the integer to a float. This prevents some, but not all, semi-spurious 246 ;; nonzero imaginary parts for (negative real)^integer. 247 ((and (eq 'maxima::mexpt (caar e)) (integerp (third e))) 248 (running-error-expt (list (running-error-eval (second e) subs bits) (list (third e) 0)))) 249 250 ;; main function dispatch. 251 ((setq f (gethash (maxima::mop e) *running-error-op*)) 252 ;(print `(e = ,e mop = ,(maxima::mop e))) 253 (setq e (mapcar #'(lambda (s) (running-error-eval s subs bits)) (maxima::margs e))) 254 (funcall f e)) 255 256 ;; f(x + ex) = f(x) + ex * f'(x) + ... Functions without bigfloat 257 ;; evaluation, for example the Bessel functions, need to be excluded. 258 ;; For now, this code rejects functions of two or more variables. 259 ((and (get (caar e) 'maxima::grad) (null (cdr (maxima::margs e))) 260 (or (gethash (caar e) maxima::*big-float-op*) (maxima::trigp (caar e)) 261 (maxima::arcp (caar e)))) 262 263 (let ((x (running-error-eval (cadr e) subs bits)) (f) (df)) 264 (setq f (maxima::take (list (caar e)) (maxima::to (first x)))) 265 (setq df (get (caar e) 'maxima::grad)) 266 (setq df (maxima::$rectform (maxima::$substitute f (caar df) (cadr df)))) 267 (setq df (bigfloat::to df)) 268 (list (bigfloat::to f) (* (second x) (abs df))))) 269 270 ;; special case hypergeometric 271 ((eq (caar e) 'maxima::$hypergeometric) 272 (running-error-hypergeometric (second e) (third e) (fourth e) subs bits)) 273 274 ;; special case sum. 275 ((or (eq (caar e) 'maxima::%sum) (eq (caar e) 'maxima::$sum)) 276 (running-error-sum (cdr e) subs bits)) 277 278 ;; special case product 279 ((or (eq (caar e) 'maxima::$product) (eq (caar e) 'maxima::%product)) 280 (running-error-product (cdr e) subs bits)) 281 282 ;; special case assignment 283 ((eq (caar e) 'maxima::msetq) 284 (maxima::mset (car e) (car (running-error-eval (cadr e) subs bits)))) 285 286 ;; Yes, via nformat, this can happen. Try, for example, nfloat('(a,b),[a=3,b=7]). 287 ((eq (caar e) 'maxima::mprogn) 288 (let ((q)) 289 (setq e (cdr e)) 290 (dolist (ek e q) 291 (setq q (running-error-eval ek subs bits))))) 292 293 (t (throw 'maxima::nfloat-nounform-return 'return-nounform))))) 294 295;; d * eps is a upper bound for how much e differs from its true value, where eps is 296;; the machine epsilon. 297 298;; First (log = natural log) 299;; 300;; log10(x) = log(x) / log(10) and log2(x) = log(x) / log(2). 301;; So 302;; log10(x) = log2(x) * (log(2) / log(10)). 303;; 304;; Second 305;; -log10(abs(d * eps / e)) = log10(abs(e)) - log10(abs(d)) - log10(eps), 306;; = (log2(abs(e)) - log2(abs(d)) - log2(eps)) * (log(2) / log(10). 307 308;; For log2 we use the binary exponent of the number. Common Lisp gives 309;; (decode-float 0.0) --> 0.0 0 1.0, by the way. 310 311(defun log10-relative-error (d e) 312 313 (if (rationalp d) (setq d (bigfloat (maxima::$bfloat (maxima::to d))))) 314 (if (rationalp e) (setq e (bigfloat (maxima::$bfloat (maxima::to e))))) 315 316 (floor (* 317 (- 318 (second (multiple-value-list (decode-float (abs e)))) 319 (+ 320 (second (multiple-value-list (decode-float (abs d)))) 321 (second (multiple-value-list (decode-float (epsilon (abs d))))))) 322 #.(/ (log 2.0) (log 10.0))))) 323 324(defun not-done (err f eps machine-eps) 325 (> (* machine-eps err) (* eps (max (abs f) 1)))) 326 327;;(defmethod epsilon ((x integer)) 0) 328 329(in-package :maxima) 330 331(defun nfloat (e subs digits max-digits) 332 (let ((z (list nil nil)) (dig digits) (eps) (machine-epsilon nil)) 333 (cond ((or (mbagp e) (mrelationp e) ($setp e)) 334 (simplify (cons (list (caar e)) 335 (mapcar #'(lambda (s) (nfloat s subs digits max-digits)) (margs e))))) 336 337 (t 338 (catch 'nfloat-nounform-return 339 (setq e (nformat e)) 340 (setq eps (expt 10.0 (- digits))) 341 (setq eps (/ eps (- 1 eps))) 342 (while (and (or (null (first z)) (bigfloat::not-done (second z) (first z) eps machine-epsilon)) 343 (< digits max-digits)) 344 (bind-fpprec digits 345 (setq z (bigfloat::running-error-eval e subs fpprec)) 346 (setq machine-epsilon 347 (cond ((not (second z)) nil) 348 ((integerp (second z)) 0) 349 (t (bigfloat::epsilon (second z))))) 350 351 (setq digits (* 2 digits)))) 352 353 (if (or (null (first z)) (>= digits max-digits)) 354 (merror "Unable to evaluate to requested number of digits") 355 (maxima::bind-fpprec dig (values (maxima::to (first z)) (maxima::to (second z)))))))))) 356 357(setf (get '$nfloat 'operators) 'simp-nfloat) 358 359(defun simp-nfloat (x yy z) 360 (declare (ignore yy)) 361 (declare (special $max_fpprec)) 362 (pop x) ;; remove ($nfloat) 363 (let* ((e (if x (simpcheck (pop x) z) (wna-err '$nfloat))) 364 (subs (if x (simpcheck (pop x) z) (take '(mlist)))) 365 (digits (if x (simpcheck (pop x) z) $fpprec)) 366 (max-digits (if x (simpcheck (pop x) z) $max_fpprec)) 367 (f)) 368 369 (cond ((and ($listp subs) 370 (every #'(lambda (s) (and (mequalp s) (symbolp ($lhs s)) (complex-number-p ($rhs s) 'mnump))) 371 (cdr subs))) 372 (cond ((or (mbagp e) (mrelationp e) ($setp e)) 373 (simplify (cons (list (caar e)) 374 (mapcar #'(lambda (s) (take '($nfloat) s subs digits max-digits)) 375 (margs e))))) 376 (t 377 (setq f (nfloat e subs digits max-digits)) 378 (if (complex-number-p f 'bigfloat-or-number-p) f 379 `(($nfloat simp) ,e ,subs ,digits ,$max_fpprec))))) 380 (t `(($nfloat simp) ,e ,subs ,digits ,$max_fpprec))))) 381