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 a40062c6860e8e778c1b15a48c2687c6cd2f5334 (commit) via a666f3923b82be0eb435ddd3e9b20697097fafe0 (commit) via a93aca9a6cfdb5b1305cb2570eea23e6f44b51c0 (commit) from c6f0a5ef00138c7f504432e0d0f485c23328274a (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 a40062c6860e8e778c1b15a48c2687c6cd2f5334 Author: Raymond Toy <toy.raymond@gmail.com> Date: Wed Mar 16 19:43:54 2011 -0400 Add gamma and log-gamma functions; work in progress. oct.asd: o Add qd-gamma.lisp. The implementations need some work. The accuracy is less than desired because gamma(2.0) /= 1. It's close but not quite right. rt-tests.lisp: o Basic tests of the gamma function. Accuracy is not as good as we would ike. qd-gamma.lisp: o New file for implementation of gamma function. diff --git a/oct.asd b/oct.asd index 89b56ee..e4011a9 100644 --- a/oct.asd +++ b/oct.asd @@ -65,6 +65,8 @@ :depends-on ("qd-methods" "qd-reader")) (:file "qd-theta" :depends-on ("qd-methods" "qd-reader")) + (:file "qd-gamma" + :depends-on ("qd-methods")) )) (defmethod perform ((op test-op) (c (eql (find-system :oct)))) diff --git a/qd-gamma.lisp b/qd-gamma.lisp new file mode 100644 index 0000000..bb9004a --- /dev/null +++ b/qd-gamma.lisp @@ -0,0 +1,170 @@ +;;;; -*- Mode: lisp -*- +;;;; +;;;; Copyright (c) 2011 Raymond Toy +;;;; Permission is hereby granted, free of charge, to any person +;;;; obtaining a copy of this software and associated documentation +;;;; files (the "Software"), to deal in the Software without +;;;; restriction, including without limitation the rights to use, +;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell +;;;; copies of the Software, and to permit persons to whom the +;;;; Software is furnished to do so, subject to the following +;;;; conditions: +;;;; +;;;; The above copyright notice and this permission notice shall be +;;;; included in all copies or substantial portions of the Software. +;;;; +;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +;;;; OTHER DEALINGS IN THE SOFTWARE. + +(in-package #:oct) + +(eval-when (:compile-toplevel :load-toplevel :execute) + (setf *readtable* *oct-readtable*)) + +;; For log-gamma we use the asymptotic formula +;; +;; log(gamma(z)) ~ (z - 1/2)*log(z) + log(2*%pi)/2 +;; + sum(bern(2*k)/(2*k)/(2*k-1)/z^(2k-1), k, 1, inf) +;; +;; = (z - 1/2)*log(z) + log(2*%pi)/2 +;; + 1/12/z*(1 - 1/30/z^2 + 1/105/z^4 + 1/140/z^6 + ... +;; + 174611/10450/z^18 + ...) +;; +;; For double-floats, let's stop the series at the power z^18. The +;; next term is 77683/483/z^20. This means that for |z| > 8.09438, +;; the series has double-float precision. +;; +;; For quad-doubles, let's stop the series at the power z^62. The +;; next term is about 6.364d37/z^64. So for |z| > 38.71, the series +;; has quad-double precision. +(defparameter *log-gamma-asymp-coef* + #(-1/30 1/105 -1/140 1/99 -691/30030 1/13 -3617/10200 43867/20349 + -174611/10450 77683/483 -236364091/125580 657931/25 -3392780147/7830 + 1723168255201/207669 -7709321041217/42160 151628697551/33 + -26315271553053477373/201514950 154210205991661/37 + -261082718496449122051/1758900 1520097643918070802691/259161 + -2530297234481911294093/9890 25932657025822267968607/2115 + -5609403368997817686249127547/8725080 19802288209643185928499101/539 + -61628132164268458257532691681/27030 29149963634884862421418123812691/190323 + -354198989901889536240773677094747/31900 + 2913228046513104891794716413587449/3363 + -1215233140483755572040304994079820246041491/16752085350 + 396793078518930920708162576045270521/61 + -106783830147866529886385444979142647942017/171360 + 133872729284212332186510857141084758385627191/2103465 + )) + +#+nil +(defun log-gamma-asymp-series (z nterms) + ;; Sum the asymptotic formula for n terms + ;; + ;; 1 + sum(c[k]/z^(2*k+2), k, 0, nterms) + (let ((z2 (* z z)) + (sum 1) + (term 1)) + (dotimes (k nterms) + (setf term (* term z2)) + (incf sum (/ (aref *log-gamma-asymp-coef* k) term))) + sum)) + +(defun log-gamma-asymp-series (z nterms) + (loop with y = (* z z) + for k from 1 to nterms + for x = 0 then + (setf x (/ (+ x (aref *log-gamma-asymp-coef* (- nterms k))) + y)) + finally (return (+ 1 x)))) + + +(defun log-gamma-asymp-principal (z nterms log2pi/2) + (+ (- (* (- z 1/2) + (log z)) + z) + log2pi/2)) + +(defun log-gamma-asymp (z nterms log2pi/2) + (+ (log-gamma-asymp-principal z nterms log2pi/2) + (* 1/12 (/ (log-gamma-asymp-series z nterms) z)))) + +(defun log2pi/2 (precision) + (ecase precision + (single-float + (coerce (/ (log (* 2 pi)) 2) 'single-float)) + (double-float + (coerce (/ (log (* 2 pi)) 2) 'double-float)) + (qd-real + (/ (log +2pi+) 2)))) + +(defun log-gamma-aux (z limit nterms) + (let ((precision (float-contagion z))) + (cond ((minusp (realpart z)) + ;; Use reflection formula if realpart(z) < 0 + ;; log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z)) + ;; Or + ;; log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z)) + (- (apply-contagion (log pi) precision) + (log (- z)) + (apply-contagion (log (sin (* pi z))) precision) + (log-gamma (- z)))) + (t + (let ((absz (abs z))) + (cond ((>= absz limit) + ;; Can use the asymptotic formula directly with 9 terms + (log-gamma-asymp z nterms (log2pi/2 precision))) + (t + ;; |z| is too small. Use the formula + ;; log(gamma(z)) = log(gamma(z+1)) - log(z) + (- (log-gamma (+ z 1)) + (log z))))))))) + +(defmethod log-gamma ((z cl:number)) + (log-gamma-aux z 9 9)) + +(defmethod log-gamma ((z qd-real)) + (log-gamma-aux z 26 26)) + +(defmethod log-gamma ((z qd-complex)) + (log-gamma-aux z 26 26)) + +(defun gamma-aux (z limit nterms) + (let ((precision (float-contagion z))) + (cond ((minusp (realpart z)) + ;; 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))) + (t + (let ((absz (abs z))) + (cond ((>= absz limit) + ;; Use log gamma directly: + ;; log(gamma(z)) = principal part + 1/12/z*(series part) + ;; so + ;; gamma(z) = exp(principal part)*exp(1/12/z*series) + (exp (log-gamma z)) + #+nil + (* (exp (log-gamma-asymp-principal z nterms + (log2pi/2 precision))) + (exp (* 1/12 (/ (log-gamma-asymp-series z nterms) z))))) + (t + ;; 1 <= |z| <= limit + ;; gamma(z) = gamma(z+1)/z + (/ (gamma-aux (+ 1 z) limit nterms) z)))))))) + +(defmethod gamma ((z cl:number)) + (gamma-aux z 9 9)) + +(defmethod gamma ((z qd-real)) + (gamma-aux z 39 32)) + +(defmethod gamma ((z qd-complex)) + (gamma-aux z 39 32)) + diff --git a/rt-tests.lisp b/rt-tests.lisp index c624f4a..59d798f 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1141,3 +1141,57 @@ when result append (list (list (list k m) result))) nil) + +(rt:deftest gamma.1.d + (let ((g (gamma 0.5d0)) + (true (sqrt pi))) + ;; This should give full accuracy but doesn't. + (check-accuracy 51 g true)) + nil) + +(rt:deftest gamma.1.q + (let ((g (gamma #q0.5)) + (true (sqrt +pi+))) + ;; This should give full accuracy but doesn't. + (check-accuracy 197 g true)) + nil) + +(rt:deftest gamma.2.d + (loop for k from 0 below 100 + 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) + when result + append (list (list (list k y) result))) + nil) + +(rt:deftest gamma.2.q + (loop for k from 0 below 100 + for y = (+ 1 (random #q100)) + for g = (abs (gamma (complex 0 y))) + for true = (sqrt (/ +pi+ y (sinh (* +pi+ y)))) + for result = (check-accuracy 196 g true) + when result + append (list (list (list k y) result))) + nil) + +(rt:deftest gamma.3.d + (loop for k from 0 below 100 + for y = (+ 1 (random 100d0)) + for g = (abs (gamma (complex 1/2 y))) + for true = (sqrt (/ pi (cosh (* pi y)))) + for result = (check-accuracy 44 g true) + when result + append (list (list (list k y) result))) + nil) + +(rt:deftest gamma.3.q + (loop for k from 0 below 100 + for y = (+ 1 (random #q100)) + for g = (abs (gamma (complex 1/2 y))) + for true = (sqrt (/ +pi+ (cosh (* +pi+ y)))) + for result = (check-accuracy 196 g true) + when result + append (list (list (list k y) result))) + nil) commit a666f3923b82be0eb435ddd3e9b20697097fafe0 Author: Raymond Toy <toy.raymond@gmail.com> Date: Wed Mar 16 19:41:33 2011 -0400 QCOMPLEX takes two required args now. Fixes issues like (complex 1/2 #q1), which was signaling an error. qd-class.lisp: o Update defgeneric for QCOMPLEX for two required args. qd-methods.lisp: o Update existing QCOMPLEX methods to take two args. o Add methods to QCOMPLEX to handle the missing cases. diff --git a/qd-class.lisp b/qd-class.lisp index 2134ef7..43cc2e9 100644 --- a/qd-class.lisp +++ b/qd-class.lisp @@ -203,8 +203,8 @@ (defgeneric qexpt (x y) (:documentation "X^Y")) -(defgeneric qcomplex (x &optional y) - (:documentation "Create a complex number with components X and Y. If Y not given, assume 0")) +(defgeneric qcomplex (x y) + (:documentation "Create a complex number with components X and Y.")) (defgeneric qinteger-decode-float (f) (:documentation "integer-decode-float")) diff --git a/qd-methods.lisp b/qd-methods.lisp index fa27ce2..9341311 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -507,13 +507,21 @@ underlying floating-point format" (return nil))) (return nil)))) -(defmethod qcomplex ((x real) &optional y) - (cl:complex x (if y y 0))) +(defmethod qcomplex ((x cl:real) (y cl:real)) + (cl:complex x y)) -(defmethod qcomplex ((x qd-real) &optional y) +(defmethod qcomplex ((x cl:real) (y qd-real)) + (qcomplex (make-qd x) y)) + +(defmethod qcomplex ((x qd-real) (y qd-real)) + (make-instance 'qd-complex + :real (qd-value x) + :imag (qd-value y))) + +(defmethod qcomplex ((x qd-real) (y cl:real)) (make-instance 'qd-complex :real (qd-value x) - :imag (if y (qd-value y) +qd-zero+))) + :imag (make-qd-d y))) (defun complex (x &optional (y 0)) (qcomplex x y)) commit a93aca9a6cfdb5b1305cb2570eea23e6f44b51c0 Author: Raymond Toy <toy.raymond@gmail.com> Date: Tue Mar 15 13:14:05 2011 -0400 Move float-contagion stuff to qd-methods. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 385f00f..475763b 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -29,59 +29,6 @@ (declaim (inline descending-transform ascending-transform)) -;; 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))))) - ;;; Jacobian elliptic functions (defun ascending-transform (u m) diff --git a/qd-methods.lisp b/qd-methods.lisp index 814672c..fa27ce2 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1048,4 +1048,58 @@ the same precision as the argument. The argument can be complex.")) (float pi (realpart z))) (defmethod float-pi ((z qd-complex)) - +pi+) \ No newline at end of file + +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: oct.asd | 2 + qd-class.lisp | 4 +- qd-elliptic.lisp | 53 ----------------- qd-gamma.lisp | 170 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ qd-methods.lisp | 72 +++++++++++++++++++++-- rt-tests.lisp | 54 +++++++++++++++++ 6 files changed, 295 insertions(+), 60 deletions(-) create mode 100644 qd-gamma.lisp hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats