source: qd-gamma.lisp @ d17228

Last change on this file since d17228 was d17228, checked in by Raymond Toy <toy.raymond@…>, 4 years ago

Add Fresnel integrals; fix issue in incomplete-gamma for large args.

o INCOMPLETE-GAMMA was returning bad values for large (complex)

arguments. Fix this by using incomplete gamma tail function since
the incomplete gamma function approaches gamma for large arguments.

o Implement Fresnel S and C functions.

  • Property mode set to 100644
File size: 11.5 KB
Line 
1;;;; -*- Mode: lisp -*-
2;;;;
3;;;; Copyright (c) 2011 Raymond Toy
4;;;; Permission is hereby granted, free of charge, to any person
5;;;; obtaining a copy of this software and associated documentation
6;;;; files (the "Software"), to deal in the Software without
7;;;; restriction, including without limitation the rights to use,
8;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
9;;;; copies of the Software, and to permit persons to whom the
10;;;; Software is furnished to do so, subject to the following
11;;;; conditions:
12;;;;
13;;;; The above copyright notice and this permission notice shall be
14;;;; included in all copies or substantial portions of the Software.
15;;;;
16;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
18;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
20;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
21;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
22;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
23;;;; OTHER DEALINGS IN THE SOFTWARE.
24
25(in-package #:oct)
26
27(eval-when (:compile-toplevel :load-toplevel :execute)
28  (setf *readtable* *oct-readtable*))
29
30;; For log-gamma we use the asymptotic formula
31;;
32;; log(gamma(z)) ~ (z - 1/2)*log(z) + log(2*%pi)/2
33;;                   + sum(bern(2*k)/(2*k)/(2*k-1)/z^(2k-1), k, 1, inf)
34;;
35;;               = (z - 1/2)*log(z) + log(2*%pi)/2
36;;                  + 1/12/z*(1 - 1/30/z^2 + 1/105/z^4 + 1/140/z^6 + ...
37;;                              + 174611/10450/z^18 + ...)
38;;
39;; For double-floats, let's stop the series at the power z^18.  The
40;; next term is 77683/483/z^20.  This means that for |z| > 8.09438,
41;; the series has double-float precision.
42;;
43;; For quad-doubles, let's stop the series at the power z^62.  The
44;; next term is about 6.364d37/z^64.  So for |z| > 38.71, the series
45;; has quad-double precision.
46(defparameter *log-gamma-asymp-coef*
47  #(-1/30 1/105 -1/140 1/99 -691/30030 1/13 -3617/10200 43867/20349
48    -174611/10450 77683/483 -236364091/125580 657931/25 -3392780147/7830
49    1723168255201/207669 -7709321041217/42160 151628697551/33
50    -26315271553053477373/201514950 154210205991661/37
51    -261082718496449122051/1758900 1520097643918070802691/259161
52    -2530297234481911294093/9890 25932657025822267968607/2115
53    -5609403368997817686249127547/8725080 19802288209643185928499101/539
54    -61628132164268458257532691681/27030 29149963634884862421418123812691/190323
55    -354198989901889536240773677094747/31900
56    2913228046513104891794716413587449/3363
57    -1215233140483755572040304994079820246041491/16752085350
58    396793078518930920708162576045270521/61
59    -106783830147866529886385444979142647942017/171360
60    133872729284212332186510857141084758385627191/2103465
61    ))
62
63#+nil
64(defun log-gamma-asymp-series (z nterms)
65  ;; Sum the asymptotic formula for n terms
66  ;;
67  ;; 1 + sum(c[k]/z^(2*k+2), k, 0, nterms)
68  (let ((z2 (* z z))
69        (sum 1)
70        (term 1))
71    (dotimes (k nterms)
72      (setf term (* term z2))
73      (incf sum (/ (aref *log-gamma-asymp-coef* k) term)))
74    sum))
75
76(defun log-gamma-asymp-series (z nterms)
77  (loop with y = (* z z)
78     for k from 1 to nterms
79     for x = 0 then
80       (setf x (/ (+ x (aref *log-gamma-asymp-coef* (- nterms k)))
81                  y))
82     finally (return (+ 1 x))))
83       
84
85(defun log-gamma-asymp-principal (z nterms log2pi/2)
86  (+ (- (* (- z 1/2)
87           (log z))
88        z)
89     log2pi/2))
90
91(defun log-gamma-asymp (z nterms log2pi/2)
92  (+ (log-gamma-asymp-principal z nterms log2pi/2)
93     (* 1/12 (/ (log-gamma-asymp-series z nterms) z))))
94
95(defun log2pi/2 (precision)
96  (ecase precision
97    (single-float
98     (coerce (/ (log (* 2 pi)) 2) 'single-float))
99    (double-float
100     (coerce (/ (log (* 2 pi)) 2) 'double-float))
101    (qd-real
102     (/ (log +2pi+) 2))))
103
104(defun log-gamma-aux (z limit nterms)
105  (let ((precision (float-contagion z)))
106    (cond ((minusp (realpart z))
107           ;; Use reflection formula if realpart(z) < 0
108           ;;   log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z))
109           ;; Or
110           ;;   log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z))
111           (- (apply-contagion (log pi) precision)
112              (log (- z))
113              (apply-contagion (log (sin (* pi z))) precision)
114              (log-gamma (- z))))
115          (t
116           (let ((absz (abs z)))
117             (cond ((>= absz limit)
118                    ;; Can use the asymptotic formula directly with 9 terms
119                    (log-gamma-asymp z nterms (log2pi/2 precision)))
120                   (t
121                    ;; |z| is too small.  Use the formula
122                    ;; log(gamma(z)) = log(gamma(z+1)) - log(z)
123                    (- (log-gamma (+ z 1))
124                       (log z)))))))))
125
126(defmethod log-gamma ((z cl:number))
127  (log-gamma-aux z 9 9))
128
129(defmethod log-gamma ((z qd-real))
130  (log-gamma-aux z 26 26))
131
132(defmethod log-gamma ((z qd-complex))
133  (log-gamma-aux z 26 26))
134
135(defun gamma-aux (z limit nterms)
136  (let ((precision (float-contagion z)))
137    (cond ((minusp (realpart z))
138           ;; Use reflection formula if realpart(z) < 0:
139           ;;  gamma(-z) = -pi*csc(pi*z)/gamma(z+1)
140           ;; or
141           ;;  gamma(z) = pi*csc(pi*z)/gamma(1-z)
142           (/ (float-pi z)
143              (sin (* (float-pi z) z))
144              (gamma-aux (- 1 z) limit nterms)))
145          (t
146           (let ((absz (abs z)))
147             (cond ((>= absz limit)
148                    ;; Use log gamma directly:
149                    ;;  log(gamma(z)) = principal part + 1/12/z*(series part)
150                    ;; so
151                    ;;  gamma(z) = exp(principal part)*exp(1/12/z*series)
152                    (exp (log-gamma z))
153                    #+nil
154                    (* (exp (log-gamma-asymp-principal z nterms
155                                                       (log2pi/2 precision)))
156                       (exp (* 1/12 (/ (log-gamma-asymp-series z nterms) z)))))
157                   (t
158                    ;; 1 <= |z| <= limit
159                    ;; gamma(z) = gamma(z+1)/z
160                    (/ (gamma-aux (+ 1 z) limit nterms) z))))))))
161                 
162(defmethod gamma ((z cl:number))
163  (gamma-aux z 9 9))
164
165(defmethod gamma ((z qd-real))
166  (gamma-aux z 39 32))
167
168(defmethod gamma ((z qd-complex))
169  (gamma-aux z 39 32))
170
171
172;; Lentz's algorithm for evaluating continued fractions.
173;;
174;; Let the continued fraction be:
175;;
176;;      a1    a2    a3
177;; b0 + ----  ----  ----
178;;      b1 +  b2 +  b3 +
179;;
180(defun lentz (bf af)
181  (flet ((value-or-tiny (v)
182           (if (zerop v)
183               (etypecase v
184                 ((or double-float cl:complex)
185                  least-positive-normalized-double-float)
186                 ((or qd-real qd-complex)
187                  (make-qd least-positive-normalized-double-float)))
188               v)))
189    (let* ((f (value-or-tiny (funcall bf 0)))
190           (c f)
191           (d 0)
192           (eps (epsilon f)))
193      (loop
194         for j from 1
195         for an = (funcall af j)
196         for bn = (funcall bf j)
197         do (progn
198              (setf d (value-or-tiny (+ bn (* an d))))
199              (setf c (value-or-tiny (+ bn (/ an c))))
200              (setf d (/ d))
201              (let ((delta (* c d)))
202                (setf f (* f delta))
203                (when (<= (abs (- delta 1)) eps)
204                  (return)))))
205      f)))
206
207;; Continued fraction for erf(b):
208;;
209;; z[n] = 1+2*n-2*z^2
210;; a[n] = 4*n*z^2
211;;
212;; This works ok, but has problems for z > 3 where sometimes the
213;; result is greater than 1.
214#+nil
215(defun erf (z)
216  (let* ((z2 (* z z))
217         (twoz2 (* 2 z2)))
218    (* (/ (* 2 z)
219          (sqrt (float-pi z)))
220       (exp (- z2))
221       (/ (lentz #'(lambda (n)
222                     (- (+ 1 (* 2 n))
223                        twoz2))
224                 #'(lambda (n)
225                     (* 4 n z2)))))))
226
227;; Tail of the incomplete gamma function:
228;; integrate(x^(a-1)*exp(-x), x, z, inf)
229;;
230;; The continued fraction, valid for all z except the negative real
231;; axis:
232;;
233;; b[n] = 1+2*n+z-a
234;; a[n] = n*(a-n)
235;;
236;; See http://functions.wolfram.com/06.06.10.0003.01
237(defun cf-incomplete-gamma-tail (a z)
238  (/ (* (expt z a)
239        (exp (- z)))
240     (let ((z-a (- z a)))
241       (lentz #'(lambda (n)
242                  (+ n n 1 z-a))
243              #'(lambda (n)
244                  (* n (- a n)))))))
245
246;; Incomplete gamma function:
247;; integrate(x^(a-1)*exp(-x), x, 0, z)
248;;
249;; The continued fraction, valid for all z except the negative real
250;; axis:
251;;
252;; b[n] = n - 1 + z + a
253;; a[n] = -z*(a + n)
254;;
255;; See http://functions.wolfram.com/06.06.10.0005.01.  We modified the
256;; continued fraction slightly and discarded the first quotient from
257;; the fraction.
258(defun cf-incomplete-gamma (a z)
259  (/ (* (expt z a)
260        (exp (- z)))
261     (let ((za1 (+ z a 1)))
262       (- a (/ (* a z)
263               (lentz #'(lambda (n)
264                          (+ n za1))
265                      #'(lambda (n)
266                          (- (* z (+ a n))))))))))
267
268;; Tail of the incomplete gamma function.
269(defun incomplete-gamma-tail (a z)
270  "Tail of the incomplete gamma function defined by:
271
272  integrate(t^(a-1)*exp(-t), t, z, inf)"
273  (let* ((prec (float-contagion a z))
274         (a (apply-contagion a prec))
275         (z (apply-contagion z prec)))
276    (if (and (realp a) (realp z))
277        ;; For real values, we split the result to compute either the
278        ;; tail directly or from gamma(a) - incomplete-gamma
279        (if (> z (- a 1))
280            (cf-incomplete-gamma-tail a z)
281            (- (gamma a) (cf-incomplete-gamma a z)))
282        (cf-incomplete-gamma-tail a z))))
283
284(defun incomplete-gamma (a z)
285  "Incomplete gamma function defined by:
286
287  integrate(t^(a-1)*exp(-t), t, 0, z)"
288  (let* ((prec (float-contagion a z))
289         (a (apply-contagion a prec))
290         (z (apply-contagion z prec)))
291    (if (and (realp a) (realp z))
292        (if (< z (- a 1))
293            (cf-incomplete-gamma a z)
294            (- (gamma a) (cf-incomplete-gamma-tail a z)))
295        (if (< (abs z) (abs a))
296            (cf-incomplete-gamma a z)
297            (- (gamma a) (cf-incomplete-gamma-tail a z))))))
298
299(defun erf (z)
300  "Error function:
301
302    erf(z) = 2/sqrt(%pi)*sum((-1)^k*z^(2*k+1)/k!/(2*k+1), k, 0, inf)
303
304  For real z, this is equivalent to
305
306    erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z) for real z."
307  ;;
308  ;; Erf is an odd function: erf(-z) = -erf(z)
309  (if (minusp (realpart z))
310      (- (erf (- z)))
311      (/ (incomplete-gamma 1/2 (* z z))
312         (sqrt (float-pi z)))))
313
314(defun erfc (z)
315  "Complementary error function:
316
317    erfc(z) = 1 - erf(z)"
318  ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
319  ;; near 1.  Wolfram says
320  ;;
321  ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
322  ;;
323  ;; For real(z) > 0, sqrt(z^2)/z is 1 so
324  ;;
325  ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
326  ;;         = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
327  ;;
328  ;; For real(z) < 0, sqrt(z^2)/z is -1 so
329  ;;
330  ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
331  ;;         = 1 + 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
332  (if (>= (realpart z) 0)
333      (/ (incomplete-gamma-tail 1/2 (* z z))
334         (sqrt (float-pi z)))
335      (+ 1
336         (/ (incomplete-gamma 1/2 (* z z))
337            (sqrt (float-pi z))))))
338
339(defun exp-integral-e (v z)
340  "Exponential integral E:
341
342   E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)"
343  ;; Wolfram gives E(v,z) = z^(v-1)*gamma_incomplete_tail(1-v,z)
344  (* (expt z (- v 1))
345     (incomplete-gamma-tail (- 1 v) z)))
346
347(defun fresnel-s (z)
348  "Fresnel S:
349
350   S(z) = integrate(sin(%pi*t^2/2), t, 0, z) "
351  (let ((sqrt-pi (sqrt (float-pi z))))
352    (flet ((fs (z)
353             ;; Wolfram gives
354             ;;
355             ;;  S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z))
356             ;;
357             ;; where c = sqrt(%pi)/2*(1+%i).
358             (* #c(1/4 1/4)
359                (- (erf (* #c(1/2 1/2) sqrt-pi z))
360                   (* #c(0 1)
361                      (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
362      (if (realp z)
363          ;; FresnelS is real for a real argument. And it is odd.
364          (if (minusp z)
365              (- (realpart (fs (- z))))
366              (realpart (fs z)))
367          (fs z)))))
368
369(defun fresnel-c (z)
370  "Fresnel C:
371
372   C(z) = integrate(cos(%pi*t^2/2), t, 0, z) "
373  (let ((sqrt-pi (sqrt (float-pi z))))
374    (flet ((fs (z)
375             ;; Wolfram gives
376             ;;
377             ;;  C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z))
378             ;;
379             ;; where c = sqrt(%pi)/2*(1+%i).
380             (* #c(1/4 -1/4)
381                (+ (erf (* #c(1/2 1/2) sqrt-pi z))
382                   (* #c(0 1)
383                      (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
384      (if (realp z)
385          ;; FresnelS is real for a real argument. And it is odd.
386          (if (minusp z)
387              (- (realpart (fs (- z))))
388              (realpart (fs z)))
389          (fs z)))))
390     
Note: See TracBrowser for help on using the repository browser.