\\ This is the file definitions.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.

\\ Number of points in RP^2 (used only in definition of T matrices).

N=7;

\\ Ultraspherical polynomials of degree k in variable x, for unit
\\ sphere in R^n:

{
P(n,k,x) = 

  local(L,alpha,p);

  \\ n=2 would lead to division by zero below, but fortunately that
  \\ case is built into PARI.  Clumsy substitutions is because the
  \\ PARI function poltchebi insists on getting an actual variable
  \\ name, so we use a dummy variable xx.
  if(n==2,return(subst(poltchebi(k,xx),xx,x)/subst(poltchebi(k,xx),xx,1)));

  alpha = n/2-1;

  \\ Polynomials of degree 0 through k, calculated recursively.
  L = vector(k+1);  
  L[1] = 1;
  L[2] = 2*alpha*x;
  \\ Now use three-term recurrence:
  for(i=2,k,L[i+1] = (1/i)*(2*x*(i+alpha-1)*L[i] - (i+2*alpha-2)*L[i-1]));

  p = L[k+1];

  \\ Recalculate value at 1:  
  \\ (Not just substituting for x in case x is a number.)
  L[1] = 1;
  L[2] = 2*alpha*1;
  for(i=2,k,L[i+1] = (1/i)*(2*1*(i+alpha-1)*L[i] - (i+2*alpha-2)*L[i-1]));

  return(p/L[k+1]);  \\ Return normalized value.
}

\\ Upper left d-by-d portion of matrix S-hat for RP^(n-1), in terms of
\\ variables u,v,t:
\\ (Uses uu,vv,tt,xx as dummy variables.  Note that the hat is omitted
\\ from the names of the functions here.)

{
S(n,k,d) =

  local(M,p);

  \\ Start with monomials in uu,vv with exponents congruent to k mod 2.
  M = matrix(d,d,i,j,uu^((k%2)+2*(i-1))*vv^((k%2)+2*(j-1)));

  p = P(n-1,k,(tt-uu*vv)/xx)*xx^k;
  \\ xx will be sqrt((1-uu^2)*(1-vv^2)), but must avoid square roots
  p = sum(i=0,k,polcoeff(p,2*i,xx)*((1-uu^2)*(1-vv^2))^i);
  \\ Replace xx^(2*i) with ((1-uu^2)*(1-vv^2))^i.
  M = M*p;

  \\ Return symmetrized matrix:

  return((subst(subst(subst(M,uu,u),vv,v),tt,t) +
          subst(subst(subst(M,uu,u),vv,t),tt,v) +
          subst(subst(subst(M,uu,v),vv,u),tt,t) +
          subst(subst(subst(M,uu,v),vv,t),tt,u) +
          subst(subst(subst(M,uu,t),vv,u),tt,v) +
          subst(subst(subst(M,uu,t),vv,v),tt,u))/6);
}

\\ Upper left d-by-d portion of matrix T-hat for RP^(n-1).

{
T(n,k,d) = (N-2)*S(n,k,d) + subst(subst(S(n,k,d),v,u),t,1) 
                          + subst(subst(S(n,k,d),u,v),t,1)
                          + subst(subst(S(n,k,d),u,t),v,1);
}