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