root/qd-gamma.lisp @ 19dd5231c42c6eec40a7876bd9ffb2661fb05012

Revision 19dd5231c42c6eec40a7876bd9ffb2661fb05012, 24.3 KB (checked in by Raymond Toy <rtoy@…>, 2 years ago)

Cleanups for psi and bug fix for exp-integral-e.

* Allow (exp-integral-e 1 z) to work.
* psi

  • Handle psi(1) specially.
  • Do a better job with cot(%pi*z) when z is an odd multiple of 1/2 where cot is 0.
  • Fib bug in computing the number of terms when we try to float a complex. Just float the realpart.
  • Property mode set to 100644
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           (- (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 (zerop a)
379        ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
380        (exp-integral-e 1 z)
381        (if (and (zerop (imagpart a))
382                 (zerop (imagpart z)))
383            ;; For real values, we split the result to compute either the
384            ;; tail directly from the continued fraction or from gamma(a)
385            ;; - incomplete-gamma.  The continued fraction doesn't
386            ;; converge on the negative real axis, so we can't use that
387            ;; there.  And accuracy appears to be better if z is "small".
388            ;; We take this to mean |z| < |a-1|.  Note that |a-1| is the
389            ;; peak of the integrand.
390            (if (and (> (abs z) (abs (- a 1)))
391                     (not (minusp (realpart z))))
392                (cf-incomplete-gamma-tail a z)
393                (- (gamma a) (cf-incomplete-gamma a z)))
394            ;; If the argument is close enough to the negative real axis,
395            ;; the continued fraction for the tail is not very accurate.
396            ;; Use the incomplete gamma function to evaluate in this
397            ;; region.  (Arbitrarily selected the region to be a sector.
398            ;; But what is the correct size of this sector?)
399            (if (<= (abs (phase z)) 3.1)
400                (cf-incomplete-gamma-tail a z)
401                (- (gamma a) (cf-incomplete-gamma a z)))))))
402
403(defun incomplete-gamma (a z)
404  "Incomplete gamma function defined by:
405
406  integrate(t^(a-1)*exp(-t), t, 0, z)"
407  (let* ((prec (float-contagion a z))
408         (a (apply-contagion a prec))
409         (z (apply-contagion z prec)))
410    (if (and (< (abs a) 1) (< (abs z) 1))
411        (s-incomplete-gamma a z)
412        (if (and (realp a) (realp z))
413            (if (< z (- a 1))
414                (cf-incomplete-gamma a z)
415                (- (gamma a) (cf-incomplete-gamma-tail a z)))
416            ;; The continued fraction doesn't converge very fast if a
417            ;; and z are small.  In this case, use the series
418            ;; expansion instead, which converges quite rapidly.
419            (if (< (abs z) (abs a))
420                (cf-incomplete-gamma a z)
421                (- (gamma a) (cf-incomplete-gamma-tail a z)))))))
422
423(defun erf (z)
424  "Error function:
425
426    erf(z) = 2/sqrt(%pi)*sum((-1)^k*z^(2*k+1)/k!/(2*k+1), k, 0, inf)
427
428  For real z, this is equivalent to
429
430    erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z) for real z."
431  ;;
432  ;; Erf is an odd function: erf(-z) = -erf(z)
433  (if (minusp (realpart z))
434      (- (erf (- z)))
435      (/ (incomplete-gamma 1/2 (* z z))
436         (sqrt (float-pi z)))))
437
438(defun erfc (z)
439  "Complementary error function:
440
441    erfc(z) = 1 - erf(z)"
442  ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
443  ;; near 1.  Wolfram says
444  ;;
445  ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
446  ;;
447  ;; For real(z) > 0, sqrt(z^2)/z is 1 so
448  ;;
449  ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
450  ;;         = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
451  ;;
452  ;; For real(z) < 0, sqrt(z^2)/z is -1 so
453  ;;
454  ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
455  ;;         = 1 + 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
456  (if (>= (realpart z) 0)
457      (/ (incomplete-gamma-tail 1/2 (* z z))
458         (sqrt (float-pi z)))
459      (+ 1
460         (/ (incomplete-gamma 1/2 (* z z))
461            (sqrt (float-pi z))))))
462
463(defun cf-exp-integral-e (v z)
464  ;; We use the continued fraction
465  ;;
466  ;; E(v,z) = exp(-z)/cf(z)
467  ;;
468  ;; where the continued fraction cf(z) is
469  ;;
470  ;; a[k] = -k*(k+v-1)
471  ;; b[k] = v + 2*k + z
472  ;;
473  ;; for k = 1, inf
474  (let ((z+v (+ z v)))
475    (/ (exp (- z))
476       (lentz #'(lambda (k)
477                  (+ z+v (* 2 k)))
478              #'(lambda (k)
479                  (* (- k)
480                     (+ k v -1)))))))
481
482
483;; For v not an integer:
484;;
485;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
486;;
487;; For v an integer:
488;;
489;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z))
490;;          - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1)
491;;
492(defun s-exp-integral-e (v z)
493  ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
494  (let ((-z (- z))
495        (-v (- v))
496        (eps (epsilon z)))
497    (if (and (realp v)
498             (= v (ftruncate v)))
499        ;; v is an integer
500        (let* ((n (truncate v))
501               (n-1 (1- n)))
502          (- (* (/ (expt -z n-1)
503                   (gamma v))
504                (- (psi v) (log z)))
505             (loop for k from 0
506                   for term = 1 then (* term (/ -z k))
507                   for sum = (if (= v 1) 0 (/ (- 1 v)))
508                     then (+ sum (let ((denom (- k n-1)))
509                                   (if (zerop denom)
510                                       0
511                                       (/ term denom))))
512                   when (< (abs term) (* (abs sum) eps))
513                     return sum)))
514        (loop for k from 0
515              for term = 1 then (* term (/ -z k))
516              for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
517              when (< (abs term) (* (abs sum) eps))
518                return (- (* (gamma (- 1 v)) (expt z (- v 1)))
519                                    sum)))))
520
521(defun exp-integral-e (v z)
522  "Exponential integral E:
523
524   E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)"
525  ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);
526  ;;
527  ;; for |arg(z)| < pi.
528  ;;
529  ;;
530  (cond ((and (realp v) (minusp v))
531         ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
532         (let ((-v (- v)))
533           (* (expt z (- v 1))
534              (incomplete-gamma-tail (+ -v 1) z))))
535        ((< (abs z) 1)
536         ;; Use series for small z
537         (s-exp-integral-e v z))
538        ((>= (abs (phase z)) 3.1)
539         ;; The continued fraction doesn't converge on the negative
540         ;; real axis, and converges very slowly near the negative
541         ;; real axis, so use the incomplete-gamma-tail function in
542         ;; this region.  "Closeness" to the negative real axis is
543         ;; teken to mean that z is in a sector near the axis.
544         ;;
545         ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
546         (* (expt z (- v 1))
547            (incomplete-gamma-tail (- 1 v) z)))
548        (t
549         ;; Use continued fraction for everything else.
550         (cf-exp-integral-e v z))))
551
552;; Series for Fresnel S
553;;
554;;   S(z) = z^3*sum((%pi/2)^(2*k+1)(-z^4)^k/(2*k+1)!/(4*k+3), k, 0, inf)
555;;
556;; Compute as
557;;
558;;   S(z) = z^3*sum(a(k)/(4*k+3), k, 0, inf)
559;;
560;; where
561;;
562;;   a(k+1) = -a(k) * (%pi/2)^2 * z^4 / (2*k+2) / (2*k+3)
563;;
564;;   a(0) = %pi/2.
565(defun fresnel-s-series (z)
566  (let* ((pi/2 (* 1/2 (float-pi z)))
567         (factor (- (* (expt z 4) pi/2 pi/2)))
568         (eps (epsilon z))
569         (sum 0)
570         (term pi/2))
571    (loop for k2 from 0 by 2
572       until (< (abs term) (* eps (abs sum)))
573       do
574       (incf sum (/ term (+ 3 k2 k2)))
575       (setf term (/ (* term factor)
576                     (* (+ k2 2)
577                        (+ k2 3)))))
578    (* sum (expt z 3))))
579   
580(defun fresnel-s (z)
581  "Fresnel S:
582
583   S(z) = integrate(sin(%pi*t^2/2), t, 0, z) "
584  (let ((prec (float-contagion z))
585        (sqrt-pi (sqrt (float-pi z))))
586    (flet ((fs (z)
587             ;; Wolfram gives
588             ;;
589             ;;  S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z))
590             ;;
591             ;; where c = sqrt(%pi)/2*(1+%i).
592             ;;
593             ;; But for large z, we should use erfc.  Then
594             ;;  S(z) = 1/2 - (1+%i)/4*(erfc(c*z) - %i*erfc(conjugate(c)*z))
595             (if (and t (> (abs z) 2))
596                 (- 1/2
597                    (* #c(1/4 1/4)
598                       (- (erfc (* #c(1/2 1/2) sqrt-pi z))
599                          (* #c(0 1)
600                             (erfc (* #c(1/2 -1/2) sqrt-pi z))))))
601                 (* #c(1/4 1/4)
602                    (- (erf (* #c(1/2 1/2) sqrt-pi z))
603                       (* #c(0 1)
604                          (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
605          (rfs (z)
606            ;; When z is real, recall that erf(conjugate(z)) =
607            ;; conjugate(erf(z)).  Then
608            ;;
609            ;;  S(z) = 1/2*(realpart(erf(c*z)) - imagpart(erf(c*z)))
610            ;;
611            ;; But for large z, we should use erfc.  Then
612            ;;
613            ;;  S(z) = 1/2 - 1/2*(realpart(erfc(c*z)) - imagpart(erf(c*z)))
614            (if (> (abs z) 2)
615                (let ((s (erfc (* #c(1/2 1/2) sqrt-pi z))))
616                  (- 1/2
617                     (* 1/2 (- (realpart s) (imagpart s)))))
618                (let ((s (erf (* #c(1/2 1/2) sqrt-pi z))))
619                  (* 1/2 (- (realpart s) (imagpart s)))))))
620      ;; For small z, the erf terms above suffer from subtractive
621      ;; cancellation.  So use the series in this case.  Some simple
622      ;; tests were done to determine that for double-floats we want
623      ;; to use the series for z < 1 to give max accuracy.  For
624      ;; qd-real, the above formula is good enough for z > 1d-5.
625      (if (< (abs z) (ecase prec
626                       (single-float 1.5f0)
627                       (double-float 1d0)
628                       (qd-real #q1)))
629          (fresnel-s-series z)
630          (if (realp z)
631              ;; FresnelS is real for a real argument. And it is odd.
632              (if (minusp z)
633                  (- (rfs (- z)))
634                  (rfs z))
635              (fs z))))))
636
637(defun fresnel-c (z)
638  "Fresnel C:
639
640   C(z) = integrate(cos(%pi*t^2/2), t, 0, z) "
641  (let ((sqrt-pi (sqrt (float-pi z))))
642    (flet ((fs (z)
643             ;; Wolfram gives
644             ;;
645             ;;  C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z))
646             ;;
647             ;; where c = sqrt(%pi)/2*(1+%i).
648             (* #c(1/4 -1/4)
649                (+ (erf (* #c(1/2 1/2) sqrt-pi z))
650                   (* #c(0 1)
651                      (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
652      (if (realp z)
653          ;; FresnelS is real for a real argument. And it is odd.
654          (if (minusp z)
655              (- (realpart (fs (- z))))
656              (realpart (fs z)))
657          (fs z)))))
658
659(defun sin-integral (z)
660  "Sin integral:
661
662  Si(z) = integrate(sin(t)/t, t, 0, z)"
663  ;; Wolfram has
664  ;;
665  ;; Si(z) = %i/2*(gamma_inc_tail(0, -%i*z) - gamma_inc_tail(0, %i*z) + log(-%i*z)-log(%i*z))
666  ;;
667  (flet ((si (z)
668           (* #c(0 1/2)
669              (let ((iz (* #c(0 1) z))
670                    (-iz (* #c(0 -1) z)))
671                (+ (- (incomplete-gamma-tail 0 -iz)
672                      (incomplete-gamma-tail 0 iz))
673                   (- (log -iz)
674                      (log iz)))))))
675    (if (realp z)
676        ;; Si is odd and real for real z.  In this case, we have
677        ;;
678        ;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi)
679        ;;       = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x))
680        ;; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so
681        ;;
682        ;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x))
683        (cond ((< z 0)
684               (- (sin-integral (- z))))
685              ((= z 0)
686               (* 0 z))
687              (t
688               (+ (* 1/2 (float-pi z))
689                  (imagpart (incomplete-gamma-tail 0 (complex 0 z))))))
690        (si z))))
691
692(defun cos-integral (z)
693  "Cos integral:
694
695    Ci(z) = integrate((cos(t) - 1)/t, t, 0, z) + log(z) + gamma
696
697  where gamma is Euler-Mascheroni constant"
698  ;; Wolfram has
699  ;;
700  ;; Ci(z) = log(z) - 1/2*(gamma_inc_tail(0, -%i*z) + gamma_inc_tail(0, %i*z) + log(-%i*z)+log(%i*z))
701  ;;
702  (flet ((ci (z)
703           (- (log z)
704              (* 1/2
705                 (let ((iz (* #c(0 1) z))
706                       (-iz (* #c(0 -1) z)))
707                   (+ (+ (incomplete-gamma-tail 0 -iz)
708                         (incomplete-gamma-tail 0 iz))
709                      (+ (log -iz)
710                         (log iz))))))))
711    (if (and (realp z) (plusp z))
712        (realpart (ci z))
713        (ci z))))
714
715;; Array of values of the Bernoulli numbers.  We only have enough for
716;; the evaluation of the psi function.
717(defconstant bern-values
718  (make-array 55
719              :initial-contents
720              '(1
721                -1/2
722                1/6
723                0
724                -1/30
725                0
726                1/42
727                0
728                -1/30
729                0
730                5/66
731                0
732                -691/2730
733                0
734                7/6
735                0
736                -3617/510
737                0
738                43867/798
739                0
740                -174611/330
741                0
742                854513/138
743                0
744                -236364091/2730
745                0
746                8553103/6
747                0
748                -23749461029/870
749                0
750                8615841276005/14322
751                0
752                -7709321041217/510
753                0
754                2577687858367/6
755                0
756                -26315271553053477373/1919190
757                0
758                2929993913841559/6
759                0
760                -261082718496449122051/13530
761                0
762                1520097643918070802691/1806
763                0
764                -27833269579301024235023/690
765                0
766                596451111593912163277961/282
767                0
768                -5609403368997817686249127547/46410
769                0
770                495057205241079648212477525/66
771                0
772                -801165718135489957347924991853/1590
773                0
774                29149963634884862421418123812691/798
775                )))
776               
777(defun bern (k)
778  (aref bern-values k))
779
780(defun psi (z)
781  "Digamma function defined by
782
783  - %gamma + sum(1/k-1/(k+z-1), k, 1, inf)
784
785  where %gamma is Euler's constant"
786
787  ;; A&S 6.3.7:  Reflection formula
788  ;;
789  ;;   psi(1-z) = psi(z) + %pi*cot(%pi*z)
790  ;;
791  ;; A&S 6.3.6:  Recurrence formula
792  ;;
793  ;;   psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z)
794  ;;
795  ;; A&S 6.3.8:  Asymptotic formula
796  ;;
797  ;;   psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf)
798  ;;
799  ;; So use reflection formula if Re(z) < 0.  For z > 0, use the recurrence
800  ;; formula to increase the argument and then apply the asymptotic formula.
801
802  (cond ((= z 1)
803         ;; psi(1) = -%gamma
804         (- (float +%gamma+ (if (integerp z) 0.0 z))))
805        ((minusp (realpart z))
806         (let ((p (float-pi z)))
807           (flet ((cot-pi (z)
808                    ;; cot(%pi*z), carefully.  If z is an odd multiple
809                    ;; of 1/2, cot is 0.
810                    (if (and (realp z)
811                             (= 1/2 (- z (ftruncate z))))
812                        (float 0 z)
813                        (/ (tan (* p z))))))
814             (- (psi (- 1 z))
815                (* p (cot-pi z))))))
816        (t
817         (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon (float (realpart z))) 10)))))))
818                (m 0)
819                (y (expt (+ z k) 2))
820                (x 0))
821           (loop for i from 1 upto (floor k 2) do
822             (progn
823               (incf m (+ (/ (+ z i i -1))
824                          (/ (+ z i i -2))))
825               (setf x (/ (+ x (/ (bern (+ k 2 (* -2 i)))
826                                  (- k i i -2)))
827                          y))))
828           (- (log (+ z k))
829              (/ (* 2 (+ z k)))
830              x
831              m)))))
832
Note: See TracBrowser for help on using the browser.