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 | |
---|