from __future__ import annotations

import math
import re
from pathlib import Path

import numpy as np


def H(N: int) -> np.ndarray:
    i = np.arange(1, N + 1)
    return np.minimum.outer(i, i).astype(float)


def Glog(N: int) -> np.ndarray:
    i = np.arange(1, N + 1)
    n = 2 * N - i
    x = np.log((2 * N) / n)
    return np.minimum.outer(x, x)


def min_eig(A: np.ndarray) -> float:
    return float(np.linalg.eigvalsh((A + A.T.conj()) / 2).min())


def test_log_time_change() -> None:
    for N in [2, 3, 5, 8, 16, 31, 64]:
        A = H(N)
        G = Glog(N)
        assert min_eig(A - N * G) > -1e-10
        assert min_eig(2 * N * G - A) > -1e-10


def test_slow_multiplier() -> None:
    for N in [2, 3, 5, 8, 16, 31, 64]:
        A = H(N)
        d = (2 * N - np.arange(1, N + 1)) / N
        D = np.diag(d)
        B = D @ A @ D
        w, V = np.linalg.eigh(A)
        Ais = V @ np.diag(1 / np.sqrt(w)) @ V.T
        M = Ais @ B @ Ais
        ev = np.linalg.eigvalsh((M + M.T) / 2)
        assert ev.min() >= 0.25 - 1e-9
        assert ev.max() <= 9 + 1e-9


def block_partition(N: int, L: int):
    out = []
    a = 0
    while a < N:
        b = min(N, a + L)
        out.append(np.arange(a, b))
        a = b
    return out


def test_block_fourier() -> None:
    for N, L in [(8, 3), (8, 4), (13, 4), (16, 5), (31, 7)]:
        C = np.tril(np.ones((N, N)))
        target = C @ C.T
        loc = np.zeros((N, N), dtype=complex)
        coarse = np.zeros((N, N), dtype=complex)
        blocks = block_partition(N, L)
        for idx in blocks:
            ell = len(idx)
            s = np.arange(ell)
            for r in range(ell):
                q = np.zeros(N, dtype=complex)
                q[idx] = np.exp(2j * np.pi * r * s / ell) / np.sqrt(ell)
                psi = C @ q
                if r == 0:
                    coarse += np.outer(psi, psi.conj())
                else:
                    # integrated zero-mean Fourier atom is block-supported
                    outside = np.ones(N, dtype=bool)
                    outside[idx] = False
                    assert np.max(np.abs(psi[outside])) < 1e-9
                    expected = 1 / (2 * np.sin(np.pi * r / ell) ** 2)
                    got = float(np.vdot(psi, psi).real)
                    assert abs(got - expected) < 1e-8 * max(1.0, expected)
                    loc += np.outer(psi, psi.conj())
        assert np.max(np.abs(target - (loc + coarse))) < 1e-8
        assert np.linalg.matrix_rank(coarse.real, tol=1e-8) <= math.ceil(N / L)


def test_source_delimiters() -> None:
    here = Path(__file__).resolve().parent
    papers = list(here.glob("*.md"))
    assert papers
    for p in papers:
        text = p.read_text(encoding="utf-8")
        for bad in [r"\\(", r"\\)", r"\\[", r"\\]"]:
            assert bad not in text, (p, bad)
        assert text.count("$$") % 2 == 0, p
        stripped = re.sub(r"\$\$.*?\$\$", "", text, flags=re.S)
        singles = re.findall(r"(?<!\\)\$", stripped)
        assert len(singles) % 2 == 0, p


def main() -> None:
    test_log_time_change()
    test_slow_multiplier()
    test_block_fourier()
    test_source_delimiters()
    print("Paper 49 deterministic validation: PASS")


if __name__ == "__main__":
    main()