program.py

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

904 bytes, as of the version from 2026-10-31 22:58 (current). Recorded here, not run.

import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.arith.misc import moebius

R = RealBallField(numberdb.bits(100, losing=128))
C = (2*R.pi()).log() - R.euler_constant()
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))
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 E(N):
    assert N >= 2
    w = [R(0) if n == N else R(moebius(n))*(1-R(n).log()/R(N).log())
         for n in range(1, N+1)]
    b = [(R(n).log()+1-R.euler_constant())/n for n in range(1, N+1)]
    return 1 + 2*sum((w[i]*b[i] for i in range(N)), 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))

print(E(10))