#!/usr/bin/env python3
"""Generate independent references for the real Legendre Pi integral."""

from __future__ import annotations

import argparse
import math
from pathlib import Path

ROOT = Path(__file__).resolve().parent.parent
OUTPUT = ROOT / "tests" / "EllipticThirdKindReference.inc"
GENERATOR_VERSION = 1
QUADRATURE_ORDER = 2048
CASES = (
    (0.7, 0.4, 0.6),
    (math.pi / 2, 0.4, 0.6),
    (1.2, -1.0, 0.3),
    (1.5, 0.99, 0.2),
    (1.2, 0.8, 0.99),
    (0.7, 1.0, 0.5),
    (1.0, -16.0, 0.95),
    (math.pi / 2, 0.8, 0.99),
    (math.pi / 2, -16.0, 0.0),
    (1.56, 0.999, 0.2),
    (0.7, 1.0, 1.0),
)


def gauss_legendre(order: int) -> tuple[list[float], list[float]]:
    nodes = [0.0] * order
    weights = [0.0] * order
    for i in range((order + 1) // 2):
        x = math.cos(math.pi * (i + 0.75) / (order + 0.5))
        for _ in range(32):
            p0, p1 = 1.0, x
            for degree in range(2, order + 1):
                p0, p1 = p1, ((2 * degree - 1) * x * p1 - (degree - 1) * p0) / degree
            derivative = order * (x * p1 - p0) / (x * x - 1.0)
            step = p1 / derivative
            x -= step
            if abs(step) <= 2e-16:
                break
        weight = 2.0 / ((1.0 - x * x) * derivative * derivative)
        nodes[i], nodes[order - i - 1] = -x, x
        weights[i] = weights[order - i - 1] = weight
    return nodes, weights


NODES, WEIGHTS = gauss_legendre(QUADRATURE_ORDER)


def integrate(phi: float, n: float, m: float) -> float:
    midpoint = phi / 2.0
    scale = phi / 2.0
    total = 0.0
    for node, weight in zip(NODES, WEIGHTS):
        theta = midpoint + scale * node
        sine = math.sin(theta)
        total += scale * weight / (
            math.sqrt(1.0 - m * sine * sine) * (1.0 - n * sine * sine)
        )
    return total


def pascal(value: float) -> str:
    return f"{value:.17g}"


def render() -> str:
    lines = [
        f"{{ Generated by tools/generate_elliptic_third_kind_data.py v{GENERATOR_VERSION}.",
        f"  {QUADRATURE_ORDER}-point Gauss-Legendre quadrature of DLMF 19.2.7. }}",
        "const",
        f"  EllipticThirdKindReferenceCount = {len(CASES)};",
        "  EllipticThirdKindReferences: array[0..EllipticThirdKindReferenceCount - 1] of record",
        "    Phi, N, M, PiValue: Double;",
        "  end = (",
    ]
    for index, (phi, n, m) in enumerate(CASES):
        comma = "," if index + 1 < len(CASES) else ""
        lines.append(
            f"    (Phi: {pascal(phi)}; N: {pascal(n)}; M: {pascal(m)}; "
            f"PiValue: {pascal(integrate(phi, n, m))}){comma}"
        )
    lines.extend(["  );", ""])
    return "\n".join(lines)


def main() -> int:
    parser = argparse.ArgumentParser()
    parser.add_argument("--check", action="store_true")
    args = parser.parse_args()
    content = render()
    if args.check:
        if not OUTPUT.exists() or OUTPUT.read_text(encoding="utf-8") != content:
            print(f"Reference data is stale: {OUTPUT.relative_to(ROOT)}")
            return 1
        print(f"Reference data is current: {OUTPUT.relative_to(ROOT)}")
        return 0
    OUTPUT.write_text(content, encoding="utf-8", newline="\n")
    print(f"Wrote {OUTPUT.relative_to(ROOT)}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())