"""Structured-matrix builders and exact rational positivity certificates.

Builders are scalar-agnostic: they only index and copy, so they work with
``Fraction``, ``arb``, ``mp.mpf``, or floats.  Feed ``arb`` matrices to
:func:`rh_verify.certificates.certify_positive_ldl`; use the ``_fraction``
certificates here when the data is exactly rational and the conclusion
should not depend on precision at all.
"""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction
from itertools import combinations

from .certificates import CertificationError

__all__ = [
    "FractionLDL",
    "certify_positive_ldl_fraction",
    "certify_totally_positive_fraction",
    "determinant_fraction",
    "gram",
    "hankel",
    "leading_principal_minors_fraction",
    "toeplitz",
]


def hankel(values, order: int):
    """Order-``order`` Hankel matrix with entry ``values[row + column]``."""
    if order < 1:
        raise CertificationError("Hankel order must be at least one")
    if len(values) < 2 * order - 1:
        raise CertificationError(
            f"Hankel order {order} needs {2 * order - 1} values, got {len(values)}"
        )
    return [
        [values[row + column] for column in range(order)]
        for row in range(order)
    ]


def toeplitz(first_row, first_column):
    """Toeplitz matrix from its first row and first column.

    ``first_row[0]`` and ``first_column[0]`` must be the same corner entry.
    """
    if not first_row or not first_column:
        raise CertificationError("Toeplitz needs a nonempty row and column")
    if first_row[0] != first_column[0]:
        raise CertificationError(
            "Toeplitz corner mismatch: "
            f"{first_row[0]!r} versus {first_column[0]!r}"
        )
    return [
        [
            first_row[column - row] if column >= row else first_column[row - column]
            for column in range(len(first_row))
        ]
        for row in range(len(first_column))
    ]


def gram(vectors, inner_product):
    """Gram matrix ``G[i][j] = inner_product(vectors[i], vectors[j])``."""
    if not vectors:
        raise CertificationError("Gram matrix needs at least one vector")
    return [
        [inner_product(left, right) for right in vectors]
        for left in vectors
    ]


def _fraction_entries(matrix, label: str) -> list[list[Fraction]]:
    entries = [[Fraction(value) for value in row] for row in matrix]
    size = len(entries)
    if size == 0 or any(len(row) != size for row in entries):
        raise CertificationError(f"{label} must be a nonempty square matrix")
    return entries


def _validated_symmetric_fraction(matrix, label: str) -> list[list[Fraction]]:
    entries = _fraction_entries(matrix, label)
    for row in range(len(entries)):
        for column in range(row + 1, len(entries)):
            if entries[row][column] != entries[column][row]:
                raise CertificationError(
                    f"{label} is not symmetric at ({row}, {column}): "
                    f"{entries[row][column]} versus {entries[column][row]}"
                )
    return entries


@dataclass(frozen=True)
class FractionLDL:
    """Exact ``LDL^T`` witness: unit lower factor and positive pivots."""

    order: int
    lower: tuple[tuple[Fraction, ...], ...]
    pivots: tuple[Fraction, ...]


def certify_positive_ldl_fraction(matrix, *, label: str = "matrix") -> FractionLDL:
    """Exact positive-definiteness certificate over the rationals."""
    work = _validated_symmetric_fraction(matrix, label)
    size = len(work)
    lower = [
        [Fraction(1) if row == column else Fraction(0) for column in range(size)]
        for row in range(size)
    ]
    pivots: list[Fraction] = []
    for pivot_index in range(size):
        pivot = work[pivot_index][pivot_index]
        if pivot <= 0:
            raise CertificationError(
                f"{label} LDL pivot {pivot_index} is not positive: {pivot}"
            )
        pivots.append(pivot)
        for row in range(pivot_index + 1, size):
            multiplier = work[row][pivot_index] / pivot
            lower[row][pivot_index] = multiplier
            for column in range(pivot_index, size):
                work[row][column] -= multiplier * work[pivot_index][column]
    return FractionLDL(
        order=size,
        lower=tuple(tuple(row) for row in lower),
        pivots=tuple(pivots),
    )


def determinant_fraction(matrix, *, label: str = "matrix") -> Fraction:
    """Exact determinant by fraction-free row elimination with pivoting."""
    work = _fraction_entries(matrix, label)
    size = len(work)
    sign = 1
    for pivot_index in range(size):
        pivot_row = next(
            (
                row
                for row in range(pivot_index, size)
                if work[row][pivot_index] != 0
            ),
            None,
        )
        if pivot_row is None:
            return Fraction(0)
        if pivot_row != pivot_index:
            work[pivot_index], work[pivot_row] = work[pivot_row], work[pivot_index]
            sign = -sign
        pivot = work[pivot_index][pivot_index]
        for row in range(pivot_index + 1, size):
            multiplier = work[row][pivot_index] / pivot
            for column in range(pivot_index, size):
                work[row][column] -= multiplier * work[pivot_index][column]
    determinant = Fraction(sign)
    for index in range(size):
        determinant *= work[index][index]
    return determinant


def leading_principal_minors_fraction(
    matrix, *, label: str = "matrix"
) -> tuple[Fraction, ...]:
    """Exact leading principal minors, smallest block first."""
    entries = _fraction_entries(matrix, label)
    return tuple(
        determinant_fraction(
            [row[: order + 1] for row in entries[: order + 1]],
            label=f"{label} leading minor {order + 1}",
        )
        for order in range(len(entries))
    )


def certify_totally_positive_fraction(
    matrix, *, max_order: int | None = None, strict: bool = True,
    label: str = "matrix",
) -> int:
    """Exact total positivity: every minor up to ``max_order`` is positive.

    With ``strict=False`` the certificate is total nonnegativity instead.
    The minor count grows combinatorially; this is meant for the small
    matrices that appear in PF and TP2/TP3 certificates, not for bulk
    numerics.  Returns the number of minors checked.
    """
    entries = _fraction_entries(matrix, label)
    size = len(entries)
    top = size if max_order is None else max_order
    if top < 1 or top > size:
        raise CertificationError(
            f"{label} minor order bound {top} is outside [1, {size}]"
        )
    checked = 0
    for order in range(1, top + 1):
        for rows in combinations(range(size), order):
            for columns in combinations(range(size), order):
                minor = determinant_fraction(
                    [[entries[row][column] for column in columns] for row in rows],
                    label=f"{label} minor {rows}x{columns}",
                )
                if minor < 0 or (strict and minor == 0):
                    bound = "positive" if strict else "nonnegative"
                    raise CertificationError(
                        f"{label} minor rows={rows} columns={columns} "
                        f"is not {bound}: {minor}"
                    )
                checked += 1
    return checked
