back to table · edit · history · where entries came from · files · download
5468 bytes, as of the version from 2026-10-29 06:09. Recorded here, not run.
"""Torsional rigidity of classical planar domains -- numberdb.org/T602.
Run with SageMath (or passagemath with the named rings installed):
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # fill draft, NUMBERDB_API_KEY set
sage -python checks.py # independent identities and evidence
The 19 domains and two invertible scale conventions are justified in
argument-selection.md. Reference dimensions are radii, semiaxes, sides or legs
as specified there, never diameters. The Poisson normalization is -Delta u=1.
The engineering constant J is four times the quantity stored here.
All arithmetic uses Arb balls. The correction to the odd zeta(5) sum is
summed through 159 (80 terms) at the default 100 output digits. If M is the
next odd integer and c=pi*r, the positive omitted correction is at most
exp(-c*M)/(M^5*(1-exp(-2*c))) for tanh, with an additional factor
1/(1-exp(-c*M)) for coth. A symmetric ball of that radius is added before
the prefactors; enclosing a positive tail symmetrically is conservative.
WORKING_GUARD=168 bits gives 501 working bits for 100-digit output.
Measured minimum relative accuracy was 495 bits (over 149 decimal digits). The separate checks measure retained accuracy and repeat at 200
working decimal digits. The term count also grows for larger requested
precision, using c>=pi. No helper source controls this enumeration.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
WORKING_GUARD = 168
RECTANGLES = ((1,1),(4,3),(3,2),(5,3),(2,1),(5,2),(3,1),(4,1),(5,1),(6,1))
ELLIPSES = ((3,2),(2,1),(3,1),(4,1))
NAMED = ('disk','equilateral-triangle','right-isosceles-triangle','half-disk','slit-disk')
DOMAINS = NAMED + tuple('rectangle-%d-%d' % x for x in RECTANGLES) + tuple('ellipse-%d-%d' % x for x in ELLIPSES)
NORMALIZATIONS = ('unit-length','unit-area')
def odd_sum(field, ratio, coth=False, terms=80):
"""Enclose the complete odd tanh/coth sum, including its infinite tail."""
c = field.pi()*ratio
correction = field(0)
for k in range(1,2*terms,2):
q = (-c*k).exp()
correction += q/(field(k)**5*(1-q if coth else 1+q))
m = 2*terms+1
q = (-c*m).exp()
tail = q/(field(m)**5*(1-(-2*c).exp()))
if coth:
tail /= 1-q
correction += field(0,tail.upper())
base = field(31)/32*field(5).zeta()
return base + 2*correction if coth else base-2*correction
def area(domain, field):
if domain in ('disk','slit-disk'):
return field.pi()
if domain == 'half-disk':
return field.pi()/2
if domain == 'equilateral-triangle':
return field(3).sqrt()/4
if domain == 'right-isosceles-triangle':
return field(1)/2
family,p,q = domain.split('-')
r = field(QQ(p)/QQ(q))
if family == 'rectangle':
return r
if family == 'ellipse':
return field.pi()*r
raise ValueError(domain)
@lru_cache(maxsize=None)
def reference_value(domain, bits, terms):
field=RealBallField(bits)
pi=field.pi()
if domain == 'disk':
return pi/8
if domain == 'equilateral-triangle':
return field(3).sqrt()/320
if domain == 'half-disk':
return pi/8-1/pi
if domain == 'slit-disk':
return pi/4-16/(9*pi)
if domain == 'right-isosceles-triangle':
return field(1)/48-4/pi**5*odd_sum(field,field(1),True,terms)
family,p,q=domain.split('-')
r=field(QQ(p)/QQ(q))
if family == 'rectangle':
return r/12-16/pi**5*odd_sum(field,r,False,terms)
if family == 'ellipse':
return pi*r**3/(4*(r**2+1))
raise ValueError(domain)
def compute(domain, normalization, bits, terms=80):
value=reference_value(domain,bits,terms)
if normalization == 'unit-area':
value /= area(domain,RealBallField(bits))**2
elif normalization != 'unit-length':
raise ValueError(normalization)
assert value.is_finite() and value>0
return value
class TorsionalRigidity(numberdb.Generator):
table='T602'
files=('generate.py','checks.py','check.py','argument-selection.md',
'derivations.md','skill-manifest.json','oeis-A180310.txt',
'oeis-A180311.txt','oeis-A180314.txt')
parameters=('domain','normalization')
type='R'
digits=100
rigour='proven'
def enumerate(self):
for domain in DOMAINS:
for normalization in NORMALIZATIONS:
yield {'domain':domain,'normalization':normalization}
def value(self,params,digits):
terms=max(80,(digits+60)//2)
return compute(params['domain'],params['normalization'],
numberdb.bits(digits,losing=WORKING_GUARD),terms)
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isidentifier():
key=key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
generator=TorsionalRigidity()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
print(generator.publish(message='Fill 38 torsional rigidities with certified 100-digit enclosures'))
else:
report=generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)