program.py

back to table · edit · history · where entries came from · files · download

945 bytes, as of the version from 2026-11-10 22:01 (current). Recorded here, not run.

import numberdb.sage as numberdb
from functools import lru_cache
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.arith.misc import moebius
R = RealBallField(numberdb.bits(100, losing=128))
C = (2*R.pi()).log() - R.euler_constant()
@lru_cache(None)
def V(p, q):
    return sum((R(QQ((m*p) % q)/q)*(R.pi()*R(QQ(m)/q)).cot()
                for m in range(1,q)), R(0))
@lru_cache(None)
def A(r):
    p, q = int(r.numerator()), int(r.denominator())
    x = R(r)
    return ((1-x)*x.log()+(1+x)*C)/2 - R.pi()*(V(p,q)+V(q,p))/(2*q)
def U(N):
    assert N >= 1
    w = [moebius(n) for n in range(1,N+1)]
    return 1 + 2*sum((w[n-1]*(R(n).log()+1-R.euler_constant())/n
                     for n in range(1,N+1)), R(0)) + sum(
        (w[j-1]*w[k-1]*A(QQ(j)/k)/j
         for j in range(1,N+1) for k in range(1,N+1)), R(0))
# The stored prefix is every N=1,...,30; U also accepts larger N.
print(U(10))