This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "OCT: A portable Lisp implementation for quad-double precision floats".

The branch, master has been updated via 46e7193fb5b2fad35e67366ff375d9ebbb524c5c (commit) via 8227b233c4ba1f58d7928cdd5020201a5465c10d (commit) via 9d3c8cf17f1d644be01749fe91484c514eb80d3c (commit) via b725d6ef804b369ede3bc4bf037bd746a5d32406 (commit) from a40062c6860e8e778c1b15a48c2687c6cd2f5334 (commit)

Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below.

- Log ----------------------------------------------------------------- commit 46e7193fb5b2fad35e67366ff375d9ebbb524c5c Author: Raymond Toy toy.raymond@gmail.com Date: Thu Mar 17 13:37:29 2011 -0400

Implement exponential integral E.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 4781167..a25bc7d 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -333,3 +333,11 @@ (+ 1 (/ (incomplete-gamma 1/2 (* z z)) (sqrt (float-pi z)))))) + +(defun exp-integral-e (v z) + "Exponential integral E: + + E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)" + ;; Wolfram gives E(v,z) = z^(v-1)*gamma_incomplete_tail(1-v,z) + (* (expt z (- v 1)) + (incomplete-gamma-tail (- 1 v) z)))

commit 8227b233c4ba1f58d7928cdd5020201a5465c10d Author: Raymond Toy toy.raymond@gmail.com Date: Thu Mar 17 13:20:43 2011 -0400

Implement erfc; fix erf(-z), add comments.

o Was returning the wrong value for erf(-z). Use erf(-z) = - erf(z). o Add implementation of erfc. o Document code and algorithms a bit better.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp index dbdd4c7..4781167 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -267,6 +267,9 @@

;; Tail of the incomplete gamma function. (defun incomplete-gamma-tail (a z) + "Tail of the incomplete gamma function defined by: + + integrate(t^(a-1)*exp(-t), t, z, inf)" (let* ((prec (float-contagion a z)) (a (apply-contagion a prec)) (z (apply-contagion z prec))) @@ -279,6 +282,9 @@ (cf-incomplete-gamma-tail a z))))

(defun incomplete-gamma (a z) + "Incomplete gamma function defined by: + + integrate(t^(a-1)*exp(-t), t, 0, z)" (let* ((prec (float-contagion a z)) (a (apply-contagion a prec)) (z (apply-contagion z prec))) @@ -289,5 +295,41 @@ (cf-incomplete-gamma a z))))

(defun erf (z) - (/ (incomplete-gamma 1/2 (* z z)) - (sqrt (float-pi z)))) + "Error function: + + erf(z) = 2/sqrt(%pi)*sum((-1)^k*z^(2*k+1)/k!/(2*k+1), k, 0, inf) + + For real z, this is equivalent to + + erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z) for real z." + ;; + ;; Erf is an odd function: erf(-z) = -erf(z) + (if (minusp (realpart z)) + (- (erf (- z))) + (/ (incomplete-gamma 1/2 (* z z)) + (sqrt (float-pi z))))) + +(defun erfc (z) + "Complementary error function: + + erfc(z) = 1 - erf(z)" + ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is + ;; near 1. Wolfram says + ;; + ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2)) + ;; + ;; For real(z) > 0, sqrt(z^2)/z is 1 so + ;; + ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)) + ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2) + ;; + ;; For real(z) < 0, sqrt(z^2)/z is -1 so + ;; + ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)) + ;; = 1 + 1/sqrt(pi)*gamma_incomplete(1/2,z^2) + (if (>= (realpart z) 0) + (/ (incomplete-gamma-tail 1/2 (* z z)) + (sqrt (float-pi z))) + (+ 1 + (/ (incomplete-gamma 1/2 (* z z)) + (sqrt (float-pi z))))))

commit 9d3c8cf17f1d644be01749fe91484c514eb80d3c Author: Raymond Toy toy.raymond@gmail.com Date: Thu Mar 17 11:41:45 2011 -0400

Add implementation for incomplete gamma and erf, with tests.

qd-gamma.lisp: o Add implementation for Lentz's algorithm for evaluating continued fractions. o Implement incomplete-gamma and incomplete-gamma-tail using continued fractions. o Implement erf

rt-tests.lisp: o Add tests

diff --git a/qd-gamma.lisp b/qd-gamma.lisp index bb9004a..dbdd4c7 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -168,3 +168,126 @@ (defmethod gamma ((z qd-complex)) (gamma-aux z 39 32))

