//! Exact verifier for the 548-cell continuum RLT certificate determining delta_3.
use std::{env,fs,process};
fn fail(s:impl AsRef<str>)->!{eprintln!("FAIL: {}",s.as_ref());process::exit(1)}
struct T<'a>{v:Vec<&'a str>,i:usize}impl<'a>T<'a>{fn new(s:&'a str)->Self{Self{v:s.split_whitespace().collect(),i:0}}fn n(&mut self)->&'a str{if self.i==self.v.len(){fail("unexpected EOF")}let z=self.v[self.i];self.i+=1;z}fn e(&mut self,w:&str){let z=self.n();if z!=w{fail(format!("expected {w}, got {z}"))}}fn z(&mut self,w:&str)->i128{self.n().parse().unwrap_or_else(|_|fail(format!("bad {w}")))}}
fn add(a:i128,b:i128,w:&str)->i128{a.checked_add(b).unwrap_or_else(||fail(format!("overflow {w}")))}fn sub(a:i128,b:i128,w:&str)->i128{a.checked_sub(b).unwrap_or_else(||fail(format!("overflow {w}")))}fn mul(a:i128,b:i128,w:&str)->i128{a.checked_mul(b).unwrap_or_else(||fail(format!("overflow {w}")))}fn abs(a:i128)->i128{a.checked_abs().unwrap_or_else(||fail("abs overflow"))}
fn main(){let path=env::args().nth(1).unwrap_or_else(||fail("usage: check-continuum-rlt CERT"));if env::args().nth(2).is_some(){fail("too many args")}let text=fs::read_to_string(path).unwrap_or_else(|e|fail(format!("read: {e}")));if text.len()>50_000_000{fail("too large")}let mut t=T::new(&text);t.e("AP_CONTINUUM_RLT_V1");t.e("L");let l=t.z("L")as usize;t.e("denom");let d=t.z("D");t.e("factor_scale");let rs=t.z("R");if l!=548||d<1||d>1_000_000_000||rs<1||rs>1_000_000_000{fail("parameters")}let n=l/2;let mut read_matrix=|label:&str|{t.e(label);let mut a=vec![0i128;n*n];for z in &mut a{*z=t.z(label);if *z<0||*z>10_i128.pow(15){fail("matrix range")}}a};let u=read_matrix("U");let v=read_matrix("V");
 fn read_sym(t:&mut T,label:&str,n:usize)->Vec<i128>{t.e(label);let mut a=vec![0;n*n];for i in 0..n{for j in 0..=i{let z=t.z(label);if z<0||z>10_i128.pow(15){fail("N range")}a[i*n+j]=z;a[j*n+i]=z}}a}
 let na=read_sym(&mut t,"NA",n);let nb=read_sym(&mut t,"NB",n);t.e("eta");let mut eta=vec![0;n];for z in &mut eta{*z=t.z("eta");if *z<0||*z>10_i128.pow(15){fail("eta range")}}
 fn read_factor(t:&mut T,label:&str,n:usize)->Vec<i128>{t.e(label);let mut a=vec![0;n*n];for i in 0..n{for j in 0..=i{let z=t.z(label);if abs(z)>10_i128.pow(12){fail("factor range")}a[i*n+j]=z}}a}
 let fs=read_factor(&mut t,"factor_sym",n);let fa=read_factor(&mut t,"factor_anti",n);if t.i!=t.v.len(){fail("trailing tokens")}
 // The published 12 runs; +1 and -1 are spin colors.
 let cuts=[0usize,28,34,62,99,158,274,390,449,486,514,520,548];let mut spin=vec![0i128;l];for b in 0..12{for z in &mut spin[cuts[b]..cuts[b+1]]{*z=if b%2==0{1}else{-1}}}if spin.iter().sum::<i128>()!=0{fail("candidate is not balanced")}
 let mut pref=vec![0i128;l+1];for i in 0..l{pref[i+1]=pref[i]+spin[i]}
 // Return 4L times the integral of the candidate on [0,k/(4L)].
 let lf4=|kk:i128|->i128{let k=kk.clamp(0,(4*l)as i128);let q=(k/4)as usize;let rem=k%4;if q==l{4*pref[l]}else{4*pref[q]+rem*spin[q]}};
 let mut m=vec![0i128;l];let mut slope=vec![0i128;l];for i in 0..l{let mut z=[0i128;3];for r in 0..=2i128{let two=8*(i as i128)+4*r;let h4=lf4(two)-lf4(two-4*l as i128)+lf4(2*i as i128+r+2*l as i128)-lf4(2*i as i128+r);z[r as usize]=-spin[i]*h4}let (m4,c4)=if 2*z[1]==z[0]+z[2]{(z[0].min(z[2]),abs(z[2]-z[0]))}else if z[0]==z[2]{(z[0].min(z[1]),abs(z[0]-z[1]))}else{fail(format!("cell {i} field is not uniform-linear"))};if m4%4!=0||c4%4!=0{fail(format!("nonintegral field distribution cell {i}"))}m[i]=m4/4;slope[i]=c4/4;if m[i]<0||!(1..=3).contains(&slope[i]){fail(format!("bad field distribution cell {i}"))}}
 for i in 0..n{if m[i]!=m[l-1-i]||slope[i]!=slope[l-1-i]{fail("field reflection")}}
 // Build 2C, where C is the conservative cell-pair lower-bound matrix.
 let mut orig=vec![0i128;l*l];for i in 0..l{for j in 0..l{let rr=2*(i as isize)-j as isize;let full=rr>=1&&rr<=l as isize-2;let empty=rr<=-2||rr>=l as isize+1;orig[i*l+j]=if full{2*spin[i]*spin[j]}else if !empty&&spin[i]*spin[j]<0{-2}else{0}}}let mut h2=vec![0i128;l*l];for i in 0..l{for j in 0..l{h2[i*l+j]=orig[i*l+j]+orig[j*l+i]}h2[i*l+i]+=slope[i]}for i in 0..l{for j in 0..l{if h2[i*l+j]!=h2[(l-1-i)*l+(l-1-j)]{fail("quadratic matrix reflection")}}}
 // Reconstruct symmetric and antisymmetric sector numerators, scaled by 2D.
 let mut ps=vec![0i128;n*n];let mut pa=vec![0i128;n*n];for i in 0..n{let mut budget=mul(n as i128,eta[i],"budget");for j in 0..n{budget=add(budget,add(u[i*n+j],v[i*n+j],"budget"),"budget")}if budget>mul(m[i],d,"budget limit"){fail(format!("row budget {i}"))}for j in 0..n{let ha=h2[i*l+j];let hb=h2[i*l+(l-1-j)];let su=add(add(u[i*n+j],u[j*n+i],"U sym"),add(v[i*n+j],v[j*n+i],"V sym"),"UV");let au=sub(add(u[i*n+j],u[j*n+i],"U sym"),add(v[i*n+j],v[j*n+i],"V sym"),"UV");ps[i*n+j]=add(sub(mul(d,ha+hb,"base sym"),2*(na[i*n+j]+nb[i*n+j]),"N sym"),add(su,2*(eta[i]+eta[j]),"RLT sym"),"P sym");pa[i*n+j]=add(sub(mul(d,ha-hb,"base anti"),2*na[i*n+j]-2*nb[i*n+j],"N anti"),au,"P anti")}}for i in 0..n{for j in 0..n{if ps[i*n+j]!=ps[j*n+i]||pa[i*n+j]!=pa[j*n+i]{fail("sector nonsymmetric")}}}
 fn check_psd(a:&[i128],f:&[i128],n:usize,d:i128,r:i128)->i128{let rr=mul(r,r,"R2");let mut least:i128=i128::MAX;for i in 0..n{let mut rad=0;let mut dia=0;for j in 0..n{let mut dot=0;for k in 0..=i.min(j){dot=add(dot,mul(f[i*n+k],f[j*n+k],"dot"),"dot")}let q=sub(mul(rr,a[i*n+j],"R2P"),mul(2*d,dot,"2DFF"),"residual");if i==j{dia=q}else{rad=add(rad,abs(q),"radius")}}if dia<rad{fail(format!("PSD residual row {i}"))}least=least.min(dia-rad)}least}
 let ms=check_psd(&ps,&fs,n,d,rs);let ma=check_psd(&pa,&fa,n,d,rs);
 // Independently enumerate the candidate's exact equal-bead upper score.
 let mut score=0u128;for i in 0..l{for j in 0..l{if (i+j)%2==0{let k=(i+j)/2;if spin[i]==spin[j]&&spin[j]==spin[k]{score+=2}}else{let k=(i+j)/2;for q in [k,k+1]{if spin[i]==spin[j]&&spin[j]==spin[q]{score+=1}}}}}if score!=128232{fail(format!("candidate score {score}"))}
 println!("PASS: all 548 field distributions and RLT row budgets exact");println!("PASS: symmetric/antisymmetric PSD residual margins {ms}, {ma}");println!("PASS: continuum lower bound equals candidate; delta_3 = 117/2192");
}
