back to table · edit · history · where entries came from · files · download
3283 bytes, as of the version from 2026-10-29 07:39 (current). Recorded here, not run.
"""Independent continuation of the defining Dirichlet hemisphere spectrum.
This module imports no generator code or combinatorial coefficients. Put
x=k+a, a=(d-1)/2, and expand binomial(k+d-2,d-1)=sum c_j x^j.
For k>=K the binomial expansion gives
Z_tail(s)=sum_j c_j sum_r (s)_r/r! a^(2r) zeta(2s+2r-j,K+a).
The r=0 derivative is 2*zeta'(-j,K+a). For r>=1 its derivative is
a^(2r)/r * zeta(2r-j,K+a),
except at 2r-j=1, where the Laurent expansion gives
a^(2r)/r * (H_(r-1)/2 - psi(K+a)).
The finite k<K part differentiates to -m_k log(k(k+d-1)).
For N>=floor(d/2), q=K+a and t=(a/q)^2<1, the omitted derivative has
absolute value at most
t^(N+1)/((N+1)(1-t)) * sum_j |c_j| q^j (1+q/(2N+1-j)).
This follows from the integral bound zeta(p,q)<=q^(-p)(1+q/(p-1)),
p>1, and a geometric sum. We include this error in the returned Arb ball.
All arithmetic and Hurwitz-zeta series evaluations are enclosing arithmetic.
"""
from math import comb, factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
def continuation(d, working=240, K=40, tail_digits=170, boundary='D'):
if d<1 or boundary not in ('D','N'):
raise ValueError('invalid dimension or boundary condition')
a=QQ(d-1)/2
q=K+a
# Expand the defining multiplicity, without the generator's Stirling or
# central factorial recurrences. For N use binomial(k+d-1,d-1).
c=[QQ(1)]
for j in range(d-1):
root=QQ(j+(boundary=='N'))-a
out=[QQ(0)]*(len(c)+1)
for i,ci in enumerate(c):
out[i]+=root*ci
out[i+1]+=ci
c=out
c=[ci/QQ(factorial(d-1)) for ci in c]
# A polynomial identity control verifies every defining multiplicity.
for k in range(1,d+2):
assert sum(ci*(k+a)**j for j,ci in enumerate(c))==comb(d+k-2+(boundary=='N'),d-1)
t=(a/q)**2
N=max(1,d//2)
def bound(N):
return t**(N+1)/((N+1)*(1-t))*sum(
abs(ci)*q**j*(1+q/QQ(2*N+1-j)) for j,ci in enumerate(c))
E=bound(N)
while E>QQ(1)/10**tail_digits:
N+=1
E=bound(N)
C=ComplexBallField(numberdb.bits(working));R=RealBallField(numberdb.bits(working))
P=PolynomialRing(C,'s');s=P.gen()
deriv=-sum((R(comb(d+k-2+(boundary=='N'),d-1))*R(k*(k+d-1)).log()
for k in range(1,K)),R(0))
zetas={}
psi=C(q).psi().real()
for j,cj in enumerate(c):
if not cj:
continue
val=2*(s-j)._zeta_series(2,a=C(q))[1].real()
if a:
for r in range(1,N+1):
p=2*r-j
if p==1:
h=sum((QQ(1)/QQ(v) for v in range(1,r)),QQ(0))
z=R(h/2)-psi
else:
if p not in zetas:
zetas[p]=C(p).zeta(C(q)).real()
z=zetas[p]
val+=R(a**(2*r)/r)*z
deriv+=R(cj)*val
result=(-deriv).add_error(R(E))
assert result.is_finite() and result.rad_as_ball()<R(10)**(-160)
return result, {'K':K,'N':N,'tail_bound':str(R(E)),
'relative_accuracy_bits':int(result.accuracy())}