#!/usr/bin/env python3
"""Verify a strict finite-degree separation from all four-factor products.

The candidate has degree 61.  Its exact value at lambda=17/4 is below the
exact minimum of the entire degree<=80 four-factor Dirichlet product
family.  The comparison uses rational numbers only.
"""

from __future__ import annotations

import argparse
import json
import math
from fractions import Fraction
from pathlib import Path

from exact_toeplitz import normalized_fourier, normalized_objective

F = Fraction
Q_TERMS = {0: 1, 9: -1, 10: -1, 31: 1, 34: 1, 55: -1, 56: -1, 65: 1}


def witness_polynomial() -> list[int]:
    coefficients = [Q_TERMS.get(j, 0) for j in range(66)]
    for _ in range(4):
        if sum(coefficients) != 0:
            raise AssertionError("The supplied numerator is not divisible by (1-z)^4")
        running = 0
        quotient = []
        for coefficient in coefficients[:-1]:
            running += coefficient
            quotient.append(running)
        coefficients = quotient
    if coefficients[-1] != 1:
        raise AssertionError("The witness must be monic")
    return coefficients


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--results-dir", type=Path, default=Path(__file__).resolve().parent / "results")
    args = parser.parse_args()
    directory = args.results_dir
    directory.mkdir(parents=True, exist_ok=True)
    lam, alpha = F(17, 4), F(1, 4)
    polynomial = witness_polynomial()
    value = normalized_objective(polynomial, normalized_fourier(lam, 61))
    q_coefficients = [Q_TERMS.get(j, 0) for j in range(66)]
    q_value = normalized_objective(q_coefficients, normalized_fourier(alpha, 65))
    c0_ratio = F(1)
    for j in range(4):
        c0_ratio *= 2*(2*(alpha+j)+1)/(alpha+j+1)
    if q_value != value*c0_ratio:
        raise AssertionError("The two exact energy calculations disagree")
    moments = [sum(coefficient*exponent**r for exponent, coefficient in Q_TERMS.items())
               for r in range(5)]
    if moments[:4] != [0, 0, 0, 0] or moments[4] == 0:
        raise AssertionError("The numerator must have order exactly four at one")
    products = json.loads((directory / "products_n80.json").read_text())["results"]
    product_row = next(row for row in products if F(row["lambda"]) == lam)
    if product_row["status"] != "complete_product_family" or product_row["degree_bound"] != 80:
        raise AssertionError("The product comparison requires a completed degree-80 product search")
    product_minimum = F(product_row["ratio_to_c0_lambda"])
    difference = product_minimum - value
    if difference <= 0:
        raise AssertionError("The witness does not strictly improve on the product-family minimum")
    c0_display = math.gamma(2*float(lam)+1)/math.gamma(float(lam)+1)**2
    report = {
        "lambda": str(lam), "status": "exact_strict_product_family_separation",
        "degree": len(polynomial)-1, "polynomial": polynomial,
        "q_terms": [[j, coefficient] for j, coefficient in Q_TERMS.items()],
        "positive_support": [j for j, coefficient in Q_TERMS.items() if coefficient > 0],
        "negative_support": [j for j, coefficient in Q_TERMS.items() if coefficient < 0],
        "moments_zero_through_three_and_nonzero_at_four": moments,
        "normalized_value": str(value), "ratio_to_c0_alpha": str(q_value),
        "integral_over_2pi_display": float(value)*c0_display,
        "product_degree_bound": 80, "product_family_tuple_count": product_row["tuples_visited"],
        "product_minimum_lengths": product_row["lengths"],
        "normalized_product_minimum": str(product_minimum),
        "product_minimum_integral_over_2pi_display": float(product_minimum)*c0_display,
        "positive_normalized_gap": str(difference),
        "positive_gap_integral_over_2pi_display": float(difference)*c0_display,
        "normalization": "normalized values divide the integral by 2*pi*c0(17/4)",
    }
    (directory / "nonproduct_17_4.json").write_text(json.dumps(report, indent=2, sort_keys=True)+"\n")
    print(json.dumps({key: report[key] for key in
                      ("status", "degree", "integral_over_2pi_display",
                       "product_minimum_integral_over_2pi_display", "positive_gap_integral_over_2pi_display")}))


if __name__ == "__main__":
    main()