//! Exact verifier for box-SDP/SOS lower-bound certificates for delta_3.
use std::{env,fs,process};
fn fail(s:impl AsRef<str>)->!{eprintln!("FAIL: {}",s.as_ref());process::exit(1)}
struct T<'a>{x:Vec<&'a str>,i:usize} impl<'a>T<'a>{fn new(s:&'a str)->Self{Self{x:s.split_whitespace().collect(),i:0}}fn n(&mut self)->&'a str{if self.i==self.x.len(){fail("unexpected EOF")}let z=self.x[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 u(&mut self,w:&str)->usize{self.n().parse().unwrap_or_else(|_|fail(format!("bad {w}")))}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 in {w}")))}fn sub(a:i128,b:i128,w:&str)->i128{a.checked_sub(b).unwrap_or_else(||fail(format!("overflow in {w}")))}fn mul(a:i128,b:i128,w:&str)->i128{a.checked_mul(b).unwrap_or_else(||fail(format!("overflow in {w}")))}fn abs(a:i128)->i128{a.checked_abs().unwrap_or_else(||fail("overflow in abs"))}fn gcd(mut a:i128,mut b:i128)->i128{a=abs(a);b=abs(b);while b!=0{let r=a%b;a=b;b=r}a}
fn ae(i0:usize,j0:usize,l:usize)->i128{let i=i0+1;let j=j0+1;((j+1<=2*i&&2*i<=j+l)as i128)+((i+1<=2*j&&2*j<=i+l)as i128)}
fn main(){let p=env::args().nth(1).unwrap_or_else(||fail("usage: check-box-sos-lower CERT"));if env::args().nth(2).is_some(){fail("too many arguments")}let s=fs::read_to_string(&p).unwrap_or_else(|e|fail(format!("read: {e}")));if s.len()>100_000_000{fail("too large")}let mut t=T::new(&s);t.e("AP_BOX_SOS_LOWER_V1");t.e("L");let l=t.u("L");t.e("denom");let d=t.z("D");t.e("factor_scale");let r=t.z("R");if l==0||l>512||l%2!=0||d<1||d>1_000_000_000||r<1||r>1_000_000_000{fail("parameter range")}
 t.e("diag");let mut av=vec![0i128;l];for x in &mut av{*x=t.z("diag");if *x<0||*x>10_i128.pow(15){fail("diag range")}}
 let mut h=vec![0i128;l*l];for i in 0..l{for j in 0..l{h[i*l+j]=mul(2*d,ae(i,j,l),"2DA")}}for i in 0..l{h[i*l+i]=add(h[i*l+i],2*av[i],"diag")}
 t.e("triangles");let nt=t.u("triangle count");if nt>2_000_000{fail("too many triangles")}let mut lsum=0i128;for _ in 0..nt{let typ=t.u("type");let i=t.u("i");let j=t.u("j");let k=t.u("k");let m=t.z("lambda");if typ>3||!(i<j&&j<k&&k<l)||m<0||m>10_i128.pow(15){fail("invalid triangle")}lsum=add(lsum,m,"lambda sum");let sg:[i128;3]=match typ{0=>[1,1,1],1=>[1,-1,-1],2=>[-1,1,-1],_=>[-1,-1,1]};for ((u,v),q)in[(i,j),(i,k),(j,k)].into_iter().zip(sg){let w=mul(m,q,"triangle");h[u*l+v]=sub(h[u*l+v],w,"triangle");h[v*l+u]=sub(h[v*l+u],w,"triangle")}}
 t.e("factor");let mut f=vec![0i128;l*l];for i in 0..l{for j in 0..=i{f[i*l+j]=t.z("factor");if abs(f[i*l+j])>10_i128.pow(12){fail("factor range")}}}
 t.e("claimed_multiplier_sum");let cs=t.z("sum");t.e("claimed_q_bound");let qn=t.z("q numerator");let qd=t.z("q denominator");t.e("claimed_delta_bound");let bn=t.z("bound numerator");let bd=t.z("bound denominator");if t.i!=t.x.len(){fail("trailing tokens")}
 let asum=av.iter().fold(0i128,|x,&y|add(x,y,"diag sum"));let total=add(asum,lsum,"multiplier sum");if total!=cs{fail(format!("sum {total} != claim {cs}"))}
 let rr=mul(r,r,"R2");let two_d=2*d;let mut mm=None;for i in 0..l{let mut rad=0i128;let mut dia=0;for j in 0..l{let mut dot=0;for k in 0..=i.min(j){dot=add(dot,mul(f[i*l+k],f[j*l+k],"dot"),"dot")}let q=sub(mul(rr,h[i*l+j],"R2H"),mul(two_d,dot,"2DPP"),"residual");if i==j{dia=q}else{rad=add(rad,abs(q),"radius")}}if dia<rad{fail(format!("residual row {i} not diagonally dominant"))}let z=dia-rad;mm=Some(mm.map_or(z,|v:i128|v.min(z)))}
 let gq=gcd(total,d);if qn!=-total/gq||qd!=d/gq{fail("wrong q bound")}
 let li=l as i128;let raw= sub(sub(mul(d,mul(li,li,"L2"),"bound"),mul(2*d,li,"bound"),"bound"),2*total,"bound");let den=mul(16*d,mul(li,li,"L2"),"bound");let gb=gcd(raw,den);if raw<=0||bn!=raw/gb||bd!=den/gb{fail("wrong delta bound")}
 if mul(bn,32768,"compare")<=mul(1675,bd,"compare"){fail("does not improve PRS08")}
 println!("PASS: exact box-SOS certificate (L={l}, triangles={nt}, min residual margin={})",mm.unwrap());println!("PASS: q_L >= {qn}/{qd} = {:.9}",qn as f64/qd as f64);println!("PASS: delta_3 >= {bn}/{bd} = {:.12}",bn as f64/bd as f64);
}
