December 2011 - oct-cvs - mailman3.common-lisp.net

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. c388f81713d7b2a483000d3cee1af030ed2c1cac
by Raymond Toy 06 Dec '11

06 Dec '11

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 c388f81713d7b2a483000d3cee1af030ed2c1cac (commit) from 26ef83a4543cdca6f752abca6f68e9d7f586c68e (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 c388f81713d7b2a483000d3cee1af030ed2c1cac Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 19:58:44 2011 -0800 Better exp-integral-e computation and fix for incomplete-gamma-tail. For exp-integral-e, use the series for small z and the incomplete-gamma-tail for near the negative real axis. Otherwise, use the continued fraction. In incomplete-gamma-tail, we were using the continued fraction instead of the incomplete-gamma function for the region just below the negative real axis. We should use the cf except in that region. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index f13b9a8..48f8345 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -396,7 +396,7 @@ ;; Use the incomplete gamma function to evaluate in this ;; region. (Arbitrarily selected the region to be a sector. ;; But what is the correct size of this sector?) - (if (<= (phase z) 3.1) + (if (<= (abs (phase z)) 3.1) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z))))))) @@ -460,15 +460,7 @@ (/ (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)" - ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); - ;; - ;; for |arg(z)| < pi. - ;; - ;; +(defun cf-exp-integral-e (v z) ;; We use the continued fraction ;; ;; E(v,z) = exp(-z)/cf(z) @@ -485,7 +477,45 @@ (+ z+v (* 2 k))) #'(lambda (k) (* (- k) - (1- (+ k v)))))))) + (+ k v -1))))))) + +(defun s-exp-integral-e (v z) + ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf) + (let ((-z (- z)) + (-v (- v)) + (eps (epsilon z))) + (loop for k from 0 + for term = 1 then (* term (/ -z k)) + for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v))) + when (< (abs term) (* (abs sum) eps)) + return (- (* (gamma (- 1 v)) (expt z (- v 1))) + sum)))) + +(defun exp-integral-e (v z) + "Exponential integral E: + + E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)" + ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); + ;; + ;; for |arg(z)| < pi. + ;; + ;; + (cond ((< (abs z) 1) + ;; Use series for small z + (s-exp-integral-e v z)) + ((>= (abs (phase z)) 3.1) + ;; The continued fraction doesn't converge on the negative + ;; real axis, and converges very slowly near the negative + ;; real axis, so use the incomplete-gamma-tail function in + ;; this region. "Closeness" to the negative real axis is + ;; teken to mean that z is in a sector near the axis. + ;; + ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z) + (* (expt z (- v 1)) + (incomplete-gamma-tail (- 1 v) z))) + (t + ;; Use continued fraction for everything else. + (cf-exp-integral-e v z)))) ;; Series for Fresnel S ;; ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 52 +++++++++++++++++++++++++++++++++++++++++----------- 1 files changed, 41 insertions(+), 11 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 26ef83a4543cdca6f752abca6f68e9d7f586c68e
by Raymond Toy 05 Dec '11

