[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 46e7193fb5b2fad35e67366ff375d9ebbb524c5c
Raymond Toy
rtoy at common-lisp.net
Thu Mar 17 17:37:55 UTC 2011
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 at 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 at 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 at 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 at 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
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list