+ +;; Lentz's algorithm for evaluating continued fractions. +;; +;; Let the continued fraction be: +;; +;; a1 a2 a3 +;; b0 + ---- ---- ---- +;; b1 + b2 + b3 + +;; +(defun lentz (bf af) + (flet ((value-or-tiny (v) + (if (zerop v) + (etypecase v + ((or double-float cl:complex) + least-positive-normalized-double-float) + ((or qd-real qd-complex) + (make-qd least-positive-normalized-double-float))) + v))) + (let* ((f (value-or-tiny (funcall bf 0))) + (c f) + (d 0) + (eps (epsilon f))) + (loop + for j from 1 + for an = (funcall af j) + for bn = (funcall bf j) + do (progn + (setf d (value-or-tiny (+ bn (* an d)))) + (setf c (value-or-tiny (+ bn (/ an c)))) + (setf d (/ d)) + (let ((delta (* c d))) + (setf f (* f delta)) + (when (<= (abs (- delta 1)) eps) + (return))))) + f))) + +;; Continued fraction for erf(b): +;; +;; z[n] = 1+2*n-2*z^2 +;; a[n] = 4*n*z^2 +;; +;; This works ok, but has problems for z > 3 where sometimes the +;; result is greater than 1. +#+nil +(defun erf (z) + (let* ((z2 (* z z)) + (twoz2 (* 2 z2))) + (* (/ (* 2 z) + (sqrt (float-pi z))) + (exp (- z2)) + (/ (lentz #'(lambda (n) + (- (+ 1 (* 2 n)) + twoz2)) + #'(lambda (n) + (* 4 n z2))))))) + +;; Tail of the incomplete gamma function: +;; integrate(x^(a-1)*exp(-x), x, z, inf) +;; +;; The continued fraction, valid for all z except the negative real +;; axis: +;; +;; b[n] = 1+2*n+z-a +;; a[n] = n*(a-n) +;; +;; See http://functions.wolfram.com/06.06.10.0003.01 +(defun cf-incomplete-gamma-tail (a z) + (/ (* (expt z a) + (exp (- z))) + (let ((z-a (- z a))) + (lentz #'(lambda (n) + (+ n n 1 z-a)) + #'(lambda (n) + (* n (- a n))))))) + +;; Incomplete gamma function: +;; integrate(x^(a-1)*exp(-x), x, 0, z) +;; +;; The continued fraction, valid for all z except the negative real +;; axis: +;; +;; b[n] = n - 1 + z + a +;; a[n] = -z*(a + n) +;; +;; See http://functions.wolfram.com/06.06.10.0005.01. We modified the +;; continued fraction slightly and discarded the first quotient from +;; the fraction. +(defun cf-incomplete-gamma (a z) + (/ (* (expt z a) + (exp (- z))) + (let ((za1 (+ z a 1))) + (- a (/ (* a z) + (lentz #'(lambda (n) + (+ n za1)) + #'(lambda (n) + (- (* z (+ a n)))))))))) + +;; Tail of the incomplete gamma function. +(defun incomplete-gamma-tail (a z) + (let* ((prec (float-contagion a z)) + (a (apply-contagion a prec)) + (z (apply-contagion z prec))) + (if (and (realp a) (realp z)) + ;; For real values, we split the result to compute either the + ;; tail directly or from gamma(a) - incomplete-gamma + (if (> z (- a 1)) + (cf-incomplete-gamma-tail a z) + (- (gamma a) (cf-incomplete-gamma a z))) + (cf-incomplete-gamma-tail a z)))) + +(defun incomplete-gamma (a z) + (let* ((prec (float-contagion a z)) + (a (apply-contagion a prec)) + (z (apply-contagion z prec))) + (if (and (realp a) (realp z)) + (if (< z (- a 1)) + (cf-incomplete-gamma a z) + (- (gamma a) (cf-incomplete-gamma-tail a z))) + (cf-incomplete-gamma a z)))) + +(defun erf (z) + (/ (incomplete-gamma 1/2 (* z z)) + (sqrt (float-pi z)))) diff --git a/rt-tests.lisp b/rt-tests.lisp index 59d798f..5f4cdd5 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1111,7 +1111,7 @@ for m = (random #q1) for t3 = (elliptic-theta-2 0 (elliptic-nome m)) for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m))) - for result = (check-accuracy 206 t3 true) + for result = (check-accuracy 205 t3 true) when result append (list (list (list k m) result))) nil) @@ -1195,3 +1195,33 @@ when result append (list (list (list k y) result))) nil) + +;; gamma_incomplete(2,z) = integrate(t*exp(-t), t, z, inf) +;; = (z+1)*exp(-z) +(rt:deftest gamma-incomplete-tail.1.d + (let* ((z 5d0) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 52 gi true)) + nil) + +(rt:deftest gamma-incomplete-tail.2.d + (let* ((z #c(1 5d0)) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 50 gi true)) + nil) + +(rt:deftest gamma-incomplete-tail.1.q + (let* ((z 5d0) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 212 gi true)) + nil) + +(rt:deftest gamma-incomplete-tail.1.q + (let* ((z #q(1 5)) + (gi (incomplete-gamma-tail 2 z)) + (true (* (+ z 1) (exp (- z))))) + (check-accuracy 206 gi true)) + nil)

commit b725d6ef804b369ede3bc4bf037bd746a5d32406 Author: Raymond Toy toy.raymond@gmail.com Date: Thu Mar 17 11:38:43 2011 -0400

Add missing methods for QEXPT.

o We weren't handling the case of (expt real qd-complex). Add a method for this and other missing methods for QEXPT. o Move float contagion stuff from the end to the beginning so we can use it in this file.

