interface(quiet=true):
#Dependent.maple
#
#  We manually solve the equations when the spheres have affinely dependent
#  centers at
#     (0,a,0,...,0), (0,-1,0,...,0), and (0,0,...,+/-1,...,0)
lprint(`    `);
lprint(`    `);
lprint(`We manually solve the equations when the spheres have centers`);
lprint(`  a e_2,   -e_2,   and    \pm e_j   for  j=3,...,n`);
lprint(`Here are the equations`);
lprint(`        V^2*(P^2 \\pm 2pj+1-r^2)-vj^2   for  j=3,...,n`);
lprint(`        V^2*(P^2 - 2ap2+a^2-r^2)-v2^2`);
lprint(`        V^2*(P^2 + 2p2+1-r^2)-v2^2`);
lprint(`        P \cdot V`);
lprint(`###################################################################`);
lprint(`We have the 2n-4 equations for the spheres with centers at \\pm e_j`);
V^2*(P^2-2*pj+1-r^2)-vj^2;
V^2*(P^2+2*pj+1-r^2)-vj^2;
lprint(`Combining these, we obtain`);
lprint(`p3 = p4 = ... = pn = 0,`);
lprint(`v3^2 = v4^2 = ... = vn^2 = 0,   and`);
lprint(`the equation from the spheres with center \pm e_3:`);
v3^2-V^2*(P^2+1-r^2);
lprint(`################################################################`);
lprint(`The equations for the spheres with centers at ae_2 and -e_2 are:`);
f:=V^2*(P^2-2*a*p2+a^2-r^2)-a^2*v2^2;
g:=V^2*(P^2+2*p2+1-r^2)-v2^2;
lprint(`Consider f+a*g:`);
factor(simplify(f+a*g));
lprint(`and also f-a^2*g:`);
factor(simplify(f-a^2*g));
lprint(`Using V^2\\neq 0 and a\\neq -1, (which allow us to cancel factors) and previous `);
lprint(`equations, (so that P^2=p1^2+p2^2) we obtain the following system of equations:`);
lprint(``);
lprint(`p3 = p4 = ... = pn = 0`);
lprint(`v3^2 = v4^2 = ... = vn^2 = 0`);
w:=p1*v1+p2*v2;
x:=(1-a)*(p1^2+p2^2-r^2) - 2*a*p2;
y:=V^2*(p1^2+p2^2-r^2+a)-a*v2^2;
z:= v3^2-V^2*(p1^2+p2^2-r^2+1);

lprint(`We combine the second and third equations to solve for p2:`);
P2:=factor(solve(simplify(((1-a)*y-V^2*x)/a)=0,p2));
lprint(`Plugging into the first equation, we solve for p1:`);
P1:=solve(subs(p2=P2,w)=0,p1);

lprint(`We substitute these into the second equation and clear the denominator`);
lprint(`to obtain a sextic.  We use V^2=v1^2+v2^2+(n-2)v3^2`);
W:=sort(collect(expand(simplify(subs(V=sqrt(v1^2+v2^2+(n-2)*v3^2),
        simplify(4*v1^2*V^4*subs(p1=P1,p2=P2,x)/(1-a)))))
        ,[v1,v2,v3], distributed ),[v1,v2,v3]);
lprint(``);
lprint(`Here is the formula in the paper:`);
X:=(1-a)^2*(v1^2+v2^2)*(V^2-v2^2)^2 -4*r^2*v1^2*V^4+4*a*v1^2*V^2*(V^2-v2^2);
lprint(`We check this:    simplify(subs(V=sqrt(v1^2+v2^2+(n-2)*v3^2),W-X))`);
simplify(subs(V=sqrt(v1^2+v2^2+(n-2)*v3^2),W-X));
lprint(`  `);
lprint(`Adding the third equation and the fourth equation, z+y, we obtain`);
Y:=collect(simplify(subs(V=sqrt(v1^2+v2^2+(n-2)*v3^2),z+y)),[v1,v2,v3]);
lprint(`We substitute for the squared variables in this equation to obtain`);
lprint(`  `);

Z:=collect(subs(v1=sqrt(alpha),v2=sqrt(beta),v3=sqrt(gamma),Y),[alpha,beta,gamma]):
Z;
lprint(`Solving for beta and plugging into the sextic we obtain the cubic`);
lprint(` `);
s:=solve(Z=0,beta):
Xx:=collect(simplify(subs(v1=sqrt(alpha),v2=sqrt(s),v3=sqrt(gamma),W)),[alpha,gamma]):
factor(coeff(Xx,gamma^3))*gamma^3+
factor(coeff(Xx,gamma^2))*gamma^2+
factor(coeff(Xx,gamma))*gamma+
factor(coeff(Xx,alpha^3))*alpha^3;

lprint(` `);
lprint(`We compute the discriminant of this cubic, discrim(subs(gamma=1,Xx),alpha)`);
Disc:=discrim(subs(gamma=1,Xx),alpha):
lprint(`This discriminant is quite large, and it does not factor`);
interface(quiet=false):
nops(Disc);
degree(Disc);
nops(factor(Disc));
interface(quiet=true):
lprint(` `);
lprint(`We consider the leading and trailing coefficients of the cubic`);
lprint(`Coefficent of gamma^3:`);
factor(coeff(Xx,gamma^3));
lprint(`Coefficent of alpha^3:`);
factor(coeff(Xx,alpha^3));
lprint(`This shows that for each n, neither alpha nor gamma vanishes for general a, r.`);
lprint(` `);
lprint(`To check for beta, we solve the linear equation for alpha and plug that \n  into the cubic `);
lprint(` `);
s:=solve(Z=0,alpha):
Xx:=collect(simplify((a-1)^3*subs(v1=sqrt(s),v2=sqrt(beta),v3=sqrt(gamma),W)),[alpha,gamma]):

factor(coeff(Xx,beta^3))*beta^3+
factor(coeff(coeff(Xx,gamma),beta^2))*beta^2*gamma+
factor(coeff(coeff(Xx,gamma^2),beta))*beta*gamma^2+
factor(coeff(Xx,gamma^3))*gamma^3;
lprint(` `);
lprint(`We consider the leading and trailing coefficients of this cubic`);
lprint(`Coefficent of gamma^3:`);
factor(coeff(Xx,gamma^3));
lprint(`Coefficent of beta^3:`);
factor(coeff(Xx,beta^3));
lprint(`This shows that for each n, neither beta nor gamma vanishes for general a, r.`);
lprint(` `);
lprint(` `);
lprint(` `);