#AffinelyInd.maple
#
# Frank Sottile
# Thorsten Theobald
# Amherst, MA, April 2001
#
#  This generates a Singular input file for the geometric configuration of
# Section 3.  
#
#   The first 4 spheres to have radius r and centers in R^3 at 
#    (1,1,1), (1,-1,-1), (-1,1,-1), and (-1,-1,1) 
#  Then the additional spheres will also have radius r, but centers at
#  +/- A * e_i, where e_i are the basis elements with a 1 in the ith slot.
#
#  Both A and r are parameters.  When n=3, 4, and 5, we computed the full
#  Groebner basis in the coordinates (affine for p(i) and homogeneous for
#  v(i)).   This requires a saturation.  In the three cases, the ideal factors
#  into three pieces, each piece has n linear and (n-1) quadratic generators.
#  Since the case n=5 takes 6500+ seconds to do, we can only test larger
#  instances if we use the (expected) linear relations to reduce the number of 
#  variables.
#
#   Here are the generators of one of the 3 components for n=3, 4, and 5.
#   (Observe that if we allow sqrt(r^2-2) and some other square roots of
#    expressions in the parameters, the ideals factor completely.)
#
# n=3:   
#     p(3), p(2), (1)
#     p(1)^2+(-r^2+2)
#     p(1)*v(2)*v(3)+(r^2-2)*v(2)^2+(r^2-2)*v(3)^2
# n=4:
#     p(4), p(3), p(2), v(1)
#     (-A^2+2)*v(2)^2+(-A^2+2)*v(3)^2+3*v(4)^2
#     (-A^2-1)*p(1)^2+(r^2*A^2+r^2-3*A^2)
#     (-A^2-1)*p(1)*v(4)^2+(-A^2+2)*v(2)*v(3)
# n=5:
#     p(5), p(4), p(3), p(2), v(1)
#     v(4)^2-v(5)^2
#     (A^2-2)*v(2)^2+(A^2-2)*v(3)^2+(A^2-6)*v(5)^2
#     (-A^2-2)*p(1)^2+(r^2*A^2+2*r^2-4*A^2)
#     (-A^2-2)*p(1)*v(5)^2+(-A^2+2)*v(2)*v(3)
#
#   From these examples, it seems that the component has
#      (*)     p(n) = p(n-1) = ... = p(2) = 0 = v(1), 
#   and there is some homogeneous quadric involving only v(2) and v(3).
#
#   We shall assume (*), and use it to halve the number of variables.   
#  It suffices to saturate wrt the ideal v(2), as we have already
#  removed the zero dimensional pieces.
#
interface(quiet=true):
with(linalg,vector):

n:=18:
dot := proc(a,b)
  linalg[dotprod](a,b,`orthogonal`)
end:

EQ := proc(X,r)
     simplify(dot(V,V)*(dot(P,P)-2*dot(P,X)+dot(X,X))-dot(V,X)^2
              -r^2*dot(V,V))
end:

########################################

V:=[]:     P:=[]:
for i from 1 to n do
 V:=[V[],v.i]:   P:=[P[],p.i]:
od:

V:=vector(V):    P:=vector(P):

P1:=[ 1, 1, 1]:
P2:=[ 1,-1,-1]:
P3:=[-1, 1,-1]:
P4:=[-1,-1, 1]:

for i from 1 to 4 do
 for j from 4 to n do
  P.i := [P.i[],0]:
 od:
 X.i := linalg[vector](P.i):
od:

for i from 1 to 4 do
 r.i:=r:
od:
for i from 5 to 2*(n-1) do
 r.i:=r:
od:

for j from 3 to n-1 do
 X.(2*j-1):=[0,0,0]:
 X.(2*j):=[0,0,0]:
 for i from 3 to n-1 do
  if i<>j then 
   X.(2*j-1):=[X.(2*j-1)[],0]:
   X.(2*j)  :=[X.(2*j)[],0]:
  else
   X.(2*j-1):=[X.(2*j-1)[],A]:
   X.(2*j)  :=[X.(2*j)[],-A]:
  fi:
 od:
 X.(2*j-1):=vector(X.(2*j-1)):
 X.(2*j):=vector(X.(2*j)):
od:

#   Here is the step where we set 
#      (*)     p(n) = p(n-1) = ... = p(2) = 0 = v(1), 
#

SUBS:=[p1=p(1),v1=0]: 
VARS:=cat(convert(p(1),string)):
IRREL:=convert(v(2),string):
for i from 2 to n do
 VARS:=cat(VARS,",",convert(v(i),string)):
 IRREL:=cat(IRREL,",",convert(v(i),string)):
 SUBS:=[SUBS[],p.i=0,v.i=v(i)]:
od:

################################################################
lprint(`LIB "elim.lib";`);
lprint(`option(redSB);`);
lprint(`int T = timer;`);
lprint(`ring R= (0,r,A), (`,convert(VARS,symbol),`),dp;`);
lprint(`ideal Irel = `,convert(IRREL,symbol),`;`);
lprint(`ideal I =`);
for i from 1 to 4 do
 lprint(subs(SUBS[],primpart(EQ(X.i,r.i))),`,`);
od:
for j from 3 to n-1 do
 lprint(subs(SUBS[],simplify((EQ(X.(2*j-1),r.(2*j-1))-EQ(X.(2*j),r.(2*j)))/A)),`,`);
 lprint(subs(SUBS[],primpart(EQ(X.(2*j-1),r.(2*j-1))+EQ(X.(2*j),r.(2*j)))),`,`);
od:
lprint(subs(SUBS[],dot(V,P)),`;`);
lprint(`ideal G = std(sat(I,v(2))[1]);`);
lprint(`G;`);
lprint(`degree(G);`);
lprint(`timer - T;`);
lprint(`quit;`);

quit;