From 69542e883f75a0152d59583920a5cddce638d50e Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 1 Sep 2016 23:05:10 -0400 Subject: [PATCH 01/10] more accurate and faster lgamma, and use a more standard branch cut for real(x)<0 --- base/special/gamma.jl | 88 ++++++++++++++++++++++++++----------------- test/math.jl | 14 ++++++- 2 files changed, 66 insertions(+), 36 deletions(-) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 083ad14c30fdb..1112a1046bd02 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -33,49 +33,67 @@ Compute the logarithmic factorial of `x` lfact(x::Real) = (x<=1 ? zero(float(x)) : lgamma(x+one(x))) @vectorize_1arg Number lfact -const clg_coeff = [76.18009172947146, - -86.50532032941677, - 24.01409824083091, - -1.231739572450155, - 0.1208650973866179e-2, - -0.5395239384953e-5] - -function clgamma_lanczos(z) - const sqrt2pi = 2.5066282746310005 - - y = x = z - temp = x + 5.5 - zz = log(temp) - zz = zz * (x+0.5) - temp -= zz - ser = complex(1.000000000190015, 0) - for j=1:6 - y += 1.0 - zz = clg_coeff[j]/y - ser += zz - end - zz = sqrt2pi*ser / x - return log(zz) - temp -end - """ lgamma(x) Compute the logarithm of the absolute value of [`gamma`](:func:`gamma`) for [`Real`](:obj:`Real`) `x`, while for [`Complex`](:obj:`Complex`) `x` it computes the -logarithm of `gamma(x)`. +principal branch cut of the logarithm of `gamma(x)` (defined for negative `real(x)` +by analytic continuation from positive `real(x)`). """ -function lgamma(z::Complex) - if real(z) <= 0.5 - a = clgamma_lanczos(1-z) - b = log(sinpi(z)) - const logpi = 1.14472988584940017 - z = logpi - b - a - else - z = clgamma_lanczos(z) +function lgamma end + +# asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)| +@inline function lgamma_asymptotic(z::Complex{Float64}) + zinv = inv(z) + t = zinv*zinv + # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1)) + return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2 + zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03,7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04,8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03,6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) +end + +function lgamma(z::Complex{Float64}) + x = real(z) + y = imag(z) + yabs = abs(y) + if !isfinite(x) || !isfinite(y) # Inf or NaN + if isinf(x) && isfinite(y) + return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y)) + elseif isfinite(x) && isinf(y) + return Complex(-Inf, y) + else + return Complex(NaN, NaN) + end + elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y| + return lgamma_asymptotic(z) + elseif x < 0.1 # use reflection formula to transform to x > 0 + if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0 + return Complex(Inf, signbit(x) ? copysign(oftype(x, pi), -y) : -y) + end + # the 2pi * floor(...) stuff is to choose the correct branch cut for log(sinpi(z)) + return Complex(1.1447298858494001741434262, # log(pi) + copysign(6.2831853071795864769252842, y) # 2pi + * floor(0.5*x+0.25)) - + log(sinpi(z)) - lgamma(1-z) + elseif abs(x - 1) + yabs < 0.1 + # taylor series around zero at z=1 + # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] + w = Complex(x - 1, y) + return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02) + end + # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series + shift = log(z) + x += 1 + while x <= 7 + shift += log(Complex(x,y)) + x += 1 end - complex(real(z), angle_restrict_symm(imag(z))) + # note: it is faster to take the product of the log arguments, and *then* + # take the log to get `shift`, but this gives the wrong branch cut + return lgamma_asymptotic(Complex(x,y)) - shift end +lgamma{T<:Union{Integer,Rational}}(z::Complex{T}) = lgamma(float(z)) +lgamma{T<:Union{Float32,Float16}}(z::Complex{T}) = Complex{T}(lgamma(Complex{Float64}(z))) gamma(z::Complex) = exp(lgamma(z)) diff --git a/test/math.jl b/test/math.jl index 26a744b3211f0..a65ec6bb8aa61 100644 --- a/test/math.jl +++ b/test/math.jl @@ -554,7 +554,7 @@ for elty in (Float32, Float64) end @test lgamma(1.4+3.7im) ≈ -3.7094025330996841898 + 2.4568090502768651184im @test lgamma(1.4+3.7im) ≈ log(gamma(1.4+3.7im)) -@test lgamma(-4.2+0im) ≈ lgamma(-4.2)-pi*im +@test lgamma(-4.2+0im) ≈ lgamma(-4.2)-5pi*im @test factorial(3.0) == gamma(4.0) == factorial(3) for x in (3.2, 2+1im, 3//2, 3.2+0.1im) @test factorial(x) == gamma(1+x) @@ -654,6 +654,18 @@ for x in -10.2:0.3456:50 @test 1e-12 > err(digamma(x+0im), digamma(x)) end +# lgamma test cases (from Wolfram Alpha) +@test lgamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im +@test lgamma(3.099) ≅ lgamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796 +@test lgamma(1.15) ≅ lgamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554 +@test lgamma(0.89) ≅ lgamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963 +@test lgamma(0.91) ≅ lgamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167 +@test lgamma(0.01) ≅ lgamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419 +@test lgamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im +@test lgamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im +@test lgamma(-13.4+0im) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im +@test lgamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im + # digamma, trigamma, polygamma & zeta test cases (compared to Wolfram Alpha) @test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232 @test digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im From f2db1caef8381a2b7ade6ab25f201e95eab73cf6 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 12:04:26 -0400 Subject: [PATCH 02/10] more tests --- test/math.jl | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/test/math.jl b/test/math.jl index a65ec6bb8aa61..ee23e420e75d3 100644 --- a/test/math.jl +++ b/test/math.jl @@ -665,6 +665,23 @@ end @test lgamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im @test lgamma(-13.4+0im) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im @test lgamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im +@test lgamma(1+exp2(-20)) ≅ lgamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7 +@test lgamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im +@test lgamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im +@test lgamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im +@test lgamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im +@test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im + +# lgamma for non-finite arguments: +@test lgamma(Inf + 0im) === Inf + 0im +@test lgamma(Inf - 0.0im) === Inf - 0.0im +@test lgamma(Inf + 1im) === Inf + Inf*im +@test lgamma(Inf - 1im) === Inf - Inf*im +@test lgamma(-Inf + 0.0im) === -Inf - Inf*im +@test lgamma(-Inf - 0.0im) === -Inf + Inf*im +@test lgamma(Inf*im) === -Inf + Inf*im +@test lgamma(-Inf*im) === -Inf - Inf*im +@test lgamma(Inf + Inf*im) === lgamma(NaN + 0im) === lgamma(NaN*im) === NaN + NaN*im # digamma, trigamma, polygamma & zeta test cases (compared to Wolfram Alpha) @test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232 From e1d72fc298e0ecf4debe3169c061724de4d591f5 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 14:17:43 -0400 Subject: [PATCH 03/10] use trick from Hare (1997) to compute log(prod of shifts) rather than sum(logs of shifts) with the correct branch cut --- base/special/gamma.jl | 27 +++++++++++++++++++++++---- test/math.jl | 2 ++ 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 1112a1046bd02..9203b6dd1545b 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -52,6 +52,13 @@ function lgamma end zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03,7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04,8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03,6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) end +# Compute the logΓ(z) function using a combination of the asymptotic series, +# the Taylor series around z=1, the reflection formula, and the shift formula. +# Many details of these techniques are discussed in D. E. G. Hare, +# "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997), +# and similar techniques are used (in a somewhat different way) by the +# SciPy loggamma function. The key identities are also described +# at http://functions.wolfram.com/GammaBetaErf/LogGamma/ function lgamma(z::Complex{Float64}) x = real(z) y = imag(z) @@ -82,14 +89,26 @@ function lgamma(z::Complex{Float64}) return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02) end # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series - shift = log(z) + shiftprod = Complex(x,yabs) x += 1 + sb = false # == signbit(imag(shiftprod)) == signbit(yabs) + # To use log(product of shifts) rather than sum(logs of shifts), + # we need to keep track of the number of + to - sign flips in + # imag(shiftprod), as described in Hare (1997), proposition 2.2. + signflips = 0 while x <= 7 - shift += log(Complex(x,y)) + shiftprod *= Complex(x,yabs) + sb′ = signbit(imag(shiftprod)) + signflips += sb′ & (sb′ != sb) + sb = sb′ x += 1 end - # note: it is faster to take the product of the log arguments, and *then* - # take the log to get `shift`, but this gives the wrong branch cut + shift = log(shiftprod) + if signbit(y) # if y is negative, conjugate the shift + shift = Complex(real(shift), signflips*-6.2831853071795864769252842 - imag(shift)) + else + shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842) + end return lgamma_asymptotic(Complex(x,y)) - shift end lgamma{T<:Union{Integer,Rational}}(z::Complex{T}) = lgamma(float(z)) diff --git a/test/math.jl b/test/math.jl index ee23e420e75d3..1d5f250ad7658 100644 --- a/test/math.jl +++ b/test/math.jl @@ -671,6 +671,8 @@ end @test lgamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im @test lgamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im @test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im +@test lgamma(1+6.5im) ≅ conj(lgamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im +@test lgamma(1+1im) ≅ conj(lgamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im # lgamma for non-finite arguments: @test lgamma(Inf + 0im) === Inf + 0im From be2084cc95c34c83aa55b01f51120c57dc32ee8d Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 14:31:42 -0400 Subject: [PATCH 04/10] even more tests --- test/math.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/math.jl b/test/math.jl index 1d5f250ad7658..7dfc5a2658c31 100644 --- a/test/math.jl +++ b/test/math.jl @@ -663,7 +663,7 @@ end @test lgamma(0.01) ≅ lgamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419 @test lgamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im @test lgamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im -@test lgamma(-13.4+0im) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im +@test lgamma(-13.4+0.0im) ≅ conj(lgamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im @test lgamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im @test lgamma(1+exp2(-20)) ≅ lgamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7 @test lgamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im @@ -673,6 +673,10 @@ end @test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im @test lgamma(1+6.5im) ≅ conj(lgamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im @test lgamma(1+1im) ≅ conj(lgamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im +@test lgamma(-pi*10^7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im +@test lgamma(-pi*10^7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im +@test lgamma(-pi*10^14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im +@test lgamma(-pi*10^14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im # lgamma for non-finite arguments: @test lgamma(Inf + 0im) === Inf + 0im From d2aa66d666ee8e6d7def88d6c6e24292a1846195 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 15:26:02 -0400 Subject: [PATCH 05/10] whoops, use 1e14 and not 10^14 to avoid integer overflow on 32-bit Windows --- test/math.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/math.jl b/test/math.jl index 7dfc5a2658c31..ad476c686d343 100644 --- a/test/math.jl +++ b/test/math.jl @@ -673,10 +673,10 @@ end @test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im @test lgamma(1+6.5im) ≅ conj(lgamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im @test lgamma(1+1im) ≅ conj(lgamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im -@test lgamma(-pi*10^7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im -@test lgamma(-pi*10^7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im -@test lgamma(-pi*10^14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im -@test lgamma(-pi*10^14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im +@test lgamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im +@test lgamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im +@test lgamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im +@test lgamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im # lgamma for non-finite arguments: @test lgamma(Inf + 0im) === Inf + 0im From afbcdbf86eea2eca52b8d124f686f02a4b23b65c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 16:45:12 -0400 Subject: [PATCH 06/10] fix accuracy near zero at z=2 --- base/special/gamma.jl | 7 ++++++- test/math.jl | 2 ++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 9203b6dd1545b..2e171c4cf7548 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -53,7 +53,7 @@ function lgamma end end # Compute the logΓ(z) function using a combination of the asymptotic series, -# the Taylor series around z=1, the reflection formula, and the shift formula. +# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula. # Many details of these techniques are discussed in D. E. G. Hare, # "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997), # and similar techniques are used (in a somewhat different way) by the @@ -87,6 +87,11 @@ function lgamma(z::Complex{Float64}) # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] w = Complex(x - 1, y) return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02) + elseif abs(x - 2) + yabs < 0.1 + # taylor series around zero at z=2 + # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:15]] + w = Complex(x - 2, y) + return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01,-6.7352301053198095133246196e-02,2.0580808427784547879000897e-02,-7.3855510286739852662729527e-03,2.8905103307415232857531201e-03,-1.1927539117032609771139825e-03,5.0966952474304242233558822e-04,-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05,-9.4394882752683959040399766e-06,4.374866789907487804177898e-06,-2.0392157538013662367869141e-06) end # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series shiftprod = Complex(x,yabs) diff --git a/test/math.jl b/test/math.jl index ad476c686d343..665dd08021b97 100644 --- a/test/math.jl +++ b/test/math.jl @@ -677,6 +677,8 @@ end @test lgamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im @test lgamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im @test lgamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im +@test lgamma(2.05 + 0.03im) ≅ conj(lgamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im +@test lgamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im # lgamma for non-finite arguments: @test lgamma(Inf + 0im) === Inf + 0im From cb8caba3113ec604f5b7947c100172a548695df4 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 16:45:57 -0400 Subject: [PATCH 07/10] update manual for lgamma --- doc/stdlib/math.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/stdlib/math.rst b/doc/stdlib/math.rst index 206f9e38733f7..26a634a0b752c 100644 --- a/doc/stdlib/math.rst +++ b/doc/stdlib/math.rst @@ -1364,7 +1364,7 @@ Mathematical Functions .. Docstring generated from Julia source - Compute the logarithm of the absolute value of :func:`gamma` for :obj:`Real` ``x``\ , while for :obj:`Complex` ``x`` it computes the logarithm of ``gamma(x)``\ . + Compute the logarithm of the absolute value of :func:`gamma` for :obj:`Real` ``x``\ , while for :obj:`Complex` ``x`` it computes the principal branch cut of the logarithm of ``gamma(x)`` (defined for negative ``real(x)`` by analytic continuation from positive ``real(x)``\ ). .. function:: lfact(x) From 7f5cede131d01f95497e7007ae1b7de92afda52b Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 16:48:01 -0400 Subject: [PATCH 08/10] news for lgamma changes --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 530777425f5b8..ca6930687b372 100644 --- a/NEWS.md +++ b/NEWS.md @@ -24,6 +24,10 @@ This section lists changes that do not have deprecation warnings. * Keyword arguments are processed left-to-right: if the same keyword is specified more than once, the rightmost occurrence takes precedence ([#17785]). + * The `lgamma(z)` function now uses a different (more standard) branch cut + for `real(z) < 0`, which differs from `log(gamma(z))` by multiples of 2π + in the imaginary part ([#18330]). + Library improvements -------------------- @@ -637,3 +641,4 @@ Language tooling improvements [#17546]: https://github.com/JuliaLang/julia/issues/17546 [#17668]: https://github.com/JuliaLang/julia/issues/17668 [#17785]: https://github.com/JuliaLang/julia/issues/17785 +[#18330]: https://github.com/JuliaLang/julia/issues/18330 From 0debff52269895ec23792bce16c58cb94dcedfde Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 2 Sep 2016 16:53:55 -0400 Subject: [PATCH 09/10] can use a lower-degree Taylor series around z=2 because the coefficients decrease faster --- base/special/gamma.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 2e171c4cf7548..83df0044853e8 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -89,9 +89,9 @@ function lgamma(z::Complex{Float64}) return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02) elseif abs(x - 2) + yabs < 0.1 # taylor series around zero at z=2 - # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:15]] + # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]] w = Complex(x - 2, y) - return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01,-6.7352301053198095133246196e-02,2.0580808427784547879000897e-02,-7.3855510286739852662729527e-03,2.8905103307415232857531201e-03,-1.1927539117032609771139825e-03,5.0966952474304242233558822e-04,-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05,-9.4394882752683959040399766e-06,4.374866789907487804177898e-06,-2.0392157538013662367869141e-06) + return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01,-6.7352301053198095133246196e-02,2.0580808427784547879000897e-02,-7.3855510286739852662729527e-03,2.8905103307415232857531201e-03,-1.1927539117032609771139825e-03,5.0966952474304242233558822e-04,-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) end # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series shiftprod = Complex(x,yabs) From 6f6ff37dce39af9ee6bd5067e577ae7cf355feee Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Sat, 3 Sep 2016 15:48:42 -0400 Subject: [PATCH 10/10] linewrap poly coefs --- base/special/gamma.jl | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 83df0044853e8..15585f2339644 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -49,7 +49,10 @@ function lgamma end t = zinv*zinv # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1)) return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2 - zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03,7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04,8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03,6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) + zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03, + 7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04, + 8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03, + 6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) end # Compute the logΓ(z) function using a combination of the asymptotic series, @@ -86,12 +89,24 @@ function lgamma(z::Complex{Float64}) # taylor series around zero at z=1 # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] w = Complex(x - 1, y) - return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02) + return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01, + -4.0068563438653142846657956e-01,2.705808084277845478790009e-01, + -2.0738555102867398526627303e-01,1.6955717699740818995241986e-01, + -1.4404989676884611811997107e-01,1.2550966952474304242233559e-01, + -1.1133426586956469049087244e-01,1.000994575127818085337147e-01, + -9.0954017145829042232609344e-02,8.3353840546109004024886499e-02, + -7.6932516411352191472827157e-02,7.1432946295361336059232779e-02, + -6.6668705882420468032903454e-02) elseif abs(x - 2) + yabs < 0.1 # taylor series around zero at z=2 # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]] w = Complex(x - 2, y) - return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01,-6.7352301053198095133246196e-02,2.0580808427784547879000897e-02,-7.3855510286739852662729527e-03,2.8905103307415232857531201e-03,-1.1927539117032609771139825e-03,5.0966952474304242233558822e-04,-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) + return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01, + -6.7352301053198095133246196e-02,2.0580808427784547879000897e-02, + -7.3855510286739852662729527e-03,2.8905103307415232857531201e-03, + -1.1927539117032609771139825e-03,5.0966952474304242233558822e-04, + -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05, + -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) end # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series shiftprod = Complex(x,yabs)