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