diff --git a/qd-methods.lisp b/qd-methods.lisp index 9341311..cbf0daa 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -29,6 +29,59 @@ ;; We should do something better than this. (make-instance 'qd-real :value (rational-to-qd x)))

+;; Determine which of x and y has the higher precision and return the +;; value of the higher precision number. If both x and y are +;; rationals, just return 1f0, for a single-float value. +(defun float-contagion-2 (x y) + (etypecase x + (cl:rational + (etypecase y + (cl:rational + 1f0) + (cl:float + y) + (qd-real + y))) + (single-float + (etypecase y + ((or cl:rational single-float) + x) + ((or double-float qd-real) + y))) + (double-float + (etypecase y + ((or cl:rational single-float double-float) + x) + (qd-real + y))) + (qd-real + x))) + +;; Return a floating point (or complex) type of the highest precision +;; among all of the given arguments. +(defun float-contagion (&rest args) + ;; It would be easy if we could just add the args together and let + ;; normal contagion do the work, but we could easily introduce + ;; overflows or other errors that way. So look at each argument and + ;; determine the precision and choose the highest precision. + (let ((complexp (some #'complexp args)) + (max-type + (etypecase (reduce #'float-contagion-2 (mapcar #'realpart (if (cdr args) + args + (list (car args) 0)))) + (single-float 'single-float) + (double-float 'double-float) + (qd-real 'qd-real)))) + max-type)) + +(defun apply-contagion (number precision) + (etypecase number + ((or cl:real qd-real) + (coerce number precision)) + ((or cl:complex qd-complex) + (complex (coerce (realpart number) precision) + (coerce (imagpart number) precision))))) +

(defmethod add1 ((a number)) (cl::1+ a)) @@ -425,10 +478,13 @@ underlying floating-point format" (defmethod qexpt ((x number) (y number)) (cl:expt x y))

-(defmethod qexpt ((x qd-real) (y real)) - (exp (* y (log x)))) +(defmethod qexpt ((x number) (y qd-real)) + (exp (* y (log (apply-contagion x 'qd-real))))) + +(defmethod qexpt ((x number) (y qd-complex)) + (exp (* y (log (apply-contagion x 'qd-real)))))

-(defmethod qexpt ((x real) (y qd-real)) +(defmethod qexpt ((x qd-real) (y real)) (exp (* y (log x))))

(defmethod qexpt ((x qd-real) (y cl:complex)) @@ -437,12 +493,6 @@ underlying floating-point format" :imag (qd-value (imagpart y))) (log x))))

-(defmethod qexpt ((x cl:complex) (y qd-real)) - (exp (* y - (log (make-instance 'qd-complex - :real (qd-value (realpart x)) - :imag (qd-value (imagpart x))))))) - (defmethod qexpt ((x qd-real) (y qd-real)) ;; x^y = exp(y*log(x)) (exp (* y (log x)))) @@ -451,6 +501,15 @@ underlying floating-point format" (make-instance 'qd-real :value (npow (qd-value x) y)))

+(defmethod qexpt ((x qd-complex) (y number)) + (exp (* y (log x)))) + +(defmethod qexpt ((x qd-complex) (y qd-real)) + (exp (* y (log x)))) + +(defmethod qexpt ((x qd-complex) (y qd-complex)) + (exp (* y (log x)))) + (declaim (inline expt)) (defun expt (x y) (qexpt x y)) @@ -1058,56 +1117,3 @@ the same precision as the argument. The argument can be complex.")) (defmethod float-pi ((z qd-complex)) +pi+)

-;; Determine which of x and y has the higher precision and return the -;; value of the higher precision number. If both x and y are -;; rationals, just return 1f0, for a single-float value. -(defun float-contagion-2 (x y) - (etypecase x - (cl:rational - (etypecase y - (cl:rational - 1f0) - (cl:float - y) - (qd-real - y))) - (single-float - (etypecase y - ((or cl:rational single-float) - x) - ((or double-float qd-real) - y))) - (double-float - (etypecase y - ((or cl:rational single-float double-float) - x) - (qd-real - y))) - (qd-real - x))) - -;; Return a floating point (or complex) type of the highest precision -;; among all of the given arguments. -(defun float-contagion (&rest args) - ;; It would be easy if we could just add the args together and let - ;; normal contagion do the work, but we could easily introduce - ;; overflows or other errors that way. So look at each argument and - ;; determine the precision and choose the highest precision. - (let ((complexp (some #'complexp args)) - (max-type - (etypecase (reduce #'float-contagion-2 (mapcar #'realpart (if (cdr args) - args - (list (car args) 0)))) - (single-float 'single-float) - (double-float 'double-float) - (qd-real 'qd-real)))) - max-type)) - -(defun apply-contagion (number precision) - (etypecase number - ((or cl:real qd-real) - (coerce number precision)) - ((or cl:complex qd-complex) - (complex (coerce (realpart number) precision) - (coerce (imagpart number) precision))))) -

-----------------------------------------------------------------------

Summary of changes: qd-gamma.lisp | 173 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ qd-methods.lisp | 130 ++++++++++++++++++++++-------------------- rt-tests.lisp | 32 ++++++++++- 3 files changed, 272 insertions(+), 63 deletions(-)

hooks/post-receive