"""Zeta polynomials of Fermat curves -- numberdb.org/T666.

Run with SageMath:
    sage -pip install numberdb             # once, install the client
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --fill          # write entries to this draft
    sage -python checks.py                   # independent exact checks

The equation is X^N + Y^N = Z^N. The value is the zeta numerator
P(T), with constant coefficient 1, not its reversed characteristic
polynomial. For p = 1 mod N use every ordered a,b with a,b,a+b nonzero
modulo N, and the positive Jacobi sum J(a,b). Frobenius has eigenvalue
-J(a,b), so P(T) is the product of (1 + J(a,b)*T), retaining multiplicities.
All arithmetic is exact. No floating-point approximation or digit cutoff occurs.
The finite extent and its rationale are in argument-selection.md.

With agents/sage.sh the key can arrive on stdin using
NUMBERDB_KEY_FROM_STDIN=1; NUMBERDB_PUBLISH=1 selects filling the draft.
That client operation writes entries; it cannot publish the table.
"""
import os
import sys
from functools import lru_cache
from math import isqrt
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.number_field.number_field import CyclotomicField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

R = PolynomialRing(ZZ, 'T')
T = R.gen()


def prime(n):
    return n >= 2 and all(n % d for d in range(2, isqrt(n) + 1))


def primitive_root(p):
    for g in range(2, p):
        if len({pow(g, k, p) for k in range(p - 1)}) == p - 1:
            return g
    raise ValueError('No primitive root')


@lru_cache(maxsize=None)
def character_sums(N, p):
    if N < 3 or not prime(p) or (p - 1) % N:
        raise ValueError('This algorithm requires N >= 3 and a prime p = 1 mod N')
    K = CyclotomicField(N, 'z')
    z = K.gen()
    g = primitive_root(p)
    logs = {pow(g, k, p): k % N for k in range(p - 1)}
    roots = [z**k for k in range(N)]
    sums = {}
    for a in range(1, N):
        for b in range(1, N):
            if a + b == N:
                continue
            histogram = [0] * N
            for x in range(2, p):
                histogram[(a*logs[x] + b*logs[(1-x) % p]) % N] += 1
            sums[a, b] = sum((ZZ(c)*roots[k] for k,c in enumerate(histogram)), K(0))
    return K, sums


@lru_cache(maxsize=None)
def zeta_polynomial(N, p):
    N, p = int(N), int(p)
    K, sums = character_sums(N, p)
    S = PolynomialRing(K, 'T')
    u = S.gen()
    polynomial = S.one()
    for J in sums.values():
        polynomial *= 1 + J*u
    # Coercion to ZZ refuses any nonrational or nonintegral coefficient.
    return R([ZZ(c) for c in polynomial.list()])


class FermatZeta(numberdb.Generator):
    table = 'T666'
    parameters = ('N', 'p')
    type = 'Z[]'
    rigour = 'exact'
    digits = 100  # exact polynomial serialization ignores this precision setting

    def enumerate(self, max_degree=8, max_prime=101):
        for N in range(3, max_degree + 1):
            for p in range(2, max_prime + 1):
                if prime(p) and (p - 1) % N == 0:
                    yield {'N': N, 'p': p}

    def value(self, params, digits):
        return zeta_polynomial(int(params['N']), int(params['p']))


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = FermatZeta()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Exact Fermat zeta numerators with independent Gauss-sum and point-count checks', assisted_by='codex-cli'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)