"""Untapered Möbius approximation errors -- numberdb.org/T706.

Run with SageMath or compatible passagemath and numberdb:
    sage -pip install numberdb
    sage -python generate.py             # verify every stored entry; no write
    sage -python generate.py --publish   # fill the authorized draft, no offering
For a private draft configure NUMBERDB_API_KEY. No credentials belong here.
The campaign wrapper may provide the key on stdin with NUMBERDB_KEY_FROM_STDIN=1.

U_N = ||chi + sum_(n=1)^N mu(n) rho_n||^2, in L^2((0,infinity),dx),
rho_n(x)={1/(nx)}, chi=1_(0,1). N is a positive integer; default N=1,...,30.
Finite cotangent sums (BBLS, math/0306251v1, Proposition 89) evaluate the
Gram kernel at exact rational ratios. All transcendental operations and norm
arithmetic use Arb balls, at numberdb.bits(digits, losing=128) working bits.
No companion-table decimals are input. See derivation.md for the formulas,
checks.py for the separate integral route, and numerical-evidence.json for
measured enclosure accuracy. The gram_system implementation is reused from
T641's attached generate.py, SHA-256 7e155af1d4a4963f874e0c536051c1cf0b7afd01e21b93d64a674c7fab8592c8,
after inspecting its formula and checking every required kernel independently.
"""
import os
import sys
from functools import lru_cache
from math import gcd

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

MAX_N = 30
WORKING_GUARD = 128

def gram_system(limit, precision):
    """Fresh finite cotangent sums, Gram entries, and target inner products."""
    R = RealBallField(precision)
    C = (2*R.pi()).log() - R.euler_constant()

    @lru_cache(None)
    def vasyunin(p, q):
        return sum((R(QQ((m*p) % q)/QQ(q)) *
                    (R.pi()*R(QQ(m)/QQ(q))).cot()
                    for m in range(1, q)), R(0))

    @lru_cache(None)
    def autocorrelation(p, q):
        x = R(QQ(p)/QQ(q))
        return (((1-x)*x.log() + (1+x)*C)/2
                - R.pi()*(vasyunin(p,q)+vasyunin(q,p))/(2*q))

    G = [[R(0) for k in range(limit)] for j in range(limit)]
    for j in range(1, limit+1):
        for k in range(1, j+1):
            g = gcd(j,k)
            G[j-1][k-1] = G[k-1][j-1] = autocorrelation(j//g,k//g)/j
    b = [(R(n).log()+1-R.euler_constant())/n for n in range(1,limit+1)]
    return R, G, b


@lru_cache(None)
def computation(limit, digits):
    R, G, b = gram_system(limit, numberdb.bits(digits, losing=WORKING_GUARD))
    weights = [moebius(n) for n in range(1, limit+1)]
    total = R(1)
    values = []
    for i, w in enumerate(weights):
        total += 2*w*b[i] + w*w*G[i][i]
        total += 2*w*sum((weights[j]*G[i][j] for j in range(i)), R(0))
        assert total.is_finite() and total > 0
        assert total.accuracy() >= numberdb.bits(digits)+16
        values.append(total)
    return R, G, b, weights, values


class UntaperedErrors(numberdb.Generator):
    table = 'T706'
    parameters = ('N',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'integral_kernel.py',
             'argument-selection.md', 'derivation.md', 'numerical-evidence.json',
             'skill-manifest.json', 'research.md')

    def enumerate(self, limit=MAX_N):
        for N in range(1, limit+1):
            yield {'N': str(N)}

    def value(self, params, digits):
        N = int(params['N'])
        if N < 1:
            raise ValueError('N must be a positive integer')
        return computation(max(MAX_N, N), digits)[4][N-1]


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator = UntaperedErrors()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message=os.environ.get('NUMBERDB_EDIT_MESSAGE',
              'Recompute untapered Mobius errors for N=1,...,30')))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)