;;;; -*- Mode: lisp -*-;;;;;;;; Copyright (c) 2011 Raymond Toy;;;; Permission is hereby granted, free of charge, to any person;;;; obtaining a copy of this software and associated documentation;;;; files (the "Software"), to deal in the Software without;;;; restriction, including without limitation the rights to use,;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell;;;; copies of the Software, and to permit persons to whom the;;;; Software is furnished to do so, subject to the following;;;; conditions:;;;;;;;; The above copyright notice and this permission notice shall be;;;; included in all copies or substantial portions of the Software.;;;;;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR;;;; OTHER DEALINGS IN THE SOFTWARE.(in-package#:oct)(eval-when(:compile-toplevel:load-toplevel:execute)(setf*readtable**oct-readtable*));; For log-gamma we use the asymptotic formula;;;; log(gamma(z)) ~ (z - 1/2)*log(z) + log(2*%pi)/2;; + sum(bern(2*k)/(2*k)/(2*k-1)/z^(2k-1), k, 1, inf);;;; = (z - 1/2)*log(z) + log(2*%pi)/2;; + 1/12/z*(1 - 1/30/z^2 + 1/105/z^4 + 1/140/z^6 + ...;; + 174611/10450/z^18 + ...);;;; For double-floats, let's stop the series at the power z^18. The;; next term is 77683/483/z^20. This means that for |z| > 8.09438,;; the series has double-float precision.;;;; For quad-doubles, let's stop the series at the power z^62. The;; next term is about 6.364d37/z^64. So for |z| > 38.71, the series;; has quad-double precision.(defparameter*log-gamma-asymp-coef*#(-1/301/105-1/1401/99-691/300301/13-3617/1020043867/20349-174611/1045077683/483-236364091/125580657931/25-3392780147/78301723168255201/207669-7709321041217/42160151628697551/33-26315271553053477373/201514950154210205991661/37-261082718496449122051/17589001520097643918070802691/259161-2530297234481911294093/989025932657025822267968607/2115-5609403368997817686249127547/872508019802288209643185928499101/539-61628132164268458257532691681/2703029149963634884862421418123812691/190323-354198989901889536240773677094747/319002913228046513104891794716413587449/3363-1215233140483755572040304994079820246041491/16752085350396793078518930920708162576045270521/61-106783830147866529886385444979142647942017/171360133872729284212332186510857141084758385627191/2103465))#+nil(defunlog-gamma-asymp-series(znterms);; Sum the asymptotic formula for n terms;;;; 1 + sum(c[k]/z^(2*k+2), k, 0, nterms)(let((z2(*zz))(sum1)(term1))(dotimes(knterms)(setfterm(*termz2))(incfsum(/(aref*log-gamma-asymp-coef*k)term)))sum))(defunlog-gamma-asymp-series(znterms)(loopwithy=(*zz)forkfrom1tontermsforx=0then(setfx(/(+x(aref*log-gamma-asymp-coef*(-ntermsk)))y))finally(return(+1x))))(defunlog-gamma-asymp-principal(zntermslog2pi/2)(+(-(*(-z1/2)(logz))z)log2pi/2))(defunlog-gamma-asymp(zntermslog2pi/2)(+(log-gamma-asymp-principalzntermslog2pi/2)(*1/12(/(log-gamma-asymp-seriesznterms)z))))(defunlog2pi/2(precision)(ecaseprecision(single-float(coerce(/(log(*2pi))2)'single-float))(double-float(coerce(/(log(*2pi))2)'double-float))(qd-real(/(log+2pi+)2))))(defunlog-gamma-aux(zlimitnterms)(let((precision(float-contagionz)))(cond((minusp(realpartz));; Use reflection formula if realpart(z) < 0;; log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z));; Or;; log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z))(let((p(float-piz)))(-(logp)(log(-z))(log(sin(*pz)))(log-gamma(-z)))))(t(let((absz(absz)))(cond((>=abszlimit);; Can use the asymptotic formula directly with 9 terms(log-gamma-asympznterms(log2pi/2precision)))(t;; |z| is too small. Use the formula;; log(gamma(z)) = log(gamma(z+1)) - log(z)(-(log-gamma(+z1))(logz)))))))))(defmethodlog-gamma((zcl:number))(log-gamma-auxz99))(defmethodlog-gamma((zqd-real))(log-gamma-auxz2626))(defmethodlog-gamma((zqd-complex))(log-gamma-auxz2626))(defungamma-aux(zlimitnterms)(let((precision(float-contagionz)))(cond((<=(realpartz)0);; Use reflection formula if realpart(z) < 0:;; gamma(-z) = -pi*csc(pi*z)/gamma(z+1);; or;; gamma(z) = pi*csc(pi*z)/gamma(1-z)(if(and(realpz)(=(truncatez)z));; Gamma of a negative integer is infinity. Signal an error(error"Gamma of non-positive integer ~S"z)(/(float-piz)(sin(*(float-piz)z))(gamma-aux(-1z)limitnterms))))((and(zerop(imagpartz))(=z(truncatez)));; We have gamma(n) where an integer value n and is small;; enough. In this case, just compute the product;; directly. We do this because our current implementation;; has some round-off for these values, and that's annoying;; and unexpected.(let((n(truncatez)))(loopforprod=(apply-contagion1precision)then(*prodk)forkfrom2belownfinally(return(apply-contagionprodprecision)))))(t(let((absz(absz)))(cond((>=abszlimit);; Use log gamma directly:;; log(gamma(z)) = principal part + 1/12/z*(series part);; so;; gamma(z) = exp(principal part)*exp(1/12/z*series)(exp(log-gammaz))#+nil(*(exp(log-gamma-asymp-principalznterms(log2pi/2precision)))(exp(*1/12(/(log-gamma-asymp-seriesznterms)z)))))(t;; 1 <= |z| <= limit;; gamma(z) = gamma(z+1)/z(/(gamma-aux(+1z)limitnterms)z))))))))(defmethodgamma((zcl:number))(gamma-auxz99))(defmethodgamma((zqd-real))(gamma-auxz3932))(defmethodgamma((zqd-complex))(gamma-auxz3932));; Lentz's algorithm for evaluating continued fractions.;;;; Let the continued fraction be:;;;; a1 a2 a3;; b0 + ---- ---- ----;; b1 + b2 + b3 +;;(defvar*debug-cf-eval*nil"When true, enable some debugging prints when evaluating a continued fraction.");; Max number of iterations allowed when evaluating the continued;; fraction. When this is reached, we assume that the continued;; fraction did not converge.(defvar*max-cf-iterations*10000"Max number of iterations allowed when evaluating the continued fraction. When this is reached, we assume that the continued fraction did not converge.")(defunlentz(bfaf)(let((tiny-value-count0))(flet((value-or-tiny(v)(if(zeropv)(progn(incftiny-value-count)(etypecasev((ordouble-floatcl:complex)(sqrtleast-positive-normalized-double-float))((orqd-realqd-complex)(make-qd(sqrtleast-positive-normalized-double-float)))))v)))(let*((f(value-or-tiny(funcallbf0)))(cf)(d0)(eps(epsilonf)))(loopforjfrom1upto*max-cf-iterations*foran=(funcallafj)forbn=(funcallbfj)do(progn(setfd(value-or-tiny(+bn(*and))))(setfc(value-or-tiny(+bn(/anc))))(when*debug-cf-eval*(formatt"~&j = ~d~%"j)(formatt" an = ~s~%"an)(formatt" bn = ~s~%"bn)(formatt" c = ~s~%"c)(formatt" d = ~s~%"d))(let((delta(/cd)))(setfd(/d))(setff(*fdelta))(when*debug-cf-eval*(formatt" dl= ~S (|dl - 1| = ~S)~%"delta(abs(1-delta)))(formatt" f = ~S~%"f))(when(<=(abs(-delta1))eps)(return-fromlentz(valuesfjtiny-value-count)))))finally(error'simple-error:format-control"~<Continued fraction failed to converge after ~D iterations.~% Delta = ~S~>":format-arguments(list*max-cf-iterations*(/cd))))))));; Continued fraction for erf(z):;;;; erf(z) = 2*z/sqrt(pi)*exp(-z^2)/K;;;; where K is the continued fraction with;;;; b[n] = 1+2*n-2*z^2;; a[n] = 4*n*z^2;;;; This works ok, but has problems for z > 3 where sometimes the;; result is greater than 1 and for larger values, negative numbers;; are returned.#+nil(defuncf-erf(z)(let*((z2(*zz))(twoz2(*2z2)))(*(/(*2z)(sqrt(float-piz)))(exp(-z2))(/(lentz#'(lambda(n)(-(+1(*2n))twoz2))#'(lambda(n)(*4nz2)))))));; From Cuyt;;;; sqrt(%pi)*z*exp(z^2)*erf(z) = K;;;; where K is the continued fraction with terms F[m]*z^2/(1 + G[m]*z^2);; with F[1] = 2 and F[m] = 4*(m-1)/(2*m-3)/(2*m-1) and G[m] = -2/(2*m-1);;(defuncf-erf(z)(let((z2(*zz)))(*(/(exp(-z2))(sqrt(float-piz))z)(lentz#'(lambda(n)(if(zeropn)(float0(realpartz))(1+(/(*-2z2)(+nn-1)))))#'(lambda(n)(if(=n1)(*2z2)(/(*z24(-n1))(*(+nn-3)(+nn-1)))))))));; From the above, we also have Dawson's integral:;;;; exp(-z^-2)*integrate(exp(t^2), t, 0, z) = %i*sqrt(%pi)/2*exp(-z^2)*erf(-%i*z);;;;; -2*z*exp(-z^2)*integrate(exp(t^2), t, 0, z) = K;;;; with K = -F[m)*z^2/(1 - G[m]*z^2), where F[m] and G[m] are as above.;;;; Also erf(-%i*z) = dawson(z) * 2*exp(-z^2)/(*%i*sqrt(%pi))(defuncf-dawson(z)(let((z2(*zz)))(/(lentz#'(lambda(n)(if(zeropn)(float0(realpartz))(-1(/(*-2z2)(+nn-1)))))#'(lambda(n)(if(=n1)(*-2z2)(/(*z2-4(-n1))(*(+nn-3)(+nn-1))))))(*-2z))));; erfc(z) = z/sqrt(%pi)*exp(-z^2)*K;;;; where K is the continued fraction with a[1] = 1, a[m] = (m-1)/2,;; for m >= 2 and b[0] = 0, b[2*m+1] = z^2, b[2*m] = 1.;;;; This is valid only if Re(z) > 0.(defuncf-erfc(z)(let((z2(*zz))(zero(float0(realpartz)))(one(float1(realpartz))))(*(exp(-z2))z(/(sqrt(float-pi(realpartz))))(lentz#'(lambda(n)(if(zeropn)zero(if(evenpn)onez2)))#'(lambda(n)(if(=n1)one(/(-n1)2)))))));; w(z) = exp(-z^2)*erfc(-%i*z);;;; = -%i*z/sqrt(%pi)*K;;;; where K is the continued fraction with a[n] the same as for erfc;; and b[0] = 0, b[2*m+1] = -z^2, b[2*m] = 1.;;;; This is valid only if Im(z) > 0. We can use the following;; identities:;;;; w(-z) = 2*exp(-z^2) - w(z);; w(conj(z)) = conj(w(-z))(defuncf-w(z)(let((z2(*zz))(zero(float0(realpartz)))(one(float1(realpartz))))(*#c(0-1)z(/(sqrt(float-pi(realpartz))))(lentz#'(lambda(n)(if(zeropn)zero(if(evenpn)one(-z2))))#'(lambda(n)(if(=n1)one(/(-n1)2)))))));; Tail of the incomplete gamma function:;; integrate(x^(a-1)*exp(-x), x, z, inf);;;; The continued fraction, valid for all z except the negative real;; axis:;;;; b[n] = 1+2*n+z-a;; a[n] = n*(a-n);;;; See http://functions.wolfram.com/06.06.10.0003.01(defuncf-incomplete-gamma-tail(az)(when(and(zerop(imagpartz))(minusp(realpartz)))(error'domain-error:function-name'cf-incomplete-gamma-tail:format-arguments(list'zz):format-control"Argument ~S should not be on the negative real axis: ~S"))(/(handler-case(*(exptza)(exp(-z)))(arithmetic-error();; z^a*exp(-z) can overflow prematurely. In this case, use;; the equivalent exp(a*log(z)-z). We don't use this latter;; form because it has more roundoff error than the former.(exp(-(*a(logz))z))))(let((z-a(-za)))(lentz#'(lambda(n)(+nn1z-a))#'(lambda(n)(*n(-an)))))));; Incomplete gamma function:;; integrate(x^(a-1)*exp(-x), x, 0, z);;;; The continued fraction, valid for all z:;;;; b[n] = n - 1 + z + a;; a[n] = -z*(a + n);;;; See http://functions.wolfram.com/06.06.10.0007.01. We modified the;; continued fraction slightly and discarded the first quotient from;; the fraction.#+nil(defuncf-incomplete-gamma(az)(/(handler-case(*(exptza)(exp(-z)))(arithmetic-error();; z^a*exp(-z) can overflow prematurely. In this case, use;; the equivalent exp(a*log(z)-z). We don't use this latter;; form because it has more roundoff error than the former.(exp(-(*a(logz))z))))(let((za1(+za1)))(-a(/(*az)(lentz#'(lambda(n)(+nza1))#'(lambda(n)(-(*z(+an))))))))));; Incomplete gamma function:;; integrate(x^(a-1)*exp(-x), x, 0, z);;;; The continued fraction, valid for all z:;;;; b[n] = a + n;; a[n] = -(a+n/2)*z if n odd;; n/2*z if n even;;;; See http://functions.wolfram.com/06.06.10.0009.01.;;;; Some experiments indicate that this converges faster than the above;; and is actually quite a bit more accurate, expecially near the;; negative real axis.(defuncf-incomplete-gamma(az)(/(handler-case(*(exptza)(exp(-z)))(arithmetic-error();; z^a*exp(-z) can overflow prematurely. In this case, use;; the equivalent exp(a*log(z)-z). We don't use this latter;; form because it has more roundoff error than the former.(exp(-(*a(logz))z))))(lentz#'(lambda(n)(+na))#'(lambda(n)(if(evenpn)(*(ashn-1)z)(-(*(+a(ashn-1))z)))))));; Series expansion for incomplete gamma. Intended for |a|<1 and;; |z|<1. The series is;;;; g(a,z) = z^a * sum((-z)^k/k!/(a+k), k, 0, inf)(defuns-incomplete-gamma(az)(let((-z(-z))(eps(epsilonz)))(loopforkfrom0forterm=1then(*term(/-zk))forsum=(/a)then(+sum(/term(+ak)))when(<(absterm)(*(abssum)eps))return(*sum(exptza)))));; Tail of the incomplete gamma function.(defunincomplete-gamma-tail(az)"Tail of the incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, z, inf)"(with-floating-point-contagion(az)(if(and(realpa)(<=a0));; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z)(*(exptza)(exp-integral-e(-1a)z))(if(and(zerop(imagparta))(zerop(imagpartz)));; For real values, we split the result to compute either the;; tail directly from the continued fraction or from gamma(a);; - incomplete-gamma. The continued fraction doesn't;; converge on the negative real axis, so we can't use that;; there. And accuracy appears to be better if z is "small".;; We take this to mean |z| < |a-1|. Note that |a-1| is the;; peak of the integrand.(if(and(>(absz)(abs(-a1)))(not(minusp(realpartz))))(cf-incomplete-gamma-tailaz)(-(gammaa)(cf-incomplete-gammaaz)));; If the argument is close enough to the negative real axis,;; the continued fraction for the tail is not very accurate.;; Use the incomplete gamma function to evaluate in this;; region. (Arbitrarily selected the region to be a sector.;; But what is the correct size of this sector?)(if(<=(abs(phasez))3.1)(cf-incomplete-gamma-tailaz)(-(gammaa)(cf-incomplete-gammaaz)))))))(defunincomplete-gamma(az)"Incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, 0, z)"(with-floating-point-contagion(az)(if(and(<(absa)1)(<(absz)1))(s-incomplete-gammaaz)(if(and(realpa)(realpz))(if(<z(-a1))(cf-incomplete-gammaaz)(-(gammaa)(cf-incomplete-gamma-tailaz)));; The continued fraction doesn't converge very fast if a;; and z are small. In this case, use the series;; expansion instead, which converges quite rapidly.(if(<(absz)(absa))(cf-incomplete-gammaaz)(-(gammaa)(cf-incomplete-gamma-tailaz)))))))(defunerf(z)"Error function: erf(z) = 2/sqrt(%pi)*sum((-1)^k*z^(2*k+1)/k!/(2*k+1), k, 0, inf) For real z, this is equivalent to erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z) for real z.";;;; Erf is an odd function: erf(-z) = -erf(z)(if(minusp(realpartz))(-(erf(-z)))(/(incomplete-gamma1/2(*zz))(sqrt(float-piz)))))(defunerfc(z)"Complementary error function: erfc(z) = 1 - erf(z)";; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is;; near 1. Wolfram says;;;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2));;;; For real(z) > 0, sqrt(z^2)/z is 1 so;;;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2));; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2);;;; For real(z) < 0, sqrt(z^2)/z is -1 so;;;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2));; = 1 + 1/sqrt(pi)*gamma_incomplete(1/2,z^2)(if(>=(realpartz)0)(/(incomplete-gamma-tail1/2(*zz))(sqrt(float-piz)))(+1(/(incomplete-gamma1/2(*zz))(sqrt(float-piz))))))(defuncf-exp-integral-e(vz);; We use the continued fraction;;;; E(v,z) = exp(-z)/cf(z);;;; where the continued fraction cf(z) is;;;; a[k] = -k*(k+v-1);; b[k] = v + 2*k + z;;;; for k = 1, inf(let((z+v(+zv)))(/(exp(-z))(lentz#'(lambda(k)(+z+v(*2k)))#'(lambda(k)(*(-k)(+kv-1)))))));; For v not an integer:;;;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf);;;; For v an integer:;;;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z));; - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1);;(defuns-exp-integral-e(vz);; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)(let((-z(-z))(-v(-v))(eps(epsilonz)))(if(and(realpv)(=v(ftruncatev)));; v is an integer(let*((n(truncatev))(n-1(1-n)))(-(*(/(expt-zn-1)(gammav))(-(psiv)(logz)))(loopforkfrom0forterm=1then(*term(/-zk))forsum=(if(zeropn-1)0(/(-1v)))then(+sum(let((denom(-kn-1)))(if(zeropdenom)0(/termdenom))))when(<(absterm)(*(abssum)eps))returnsum)))(loopforkfrom0forterm=1then(*term(/-zk))forsum=(/(-1v))then(+sum(/term(+k1-v)))when(<(absterm)(*(abssum)eps))return(-(*(gamma(-1v))(exptz(-v1)))sum)))))(defunexp-integral-e(vz)"Exponential integral E: E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)";; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);;;;; for |arg(z)| < pi.;;;;(with-floating-point-contagion(vz)(cond((and(realpv)(minuspv));; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)(let((-v(-v)))(*(exptz(-v1))(incomplete-gamma-tail(+-v1)z))))((or(<(absz)1)(>=(abs(phasez))3.1));; Use series for small z or if z is near the negative real;; axis because the continued fraction does not converge on;; the negative axis and converges slowly near the negative;; axis.(s-exp-integral-evz))(t;; Use continued fraction for everything else.(cf-exp-integral-evz)))));; Series for Fresnel S;;;; S(z) = z^3*sum((%pi/2)^(2*k+1)(-z^4)^k/(2*k+1)!/(4*k+3), k, 0, inf);;;; Compute as;;;; S(z) = z^3*sum(a(k)/(4*k+3), k, 0, inf);;;; where;;;; a(k+1) = -a(k) * (%pi/2)^2 * z^4 / (2*k+2) / (2*k+3);;;; a(0) = %pi/2.(defunfresnel-s-series(z)(let*((pi/2(*1/2(float-piz)))(factor(-(*(exptz4)pi/2pi/2)))(eps(epsilonz))(sum0)(termpi/2))(loopfork2from0by2until(<(absterm)(*eps(abssum)))do(incfsum(/term(+3k2k2)))(setfterm(/(*termfactor)(*(+k22)(+k23)))))(*sum(exptz3))))(defunfresnel-s(z)"Fresnel S: S(z) = integrate(sin(%pi*t^2/2), t, 0, z) "(let((prec(float-contagionz))(sqrt-pi(sqrt(float-piz))))(flet((fs(z);; Wolfram gives;;;; S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z));;;; where c = sqrt(%pi)/2*(1+%i).;;;; But for large z, we should use erfc. Then;; S(z) = 1/2 - (1+%i)/4*(erfc(c*z) - %i*erfc(conjugate(c)*z))(if(andt(>(absz)2))(-1/2(*#c(1/41/4)(-(erfc(*#c(1/21/2)sqrt-piz))(*#c(01)(erfc(*#c(1/2-1/2)sqrt-piz))))))(*#c(1/41/4)(-(erf(*#c(1/21/2)sqrt-piz))(*#c(01)(erf(*#c(1/2-1/2)sqrt-piz)))))))(rfs(z);; When z is real, recall that erf(conjugate(z)) =;; conjugate(erf(z)). Then;;;; S(z) = 1/2*(realpart(erf(c*z)) - imagpart(erf(c*z)));;;; But for large z, we should use erfc. Then;;;; S(z) = 1/2 - 1/2*(realpart(erfc(c*z)) - imagpart(erf(c*z)))(if(>(absz)2)(let((s(erfc(*#c(1/21/2)sqrt-piz))))(-1/2(*1/2(-(realparts)(imagparts)))))(let((s(erf(*#c(1/21/2)sqrt-piz))))(*1/2(-(realparts)(imagparts)))))));; For small z, the erf terms above suffer from subtractive;; cancellation. So use the series in this case. Some simple;; tests were done to determine that for double-floats we want;; to use the series for z < 1 to give max accuracy. For;; qd-real, the above formula is good enough for z > 1d-5.(if(<(absz)(ecaseprec(single-float1.5f0)(double-float1d0)(qd-real#q1)))(fresnel-s-seriesz)(if(realpz);; FresnelS is real for a real argument. And it is odd.(if(minuspz)(-(rfs(-z)))(rfsz))(fsz))))))(defunfresnel-c(z)"Fresnel C: C(z) = integrate(cos(%pi*t^2/2), t, 0, z) "(let((sqrt-pi(sqrt(float-piz))))(flet((fs(z);; Wolfram gives;;;; C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z));;;; where c = sqrt(%pi)/2*(1+%i).(*#c(1/4-1/4)(+(erf(*#c(1/21/2)sqrt-piz))(*#c(01)(erf(*#c(1/2-1/2)sqrt-piz)))))))(if(realpz);; FresnelS is real for a real argument. And it is odd.(if(minuspz)(-(realpart(fs(-z))))(realpart(fsz)))(fsz)))))(defunsin-integral(z)"Sin integral: Si(z) = integrate(sin(t)/t, t, 0, z)";; Wolfram has;;;; Si(z) = %i/2*(gamma_inc_tail(0, -%i*z) - gamma_inc_tail(0, %i*z) + log(-%i*z)-log(%i*z));;(flet((si(z)(*#c(01/2)(let((iz(*#c(01)z))(-iz(*#c(0-1)z)))(+(-(incomplete-gamma-tail0-iz)(incomplete-gamma-tail0iz))(-(log-iz)(logiz)))))))(if(realpz);; Si is odd and real for real z. In this case, we have;;;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi);; = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x));; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so;;;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x))(cond((<z0)(-(sin-integral(-z))))((=z0)(*0z))(t(+(*1/2(float-piz))(imagpart(incomplete-gamma-tail0(complex0z))))))(siz))))(defuncos-integral(z)"Cos integral: Ci(z) = integrate((cos(t) - 1)/t, t, 0, z) + log(z) + gamma where gamma is Euler-Mascheroni constant";; Wolfram has;;;; Ci(z) = log(z) - 1/2*(gamma_inc_tail(0, -%i*z) + gamma_inc_tail(0, %i*z) + log(-%i*z)+log(%i*z));;(flet((ci(z)(-(logz)(*1/2(let((iz(*#c(01)z))(-iz(*#c(0-1)z)))(+(+(incomplete-gamma-tail0-iz)(incomplete-gamma-tail0iz))(+(log-iz)(logiz))))))))(if(and(realpz)(pluspz))(realpart(ciz))(ciz))));; Array of values of the Bernoulli numbers. We only have enough for;; the evaluation of the psi function.(defconstantbern-values(make-array55:initial-contents'(1-1/21/60-1/3001/420-1/3005/660-691/273007/60-3617/510043867/7980-174611/3300854513/1380-236364091/273008553103/60-23749461029/87008615841276005/143220-7709321041217/51002577687858367/60-26315271553053477373/191919002929993913841559/60-261082718496449122051/1353001520097643918070802691/18060-27833269579301024235023/6900596451111593912163277961/2820-5609403368997817686249127547/464100495057205241079648212477525/660-801165718135489957347924991853/1590029149963634884862421418123812691/798)))(defunbern(k)(arefbern-valuesk))(defunpsi(z)"Digamma function defined by - %gamma + sum(1/k-1/(k+z-1), k, 1, inf) where %gamma is Euler's constant";; A&S 6.3.7: Reflection formula;;;; psi(1-z) = psi(z) + %pi*cot(%pi*z);;;; A&S 6.3.6: Recurrence formula;;;; psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z);;;; A&S 6.3.8: Asymptotic formula;;;; psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf);;;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence;; formula to increase the argument and then apply the asymptotic formula.(cond((=z1);; psi(1) = -%gamma(-(float+%gamma+(if(integerpz)0.0z))))((minusp(realpartz))(let((p(float-piz)))(flet((cot-pi(z);; cot(%pi*z), carefully. If z is an odd multiple;; of 1/2, cot is 0.(if(and(realpz)(=1/2(abs(-z(ftruncatez)))))(float0z)(/(tan(*pz))))))(-(psi(-1z))(*p(cot-piz))))))(t(let*((k(*2(1+(floor(*.41(-(log(epsilon(float(realpartz)))10)))))))(m0)(y(expt(+zk)2))(x0))(loopforifrom1upto(floork2)do(progn(incfm(+(/(+zii-1))(/(+zii-2))))(setfx(/(+x(/(bern(+k2(*-2i)))(-kii-2)))y))))(-(log(+zk))(/(*2(+zk)))xm)))))