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 | (let ((p (float-pi z))) |
---|
112 | (- (log p) |
---|
113 | (log (- z)) |
---|
114 | (log (sin (* p z))) |
---|
115 | (log-gamma (- z))))) |
---|
116 | (t |
---|
117 | (let ((absz (abs z))) |
---|
118 | (cond ((>= absz limit) |
---|
119 | ;; Can use the asymptotic formula directly with 9 terms |
---|
120 | (log-gamma-asymp z nterms (log2pi/2 precision))) |
---|
121 | (t |
---|
122 | ;; |z| is too small. Use the formula |
---|
123 | ;; log(gamma(z)) = log(gamma(z+1)) - log(z) |
---|
124 | (- (log-gamma (+ z 1)) |
---|
125 | (log z))))))))) |
---|
126 | |
---|
127 | (defmethod log-gamma ((z cl:number)) |
---|
128 | (log-gamma-aux z 9 9)) |
---|
129 | |
---|
130 | (defmethod log-gamma ((z qd-real)) |
---|
131 | (log-gamma-aux z 26 26)) |
---|
132 | |
---|
133 | (defmethod log-gamma ((z qd-complex)) |
---|
134 | (log-gamma-aux z 26 26)) |
---|
135 | |
---|
136 | (defun gamma-aux (z limit nterms) |
---|
137 | (let ((precision (float-contagion z))) |
---|
138 | (cond ((<= (realpart z) 0) |
---|
139 | ;; Use reflection formula if realpart(z) < 0: |
---|
140 | ;; gamma(-z) = -pi*csc(pi*z)/gamma(z+1) |
---|
141 | ;; or |
---|
142 | ;; gamma(z) = pi*csc(pi*z)/gamma(1-z) |
---|
143 | (if (and (realp z) |
---|
144 | (= (truncate z) z)) |
---|
145 | ;; Gamma of a negative integer is infinity. Signal an error |
---|
146 | (error "Gamma of non-positive integer ~S" z) |
---|
147 | (/ (float-pi z) |
---|
148 | (sin (* (float-pi z) z)) |
---|
149 | (gamma-aux (- 1 z) limit nterms)))) |
---|
150 | ((and (zerop (imagpart z)) |
---|
151 | (= z (truncate z))) |
---|
152 | ;; We have gamma(n) where an integer value n and is small |
---|
153 | ;; enough. In this case, just compute the product |
---|
154 | ;; directly. We do this because our current implementation |
---|
155 | ;; has some round-off for these values, and that's annoying |
---|
156 | ;; and unexpected. |
---|
157 | (let ((n (truncate z))) |
---|
158 | (loop |
---|
159 | for prod = (apply-contagion 1 precision) then (* prod k) |
---|
160 | for k from 2 below n |
---|
161 | finally (return (apply-contagion prod precision))))) |
---|
162 | (t |
---|
163 | (let ((absz (abs z))) |
---|
164 | (cond ((>= absz limit) |
---|
165 | ;; Use log gamma directly: |
---|
166 | ;; log(gamma(z)) = principal part + 1/12/z*(series part) |
---|
167 | ;; so |
---|
168 | ;; gamma(z) = exp(principal part)*exp(1/12/z*series) |
---|
169 | (exp (log-gamma z)) |
---|
170 | #+nil |
---|
171 | (* (exp (log-gamma-asymp-principal z nterms |
---|
172 | (log2pi/2 precision))) |
---|
173 | (exp (* 1/12 (/ (log-gamma-asymp-series z nterms) z))))) |
---|
174 | (t |
---|
175 | ;; 1 <= |z| <= limit |
---|
176 | ;; gamma(z) = gamma(z+1)/z |
---|
177 | (/ (gamma-aux (+ 1 z) limit nterms) z)))))))) |
---|
178 | |
---|
179 | (defmethod gamma ((z cl:number)) |
---|
180 | (gamma-aux z 9 9)) |
---|
181 | |
---|
182 | (defmethod gamma ((z qd-real)) |
---|
183 | (gamma-aux z 39 32)) |
---|
184 | |
---|
185 | (defmethod gamma ((z qd-complex)) |
---|
186 | (gamma-aux z 39 32)) |
---|
187 | |
---|
188 | ;; Lentz's algorithm for evaluating continued fractions. |
---|
189 | ;; |
---|
190 | ;; Let the continued fraction be: |
---|
191 | ;; |
---|
192 | ;; a1 a2 a3 |
---|
193 | ;; b0 + ---- ---- ---- |
---|
194 | ;; b1 + b2 + b3 + |
---|
195 | ;; |
---|
196 | |
---|
197 | (defvar *debug-cf-eval* |
---|
198 | nil |
---|
199 | "When true, enable some debugging prints when evaluating a |
---|
200 | continued fraction.") |
---|
201 | |
---|
202 | ;; Max number of iterations allowed when evaluating the continued |
---|
203 | ;; fraction. When this is reached, we assume that the continued |
---|
204 | ;; fraction did not converge. |
---|
205 | (defvar *max-cf-iterations* |
---|
206 | 10000 |
---|
207 | "Max number of iterations allowed when evaluating the continued |
---|
208 | fraction. When this is reached, we assume that the continued |
---|
209 | fraction did not converge.") |
---|
210 | |
---|
211 | (defun lentz (bf af) |
---|
212 | (let ((tiny-value-count 0)) |
---|
213 | (flet ((value-or-tiny (v) |
---|
214 | (if (zerop v) |
---|
215 | (progn |
---|
216 | (incf tiny-value-count) |
---|
217 | (etypecase v |
---|
218 | ((or double-float cl:complex) |
---|
219 | (sqrt least-positive-normalized-double-float)) |
---|
220 | ((or qd-real qd-complex) |
---|
221 | (make-qd (sqrt least-positive-normalized-double-float))))) |
---|
222 | v))) |
---|
223 | (let* ((f (value-or-tiny (funcall bf 0))) |
---|
224 | (c f) |
---|
225 | (d 0) |
---|
226 | (eps (epsilon f))) |
---|
227 | (loop |
---|
228 | for j from 1 upto *max-cf-iterations* |
---|
229 | for an = (funcall af j) |
---|
230 | for bn = (funcall bf j) |
---|
231 | do (progn |
---|
232 | (setf d (value-or-tiny (+ bn (* an d)))) |
---|
233 | (setf c (value-or-tiny (+ bn (/ an c)))) |
---|
234 | (when *debug-cf-eval* |
---|
235 | (format t "~&j = ~d~%" j) |
---|
236 | (format t " an = ~s~%" an) |
---|
237 | (format t " bn = ~s~%" bn) |
---|
238 | (format t " c = ~s~%" c) |
---|
239 | (format t " d = ~s~%" d)) |
---|
240 | (let ((delta (/ c d))) |
---|
241 | (setf d (/ d)) |
---|
242 | (setf f (* f delta)) |
---|
243 | (when *debug-cf-eval* |
---|
244 | (format t " dl= ~S (|dl - 1| = ~S)~%" delta (abs (1- delta))) |
---|
245 | (format t " f = ~S~%" f)) |
---|
246 | (when (<= (abs (- delta 1)) eps) |
---|
247 | (return-from lentz (values f j tiny-value-count))))) |
---|
248 | finally |
---|
249 | (error 'simple-error |
---|
250 | :format-control "~<Continued fraction failed to converge after ~D iterations.~% Delta = ~S~>" |
---|
251 | :format-arguments (list *max-cf-iterations* (/ c d)))))))) |
---|
252 | |
---|
253 | ;; Continued fraction for erf(z): |
---|
254 | ;; |
---|
255 | ;; erf(z) = 2*z/sqrt(pi)*exp(-z^2)/K |
---|
256 | ;; |
---|
257 | ;; where K is the continued fraction with |
---|
258 | ;; |
---|
259 | ;; b[n] = 1+2*n-2*z^2 |
---|
260 | ;; a[n] = 4*n*z^2 |
---|
261 | ;; |
---|
262 | ;; This works ok, but has problems for z > 3 where sometimes the |
---|
263 | ;; result is greater than 1 and for larger values, negative numbers |
---|
264 | ;; are returned. |
---|
265 | #+nil |
---|
266 | (defun cf-erf (z) |
---|
267 | (let* ((z2 (* z z)) |
---|
268 | (twoz2 (* 2 z2))) |
---|
269 | (* (/ (* 2 z) |
---|
270 | (sqrt (float-pi z))) |
---|
271 | (exp (- z2)) |
---|
272 | (/ (lentz #'(lambda (n) |
---|
273 | (- (+ 1 (* 2 n)) |
---|
274 | twoz2)) |
---|
275 | #'(lambda (n) |
---|
276 | (* 4 n z2))))))) |
---|
277 | |
---|
278 | ;; Tail of the incomplete gamma function: |
---|
279 | ;; integrate(x^(a-1)*exp(-x), x, z, inf) |
---|
280 | ;; |
---|
281 | ;; The continued fraction, valid for all z except the negative real |
---|
282 | ;; axis: |
---|
283 | ;; |
---|
284 | ;; b[n] = 1+2*n+z-a |
---|
285 | ;; a[n] = n*(a-n) |
---|
286 | ;; |
---|
287 | ;; See http://functions.wolfram.com/06.06.10.0003.01 |
---|
288 | (defun cf-incomplete-gamma-tail (a z) |
---|
289 | (when (and (zerop (imagpart z)) (minusp (realpart z))) |
---|
290 | (error 'domain-error |
---|
291 | :function-name 'cf-incomplete-gamma-tail |
---|
292 | :format-arguments (list 'z z) |
---|
293 | :format-control "Argument ~S should not be on the negative real axis: ~S")) |
---|
294 | (/ (handler-case (* (expt z a) |
---|
295 | (exp (- z))) |
---|
296 | (arithmetic-error () |
---|
297 | ;; z^a*exp(-z) can overflow prematurely. In this case, use |
---|
298 | ;; the equivalent exp(a*log(z)-z). We don't use this latter |
---|
299 | ;; form because it has more roundoff error than the former. |
---|
300 | (exp (- (* a (log z)) z)))) |
---|
301 | (let ((z-a (- z a))) |
---|
302 | (lentz #'(lambda (n) |
---|
303 | (+ n n 1 z-a)) |
---|
304 | #'(lambda (n) |
---|
305 | (* n (- a n))))))) |
---|
306 | |
---|
307 | ;; Incomplete gamma function: |
---|
308 | ;; integrate(x^(a-1)*exp(-x), x, 0, z) |
---|
309 | ;; |
---|
310 | ;; The continued fraction, valid for all z: |
---|
311 | ;; |
---|
312 | ;; b[n] = n - 1 + z + a |
---|
313 | ;; a[n] = -z*(a + n) |
---|
314 | ;; |
---|
315 | ;; See http://functions.wolfram.com/06.06.10.0007.01. We modified the |
---|
316 | ;; continued fraction slightly and discarded the first quotient from |
---|
317 | ;; the fraction. |
---|
318 | #+nil |
---|
319 | (defun cf-incomplete-gamma (a z) |
---|
320 | (/ (handler-case (* (expt z a) |
---|
321 | (exp (- z))) |
---|
322 | (arithmetic-error () |
---|
323 | ;; z^a*exp(-z) can overflow prematurely. In this case, use |
---|
324 | ;; the equivalent exp(a*log(z)-z). We don't use this latter |
---|
325 | ;; form because it has more roundoff error than the former. |
---|
326 | (exp (- (* a (log z)) z)))) |
---|
327 | (let ((za1 (+ z a 1))) |
---|
328 | (- a (/ (* a z) |
---|
329 | (lentz #'(lambda (n) |
---|
330 | (+ n za1)) |
---|
331 | #'(lambda (n) |
---|
332 | (- (* z (+ a n)))))))))) |
---|
333 | |
---|
334 | ;; Incomplete gamma function: |
---|
335 | ;; integrate(x^(a-1)*exp(-x), x, 0, z) |
---|
336 | ;; |
---|
337 | ;; The continued fraction, valid for all z: |
---|
338 | ;; |
---|
339 | ;; b[n] = a + n |
---|
340 | ;; a[n] = -(a+n/2)*z if n odd |
---|
341 | ;; n/2*z if n even |
---|
342 | ;; |
---|
343 | ;; See http://functions.wolfram.com/06.06.10.0009.01. |
---|
344 | ;; |
---|
345 | ;; Some experiments indicate that this converges faster than the above |
---|
346 | ;; and is actually quite a bit more accurate, expecially near the |
---|
347 | ;; negative real axis. |
---|
348 | (defun cf-incomplete-gamma (a z) |
---|
349 | (/ (handler-case (* (expt z a) |
---|
350 | (exp (- z))) |
---|
351 | (arithmetic-error () |
---|
352 | ;; z^a*exp(-z) can overflow prematurely. In this case, use |
---|
353 | ;; the equivalent exp(a*log(z)-z). We don't use this latter |
---|
354 | ;; form because it has more roundoff error than the former. |
---|
355 | (exp (- (* a (log z)) z)))) |
---|
356 | (lentz #'(lambda (n) |
---|
357 | (+ n a)) |
---|
358 | #'(lambda (n) |
---|
359 | (if (evenp n) |
---|
360 | (* (ash n -1) z) |
---|
361 | (- (* (+ a (ash n -1)) z))))))) |
---|
362 | |
---|
363 | ;; Series expansion for incomplete gamma. Intended for |a|<1 and |
---|
364 | ;; |z|<1. The series is |
---|
365 | ;; |
---|
366 | ;; g(a,z) = z^a * sum((-z)^k/k!/(a+k), k, 0, inf) |
---|
367 | (defun s-incomplete-gamma (a z) |
---|
368 | (let ((-z (- z)) |
---|
369 | (eps (epsilon z))) |
---|
370 | (loop for k from 0 |
---|
371 | for term = 1 then (* term (/ -z k)) |
---|
372 | for sum = (/ a) then (+ sum (/ term (+ a k))) |
---|
373 | when (< (abs term) (* (abs sum) eps)) |
---|
374 | return (* sum (expt z a))))) |
---|
375 | |
---|
376 | ;; Tail of the incomplete gamma function. |
---|
377 | (defun incomplete-gamma-tail (a z) |
---|
378 | "Tail of the incomplete gamma function defined by: |
---|
379 | |
---|
380 | integrate(t^(a-1)*exp(-t), t, z, inf)" |
---|
381 | (with-floating-point-contagion (a z) |
---|
382 | (if (and (realp a) (<= a 0)) |
---|
383 | ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z) |
---|
384 | (* (expt z a) |
---|
385 | (exp-integral-e (- 1 a) z)) |
---|
386 | (if (and (zerop (imagpart a)) |
---|
387 | (zerop (imagpart z))) |
---|
388 | ;; For real values, we split the result to compute either the |
---|
389 | ;; tail directly from the continued fraction or from gamma(a) |
---|
390 | ;; - incomplete-gamma. The continued fraction doesn't |
---|
391 | ;; converge on the negative real axis, so we can't use that |
---|
392 | ;; there. And accuracy appears to be better if z is "small". |
---|
393 | ;; We take this to mean |z| < |a-1|. Note that |a-1| is the |
---|
394 | ;; peak of the integrand. |
---|
395 | (if (and (> (abs z) (abs (- a 1))) |
---|
396 | (not (minusp (realpart z)))) |
---|
397 | (cf-incomplete-gamma-tail a z) |
---|
398 | (- (gamma a) (cf-incomplete-gamma a z))) |
---|
399 | ;; If the argument is close enough to the negative real axis, |
---|
400 | ;; the continued fraction for the tail is not very accurate. |
---|
401 | ;; Use the incomplete gamma function to evaluate in this |
---|
402 | ;; region. (Arbitrarily selected the region to be a sector. |
---|
403 | ;; But what is the correct size of this sector?) |
---|
404 | (if (<= (abs (phase z)) 3.1) |
---|
405 | (cf-incomplete-gamma-tail a z) |
---|
406 | (- (gamma a) (cf-incomplete-gamma a z))))))) |
---|
407 | |
---|
408 | (defun incomplete-gamma (a z) |
---|
409 | "Incomplete gamma function defined by: |
---|
410 | |
---|
411 | integrate(t^(a-1)*exp(-t), t, 0, z)" |
---|
412 | (with-floating-point-contagion (a z) |
---|
413 | (if (and (< (abs a) 1) (< (abs z) 1)) |
---|
414 | (s-incomplete-gamma a z) |
---|
415 | (if (and (realp a) (realp z)) |
---|
416 | (if (< z (- a 1)) |
---|
417 | (cf-incomplete-gamma a z) |
---|
418 | (- (gamma a) (cf-incomplete-gamma-tail a z))) |
---|
419 | ;; The continued fraction doesn't converge very fast if a |
---|
420 | ;; and z are small. In this case, use the series |
---|
421 | ;; expansion instead, which converges quite rapidly. |
---|
422 | (if (< (abs z) (abs a)) |
---|
423 | (cf-incomplete-gamma a z) |
---|
424 | (- (gamma a) (cf-incomplete-gamma-tail a z))))))) |
---|
425 | |
---|
426 | (defun erf (z) |
---|
427 | "Error function: |
---|
428 | |
---|
429 | erf(z) = 2/sqrt(%pi)*sum((-1)^k*z^(2*k+1)/k!/(2*k+1), k, 0, inf) |
---|
430 | |
---|
431 | For real z, this is equivalent to |
---|
432 | |
---|
433 | erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z) for real z." |
---|
434 | ;; |
---|
435 | ;; Erf is an odd function: erf(-z) = -erf(z) |
---|
436 | (if (minusp (realpart z)) |
---|
437 | (- (erf (- z))) |
---|
438 | (/ (incomplete-gamma 1/2 (* z z)) |
---|
439 | (sqrt (float-pi z))))) |
---|
440 | |
---|
441 | (defun erfc (z) |
---|
442 | "Complementary error function: |
---|
443 | |
---|
444 | erfc(z) = 1 - erf(z)" |
---|
445 | ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is |
---|
446 | ;; near 1. Wolfram says |
---|
447 | ;; |
---|
448 | ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2)) |
---|
449 | ;; |
---|
450 | ;; For real(z) > 0, sqrt(z^2)/z is 1 so |
---|
451 | ;; |
---|
452 | ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)) |
---|
453 | ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2) |
---|
454 | ;; |
---|
455 | ;; For real(z) < 0, sqrt(z^2)/z is -1 so |
---|
456 | ;; |
---|
457 | ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)) |
---|
458 | ;; = 1 + 1/sqrt(pi)*gamma_incomplete(1/2,z^2) |
---|
459 | (if (>= (realpart z) 0) |
---|
460 | (/ (incomplete-gamma-tail 1/2 (* z z)) |
---|
461 | (sqrt (float-pi z))) |
---|
462 | (+ 1 |
---|
463 | (/ (incomplete-gamma 1/2 (* z z)) |
---|
464 | (sqrt (float-pi z)))))) |
---|
465 | |
---|
466 | (defun cf-exp-integral-e (v z) |
---|
467 | ;; We use the continued fraction |
---|
468 | ;; |
---|
469 | ;; E(v,z) = exp(-z)/cf(z) |
---|
470 | ;; |
---|
471 | ;; where the continued fraction cf(z) is |
---|
472 | ;; |
---|
473 | ;; a[k] = -k*(k+v-1) |
---|
474 | ;; b[k] = v + 2*k + z |
---|
475 | ;; |
---|
476 | ;; for k = 1, inf |
---|
477 | (let ((z+v (+ z v))) |
---|
478 | (/ (exp (- z)) |
---|
479 | (lentz #'(lambda (k) |
---|
480 | (+ z+v (* 2 k))) |
---|
481 | #'(lambda (k) |
---|
482 | (* (- k) |
---|
483 | (+ k v -1))))))) |
---|
484 | |
---|
485 | |
---|
486 | ;; For v not an integer: |
---|
487 | ;; |
---|
488 | ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf) |
---|
489 | ;; |
---|
490 | ;; For v an integer: |
---|
491 | ;; |
---|
492 | ;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z)) |
---|
493 | ;; - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1) |
---|
494 | ;; |
---|
495 | (defun s-exp-integral-e (v z) |
---|
496 | ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf) |
---|
497 | (let ((-z (- z)) |
---|
498 | (-v (- v)) |
---|
499 | (eps (epsilon z))) |
---|
500 | (if (and (realp v) |
---|
501 | (= v (ftruncate v))) |
---|
502 | ;; v is an integer |
---|
503 | (let* ((n (truncate v)) |
---|
504 | (n-1 (1- n))) |
---|
505 | (- (* (/ (expt -z n-1) |
---|
506 | (gamma v)) |
---|
507 | (- (psi v) (log z))) |
---|
508 | (loop for k from 0 |
---|
509 | for term = 1 then (* term (/ -z k)) |
---|
510 | for sum = (if (zerop n-1) |
---|
511 | 0 |
---|
512 | (/ (- 1 v))) |
---|
513 | then (+ sum (let ((denom (- k n-1))) |
---|
514 | (if (zerop denom) |
---|
515 | 0 |
---|
516 | (/ term denom)))) |
---|
517 | when (< (abs term) (* (abs sum) eps)) |
---|
518 | return sum))) |
---|
519 | (loop for k from 0 |
---|
520 | for term = 1 then (* term (/ -z k)) |
---|
521 | for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v))) |
---|
522 | when (< (abs term) (* (abs sum) eps)) |
---|
523 | return (- (* (gamma (- 1 v)) (expt z (- v 1))) |
---|
524 | sum))))) |
---|
525 | |
---|
526 | (defun exp-integral-e (v z) |
---|
527 | "Exponential integral E: |
---|
528 | |
---|
529 | E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)" |
---|
530 | ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); |
---|
531 | ;; |
---|
532 | ;; for |arg(z)| < pi. |
---|
533 | ;; |
---|
534 | ;; |
---|
535 | (with-floating-point-contagion (v z) |
---|
536 | (cond ((and (realp v) (minusp v)) |
---|
537 | ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z) |
---|
538 | (let ((-v (- v))) |
---|
539 | (* (expt z (- v 1)) |
---|
540 | (incomplete-gamma-tail (+ -v 1) z)))) |
---|
541 | ((or (< (abs z) 1) (>= (abs (phase z)) 3.1)) |
---|
542 | ;; Use series for small z or if z is near the negative real |
---|
543 | ;; axis because the continued fraction does not converge on |
---|
544 | ;; the negative axis and converges slowly near the negative |
---|
545 | ;; axis. |
---|
546 | (s-exp-integral-e v z)) |
---|
547 | (t |
---|
548 | ;; Use continued fraction for everything else. |
---|
549 | (cf-exp-integral-e v z))))) |
---|
550 | |
---|
551 | ;; Series for Fresnel S |
---|
552 | ;; |
---|
553 | ;; S(z) = z^3*sum((%pi/2)^(2*k+1)(-z^4)^k/(2*k+1)!/(4*k+3), k, 0, inf) |
---|
554 | ;; |
---|
555 | ;; Compute as |
---|
556 | ;; |
---|
557 | ;; S(z) = z^3*sum(a(k)/(4*k+3), k, 0, inf) |
---|
558 | ;; |
---|
559 | ;; where |
---|
560 | ;; |
---|
561 | ;; a(k+1) = -a(k) * (%pi/2)^2 * z^4 / (2*k+2) / (2*k+3) |
---|
562 | ;; |
---|
563 | ;; a(0) = %pi/2. |
---|
564 | (defun fresnel-s-series (z) |
---|
565 | (let* ((pi/2 (* 1/2 (float-pi z))) |
---|
566 | (factor (- (* (expt z 4) pi/2 pi/2))) |
---|
567 | (eps (epsilon z)) |
---|
568 | (sum 0) |
---|
569 | (term pi/2)) |
---|
570 | (loop for k2 from 0 by 2 |
---|
571 | until (< (abs term) (* eps (abs sum))) |
---|
572 | do |
---|
573 | (incf sum (/ term (+ 3 k2 k2))) |
---|
574 | (setf term (/ (* term factor) |
---|
575 | (* (+ k2 2) |
---|
576 | (+ k2 3))))) |
---|
577 | (* sum (expt z 3)))) |
---|
578 | |
---|
579 | (defun fresnel-s (z) |
---|
580 | "Fresnel S: |
---|
581 | |
---|
582 | S(z) = integrate(sin(%pi*t^2/2), t, 0, z) " |
---|
583 | (let ((prec (float-contagion z)) |
---|
584 | (sqrt-pi (sqrt (float-pi z)))) |
---|
585 | (flet ((fs (z) |
---|
586 | ;; Wolfram gives |
---|
587 | ;; |
---|
588 | ;; S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z)) |
---|
589 | ;; |
---|
590 | ;; where c = sqrt(%pi)/2*(1+%i). |
---|
591 | ;; |
---|
592 | ;; But for large z, we should use erfc. Then |
---|
593 | ;; S(z) = 1/2 - (1+%i)/4*(erfc(c*z) - %i*erfc(conjugate(c)*z)) |
---|
594 | (if (and t (> (abs z) 2)) |
---|
595 | (- 1/2 |
---|
596 | (* #c(1/4 1/4) |
---|
597 | (- (erfc (* #c(1/2 1/2) sqrt-pi z)) |
---|
598 | (* #c(0 1) |
---|
599 | (erfc (* #c(1/2 -1/2) sqrt-pi z)))))) |
---|
600 | (* #c(1/4 1/4) |
---|
601 | (- (erf (* #c(1/2 1/2) sqrt-pi z)) |
---|
602 | (* #c(0 1) |
---|
603 | (erf (* #c(1/2 -1/2) sqrt-pi z))))))) |
---|
604 | (rfs (z) |
---|
605 | ;; When z is real, recall that erf(conjugate(z)) = |
---|
606 | ;; conjugate(erf(z)). Then |
---|
607 | ;; |
---|
608 | ;; S(z) = 1/2*(realpart(erf(c*z)) - imagpart(erf(c*z))) |
---|
609 | ;; |
---|
610 | ;; But for large z, we should use erfc. Then |
---|
611 | ;; |
---|
612 | ;; S(z) = 1/2 - 1/2*(realpart(erfc(c*z)) - imagpart(erf(c*z))) |
---|
613 | (if (> (abs z) 2) |
---|
614 | (let ((s (erfc (* #c(1/2 1/2) sqrt-pi z)))) |
---|
615 | (- 1/2 |
---|
616 | (* 1/2 (- (realpart s) (imagpart s))))) |
---|
617 | (let ((s (erf (* #c(1/2 1/2) sqrt-pi z)))) |
---|
618 | (* 1/2 (- (realpart s) (imagpart s))))))) |
---|
619 | ;; For small z, the erf terms above suffer from subtractive |
---|
620 | ;; cancellation. So use the series in this case. Some simple |
---|
621 | ;; tests were done to determine that for double-floats we want |
---|
622 | ;; to use the series for z < 1 to give max accuracy. For |
---|
623 | ;; qd-real, the above formula is good enough for z > 1d-5. |
---|
624 | (if (< (abs z) (ecase prec |
---|
625 | (single-float 1.5f0) |
---|
626 | (double-float 1d0) |
---|
627 | (qd-real #q1))) |
---|
628 | (fresnel-s-series z) |
---|
629 | (if (realp z) |
---|
630 | ;; FresnelS is real for a real argument. And it is odd. |
---|
631 | (if (minusp z) |
---|
632 | (- (rfs (- z))) |
---|
633 | (rfs z)) |
---|
634 | (fs z)))))) |
---|
635 | |
---|
636 | (defun fresnel-c (z) |
---|
637 | "Fresnel C: |
---|
638 | |
---|
639 | C(z) = integrate(cos(%pi*t^2/2), t, 0, z) " |
---|
640 | (let ((sqrt-pi (sqrt (float-pi z)))) |
---|
641 | (flet ((fs (z) |
---|
642 | ;; Wolfram gives |
---|
643 | ;; |
---|
644 | ;; C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z)) |
---|
645 | ;; |
---|
646 | ;; where c = sqrt(%pi)/2*(1+%i). |
---|
647 | (* #c(1/4 -1/4) |
---|
648 | (+ (erf (* #c(1/2 1/2) sqrt-pi z)) |
---|
649 | (* #c(0 1) |
---|
650 | (erf (* #c(1/2 -1/2) sqrt-pi z))))))) |
---|
651 | (if (realp z) |
---|
652 | ;; FresnelS is real for a real argument. And it is odd. |
---|
653 | (if (minusp z) |
---|
654 | (- (realpart (fs (- z)))) |
---|
655 | (realpart (fs z))) |
---|
656 | (fs z))))) |
---|
657 | |
---|
658 | (defun sin-integral (z) |
---|
659 | "Sin integral: |
---|
660 | |
---|
661 | Si(z) = integrate(sin(t)/t, t, 0, z)" |
---|
662 | ;; Wolfram has |
---|
663 | ;; |
---|
664 | ;; Si(z) = %i/2*(gamma_inc_tail(0, -%i*z) - gamma_inc_tail(0, %i*z) + log(-%i*z)-log(%i*z)) |
---|
665 | ;; |
---|
666 | (flet ((si (z) |
---|
667 | (* #c(0 1/2) |
---|
668 | (let ((iz (* #c(0 1) z)) |
---|
669 | (-iz (* #c(0 -1) z))) |
---|
670 | (+ (- (incomplete-gamma-tail 0 -iz) |
---|
671 | (incomplete-gamma-tail 0 iz)) |
---|
672 | (- (log -iz) |
---|
673 | (log iz))))))) |
---|
674 | (if (realp z) |
---|
675 | ;; Si is odd and real for real z. In this case, we have |
---|
676 | ;; |
---|
677 | ;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi) |
---|
678 | ;; = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x)) |
---|
679 | ;; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so |
---|
680 | ;; |
---|
681 | ;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x)) |
---|
682 | (cond ((< z 0) |
---|
683 | (- (sin-integral (- z)))) |
---|
684 | ((= z 0) |
---|
685 | (* 0 z)) |
---|
686 | (t |
---|
687 | (+ (* 1/2 (float-pi z)) |
---|
688 | (imagpart (incomplete-gamma-tail 0 (complex 0 z)))))) |
---|
689 | (si z)))) |
---|
690 | |
---|
691 | (defun cos-integral (z) |
---|
692 | "Cos integral: |
---|
693 | |
---|
694 | Ci(z) = integrate((cos(t) - 1)/t, t, 0, z) + log(z) + gamma |
---|
695 | |
---|
696 | where gamma is Euler-Mascheroni constant" |
---|
697 | ;; Wolfram has |
---|
698 | ;; |
---|
699 | ;; Ci(z) = log(z) - 1/2*(gamma_inc_tail(0, -%i*z) + gamma_inc_tail(0, %i*z) + log(-%i*z)+log(%i*z)) |
---|
700 | ;; |
---|
701 | (flet ((ci (z) |
---|
702 | (- (log z) |
---|
703 | (* 1/2 |
---|
704 | (let ((iz (* #c(0 1) z)) |
---|
705 | (-iz (* #c(0 -1) z))) |
---|
706 | (+ (+ (incomplete-gamma-tail 0 -iz) |
---|
707 | (incomplete-gamma-tail 0 iz)) |
---|
708 | (+ (log -iz) |
---|
709 | (log iz)))))))) |
---|
710 | (if (and (realp z) (plusp z)) |
---|
711 | (realpart (ci z)) |
---|
712 | (ci z)))) |
---|
713 | |
---|
714 | ;; Array of values of the Bernoulli numbers. We only have enough for |
---|
715 | ;; the evaluation of the psi function. |
---|
716 | (defconstant bern-values |
---|
717 | (make-array 55 |
---|
718 | :initial-contents |
---|
719 | '(1 |
---|
720 | -1/2 |
---|
721 | 1/6 |
---|
722 | 0 |
---|
723 | -1/30 |
---|
724 | 0 |
---|
725 | 1/42 |
---|
726 | 0 |
---|
727 | -1/30 |
---|
728 | 0 |
---|
729 | 5/66 |
---|
730 | 0 |
---|
731 | -691/2730 |
---|
732 | 0 |
---|
733 | 7/6 |
---|
734 | 0 |
---|
735 | -3617/510 |
---|
736 | 0 |
---|
737 | 43867/798 |
---|
738 | 0 |
---|
739 | -174611/330 |
---|
740 | 0 |
---|
741 | 854513/138 |
---|
742 | 0 |
---|
743 | -236364091/2730 |
---|
744 | 0 |
---|
745 | 8553103/6 |
---|
746 | 0 |
---|
747 | -23749461029/870 |
---|
748 | 0 |
---|
749 | 8615841276005/14322 |
---|
750 | 0 |
---|
751 | -7709321041217/510 |
---|
752 | 0 |
---|
753 | 2577687858367/6 |
---|
754 | 0 |
---|
755 | -26315271553053477373/1919190 |
---|
756 | 0 |
---|
757 | 2929993913841559/6 |
---|
758 | 0 |
---|
759 | -261082718496449122051/13530 |
---|
760 | 0 |
---|
761 | 1520097643918070802691/1806 |
---|
762 | 0 |
---|
763 | -27833269579301024235023/690 |
---|
764 | 0 |
---|
765 | 596451111593912163277961/282 |
---|
766 | 0 |
---|
767 | -5609403368997817686249127547/46410 |
---|
768 | 0 |
---|
769 | 495057205241079648212477525/66 |
---|
770 | 0 |
---|
771 | -801165718135489957347924991853/1590 |
---|
772 | 0 |
---|
773 | 29149963634884862421418123812691/798 |
---|
774 | ))) |
---|
775 | |
---|
776 | (defun bern (k) |
---|
777 | (aref bern-values k)) |
---|
778 | |
---|
779 | (defun psi (z) |
---|
780 | "Digamma function defined by |
---|
781 | |
---|
782 | - %gamma + sum(1/k-1/(k+z-1), k, 1, inf) |
---|
783 | |
---|
784 | where %gamma is Euler's constant" |
---|
785 | |
---|
786 | ;; A&S 6.3.7: Reflection formula |
---|
787 | ;; |
---|
788 | ;; psi(1-z) = psi(z) + %pi*cot(%pi*z) |
---|
789 | ;; |
---|
790 | ;; A&S 6.3.6: Recurrence formula |
---|
791 | ;; |
---|
792 | ;; psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z) |
---|
793 | ;; |
---|
794 | ;; A&S 6.3.8: Asymptotic formula |
---|
795 | ;; |
---|
796 | ;; psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf) |
---|
797 | ;; |
---|
798 | ;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence |
---|
799 | ;; formula to increase the argument and then apply the asymptotic formula. |
---|
800 | |
---|
801 | (cond ((= z 1) |
---|
802 | ;; psi(1) = -%gamma |
---|
803 | (- (float +%gamma+ (if (integerp z) 0.0 z)))) |
---|
804 | ((minusp (realpart z)) |
---|
805 | (let ((p (float-pi z))) |
---|
806 | (flet ((cot-pi (z) |
---|
807 | ;; cot(%pi*z), carefully. If z is an odd multiple |
---|
808 | ;; of 1/2, cot is 0. |
---|
809 | (if (and (realp z) |
---|
810 | (= 1/2 (abs (- z (ftruncate z))))) |
---|
811 | (float 0 z) |
---|
812 | (/ (tan (* p z)))))) |
---|
813 | (- (psi (- 1 z)) |
---|
814 | (* p (cot-pi z)))))) |
---|
815 | (t |
---|
816 | (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon (float (realpart z))) 10))))))) |
---|
817 | (m 0) |
---|
818 | (y (expt (+ z k) 2)) |
---|
819 | (x 0)) |
---|
820 | (loop for i from 1 upto (floor k 2) do |
---|
821 | (progn |
---|
822 | (incf m (+ (/ (+ z i i -1)) |
---|
823 | (/ (+ z i i -2)))) |
---|
824 | (setf x (/ (+ x (/ (bern (+ k 2 (* -2 i))) |
---|
825 | (- k i i -2))) |
---|
826 | y)))) |
---|
827 | (- (log (+ z k)) |
---|
828 | (/ (* 2 (+ z k))) |
---|
829 | x |
---|
830 | m))))) |
---|
831 | |
---|