[a40062c] | 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 | |
---|
[9d3c8c] | 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) |
---|
[8227b2] | 270 | "Tail of the incomplete gamma function defined by: |
---|
| 271 | |
---|
| 272 | integrate(t^(a-1)*exp(-t), t, z, inf)" |
---|
[9d3c8c] | 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) |
---|
[8227b2] | 285 | "Incomplete gamma function defined by: |
---|
| 286 | |
---|
| 287 | integrate(t^(a-1)*exp(-t), t, 0, z)" |
---|
[9d3c8c] | 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))) |
---|
[d17228] | 295 | (if (< (abs z) (abs a)) |
---|
| 296 | (cf-incomplete-gamma a z) |
---|
| 297 | (- (gamma a) (cf-incomplete-gamma-tail a z)))))) |
---|
[9d3c8c] | 298 | |
---|
| 299 | (defun erf (z) |
---|
[8227b2] | 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)))))) |
---|
[46e719] | 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))) |
---|
[d17228] | 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 | |
---|