share/hypergeometric/nfloat.lisp

   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)))))