#!/usr/bin/env python3
"""Turn a numerical diagonal SDP solution into an exact PSD certificate.

This is a generator only.  The independent Rust checker reconstructs every
integer matrix and checks the certificate using exact i128 arithmetic.
"""
from pathlib import Path
import numpy as np

HERE = Path(__file__).resolve().parent
L = 256
D = 10**6
R = 10**8
FACTOR_SLACK = 100
TARGET_SLACK = 200

d = np.loadtxt(HERE / "d_256_numerical.txt")
assert d.shape == (L,)
base = np.ceil(d * D).astype(np.int64)
indices = np.arange(L)
i = indices[:, None] + 1
j = indices[None, :] + 1
b = ((j + 1 <= 2*i) & (2*i <= j + L)).astype(np.int64)
a = b + b.T
factor_matrix = a + np.diag((base + FACTOR_SLACK) / D)
p = np.rint(np.linalg.cholesky(factor_matrix) * R).astype(np.int64)
diag = base + TARGET_SLACK
s = int(diag.sum())

from fractions import Fraction
bound = Fraction(1,16) - Fraction(1,8*L) - Fraction(s,8*D*L*L)
out = HERE / "l256_certificate.txt"
with out.open("w") as f:
    print("AP_SDP_LOWER_V1", file=f)
    print("L", L, file=f)
    print("denom", D, file=f)
    print("factor_scale", R, file=f)
    print("diag", *map(int, diag), file=f)
    print("factor", file=f)
    for row in range(L):
        print(*map(int, p[row,:row+1]), file=f)
    print("claimed_diag_sum", s, file=f)
    print("claimed_bound", bound.numerator, bound.denominator, file=f)
print(f"wrote {out}; sum={s}; bound={bound}={float(bound):.12f}")