05 Dec '11

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 26ef83a4543cdca6f752abca6f68e9d7f586c68e (commit) from 06f8a53123c7c6e204099fe81e05b66611ab9ca5 (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 26ef83a4543cdca6f752abca6f68e9d7f586c68e Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 14:06:39 2011 -0800 Oops. Put back QEXPT method for QD-COMPLEX, QD-REAL. diff --git a/qd-complex.lisp b/qd-complex.lisp index 409e94f..6be9983 100644 --- a/qd-complex.lisp +++ b/qd-complex.lisp @@ -742,6 +742,9 @@ Z may be any number, but the result is always a complex." (defmethod qexpt ((x qd-complex) (y number)) (exp (* (float y #q0) (log x)))) +(defmethod qexpt ((x qd-complex) (y qd-real)) + (exp (* y (log x)))) + (defmethod qexpt ((x number) (y qd-complex)) (exp (* y (log (float x #q0))))) ----------------------------------------------------------------------- Summary of changes: qd-complex.lisp | 3 +++ 1 files changed, 3 insertions(+), 0 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 06f8a53123c7c6e204099fe81e05b66611ab9ca5
by Raymond Toy 05 Dec '11

05 Dec '11

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 06f8a53123c7c6e204099fe81e05b66611ab9ca5 (commit) via 2aece6ff3951520cf025b902045a012de32b52ff (commit) from 612465af4decd8a0dfd4dbc199f0d48503e4b31a (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 06f8a53123c7c6e204099fe81e05b66611ab9ca5 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 12:48:30 2011 -0800 Clean up code: o Get rid of unused var in FLOAT-CONTAGION. o QEXPT for qd-complex first arg should be defined in qd-complex.lisp. diff --git a/qd-methods.lisp b/qd-methods.lisp index bb851c4..0eef82d 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -64,15 +64,12 @@ ;; 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)) + (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))) (defun apply-contagion (number precision) (etypecase number @@ -501,15 +498,6 @@ 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)) commit 2aece6ff3951520cf025b902045a012de32b52ff Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 12:47:02 2011 -0800 Get rid of compiler warning about unused variable in theta-1 and theta-3. diff --git a/qd-theta.lisp b/qd-theta.lisp index fb41b82..8b160c1 100644 --- a/qd-theta.lisp +++ b/qd-theta.lisp @@ -84,6 +84,7 @@ (2q^1/4 (* 2 (sqrt (sqrt q))))) (3by3rec s s 0 #'(lambda (k matfun) + (declare (ignore k)) (funcall matfun (setf -2q^2ncos (* q^2 -2q^2ncos)) 1 @@ -111,6 +112,7 @@ (cos (cos (* 2 z)))) (3by3rec (* q cos) 1 1 #'(lambda (k matfun) + (declare (ignore k)) (funcall matfun (* 2 (* (setf q^2k (* q^2 q^2k)) q cos)) 1 ----------------------------------------------------------------------- Summary of changes: qd-methods.lisp | 24 ++++++------------------ qd-theta.lisp | 2 ++ 2 files changed, 8 insertions(+), 18 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 612465af4decd8a0dfd4dbc199f0d48503e4b31a
by Raymond Toy 05 Dec '11

05 Dec '11

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 612465af4decd8a0dfd4dbc199f0d48503e4b31a (commit) from 1efd5f22d66198240c6952ab8b59b0d5436e8519 (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 612465af4decd8a0dfd4dbc199f0d48503e4b31a Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 11:29:30 2011 -0800 Reduce required accuracy so CCL can pass oct.theta2.1.d and oct.elliptic-pi.n1.d diff --git a/rt-tests.lisp b/rt-tests.lisp index 736f7f8..33b2444 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -816,7 +816,7 @@ (check-signs #'atan #q(-0d0 2d0) true)) t) -;; Test x < -1 +;; Test x < -1. CLHS says for x < -1, atanh is continuous with quadrant III. (rt:deftest oct.atanh-branch-neg.1 (let ((true (cl:atanh #c(-2d0 -1d-20)))) (check-signs #'atanh -2d0 true)) @@ -827,7 +827,7 @@ (check-signs #'atanh #q-2 true)) t) -;; Test x > 1 +;; Test x > 1. CLHS says for x > 1, atanh is continus with quadrant I. (rt:deftest oct.atanh-branch-pos.1 (let ((true (cl:atanh #c(2d0 1d-20)))) (check-signs #'atanh 2d0 true)) @@ -1078,7 +1078,7 @@ for phi = (random (/ pi 2)) for epi = (elliptic-pi 1 phi 0) for true = (tan phi) - for result = (check-accuracy 36 epi true) + for result = (check-accuracy 34.5 epi true) unless (eq nil result) append (list (list (list k phi) result))) nil) @@ -1167,7 +1167,7 @@ for m = (random 1d0) 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 49 t3 true) + for result = (check-accuracy 43.5 t3 true) when result append (list (list (list k m) result))) nil) ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 8 ++++---- 1 files changed, 4 insertions(+), 4 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 1efd5f22d66198240c6952ab8b59b0d5436e8519
by Raymond Toy 05 Dec '11

05 Dec '11

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 1efd5f22d66198240c6952ab8b59b0d5436e8519 (commit) from a0b50beda0b7bee4268de6c4a8ef8ff4136f33df (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 1efd5f22d66198240c6952ab8b59b0d5436e8519 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Dec 5 09:45:54 2011 -0800 Update version (and make it compatible with latest asdf2). diff --git a/oct.asd b/oct.asd index a6111ca..ed7a586 100644 --- a/oct.asd +++ b/oct.asd @@ -36,7 +36,7 @@ :author "Raymond Toy" :maintainer "See <http://www.common-lisp.net/project/oct>" :licence "MIT" - :version "2011-02-09" ; Just use the date + :version "2011.12.05" ; Just use the date :components ((:file "qd-package") (:file "qd-rep" :depends-on ("qd-package")) @@ -69,12 +69,12 @@ :depends-on ("qd-methods" "qd-reader")) )) -(defmethod perform ((op test-op) (c (eql (find-system :oct)))) +(defmethod perform ((op test-op) (c (eql (asdf:find-system :oct)))) (oos 'test-op 'oct-tests)) (asdf:defsystem oct-tests :depends-on (oct) - :version "2011-02-09" ; Just use the date + :version "2011.12.05" ; Just use the date :in-order-to ((compile-op (load-op :rt)) (test-op (load-op :rt :oct))) :components @@ -82,7 +82,7 @@ (:file "qd-test") (:file "rt-tests"))) -(defmethod perform ((op test-op) (c (eql (find-system :oct-tests)))) +(defmethod perform ((op test-op) (c (eql (asdf:find-system :oct-tests)))) (or (funcall (intern "DO-TESTS" (find-package "RT"))) (error "TEST-OP failed for OCT-TESTS"))) ----------------------------------------------------------------------- Summary of changes: oct.asd | 8 ++++---- 1 files changed, 4 insertions(+), 4 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. a0b50beda0b7bee4268de6c4a8ef8ff4136f33df
by Raymond Toy 05 Dec '11

05 Dec '11

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 a0b50beda0b7bee4268de6c4a8ef8ff4136f33df (commit) from 6b8a8a6fe864050f9c6371150d8070f9b38fe76e (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 a0b50beda0b7bee4268de6c4a8ef8ff4136f33df Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Dec 4 23:24:06 2011 -0800 Remove some random junk that we left in. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 8e523ed..f13b9a8 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -649,23 +649,3 @@ (if (and (realp z) (plusp z)) (realpart (ci z)) (ci z)))) - -;;; Exponential integral e defined by -;;; -;;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); -;;; -;;; for |arg(z)| < pi. -;;; -;;; -;;; We use the continued fraction -;;; -;;; E(v,z) = exp(-z)/cf(z) -;;; -;; where the continued fraction cf(z) is -;; -;; a[k] = -k*(k+v-1) -;; b[k] = v + 2*k + z -;; -;; for k = 1, inf - -(defun expintegral-e (v z) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 20 -------------------- 1 files changed, 0 insertions(+), 20 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 6b8a8a6fe864050f9c6371150d8070f9b38fe76e
by Raymond Toy 05 Dec '11

05 Dec '11

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 6b8a8a6fe864050f9c6371150d8070f9b38fe76e (commit) from b609574d330c7b0c9a6de588d3884a00f9aad13b (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 6b8a8a6fe864050f9c6371150d8070f9b38fe76e Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Dec 4 23:22:01 2011 -0800 Use exp-integral-e for (incomplete-gamma-tail 0 z). We already have exp-integral-e function so move expintegral-e implementation to exp-integral-e. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 05f85d4..8e523ed 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -375,27 +375,30 @@ (let* ((prec (float-contagion a z)) (a (apply-contagion a prec)) (z (apply-contagion z prec))) - (if (and (zerop (imagpart a)) - (zerop (imagpart z))) - ;; For real values, we split the result to compute either the - ;; tail directly from the continued fraction or from gamma(a) - ;; - incomplete-gamma. The continued fraction doesn't - ;; converge on the negative real axis, so we can't use that - ;; there. And accuracy appears to be better if z is "small". - ;; We take this to mean |z| < |a-1|. Note that |a-1| is the - ;; peak of the integrand. - (if (and (> (abs z) (abs (- a 1))) - (not (minusp (realpart z)))) - (cf-incomplete-gamma-tail a z) - (- (gamma a) (cf-incomplete-gamma a z))) - ;; If the argument is close enough to the negative real axis, - ;; the continued fraction for the tail is not very accurate. - ;; Use the incomplete gamma function to evaluate in this - ;; region. (Arbitrarily selected the region to be a sector. - ;; But what is the correct size of this sector?) - (if (<= (phase z) 3.1) - (cf-incomplete-gamma-tail a z) - (- (gamma a) (cf-incomplete-gamma a z)))))) + (if (zerop a) + ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z) + (exp-integral-e 1 z) + (if (and (zerop (imagpart a)) + (zerop (imagpart z))) + ;; For real values, we split the result to compute either the + ;; tail directly from the continued fraction or from gamma(a) + ;; - incomplete-gamma. The continued fraction doesn't + ;; converge on the negative real axis, so we can't use that + ;; there. And accuracy appears to be better if z is "small". + ;; We take this to mean |z| < |a-1|. Note that |a-1| is the + ;; peak of the integrand. + (if (and (> (abs z) (abs (- a 1))) + (not (minusp (realpart z)))) + (cf-incomplete-gamma-tail a z) + (- (gamma a) (cf-incomplete-gamma a z))) + ;; If the argument is close enough to the negative real axis, + ;; the continued fraction for the tail is not very accurate. + ;; Use the incomplete gamma function to evaluate in this + ;; region. (Arbitrarily selected the region to be a sector. + ;; But what is the correct size of this sector?) + (if (<= (phase z) 3.1) + (cf-incomplete-gamma-tail a z) + (- (gamma a) (cf-incomplete-gamma a z))))))) (defun incomplete-gamma (a z) "Incomplete gamma function defined by: @@ -461,9 +464,28 @@ "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))) + ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); + ;; + ;; for |arg(z)| < pi. + ;; + ;; + ;; We use the continued fraction + ;; + ;; E(v,z) = exp(-z)/cf(z) + ;; + ;; where the continued fraction cf(z) is + ;; + ;; a[k] = -k*(k+v-1) + ;; b[k] = v + 2*k + z + ;; + ;; for k = 1, inf + (let ((z+v (+ z v))) + (/ (exp (- z)) + (lentz #'(lambda (k) + (+ z+v (* 2 k))) + #'(lambda (k) + (* (- k) + (1- (+ k v)))))))) ;; Series for Fresnel S ;; @@ -647,10 +669,3 @@ ;; for k = 1, inf (defun expintegral-e (v z) - (let ((z+v (+ z v))) - (/ (exp (- z)) - (lentz #'(lambda (k) - (+ z+v (* 2 k))) - #'(lambda (k) - (* (- k) - (1- (+ k v)))))))) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 77 ++++++++++++++++++++++++++++++++++----------------------- 1 files changed, 46 insertions(+), 31 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. b609574d330c7b0c9a6de588d3884a00f9aad13b
by Raymond Toy 05 Dec '11

05 Dec '11

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 b609574d330c7b0c9a6de588d3884a00f9aad13b (commit) via 5079d741cabce7dbaa84410fdfe565dfff03a51c (commit) from 5da3f9d5ccf7ddc9fe3f9334391984426acaeb9d (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 b609574d330c7b0c9a6de588d3884a00f9aad13b Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Dec 4 21:29:17 2011 -0800 Signal error for gamma of negative integers. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 8cab1bb..05f85d4 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -134,14 +134,18 @@ (defun gamma-aux (z limit nterms) (let ((precision (float-contagion z))) - (cond ((minusp (realpart z)) + (cond ((<= (realpart z) 0) ;; Use reflection formula if realpart(z) < 0: ;; gamma(-z) = -pi*csc(pi*z)/gamma(z+1) ;; or ;; gamma(z) = pi*csc(pi*z)/gamma(1-z) - (/ (float-pi z) - (sin (* (float-pi z) z)) - (gamma-aux (- 1 z) limit nterms))) + (if (and (realp z) + (= (truncate z) z)) + ;; Gamma of a negative integer is infinity. Signal an error + (error "Gamma of non-positive integer ~S" z) + (/ (float-pi z) + (sin (* (float-pi z) z)) + (gamma-aux (- 1 z) limit nterms)))) ((and (zerop (imagpart z)) (= z (truncate z))) ;; We have gamma(n) where an integer value n and is small @@ -623,4 +627,30 @@ (if (and (realp z) (plusp z)) (realpart (ci z)) (ci z)))) - \ No newline at end of file + +;;; Exponential integral e defined by +;;; +;;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); +;;; +;;; for |arg(z)| < pi. +;;; +;;; +;;; We use the continued fraction +;;; +;;; E(v,z) = exp(-z)/cf(z) +;;; +;; where the continued fraction cf(z) is +;; +;; a[k] = -k*(k+v-1) +;; b[k] = v + 2*k + z +;; +;; for k = 1, inf + +(defun expintegral-e (v z) + (let ((z+v (+ z v))) + (/ (exp (- z)) + (lentz #'(lambda (k) + (+ z+v (* 2 k))) + #'(lambda (k) + (* (- k) + (1- (+ k v)))))))) commit 5079d741cabce7dbaa84410fdfe565dfff03a51c Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Dec 4 21:28:27 2011 -0800 Fix bug in CEILING and FCEILING; add tests. diff --git a/qd-methods.lisp b/qd-methods.lisp index e4df131..bb851c4 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -652,7 +652,7 @@ underlying floating-point format" (multiple-value-bind (f rem) (floor x y) (if (zerop rem) - (values (+ f 1) + (values f rem) (values (+ f 1) (- rem 1))))) @@ -661,7 +661,7 @@ underlying floating-point format" (multiple-value-bind (f rem) (ffloor x y) (if (zerop rem) - (values (+ f 1) + (values f rem) (values (+ f 1) (- rem 1))))) diff --git a/rt-tests.lisp b/rt-tests.lisp index acc8b74..736f7f8 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -55,6 +55,74 @@ ;;; Some simple tests from the Yozo Hida's qd package. +(rt:deftest ceiling-d.1 + (multiple-value-list (ceiling -50d0)) + (-50 0d0)) + +(rt:deftest ceiling-d.2 + (let ((z -50.1d0)) + (multiple-value-bind (res rem) + (ceiling -50.1d0) + (list res (= z (+ res rem))))) + (-50 t)) + +(rt:deftest ceiling-q.1 + (multiple-value-bind (res rem) + (ceiling #q-50q0) + (list res (zerop rem))) + (-50 t)) + +(rt:deftest ceiling-q.2 + (let ((z #q-50.1q0)) + (multiple-value-bind (res rem) + (ceiling z) + (list res (= z (+ res rem))))) + (-50 t)) + +(rt:deftest truncate-d.1 + (multiple-value-list (truncate -50d0)) + (-50 0d0)) + +(rt:deftest truncate-q.1 + (multiple-value-bind (res rem) + (truncate #q-50q0) + (list res (zerop rem))) + (-50 t)) + +(rt:deftest fceiling-d.1 + (multiple-value-list (fceiling -50d0)) + (-50d0 0d0)) + +(rt:deftest fceiling-d.2 + (let ((z -50.1d0)) + (multiple-value-bind (res rem) + (fceiling -50.1d0) + (list res (= z (+ res rem))))) + (-50d0 t)) + +(rt:deftest fceiling-q.1 + (multiple-value-bind (res rem) + (fceiling #q-50q0) + (list (= res -50) (zerop rem))) + (t t)) + +(rt:deftest fceiling-q.2 + (let ((z #q-50.1q0)) + (multiple-value-bind (res rem) + (fceiling z) + (list (= res -50) (= z (+ res rem))))) + (t t)) + +(rt:deftest ftruncate-d.1 + (multiple-value-list (ftruncate -50d0)) + (-50d0 0d0)) + +(rt:deftest ftruncate-q.1 + (multiple-value-bind (res rem) + (ftruncate #q-50q0) + (list (= res -50) (zerop rem))) + (t t)) + ;; Pi via Machin's formula (rt:deftest oct.pi.machin (let* ((*standard-output* *null*) @@ -1161,7 +1229,7 @@ for y = (+ 1 (random 100d0)) for g = (abs (gamma (complex 0 y))) for true = (sqrt (/ pi y (sinh (* pi y)))) - for result = (check-accuracy 45 g true) + for result = (check-accuracy 44 g true) when result append (list (list (list k y) result))) nil) @@ -1294,3 +1362,4 @@ (true (fresnel-s-series z))) (check-accuracy 212 s true)) nil) + ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 40 +++++++++++++++++++++++++++---- qd-methods.lisp | 4 +- rt-tests.lisp | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 107 insertions(+), 8 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0