"""Generalized log-sine integrals -- numberdb.org/T697.

Run with SageMath (or passagemath-symbolics):
    sage -pip install numberdb mpmath   # once
    sage -python generate.py            # verify every stored entry
    sage -python generate.py --submit   # fill the existing draft, with API key
    sage -python checks.py              # independent numerical checks

The Euler sine product gives log(2 sin(v/2)) = log(v) -
sum(zeta(2m)/m * (v/(2*pi))**(2m), m>=1).  Integrate its integer powers
coefficient by coefficient, including the logarithmic endpoint analytically.
certification.md proves the uniform tail bound used here. No quadrature,
fitted radius, or precomputed numerical values enter this generator.

Selection: 2<=n<=10, 0<=k<=n-2, t=1/3,1/2,2/3,1; also (4,0,5/6).
argument-selection.md explains every parameter and conditional slice.
Instruction provenance: skill-manifest.json, archived in this run at
agents/runs/20261109T220207Z-build-skills/manifest.json. Both skill copies:
886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b.
"""
import os
import sys
from functools import lru_cache
from math import comb, factorial
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

WORKING_GUARD_BITS = 192
TAIL_GUARD_DIGITS = 60
MAX_ORDER = 10
ANGLES = ('1/3', '1/2', '2/3', '1')

@lru_cache(maxsize=16)
def workspace(t, digits):
    R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD_BITS))
    P = PolynomialRing(R, 'x')
    N = numberdb.bits(digits + TAIL_GUARD_DIGITS, losing=0)
    theta = R.pi() * R(QQ(t))
    q = (R(QQ(t)) / 2)**2
    A = theta.log()
    h = P([R(0)] + [-R(2*m).zeta() * q**m / m for m in range(1, N+1)])
    powers = [P.one()]
    for j in range(1, MAX_ORDER):
        powers.append((powers[-1] * h).truncate(N+1))
    return R, N, theta, q, A, powers

@lru_cache(maxsize=128)
def integrals(t, digits, k):
    R, N, theta, q, A, powers = workspace(t, digits)
    # J[d][m] = integral_0^1 x**(k+2*m) (A+log(x))**d dx.
    J = [[R(1)/(k+2*m+1) for m in range(N+1)]]
    for d in range(1, MAX_ORDER):
        J.append([(A**d - d*J[-1][m])/(k+2*m+1) for m in range(N+1)])
    return J

def evaluate(n, k, t, digits=100, exact_zero=True):
    n, k = int(n), int(k)
    assert 1 <= n <= MAX_ORDER and 0 <= k < n and 0 < QQ(t) <= 1
    if exact_zero and (n, k, str(QQ(t))) == (2, 0, '1'):
        # Integral_0^pi log(2 sin(v/2)) dv=0, from the Gamma ratio.
        return ZZ(0)
    R, N, theta, q, A, powers = workspace(str(QQ(t)), digits)
    J = integrals(str(QQ(t)), digits, k)
    p = n-1-k
    total = R(0)
    for j in range(p+1):
        total += comb(p,j)*sum((c*J[p-j][m] for m,c in enumerate(powers[j].list())), R(0))
    # Positivity of zeta(2m)/m, bounded by 2/m, gives the omitted tails.
    C = 2*R(2).log()
    error = R(0)
    for j in range(1,p+1):
        d = p-j
        absolute_moment = sum((R(comb(d,l)*factorial(l))*abs(A)**(d-l)/R(k+1)**(l+1)
                              for l in range(d+1)),R(0))
        error += comb(p,j)*C**j * (2*q)**(N+1)*absolute_moment
    answer = -theta**(k+1)*total
    answer = answer.add_error((theta**(k+1)*error).upper())
    assert answer.is_finite()
    if exact_zero:
        assert answer.accuracy() >= numberdb.bits(digits+20, losing=0)
    return answer

class LogSine(numberdb.Generator):
    table = 'T697'
    parameters = ('n','k','t')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py','checks.py','argument-selection.md','certification.md',
             'skill-manifest.json','numerical-evidence.json')

    def enumerate(self):
        for n in range(2,MAX_ORDER+1):
            for k in range(n-1):
                angles = ('1/3','1/2','2/3','5/6','1') if (n,k)==(4,0) else ANGLES
                for t in angles:
                    yield {'n':str(n),'k':str(k),'t':t}

    def value(self,params,digits):
        v=evaluate(int(params['n']),int(params['k']),params['t'],digits)
        if params=={'n':'2','k':'0','t':'1'}:
            return {'number':v,'comment':r'The integral is exactly zero, as follows from the Gamma-ratio generating function CITE{gamma}.'}
        return v

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    g=LogSine()
    if '--submit' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(g.publish(message='Certified log-sine moments with explicit series tails and independent integral checks.'))
    else:
        report=g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)