program.py

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

1079 bytes, as of the version from 2026-11-10 20:59 (current). Recorded here, not run.

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.arith.misc import moebius
from math import gcd
from functools import lru_cache

def vasyunin_error(N, digits=100):
    N=ZZ(N)
    if N<1:raise ValueError('N must be positive')
    if N==1:return ZZ(1)
    R=RealBallField(numberdb.bits(digits, 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))
    def gram(j,k):
        d=gcd(j,k);p=j//d;q=k//d;x=R(QQ(p)/q)
        return ((1-x)*x.log()+(1+x)*C-R.pi()*(v(p,q)+v(q,p))/q)/(2*j)
    w=[QQ(moebius(n)) for n in range(1,N+1)]
    w[0]-=sum((w[n-1]/n for n in range(1,N+1)),QQ(0))
    return (1+2*sum(w[n-1]*(R(n).log()+1-R.euler_constant())/n
                    for n in range(1,N+1))
            +sum(w[j-1]*w[k-1]*gram(j,k)
                 for j in range(1,N+1) for k in range(1,N+1)))

vasyunin_error(30)