\\ PARI/GP code to implement the rigorous verification as
\\ described in readme.pdf.
\\ See https://pari.math.u-bordeaux.fr to download PARI/GP.

\\ The beginning of this file is a crude but rigorous implementation
\\ of interval arithmetic, since it is not built into PARI/GP.

israt(x) = (type(x) == "t_FRAC") || (type(x) == "t_INT");
\\ Is x rational?

isvalid(x) = (length(x)==2) && israt(x[1]) && israt(x[2]) && (x[1] <= x[2]);
\\ Check that x is a valid rational interval, for use in interval arithmetic.

times(x,y) = [vecmin([x[1]*y[1],x[1]*y[2],x[2]*y[1],x[2]*y[2]]),vecmax([x[1]*y[1],x[1]*y[2],x[2]*y[1],x[2]*y[2]])];
\\ Multiplication.

contained(x,y) = isvalid(x) && isvalid(y) && (x[1] >= y[1]) && (x[2] <= y[2]);
\\ Interval containment.

square(x) = if(x[1]>=0,return([x[1]^2,x[2]^2])); if(x[2]<=0,return([x[2]^2,x[1]^2])); return([0,max(x[1],x[2])^2]);
\\ Squaring.

digs = 36;
\\ How many digits to use in creating intervals.

standardize(x) = [floor(x[1]*10^digs)/10^digs,ceil(x[2]*10^digs)/10^digs];
\\ Turn a number x into an interval based on digs.

standardizevec(x) = vector(length(x),i,standardize(x[i]));
\\ Standardize a vector.

squareroot(x) = local(y); if(x[1]<0,error("Complex number.")); y=standardize([sqrt(x[1]),sqrt(x[2])]); if(!contained(x,square(y)),error("Failure in computing square root.")); return(y);
\\ Compute a square root. Note that this function uses floating-point
\\ arithmetic, but it rigorously verifies the result and raises an error
\\ if it is unable to verify it.

multby(x,y) = vector(length(x),i,times(x[i],y));
\\ Scalar multiplication of a vector x by y.

intIP(x,y) = sum(i=1,length(x),times(x[i],y[i]));
\\ Interval inner product of two vectors.

reciprocal(x) = if(x[1]*x[2]<=0,error("Division by 0.")); if(x[1]>0,return([1/x[2],1/x[1]]),return([1/x[1],1/x[2]]));
\\ Reciprocal of an interval.

intnormalize(x) = standardizevec(multby(x,reciprocal(squareroot(intIP(x,x)))));
\\ Rescale an interval vector to have length 1.

isless(x,y) = (x[2] < y[1]);
isgreater(x,y) = (x[1] > y[2]);
\\ Interval inequalities.

