back to table · edit · history · where entries came from · files · download
4374 bytes, as of the version from 2026-10-29 08:32 (current). Recorded here, not run.
"""Exact finite determinant formulas used by numberdb.org/T604.
Cunha and Freitas, arXiv:2404.12114v1, Theorems 2.7, 2.8 and 3.2.
Adapted from the attached generators of NumberDB T595 and T603, downloaded
in this run and checked against the cited formulas. These compute afresh;
no rounded database values enter the computation.
"""
from functools import lru_cache
from math import factorial
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
WORKING_DIGITS = 150
@lru_cache(None)
def stirling_row(n):
"""Signed first-kind Stirling numbers from an integer recurrence."""
if n == 0:
return (ZZ(1),)
previous = stirling_row(n - 1)
return tuple((previous[k-1] if k else ZZ(0))
- (n-1)*(previous[k] if k < n else ZZ(0))
for k in range(n+1))
def harmonic(n):
return sum((QQ(1)/k for k in range(1, n+1)), QQ(0))
@lru_cache(None)
def even_correction(n):
# Coefficients of x product_{j=1}^{n-1}(x-(2j-1)^2).
v = [ZZ(0), ZZ(1)]
for j in range(1, n):
new = [ZZ(0)]*(len(v)+1)
for k, coefficient in enumerate(v):
new[k] -= (2*j-1)**2*coefficient
new[k+1] += coefficient
v = new
scale = QQ(1)/(ZZ(4)**(n-1)*factorial(2*n-1))
return sum((scale*v[i]*ZZ(2*n-1)**(2*i)/QQ(2*i)
* (harmonic(i-1)/2-harmonic(2*i-1))
for i in range(1, n+1)), QQ(0))
@lru_cache(None)
def sphere_logdet(d, working_digits=150):
if d < 1:
raise ValueError('The dimension must be a positive integer.')
C = ComplexBallField(numberdb.bits(working_digits))
R = RealBallField(numberdb.bits(working_digits))
if d == 1:
return 2*(2*R.pi()).log()
n = d//2
if d % 2:
a, b = stirling_row(2*n), stirling_row(2*n+1)
derivative = sum((C(-2*k).zetaderiv(1)*(a[2*k]+b[2*k+1])
for k in range(1, n+1)), C(0))
derivative *= QQ(2)/factorial(2*n)
derivative += (R(n)/R.pi()).log()
else:
a, b = stirling_row(2*n-1), stirling_row(2*n)
derivative = sum((C(1-2*k).zetaderiv(1)*(a[2*k-1]+b[2*k])
for k in range(1, n+1)), C(0))
derivative *= QQ(2)/factorial(2*n-1)
derivative += R(even_correction(n)) + R(2*n-1).log()
assert derivative.real().is_finite() and derivative.imag().is_finite() and derivative.imag().contains_zero()
return -derivative.real()
@lru_cache(None)
def stirling(n,k):
if n == 0:
return int(k == 0)
if k < 1 or k > n:
return 0
return stirling(n-1,k-1)-(n-1)*stirling(n-1,k)
def roots_polynomial(roots):
"""Exact coefficient list of product (X-r); no symbolic ring needed."""
c=[1]
for r in roots:
out=[0]*(len(c)+1)
for k,v in enumerate(c):
out[k]-=r*v
out[k+1]+=v
c=out
return c
@lru_cache(None)
def dirichlet_logdet(d, working=WORKING_DIGITS):
if not isinstance(d,int) or d<1:
raise ValueError('dimension must be a positive integer')
C=ComplexBallField(numberdb.bits(working))
R=RealBallField(numberdb.bits(working))
if d == 1:
return (2*R.pi()).log()
# The signed Stirling combination is valid for either parity.
zp=sum((QQ(stirling(d,k+1)+(-1)**(d-1+k)*stirling(d-1,k))
/QQ(factorial(d-1))*C(-k).zetaderiv(1)
for k in range(d)), C(0))
if d%2:
n=(d+1)//2
u=roots_polynomial([j*j for j in range(n-1)])
correction=-sum((QQ(2*u[i])/QQ(factorial(2*n-2))
*QQ((n-1)**(2*i+1))/QQ(2*i)
*(harmonic(i-1)/2-harmonic(2*i-1))
for i in range(1,n)),QQ(0))
else:
n=d//2
v=roots_polynomial([0]+[(2*j-1)**2 for j in range(1,n)])
correction=sum((QQ(v[i])/QQ(4**(n-1)*factorial(2*n-1))
*QQ((2*n-1)**(2*i))/QQ(4*i)
*(harmonic(i-1)/2-harmonic(2*i-1))
for i in range(1,n+1)),QQ(0))
result=-(zp+correction)
assert result.real().is_finite() and result.imag().is_finite()
assert result.real().accuracy() > 400 and result.imag().contains_zero()
return result.real()