oct-cvs December 2011

oct-cvs@common-lisp.net

1 participants
8 discussions

[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

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

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

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

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

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

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

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

oct-cvs December 2011