[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. cc23f1098858d652c66728ea85c1c331bfb0fbb8
Raymond Toy
rtoy at common-lisp.net
Tue Mar 22 13:13:47 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 cc23f1098858d652c66728ea85c1c331bfb0fbb8 (commit)
via 6b8dc51d96004406780ec1faa2cd4ce3b127c287 (commit)
via 97bc71a13186735474eba77043d685b0bd0a00d6 (commit)
via 50fc6412a3effd9bea8f3e481f7c35c3c911f856 (commit)
from d17228b86105b8a6ac652459b1100266ad0311b3 (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 cc23f1098858d652c66728ea85c1c331bfb0fbb8
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Mon Mar 21 19:39:24 2011 -0400
Move the uses of foo-t before the use of foo-t.
diff --git a/qd.lisp b/qd.lisp
index fb3e4e1..9f1ec3d 100644
--- a/qd.lisp
+++ b/qd.lisp
@@ -409,11 +409,6 @@ If TARGET is given, TARGET is destructively modified to contain the result."
;; which don't do a very good job with dataflow. CMUCL is one of
;; those compilers.
-(defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
- "Return the sum of the %QUAD-DOUBLE numbers A and B.
-If TARGET is given, TARGET is destructively modified to contain the result."
- (add-qd-t a b target))
-
(defun add-qd-t (a b target)
(declare (type %quad-double a b #+oct-array target)
@@ -476,6 +471,12 @@ If TARGET is given, TARGET is destructively modified to contain the result."
(%store-qd-d target (+ a0 b0) 0d0 0d0 0d0)
(%store-qd-d target s0 s1 s2 s3)))))))))))
+(defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
+ "Return the sum of the %QUAD-DOUBLE numbers A and B.
+If TARGET is given, TARGET is destructively modified to contain the result."
+ (add-qd-t a b target))
+
+
(defun neg-qd-t (a target)
(declare (type %quad-double a #+oct-array target)
#+(and cmu (not oct-array)) (ignore target))
@@ -666,11 +667,6 @@ If TARGET is given, TARGET is destructively modified to contain the result."
;; Clisp says
;; 14.142135623730950488016887242096980785696718753769480731766797379908L0
-(defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
- "Returns the product of the %QUAD-DOUBLE numbers A and B.
-If TARGET is given, TARGET is destructively modified to contain the result."
- (mul-qd-t a b target))
-
(defun mul-qd-t (a b target)
(declare (type %quad-double a b #+oct-array target)
(optimize (speed 3)
@@ -736,6 +732,10 @@ If TARGET is given, TARGET is destructively modified to contain the result."
(%store-qd-d target p0 0d0 0d0 0d0)
(%store-qd-d target r0 r1 s0 s1))))))))))))))
+(defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
+ "Returns the product of the %QUAD-DOUBLE numbers A and B.
+If TARGET is given, TARGET is destructively modified to contain the result."
+ (mul-qd-t a b target))
;; This is the non-sloppy version. I think this works just fine, but
;; since qd defaults to the sloppy multiplication version, we do the
@@ -838,11 +838,6 @@ If TARGET is given, TARGET is destructively modified to contain the result."
(multiple-value-call #'%make-qd-d
(renorm-5 p0 p1 s0 t0 t1))))))))))))))))))))
-(defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
- "Return the square of the %QUAD-DOUBLE number A. If TARGET is also given,
-it is destructively modified with the result."
- (sqr-qd-t a target))
-
(defun sqr-qd-t (a target)
"Square A"
(declare (type %quad-double a #+oct-array target)
@@ -890,12 +885,11 @@ it is destructively modified with the result."
(multiple-value-bind (a0 a1 a2 a3)
(renorm-5 p0 p1 p2 p3 p4)
(%store-qd-d target a0 a1 a2 a3)))))))))
-
-(defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
- "Return the quotient of the two %QUAD-DOUBLE numbers A and B.
-If TARGET is given, it destrutively modified with the result."
- (div-qd-t a b target))
+(defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
+ "Return the square of the %QUAD-DOUBLE number A. If TARGET is also given,
+it is destructively modified with the result."
+ (sqr-qd-t a target))
#+nil
(defun div-qd-t (a b target)
@@ -948,6 +942,11 @@ If TARGET is given, it destrutively modified with the result."
(renorm-4 q0 q1 q2 q3)
(%store-qd-d target q0 q1 q2 q3))))))))
+(defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0)))
+ "Return the quotient of the two %QUAD-DOUBLE numbers A and B.
+If TARGET is given, it destrutively modified with the result."
+ (div-qd-t a b target))
+
(declaim (inline invert-qd))
(defun invert-qd (v)
commit 6b8dc51d96004406780ec1faa2cd4ce3b127c287
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Sat Mar 19 11:38:43 2011 -0400
Add rem and mod functions that were left out.
diff --git a/qd-methods.lisp b/qd-methods.lisp
index cbf0daa..29b296d 100644
--- a/qd-methods.lisp
+++ b/qd-methods.lisp
@@ -671,6 +671,12 @@ underlying floating-point format"
(ceiling x y)
(floor x y)))
+(defun rem (x y)
+ (nth-value 1 (truncate x y)))
+
+(defun mod (x y)
+ (nth-value 1 (floor x y)))
+
(defun ftruncate (x &optional (y 1))
(if (minusp x)
(fceiling x y)
diff --git a/qd-package.lisp b/qd-package.lisp
index b102435..7db5cd3 100644
--- a/qd-package.lisp
+++ b/qd-package.lisp
@@ -188,6 +188,8 @@
#:ftruncate
#:round
#:fround
+ #:rem
+ #:mod
#:realpart
#:imagpart
#:conjugate
@@ -258,6 +260,8 @@
#:ftruncate
#:round
#:fround
+ #:rem
+ #:mod
#:realpart
#:imagpart
#:conjugate
commit 97bc71a13186735474eba77043d685b0bd0a00d6
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Mar 18 09:06:10 2011 -0400
Add series for incomplete-gamma for when the fraction is slow.
o Add series for incomplete gamma function for small a and z. Needed
because the continued fraction is slow in this range.
o In INCOMPLETE-GAMMA-TAIL, call INCOMPLETE-GAMMA instead of
CF-INCOMPLETE-GAMMA just in case a and z are small.
o In INCOMPLETE-GAMMA, use the series for small a and z.
o Simplify evaluation of Si(z) when z is real.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index a6fc7cf..c0cacaa 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -265,6 +265,21 @@
#'(lambda (n)
(- (* z (+ a n))))))))))
+;; Series expansion for incomplete gamma. Intended for |a|<1 and
+;; |z|<1. The series is
+;;
+;; g(a,z) = z^a * sum((-z)^k/k!/(a+k), k, 0, inf)
+(defun s-incomplete-gamma (a z)
+ (let ((-z (- z))
+ (eps (epsilon z)))
+ (loop for k from 0
+ for term = 1 then (* term (/ -z k))
+ for sum = (/ a) then (+ sum (/ term (+ a k)))
+ when (< (abs term) (* (abs sum) eps))
+ return (* sum (expt z a)))))
+
+
+
;; Tail of the incomplete gamma function.
(defun incomplete-gamma-tail (a z)
"Tail of the incomplete gamma function defined by:
@@ -276,9 +291,9 @@
(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))
+ (if (> (abs z) (abs (- a 1)))
(cf-incomplete-gamma-tail a z)
- (- (gamma a) (cf-incomplete-gamma a z)))
+ (- (gamma a) (incomplete-gamma a z)))
(cf-incomplete-gamma-tail a z))))
(defun incomplete-gamma (a z)
@@ -288,13 +303,18 @@
(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)))
- (if (< (abs z) (abs a))
- (cf-incomplete-gamma a z)
- (- (gamma a) (cf-incomplete-gamma-tail a z))))))
+ (if (and (< (abs a) 1) (< (abs z) 1))
+ (s-incomplete-gamma a z)
+ (if (and (realp a) (realp z))
+ (if (< z (- a 1))
+ (cf-incomplete-gamma a z)
+ (- (gamma a) (cf-incomplete-gamma-tail a z)))
+ ;; The continued fraction doesn't converge very fast if a
+ ;; and z are small. In this case, use the series
+ ;; expansion instead, which converges quite rapidly.
+ (if (< (abs z) (abs a))
+ (cf-incomplete-gamma a z)
+ (- (gamma a) (cf-incomplete-gamma-tail a z)))))))
(defun erf (z)
"Error function:
@@ -405,9 +425,20 @@
(- (log -iz)
(log iz)))))))
(if (realp z)
- (if (< z 0)
- (- (sin-integral (- z)))
- (si z))
+ ;; Si is odd and real for real z. In this case, we have
+ ;;
+ ;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi)
+ ;; = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x))
+ ;; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so
+ ;;
+ ;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x))
+ (cond ((< z 0)
+ (- (sin-integral (- z))))
+ ((= z 0)
+ (* 0 z))
+ (t
+ (+ (* 1/2 (float-pi z))
+ (imagpart (incomplete-gamma-tail 0 (complex 0 z))))))
(si z))))
(defun cos-integral (z)
commit 50fc6412a3effd9bea8f3e481f7c35c3c911f856
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Thu Mar 17 23:12:30 2011 -0400
Implement Si and Ci, sin and cos integrals.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index dfc87a6..a6fc7cf 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -387,4 +387,49 @@
(- (realpart (fs (- z))))
(realpart (fs z)))
(fs z)))))
-
\ No newline at end of file
+
+(defun sin-integral (z)
+ "Sin integral:
+
+ Si(z) = integrate(sin(t)/t, t, 0, z)"
+ ;; Wolfram has
+ ;;
+ ;; Si(z) = %i/2*(gamma_inc_tail(0, -%i*z) - gamma_inc_tail(0, %i*z) + log(-%i*z)-log(%i*z))
+ ;;
+ (flet ((si (z)
+ (* #c(0 1/2)
+ (let ((iz (* #c(0 1) z))
+ (-iz (* #c(0 -1) z)))
+ (+ (- (incomplete-gamma-tail 0 -iz)
+ (incomplete-gamma-tail 0 iz))
+ (- (log -iz)
+ (log iz)))))))
+ (if (realp z)
+ (if (< z 0)
+ (- (sin-integral (- z)))
+ (si z))
+ (si z))))
+
+(defun cos-integral (z)
+ "Cos integral:
+
+ Ci(z) = integrate((cos(t) - 1)/t, t, 0, z) + log(z) + gamma
+
+ where gamma is Euler-Mascheroni constant"
+ ;; Wolfram has
+ ;;
+ ;; Ci(z) = log(z) - 1/2*(gamma_inc_tail(0, -%i*z) + gamma_inc_tail(0, %i*z) + log(-%i*z)+log(%i*z))
+ ;;
+ (flet ((ci (z)
+ (- (log z)
+ (* 1/2
+ (let ((iz (* #c(0 1) z))
+ (-iz (* #c(0 -1) z)))
+ (+ (+ (incomplete-gamma-tail 0 -iz)
+ (incomplete-gamma-tail 0 iz))
+ (+ (log -iz)
+ (log iz))))))))
+ (if (and (realp z) (plusp z))
+ (realpart (ci z))
+ (ci z))))
+
\ No newline at end of file
-----------------------------------------------------------------------
Summary of changes:
qd-gamma.lisp | 96 +++++++++++++++++++++++++++++++++++++++++++++++++------
qd-methods.lisp | 6 +++
qd-package.lisp | 4 ++
qd.lisp | 39 +++++++++++-----------
4 files changed, 115 insertions(+), 30 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list