Commit 481fd7e3 authored by Raymond Toy's avatar Raymond Toy
Browse files

Rename some incomplete gamma functions

Rename the functions to make it clearer what's being computed since
we're using -upper and -lower for the different parts of the
incomplete gamma function.  So:

- `cf-incomplete-gamma-tail` -> `cf-incomplete-gamma-upper`
- `cf-incomplete-gamma` -> `cf-incomplete-gamma-lower`
- `s-incomplete-gamma` -> `s-incomplete-gamma-lower`
parent a26ad0bb
Pipeline #2434 passed with stage
in 4 minutes and 56 seconds
......@@ -385,10 +385,10 @@
;; a[n] = n*(a-n)
;;
;; See http://functions.wolfram.com/06.06.10.0003.01
(defun cf-incomplete-gamma-tail (a z)
(defun cf-incomplete-gamma-upper (a z)
(when (and (zerop (imagpart z)) (minusp (realpart z)))
(error 'domain-error
:function-name 'cf-incomplete-gamma-tail
:function-name 'cf-incomplete-gamma-upper
:format-arguments (list 'z z)
:format-control "Argument ~S should not be on the negative real axis: ~S"))
(/ (handler-case (* (expt z a)
......@@ -416,7 +416,7 @@
;; continued fraction slightly and discarded the first quotient from
;; the fraction.
#+nil
(defun cf-incomplete-gamma (a z)
(defun cf-incomplete-gamma-lower (a z)
(/ (handler-case (* (expt z a)
(exp (- z)))
(arithmetic-error ()
......@@ -445,7 +445,7 @@
;; Some experiments indicate that this converges faster than the above
;; and is actually quite a bit more accurate, expecially near the
;; negative real axis.
(defun cf-incomplete-gamma (a z)
(defun cf-incomplete-gamma-lower (a z)
(/ (handler-case (* (expt z a)
(exp (- z)))
(arithmetic-error ()
......@@ -464,7 +464,7 @@
;; |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)
(defun s-incomplete-gamma-lower (a z)
(let ((-z (- z))
(eps (epsilon z)))
(loop for k from 0
......@@ -501,16 +501,16 @@
;; 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)))
(cf-incomplete-gamma-upper a z)
(- (gamma a) (cf-incomplete-gamma-lower 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 (<= (abs (phase z)) 3.1)
(cf-incomplete-gamma-tail a z)
(- (gamma a) (cf-incomplete-gamma a z))))))))
(cf-incomplete-gamma-upper a z)
(- (gamma a) (cf-incomplete-gamma-lower a z))))))))
(defun incomplete-gamma-lower (a z &key (normalized-p nil))
"Incomplete gamma function defined by:
......@@ -524,17 +524,17 @@
(incomplete-gamma-lower-normalized a z)
(with-floating-point-contagion (a z)
(if (and (< (abs a) 1) (< (abs z) 1))
(s-incomplete-gamma a z)
(s-incomplete-gamma-lower 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)))
(cf-incomplete-gamma-lower a z)
(- (gamma a) (cf-incomplete-gamma-upper 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))))))))
(cf-incomplete-gamma-lower a z)
(- (gamma a) (cf-incomplete-gamma-upper a z))))))))
;; The continued fraction for the normalized incomplete gamma upper
;; function, Q(a, z), valid for all z except for the negative real
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment