From d198ce8c5fabb3ea526e2db287053d73601dff8c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 20 Apr 2016 18:03:39 -0400 Subject: [PATCH] fixes to generalized zeta function zeta(s, z) (#15945) * fixes to generalized zeta function zeta(s, z) * re-run doc/genstdlib.jl --- base/docs/helpdb/Base.jl | 8 --- base/special/gamma.jl | 115 +++++++++++++++++++++++++++------------ doc/stdlib/math.rst | 2 +- test/math.jl | 71 +++++++++++++++--------- 4 files changed, 126 insertions(+), 70 deletions(-) diff --git a/base/docs/helpdb/Base.jl b/base/docs/helpdb/Base.jl index 487488f77dd59..83eeb1b5c9fc3 100644 --- a/base/docs/helpdb/Base.jl +++ b/base/docs/helpdb/Base.jl @@ -9579,14 +9579,6 @@ Riemann zeta function ``\\zeta(s)``. """ zeta(s) -""" - zeta(s, z) - -Hurwitz zeta function ``\\zeta(s, z)``. (This is equivalent to the Riemann zeta function -``\\zeta(s)`` for the case of `z=1`.) -""" -zeta(s,z) - """ A_mul_Bt(A, B) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index 67a8741bf6f2f..6cdf38412c569 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -220,15 +220,21 @@ macro pg_horner(x, m, p...) :(($me + 1) * ($(p[1]) + $xe * $ex)) end -# compute inv(oftype(x, y)) efficiently, choosing the correct branch cut -inv_oftype(x::Complex, y::Complex) = oftype(x, inv(y)) -function inv_oftype(x::Complex, y::Real) - yi = inv(y) # using real arithmetic for efficiency - oftype(x, Complex(yi, -zero(yi))) # get correct sign of zero! +# compute oftype(x, y)^p efficiently, choosing the correct branch cut +pow_oftype(x, y, p) = oftype(x, y)^p +pow_oftype(x::Complex, y::Real, p::Complex) = oftype(x, y^p) +function pow_oftype(x::Complex, y::Real, p::Real) + if p >= 0 + # note: this will never be called for y < 0, + # which would throw an error for non-integer p here + return oftype(x, y^p) + else + yp = y^-p # use real power for efficiency + return oftype(x, Complex(yp, -zero(yp))) # get correct sign of zero! + end end -inv_oftype(x::Real, y::Real) = oftype(x, inv(y)) -# Hurwitz zeta function, which is related to polygamma +# Generalized zeta function, which is related to polygamma # (at least for integer m > 0 and real(z) > 0) by: # polygamma(m, z) = (-1)^(m+1) * gamma(m+1) * zeta(m+1, z). # Our algorithm for the polygamma is just the m-th derivative @@ -238,15 +244,27 @@ inv_oftype(x::Real, y::Real) = oftype(x, inv(y)) # So identifying the (something) with the -zeta function, we get # the zeta function for free and might as well export it, especially # since this is a common generalization of the Riemann zeta function -# (which Julia already exports). +# (which Julia already exports). Note that this geneneralization +# is equivalent to Mathematica's Zeta[s,z], and is equivalent to the +# Hurwitz zeta function for real(z) > 0. + +""" + zeta(s, z) + +Generalized zeta function ``\\zeta(s, z)``, defined +by the sum ``\\sum_{k=0}^\\infty ((k+z)^2)^{-s/2}``, where +any term with ``k+z=0`` is excluded. For ``\\Re z > 0``, +this definition is equivalent to the Hurwitz zeta function +``\\sum_{k=0}^\\infty (k+z)^{-s}``. For ``z=1``, it yields +the Riemann zeta function ``\\zeta(s)``. +""" +zeta(s,z) + function zeta(s::Union{Int,Float64,Complex{Float64}}, z::Union{Float64,Complex{Float64}}) ζ = zero(promote_type(typeof(s), typeof(z))) - # like sqrt, require complex inputs to get complex outputs - !isa(s,Integer) && isa(ζ, Real) && z < 0 && throw(DomainError()) - - z == 1 && return oftype(ζ, zeta(s)) + (z == 1 || z == 0) && return oftype(ζ, zeta(s)) s == 2 && return oftype(ζ, trigamma(z)) x = real(z) @@ -263,9 +281,13 @@ function zeta(s::Union{Int,Float64,Complex{Float64}}, end throw(DomainError()) # nothing clever to return end - - # We need a different algorithm for the real(s) < 1 domain - real(s) < 1 && throw(ArgumentError("order $s < 1 is not implemented (issue #7228)")) + if isnan(x) + if imag(z)==0 && imag(s)==0 + return oftype(ζ, x) + else + return oftype(ζ, Complex(x,x)) + end + end m = s - 1 @@ -275,32 +297,55 @@ function zeta(s::Union{Int,Float64,Complex{Float64}}, # Note: we multiply by -(-1)^m m! in polygamma below, so this factor is # pulled out of all of our derivatives. - isnan(x) && return oftype(ζ, imag(z)==0 && isa(s,Int) ? x : Complex(x,x)) - - cutoff = 7 + real(m) + imag(m) # TODO: this cutoff is too conservative? + cutoff = 7 + real(m) + abs(imag(m)) # TODO: this cutoff is too conservative? if x < cutoff # shift using recurrence formula xf = floor(x) - if x <= 0 && xf == z - if isa(s, Int) - iseven(s) && return oftype(ζ, Inf) - x == 0 && return oftype(ζ, inv(x)) - end - throw(DomainError()) # or return NaN? - end nx = Int(xf) n = ceil(Int,cutoff - nx) - ζ += inv_oftype(ζ, z)^s - for ν = -nx:-1:1 - ζₒ= ζ - ζ += inv_oftype(ζ, z + ν)^s - ζ == ζₒ && break # prevent long loop for large -x > 0 - # FIXME: still slow for small m, large Im(z) + minus_s = -s + if nx < 0 # x < 0 + # need to use (-z)^(-s) recurrence to be correct for real z < 0 + # [the general form of the recurrence term is (z^2)^(-s/2)] + minus_z = -z + ζ += pow_oftype(ζ, minus_z, minus_s) # ν = 0 term + if xf != z + ζ += pow_oftype(ζ, z - nx, minus_s) # real(z - nx) > 0, so use correct branch cut + # otherwise, if xf==z, then the definition skips this term + end + # do loop in different order, depending on the sign of s, + # so that we are looping from largest to smallest summands and + # can halt the loop early if possible; see issue #15946 + # FIXME: still slow for small m, large Im(z) + if real(s) > 0 + for ν in -nx-1:-1:1 + ζₒ= ζ + ζ += pow_oftype(ζ, minus_z - ν, minus_s) + ζ == ζₒ && break # prevent long loop for large -x > 0 + end + else + for ν in 1:-nx-1 + ζₒ= ζ + ζ += pow_oftype(ζ, minus_z - ν, minus_s) + ζ == ζₒ && break # prevent long loop for large -x > 0 + end + end + else # x ≥ 0 && z != 0 + ζ += pow_oftype(ζ, z, minus_s) end - for ν = max(1,1-nx):n-1 - ζₒ= ζ - ζ += inv_oftype(ζ, z + ν)^s - ζ == ζₒ && break # prevent long loop for large m + # loop order depends on sign of s, as above + if real(s) > 0 + for ν in max(1,1-nx):n-1 + ζₒ= ζ + ζ += pow_oftype(ζ, z + ν, minus_s) + ζ == ζₒ && break # prevent long loop for large m + end + else + for ν in n-1:-1:max(1,1-nx) + ζₒ= ζ + ζ += pow_oftype(ζ, z + ν, minus_s) + ζ == ζₒ && break # prevent long loop for large m + end end z += n end diff --git a/doc/stdlib/math.rst b/doc/stdlib/math.rst index a6353a50d9508..4d3be3e2b70d7 100644 --- a/doc/stdlib/math.rst +++ b/doc/stdlib/math.rst @@ -1674,7 +1674,7 @@ Mathematical Functions .. Docstring generated from Julia source - Hurwitz zeta function :math:`\zeta(s, z)`\ . (This is equivalent to the Riemann zeta function :math:`\zeta(s)` for the case of ``z=1``\ .) + Generalized zeta function :math:`\zeta(s, z)`\ , defined by the sum :math:`\sum_{k=0}^\infty ((k+z)^2)^{-s/2}`\ , where any term with :math:`k+z=0` is excluded. For :math:`\Re z > 0`\ , this definition is equivalent to the Hurwitz zeta function :math:`\sum_{k=0}^\infty (k+z)^{-s}`\ . For :math:`z=1`\ , it yields the Riemann zeta function :math:`\zeta(s)`\ . .. function:: ndigits(n, b) diff --git a/test/math.jl b/test/math.jl index 1c2c89e9d45c1..51bcd915edddb 100644 --- a/test/math.jl +++ b/test/math.jl @@ -615,52 +615,71 @@ end @test set_zero_subnormals(false) @test !get_zero_subnormals() -# useful test functions for relative error -err(z, x) = abs(z - x) / abs(x) +# useful test functions for relative error, which differ from isapprox +# in that errc separately looks at the real and imaginay parts +err(z, x) = z == x ? 0.0 : abs(z - x) / abs(x) errc(z, x) = max(err(real(z),real(x)), err(imag(z),imag(x))) +≅(a,b) = errc(a,b) ≤ 1e-13 for x in -10.2:0.3456:50 @test 1e-12 > err(digamma(x+0im), digamma(x)) end # digamma, trigamma, polygamma & zeta test cases (compared to Wolfram Alpha) -@test 1e-13 > err(digamma(7+0im), 1.872784335098467139393487909917597568957840664060076401194232) -@test 1e-13 > errc(digamma(7im), 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im) -@test 1e-13 > errc(digamma(-3.2+0.1im), 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im) -@test 1e-13 > err(trigamma(8+0im), 0.133137014694031425134546685920401606452509991909746283540546) -@test 1e-13 > errc(trigamma(8im), -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im) -@test 1e-13 > errc(trigamma(-3.2+0.1im), 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im) -@test 1e-13 > err(polygamma(2, 8.1+0im), -0.01723882695611191078960494454602091934457319791968308929600) -@test 1e-13 > errc(polygamma(30, 8.1+2im), -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im) -@test 1e-13 > errc(polygamma(3, 2.1+1im), 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im) +@test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232 +@test digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im +@test digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im +@test trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546 +@test trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im +@test trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im +@test polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600 +@test polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im +@test polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im @test 1e-11 > err(polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im) -@test 1e-13 > err(polygamma(13, 5.2 - 2im), 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im) +@test polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im @test 1e-11 > err(polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291) -@test 1e-13 > errc(zeta(4.1+0.3im, -3.2+0.1im), -461.95403678374488506025596495576748255121001107881278765917+926.02552636148651929560277856510991293536052745360005500774im) -@test 1e-13 > errc(zeta(4.1+0.3im, 3.2+0.1im), 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im) -@test 1e-13 > errc(zeta(4.1, 3.2+0.1im),0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im) -@test 1e-12 > errc(zeta(1.0001, -4.5e2+3.2im), 9993.89099199843392251301993718413132850540848778561412270571-3.13257480938495907945892330398176989805350557816701044268548im) -@test_throws DomainError zeta(3.1,-4.2) -@test 1e-13 > errc(zeta(3.1,-4.2+0im), -138.06320182025311080661516120845508778572835942189570145952+45.586579397698817209431034568162819207622092308850063038062im) +@test zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im +@test zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im +@test zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im +@test 1e-12 > errc(zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im) +@test zeta(3.1,-4.2) ≅ zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295 @test 1e-15 > errc(zeta(3.1+0im,-4.2), zeta(3.1,-4.2+0im)) -@test 1e-13 > errc(zeta(3.1,4.2), 0.029938344862645948405021260567725078588893266227472565010234) -@test 1e-13 > err(zeta(27, 3.1), 5.413318813037879056337862215066960774064332961282599376e-14) -@test 1e-13 > err(zeta(27, 2), 7.4507117898354294919810041706041194547190318825658299932e-9) -@test 1e-12 > err(zeta(27, -105.3), -1.311372652244914148556295810515903234635727465138859603e14) +@test zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234 +@test zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14 +@test zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9 +@test 1e-12 > err(zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14) @test polygamma(4, -3.1+Inf*im) == polygamma(4, 3.1+Inf*im) == 0 @test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0) -@test zeta(4, +0.0) == Inf == zeta(4, -0.0) -@test zeta(5, +0.0) == Inf == -zeta(5, -0.0) +@test zeta(4, +0.0) == zeta(4, -0.0) ≅ pi^4 / 90 +@test zeta(5, +0.0) == zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974 @test zeta(Inf, 1.) == 1 @test zeta(Inf, 2.) == 0 @test isnan(zeta(NaN, 1.)) @test isa([digamma(x) for x in [1.0]], Vector{Float64}) @test isa([trigamma(x) for x in [1.0]], Vector{Float64}) @test isa([polygamma(3,x) for x in [1.0]], Vector{Float64}) -@test 1e-13 > errc(zeta(2 + 1im, -1.1), zeta(2 + 1im, -1.1+0im)) -@test 1e-13 > errc(zeta(2 + 1im, -1.1), -1525.8095173321060982383023516086563741006869909580583246557 + 1719.4753293650912305811325486980742946107143330321249869576im) +@test zeta(2 + 1im, -1.1) ≅ zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im @test_approx_eq polygamma(3,5) polygamma(3,5.) +@test zeta(-3.0, 7.0) ≅ -52919/120 +@test zeta(-3.0, -7.0) ≅ 94081/120 +@test zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719 +@test zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576 +@test zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748 +@test zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161 +@test zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im +@test zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im +@test zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im +@test zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im +@test zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im +@test zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im +@test zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im +@test zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6 +@test zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078 +@test zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977 +@test zeta(-3.2+0.1im, 0.0) ≅ zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im +@test zeta(-3.2, 0.0) ≅ zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378 + @test @evalpoly(2,3,4,5,6) == 3+2*(4+2*(5+2*6)) == @evalpoly(2+0im,3,4,5,6) @test let evalcounts=0 @evalpoly(begin