source: qd-gamma.lisp @ 104efd

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

VALUE-OR-TINY returned a value that was too tiny.

  • qd-gamma.lisp::
    • Return sqrt(least-positive-normalized-double) instead of least-positive-normalized-double.
  • rt-tests.lisp::
    • Add test for this case.
  • Property mode set to 100644
File size: 24.2 KB
Line 
1;;;; -*- Mode: lisp -*-
2;;;;
3;;;; Copyright (c) 2011 Raymond Toy
4;;;; Permission is hereby granted, free of charge, to any person
5;;;; obtaining a copy of this software and associated documentation
6;;;; files (the "Software"), to deal in the Software without
7;;;; restriction, including without limitation the rights to use,
8;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
9;;;; copies of the Software, and to permit persons to whom the
10;;;; Software is furnished to do so, subject to the following
11;;;; conditions:
12;;;;
13;;;; The above copyright notice and this permission notice shall be
14;;;; included in all copies or substantial portions of the Software.
15;;;;
16;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
18;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
20;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
21;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
22;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
23;;;; OTHER DEALINGS IN THE SOFTWARE.
24
25(in-package #:oct)
26
27(eval-when (:compile-toplevel :load-toplevel :execute)
28  (setf *readtable* *oct-readtable*))
29
30;; For log-gamma we use the asymptotic formula
31;;
32;; log(gamma(z)) ~ (z - 1/2)*log(z) + log(2*%pi)/2
33;;                   + sum(bern(2*k)/(2*k)/(2*k-1)/z^(2k-1), k, 1, inf)
34;;
35;;               = (z - 1/2)*log(z) + log(2*%pi)/2
36;;                  + 1/12/z*(1 - 1/30/z^2 + 1/105/z^4 + 1/140/z^6 + ...
37;;                              + 174611/10450/z^18 + ...)
38;;
39;; For double-floats, let's stop the series at the power z^18.  The
40;; next term is 77683/483/z^20.  This means that for |z| > 8.09438,
41;; the series has double-float precision.
42;;
43;; For quad-doubles, let's stop the series at the power z^62.  The
44;; next term is about 6.364d37/z^64.  So for |z| > 38.71, the series
45;; has quad-double precision.
46(defparameter *log-gamma-asymp-coef*
47  #(-1/30 1/105 -1/140 1/99 -691/30030 1/13 -3617/10200 43867/20349
48    -174611/10450 77683/483 -236364091/125580 657931/25 -3392780147/7830
49    1723168255201/207669 -7709321041217/42160 151628697551/33
50    -26315271553053477373/201514950 154210205991661/37
51    -261082718496449122051/1758900 1520097643918070802691/259161
52    -2530297234481911294093/9890 25932657025822267968607/2115
53    -5609403368997817686249127547/8725080 19802288209643185928499101/539
54    -61628132164268458257532691681/27030 29149963634884862421418123812691/190323
55    -354198989901889536240773677094747/31900
56    2913228046513104891794716413587449/3363
57    -1215233140483755572040304994079820246041491/16752085350
58    396793078518930920708162576045270521/61
59    -106783830147866529886385444979142647942017/171360
60    133872729284212332186510857141084758385627191/2103465
61    ))
62
63#+nil
64(defun log-gamma-asymp-series (z nterms)
65  ;; Sum the asymptotic formula for n terms
66  ;;
67  ;; 1 + sum(c[k]/z^(2*k+2), k, 0, nterms)
68  (let ((z2 (* z z))
69        (sum 1)
70        (term 1))
71    (dotimes (k nterms)
72      (setf term (* term z2))
73      (incf sum (/ (aref *log-gamma-asymp-coef* k) term)))
74    sum))
75
76(defun log-gamma-asymp-series (z nterms)
77  (loop with y = (* z z)
78     for k from 1 to nterms
79     for x = 0 then
80       (setf x (/ (+ x (aref *log-gamma-asymp-coef* (- nterms k)))
81                  y))
82     finally (return (+ 1 x))))
83       
84
85(defun log-gamma-asymp-principal (z nterms log2pi/2)
86  (+ (- (* (- z 1/2)
87           (log z))
88        z)
89     log2pi/2))
90
91(defun log-gamma-asymp (z nterms log2pi/2)
92  (+ (log-gamma-asymp-principal z nterms log2pi/2)
93     (* 1/12 (/ (log-gamma-asymp-series z nterms) z))))
94
95(defun log2pi/2 (precision)
96  (ecase precision
97    (single-float
98     (coerce (/ (log (* 2 pi)) 2) 'single-float))
99    (double-float
100     (coerce (/ (log (* 2 pi)) 2) 'double-float))
101    (qd-real
102     (/ (log +2pi+) 2))))
103
104(defun log-gamma-aux (z limit nterms)
105  (let ((precision (float-contagion z)))
106    (cond ((minusp (realpart z))
107           ;; Use reflection formula if realpart(z) < 0
108           ;;   log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z))
109           ;; Or
110           ;;   log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z))
111           (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
Note: See TracBrowser for help on using the repository browser.