"""Multiplicative fractional-part autocorrelation -- numberdb.org/T640.

Run with SageMath (or compatible passagemath):
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored entries
    sage -python checks.py                   # independent offline checks
    sage -python generate.py --fill          # write draft entries, API key required

The --fill option sends entries to an existing table; it does not offer or
publish the draft. Set NUMBERDB_API_KEY in the environment for authenticated
access. The campaign wrapper may instead pass a key on stdin with
NUMBERDB_KEY_FROM_STDIN=1 and request filling with NUMBERDB_PUBLISH=1.

BBLS, arXiv:math/0306251, Proposition 89, evaluates A(p/q) using two finite
Vasyunin sums. All input fractions and modular residues are exact. Arb encloses
pi, Euler's constant, logarithms and cotangents and propagates their errors.
The 128-bit working guard was checked over every selected rational; see the
attached numerical evidence for measured accuracy. No stored decimal is used
as a computational input. The default bound 30 comes from the shared finite
Nyman--Beurling approximation problem; argument-selection.md explains it.
"""
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

WORKING_GUARD = 128
DEFAULT_BOUND = 30


@lru_cache(maxsize=None)
def vasyunin(p, q, bits):
    """V(p,q) = sum {mp/q} cot(pi*m/q); the empty sum is zero."""
    R = RealBallField(bits)
    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))


def autocorrelation(lam, digits=100):
    lam = QQ(lam)
    if lam <= 0:
        raise ValueError('The dilation must be positive.')
    p, q = int(lam.numerator()), int(lam.denominator())
    bits = numberdb.bits(digits, losing=WORKING_GUARD)
    R = RealBallField(bits)
    x = R(lam)
    C = (2*R.pi()).log() - R.euler_constant()
    return ((1-x)*x.log()+(1+x)*C)/2 - R.pi()*(
        vasyunin(p,q,bits)+vasyunin(q,p,bits))/(2*q)


class FractionalPartAutocorrelation(numberdb.Generator):
    table = 'T640'
    parameters = ('lambda',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md', 'derivation.md',
             'numerical-evidence.json', 'skill-evidence.json')

    def enumerate(self, bound=DEFAULT_BOUND):
        for height in range(1, bound+1):
            for p in range(1, height+1):
                for q in range(1, height+1):
                    if max(p,q) == height and gcd(p,q) == 1:
                        yield {'lambda': str(QQ(p)/QQ(q))}

    def value(self, params, digits):
        value = autocorrelation(params['lambda'], digits)
        if str(params['lambda']) == '1':
            return {'number': value, 'comment': r'The squared norm of the fractional-part function in $L^2((0,\infty),t^{-2}dt)$ is $A(1)=\log(2\pi)-\gamma$.'}
        return value


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        token = sys.stdin.read().strip()
        if '=' in token and token.split('=',1)[0].isupper():
            token = token.split('=',1)[1].strip().strip('\'"')
        os.environ['NUMBERDB_API_KEY'] = token
    generator = FractionalPartAutocorrelation()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        outcome = generator.publish(message='Add 555 certified autocorrelations with independent period-integral checks')
        print(outcome)
    report = generator.verify(sample=None)
    print(report)
    sys.exit(0 if report.ok else 1)