\\ This is the file optimal.txt from the paper "Three-point bounds for
\\ energy minimization" by Henry Cohn and Jeechul Woo.  It consists of
\\ PARI/GP code (see http://pari.math.u-bordeaux.fr/), but it does not do
\\ anything sophisticated and can easily be translated to other computer
\\ algebra systems.

\\ Read in the definitions from definitions.txt:

\r definitions.txt

\\ Keep track of whether anything has failed so far:

failure = 0;

\\ Take a list of rational matrices as input and determine whether each
\\ of them is symmetric and positive semidefinite.

{
Check_Semidefinite(L) =
  local(M, p, failed, i, j);

  print(" There are ",#L," matrices to test.");
  for(i=1,#L,
    M = L[i]*denominator(L[i]);  \\ Clear denominators for efficiency.

    if(M != mattranspose(M),
      print("  Matrix ",i," is NOT symmetric!");failure=1,
      print("  Matrix ",i," is symmetric."));

    p = charpoly(M,x);  \\ Get characteristic polynomial.
    p = subst(p,x,-x)*(-1)^poldegree(p);
    \\ Will check this has no negative coefficients.

    failed = 0;  \\ Haven't failed so far in this test.
    j=0;
    while(!failed && (j < poldegree(p)),
      if(polcoeff(p,j)<0,print("  Matrix ",i," is NOT positive semidefinite!");
                         failed=1; failure=1);
      j=j+1;
    );
    if(!failed,print("  Matrix ",i," is positive semidefinite."));
  );
}

\\ Check that the list of matrices in fact defines the correct potential
\\ function, as described in Appendix D, and that the bound is sharp.
\\ Here the potential function f is given as a polynomial in x
\\ (where x represents the inner product squared).

{
Check_Equality(c, L, Degree_Bound, f) =
  local(Index, Potential_Function, Auxiliary_Function, Sum_of_Squares,
  failed, Bound, Energy, J);

  failed = 0; \\ Haven't failed so far.

  \\ Build the list of monomials of degree up to Degree_Bound in u,v,t:
  Index = [];
  for(i=0,Degree_Bound,
    for(j=0,Degree_Bound-i,
      for(k=0,Degree_Bound-i-j,
        Index=concat(Index,[u^i*v^j*t^k]))));

  Auxiliary_Function = c  \\ constant term
                       + sum(i=1,#L-1,
                             trace(L[i]*T(3,i-1,matsize(L[i])[1])));
  \\ Note that matsize(L[i])[1] ensures that we take the inner product of
  \\ L[i] with a T-hat matrix of matching size.

  Sum_of_Squares = Index*(L[#L])*mattranspose(Index);
  \\ Build a sum of squares out of the index and the last matrix in L.
  \\ Note that Index is a row vector.

  Potential_Function = (subst(f,x,u^2) + subst(f,x,v^2) + subst(f,x,t^2))/3;

  if(Potential_Function != (Auxiliary_Function + Sum_of_Squares),failed=1);

  J = matrix(matsize(L[1])[1],matsize(L[1])[1],i,j,1);
  Bound = N/2*((N-1)*c - trace(L[1]*J));
  \\ The SDP bound.

  Energy = 3*subst(f,x,0) + 6*subst(f,x,1/9) + 12*subst(f,x,1/3);
  \\ Actual energy of rhombic dodecahedron code.

  if(Energy != Bound,failed=1);

  if(failed,
    print(" Equality test FAILED!");failure=1,
    print(" Equality test succeeded."));
}

\\ Read in data files and verify calculations:

\r data1.txt;
print("Checking first potential function:");
Check_Semidefinite(L);
Check_Equality(c,L,7,x^2-1/1000*x^3*(x-1/9)^2*(x-1/3)^2);

\r data2.txt;
print("Checking second potential function:");
Check_Semidefinite(L);
Check_Equality(c,L,7,x^3-1/1000*x^3*(x-1/9)^2*(x-1/3)^2);

\r data3.txt;
print("Checking third potential function:");
Check_Semidefinite(L);
Check_Equality(c,L,7,x^3*(x-1/9)-1/1000*x^3*(x-1/9)^2*(x-1/3)^2);

\r data4.txt;
print("Checking fourth potential function:");
Check_Semidefinite(L);
Check_Equality(c,L,7,x^3*(x-1/9)^2-1/1000*x^3*(x-1/9)^2*(x-1/3)^2);

\r data5.txt;
print("Checking fifth potential function:");
Check_Semidefinite(L);
Check_Equality(c,L,7,x^3*(x-1/9)^2*(x-1/3)-1/1000*x^3*(x-1/9)^2*(x-1/3)^2);

if(!failure,print("Universal optimality has been verified."));