\\ Check a proof X of the existence of a spherical code.
\\ See readme.pdf for the input format.
{
verifyproofspherical(X) =
  local(n,G,xx,Gnum,p,a,b,q,c,d,L,N,m,eps,u1,u2,v,w,Gfull,r,aa,bb);
  if(length(X) != 9,print("Not list of nine elements.");return(0));
  n = X[1];
  G = X[2];
  p = X[3];
  a = X[4];
  b = X[5];
  q = X[6];
  c = X[7];
  d = X[8];
  L = X[9];
  if(type(G) != "t_MAT",print("Not matrix.");return(0));
  N = matsize(G)[1];
  if(N != matsize(G)[2],print("Not square matrix.");return(0));
  for(i=1,N,if(G[i,i]!=1,print("Not unit vectors.");return(0)));
  for(i=1,N,for(j=1,N,
    if((type(G[i,j])!="t_POL") && (type(G[i,j]) != "t_INT") && (type(G[i,j]) != "t_FRAC"),print("Matrix entries have wrong type.");return(0));
    for(k=0,poldegree(G[i,j]),if((type(polcoef(G[i,j],k,x)) != "t_INT")&&(type(polcoef(G[i,j],k,x)) != "t_FRAC"),print("Matrix entries have wrong type.");rreturn(0)));
  ));
  if(type(n) != "t_INT",print("Dimension is not integer.");return(0));
  if(n<=0,print("Dimension not positive.");return(0));
  if((type(a) != "t_FRAC") && (type(a) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if((type(b) != "t_FRAC") && (type(b) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if(a>=b,print("Interval is wrong.");return(0));
  if((type(c) != "t_FRAC") && (type(c) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if((type(d) != "t_FRAC") && (type(d) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if(c>=d,print("Interval is wrong.");return(0));
  if(type(x) != "t_POL",print("Variable name already used.");return(0));
  if(type(t) != "t_POL",print("Variable name already used.");return(0));
  if(type(p) != "t_POL",print("Not a polynomial.");return(0));
  if(type(q) != "t_POL",print("Not a polynomial.");return(0));
  for(i=0,poldegree(p),if(type(polcoef(p,i,x)) != "t_INT",print("Polynomial coefficient is not an integer.");return(0)));
  for(i=0,poldegree(q),if(type(polcoef(q,i,x)) != "t_INT",print("Polynomial coefficient is not an integer.");return(0)));
  if(!polisirreducible(p),print("Polynomial is not irreducible.");return(0));
  if(!polisirreducible(q),print("Polynomial is not irreducible.");return(0));
  if(subst(p,x,b)==0,print("Upper endpoint should not be the root.");return(0));
  if(polsturm(p,a,b)!=1,print("Root not isolated.);return(0));
  if(subst(q,x,d)==0,print("Upper endpoint should not be the root.");return(0));
  if(polsturm(q,c,d)!=1,print("Root not isolated.);return(0));
  f = lift(charpoly(Mod(G,p),t));
  if(N<n,print("Too few points.");return(0));
  for(i=0,N-n-1,if(polcoef(f,i,t)!=0,print("Gram matrix has incorrect rank.");return(0)));
  for(i=N-n,N,if(polsturm(polcoef(f,i,t),a,b)!=0,print("Intervals are too large to pin down signs of coefficients.");return(0)));
  for(i=N-n,N,if(subst(polcoef(f,i,t),x,(a+b)/2)*(-1)^(N-i)<=0,print("Gram matrix is not positive definite.");return(0)));
  m = -1;
  for(i=1,N,for(j=i+1,N,if(m!=G[i,j],
    if(polsturm(G[i,j]-m,a,b) != 0,print("Cannot determine minimal angle.");return(0));
    if(subst(G[i,j]-m,x,(a+b)/2)>0,m=G[i,j]);
  )));
  if(subst(q,x,m) % p != 0,print("Minimal polynomial of cosine is not correct.");return(0));
  if(m!=c,if(polsturm(m-c,a,b)!=0,print("Interval for cosine failed to check.");return(0)));
  if(subst(m-c,x,(a+b)/2)<0,print("Interval for cosine failed to check.");return(0));
  if(polsturm(m-d,a,b)!=0,print("Interval for cosine failed to check.");return(0));
  if(subst(m-d,x,(a+b)/2)>=0,print("Interval for cosine failed to check.");return(0));
  if(type(L) != "t_VEC",print("Incorrect list of rattlers.");return(0));
  for(i=1,length(L),if(type(L[i])!="t_VEC",print("Incorrect rattler.");return(0)));
  for(i=1,length(L),if(length(L[i])!=N,print("Incorrect number of coefficients for rattler.");return(0)));
  for(i=1,length(L),for(j=1,N,if((type(L[i][j]) != "t_FRAC") && (type(L[i][j]) != "t_INT"),print("Coefficient of rattler is not exact.");return(0))));
  default(realprecision,4096);
  if(poldegree(p)==1,xx=polrootsreal(p)[1],xx = solve(y=a,b,subst(p,x,y)));
  Gnum = subst(G,x,xx);
  for(i=1,N,for(j=1,N,if(poldegree(G[i,j])==0,Gnum[i,j]=polcoeff(G[i,j],0))));
  aa = floor(10^1000*xx)/10^1000;
  bb = ceil(10^1000*xx)/10^1000;
  if(!contained([aa,bb],[a,b]),print("Root refinement failed.");return(0));
  if(polsturm(p,aa,bb)!=1,print("Root refinement missed the root.");return(0));
  for(i=1,N,for(j=i+1,N,
    eps = 10^(-100);
    u1 = eps*floor(Gnum[i,j]/eps);
    u2 = eps*ceil(Gnum[i,j]/eps);
    v = subst(G[i,j],x,(aa+bb)/2);
    if(!contained([v,v],[u1,u2]),print("Rattler isolation failed.");return(0));
    if((poldegree(G[i,j])>0) && ((polsturm(G[i,j]-u1,aa,bb)>0) || (polsturm(G[i,j]-u2,aa,bb)>0)),print("Rattler intervals too large.");return(0));
    Gnum[i,j] = [u1,u2];
    Gnum[j,i] = [u1,u2];
  ));
  for(i=1,N,Gnum[i,i]=[1,1]);
  for(i=1,N,for(j=1,N,if(abs((Gnum[i,j][1]+Gnum[i,j][2])/2 - subst(G[i,j],x,xx))>10^(-50),print("Numerics failure.");return(0))));
  for(i=1,length(L),
    v = vector(length(L[i]),j,[L[i][j],L[i][j]]);
    w = sum(j=1,N,sum(k=1,N,times([L[i][j]*L[i][k],L[i][j]*L[i][k]],Gnum[j,k])));
    w = reciprocal(squareroot(w));
    L[i] = multby(v,w);
  );
  Gfull = matrix(N+length(L),N+length(L));
  for(i=1,N+length(L),Gfull[i,i]=[1,1]);
  for(i=1,N,for(j=1,N,Gfull[i,j]=Gnum[i,j]));
  for(i=1,N,for(j=1,length(L),
    v = L[j];
    w = sum(k=1,N,times(Gnum[i,k],v[k]));
    Gfull[i,N+j] = w;
    Gfull[N+j,i] = w;
    if(w[2] >= c,print("Rattlers come too close.");return(0));
  ));
  for(i=1,length(L),for(j=i+1,length(L),
    v = L[i];
    w = L[j];
    r = sum(ii=1,N,sum(jj=1,N,times(times(v[ii],w[jj]),Gnum[ii,jj])));
    Gfull[N+i,N+j] = r;
    Gfull[N+j,N+i] = r;
    if(r[2] >= c,print("Rattlers come too close to each other.");return(0));
  ));
  return(1);
}

\\ Check a proof X of the existence of a real projective code.
\\ See readme.pdf for the input format.
{
verifyproofprojective(X) =
  local(n,G,xx,Gnum,p,a,b,q,c,d,L,N,m,eps,u1,u2,v,w,Gfull,r,aa,bb);
  if(length(X) != 9,print("Not list of nine elements.");return(0));
  n = X[1];
  G = X[2];
  p = X[3];
  a = X[4];
  b = X[5];
  q = X[6];
  c = X[7];
  d = X[8];
  L = X[9];
  if(type(G) != "t_MAT",print("Not matrix.");return(0));
  N = matsize(G)[1];
  if(N != matsize(G)[2],print("Not square matrix.");return(0));
  for(i=1,N,if(G[i,i]!=1,print("Not unit vectors.");return(0)));
  for(i=1,N,for(j=1,N,
    if((type(G[i,j])!="t_POL") && (type(G[i,j]) != "t_INT") && (type(G[i,j]) != "t_FRAC"),print("Matrix entries have wrong type.");return(0));
    for(k=0,poldegree(G[i,j]),if((type(polcoef(G[i,j],k,x)) != "t_INT")&&(type(polcoef(G[i,j],k,x)) != "t_FRAC"),print("Matrix entries have wrong type.");rreturn(0)));
  ));
  if(type(n) != "t_INT",print("Dimension is not integer.");return(0));
  if(n<=0,print("Dimension not positive.");return(0));
  if((type(a) != "t_FRAC") && (type(a) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if((type(b) != "t_FRAC") && (type(b) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if(a>=b,print("Interval is wrong.");return(0));
  if((type(c) != "t_FRAC") && (type(c) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if((type(d) != "t_FRAC") && (type(d) != "t_INT"),print("Interval endpoint is not exact.");return(0));
  if(c>=d,print("Interval is wrong.");return(0));
  if(type(x) != "t_POL",print("Variable name already used.");return(0));
  if(type(t) != "t_POL",print("Variable name already used.");return(0));
  if(type(p) != "t_POL",print("Not a polynomial.");return(0));
  if(type(q) != "t_POL",print("Not a polynomial.");return(0));
  for(i=0,poldegree(p),if(type(polcoef(p,i,x)) != "t_INT",print("Polynomial coefficient is not an integer.");return(0)));
  for(i=0,poldegree(q),if(type(polcoef(q,i,x)) != "t_INT",print("Polynomial coefficient is not an integer.");return(0)));
  if(!polisirreducible(p),print("Polynomial is not irreducible.");return(0));
  if(!polisirreducible(q),print("Polynomial is not irreducible.");return(0));
  if(subst(p,x,b)==0,print("Upper endpoint should not be the root.");return(0));
  if(polsturm(p,a,b)!=1,print("Root not isolated.);return(0));
  if(subst(q,x,d)==0,print("Upper endpoint should not be the root.");return(0));
  if(polsturm(q,c,d)!=1,print("Root not isolated.);return(0));
  f = lift(charpoly(Mod(G,p),t));
  if(N<n,print("Too few points.");return(0));
  for(i=0,N-n-1,if(polcoef(f,i,t)!=0,print("Gram matrix has incorrect rank.");return(0)));
  for(i=N-n,N,if(polsturm(polcoef(f,i,t),a,b)!=0,print("Intervals are too large to pin down signs of coefficients.");return(0)));
  for(i=N-n,N,if(subst(polcoef(f,i,t),x,(a+b)/2)*(-1)^(N-i)<=0,print("Gram matrix is not positive definite.");return(0)));
  m = 0;
  for(i=1,N,for(j=i+1,N,if((m!=G[i,j]) && (m!=-G[i,j]),
    if(polsturm(G[i,j]-m,a,b) != 0,print("Cannot determine minimal angle.");return(0));
    if(polsturm(-G[i,j]-m,a,b) != 0,print("Cannot determine minimal angle.");return(0));
    if(subst(G[i,j]-m,x,(a+b)/2)>0,m=G[i,j]);
    if(subst(-G[i,j]-m,x,(a+b)/2)>0,m=-G[i,j]);
  )));
  if(subst(q,x,m) % p != 0,print("Minimal polynomial of cosine is not correct.");return(0));
  if(m!=c,if(polsturm(m-c,a,b)!=0,print("Interval for cosine failed to check.");return(0)));
  if(subst(m-c,x,(a+b)/2)<0,print("Interval for cosine failed to check.");return(0));
  if(polsturm(m-d,a,b)!=0,print("Interval for cosine failed to check.");return(0));
  if(subst(m-d,x,(a+b)/2)>=0,print("Interval for cosine failed to check.");return(0));
  if(type(L) != "t_VEC",print("Incorrect list of rattlers.");return(0));
  for(i=1,length(L),if(type(L[i])!="t_VEC",print("Incorrect rattler.");return(0)));
  for(i=1,length(L),if(length(L[i])!=N,print("Incorrect number of coefficients for rattler.");return(0)));
  for(i=1,length(L),for(j=1,N,if((type(L[i][j]) != "t_FRAC") && (type(L[i][j]) != "t_INT"),print("Coefficient of rattler is not exact.");return(0))));
  default(realprecision,4096);
  if(poldegree(p)==1,xx=polrootsreal(p)[1],xx = solve(y=a,b,subst(p,x,y)));
  Gnum = subst(G,x,xx);
  for(i=1,N,for(j=1,N,if(poldegree(G[i,j])==0,Gnum[i,j]=polcoeff(G[i,j],0))));
  aa = floor(10^1000*xx)/10^1000;
  bb = ceil(10^1000*xx)/10^1000;
  if(!contained([aa,bb],[a,b]),print("Root refinement failed.");return(0));
  if(polsturm(p,aa,bb)!=1,print("Root refinement missed the root.");return(0));
  for(i=1,N,for(j=i+1,N,
    eps = 10^(-100);
    u1 = eps*floor(Gnum[i,j]/eps);
    u2 = eps*ceil(Gnum[i,j]/eps);
    v = subst(G[i,j],x,(aa+bb)/2);
    if(!contained([v,v],[u1,u2]),print("Rattler isolation failed.");return(0));
    if((poldegree(G[i,j])>0) && ((polsturm(G[i,j]-u1,aa,bb)>0) || (polsturm(G[i,j]-u2,aa,bb)>0)),print("Rattler intervals too large.");return(0));
    Gnum[i,j] = [u1,u2];
    Gnum[j,i] = [u1,u2];
  ));
  for(i=1,N,Gnum[i,i]=[1,1]);
  for(i=1,N,for(j=1,N,if(abs((Gnum[i,j][1]+Gnum[i,j][2])/2 - subst(G[i,j],x,xx))>10^(-50),print("Numerics failure.");return(0))));
  for(i=1,length(L),
    v = vector(length(L[i]),j,[L[i][j],L[i][j]]);
    w = sum(j=1,N,sum(k=1,N,times([L[i][j]*L[i][k],L[i][j]*L[i][k]],Gnum[j,k])));
    w = reciprocal(squareroot(w));
    L[i] = multby(v,w);
  );
  Gfull = matrix(N+length(L),N+length(L));
  for(i=1,N+length(L),Gfull[i,i]=[1,1]);
  for(i=1,N,for(j=1,N,Gfull[i,j]=Gnum[i,j]));
  for(i=1,N,for(j=1,length(L),
    v = L[j];
    w = sum(k=1,N,times(Gnum[i,k],v[k]));
    Gfull[i,N+j] = w;
    Gfull[N+j,i] = w;
    w = max(abs(w[1]),abs(w[2]));
    if(w >= c,print("Rattlers come too close.");return(0));
  ));
  for(i=1,length(L),for(j=i+1,length(L),
    v = L[i];
    w = L[j];
    r = sum(ii=1,N,sum(jj=1,N,times(times(v[ii],w[jj]),Gnum[ii,jj])));
    Gfull[N+i,N+j] = r;
    Gfull[N+j,N+i] = r;
    r = max(abs(r[1]),abs(r[2]));
    if(r >= c,print("Rattlers come too close to each other.");return(0));
  ));
  return(1);
}