#GenQuad4.maple
#
# Frank Sottile
# Thorsten Theobald
# Amherst, MASS
# May 17, 2001
#
########################################################################
#
# This Maple program generates a Singular input file for the lines
# in P4 tangent to 6 random quadrics.
#
#  Adapted from code written by Bernd Sturmfels
#
interface(quiet=true):

die := rand(-7..7):
# range for random cordinates of symmetric form

circle := proc()
local A,X,M,Y,i,i1,i2,j1,j2,z1,z2:
Y:=[die()]:
for i from 2 to 15 do Y:=[Y[],die()]: od:
A := array([
   [Y[1], Y[2], Y[3], Y[4], Y[5] ],
   [Y[2], Y[6], Y[7], Y[8], Y[9] ],
   [Y[3], Y[7], Y[10], Y[11], Y[12] ],
   [Y[4], Y[8], Y[11], Y[13], Y[14] ],
   [Y[5], Y[9], Y[12], Y[14], Y[15] ]]):
#
#  A is a random symmetric matrix
#
M := array([],1..10,1..10):
z1 := 0:
for i1 from 1 to 4 do
 for j1 from i1+1 to 5 do
 z1 := z1+1:
 z2 := 0:
 for i2 from 1 to 4 do
  for j2 from i2+1 to 5 do
   z2 := z2+1:
    M[z1,z2] := 
     linalg[det](linalg[submatrix](A,[i1,j1],[i2,j2])):
od:od:od:od:
X := array([[p12,p13,p14,p15,p23,p24,p25,p34,p35,p45]]):
expand(linalg[multiply](linalg[multiply](X,M),linalg[transpose](X))[1,1])
end:

##########################################
#
#  Singular input file in local coordinates
#
lprint(`int t = timer; `);
lprint(`option(redSB); `);
lprint(`print("/////////////////////////////////////////////////////////////////////////");`);
lprint(`print("//");`);
lprint(`print("//  This calculation proceeds in positive characteristic, rather than");`);
lprint(`print("//  characteristic zero.  Since we do find that the eliminant has degree");`);
lprint(`print("//  320, the Shape Lemma holds and so the ideal is reduced if and only");`);
lprint(`print("//  if the eliminant has no repeated factors.  We show that it does not.");`);
lprint(`print("//");`);
lprint(`print("//  Our conclusion for characteristic zero uses the contrapositive of");`);
lprint(`print("//    {Multiple solutions of an ideal in characteristic zero}");`);
lprint(`print("//             ====>");`);
lprint(`print("//    {Multiple solutions in any reduction of that ideal}");`);
lprint(`print("//");`);
lprint(`print("/////////////////////////////////////////////////////////////////////////");`);
lprint(`print("//");`);
lprint(`print("//  We set up the rings and the ideal");`);
lprint(`print("//");`);
lprint(`ring S=31991, (a,b,c,d,e,f), (dp(5),dp(1));`);
lprint(`ring R=31991, (a,b,c,d,e,f), dp;`);
lprint(`ideal I = `);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`,`);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`,`);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`,`);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`,`);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`,`);
lprint(primpart(subs(p12=1,p13=d,p14=e,p15=f,p23=-a,p24=-b,p25=-c,
             p34=a*e-b*d,p35=a*f-d*c,p45=b*f-e*c,
             circle())),`;`);

lprint(`print("/////////////////////////////////////////////////////////////////////////");`);
lprint(`print("//");`);
lprint(`print("//  We compute the Groebner basis of the this ideal");`);
lprint(`print("//");`);
lprint(`ideal G=std(I);`);
lprint(`print("//  The ideal has degree 320 and dimension 0");`);
lprint(`degree(G);`);
lprint(`print("//  This is the time for the first Groebner basis calculation");`);
lprint(`timer - t;`);
lprint(`setring S;`);
lprint(`print("//  Now we compute an eliminant using fglm, it will be H[1]");`);
lprint(`ideal H = fglm(R,G);`);
lprint(`print("//  The eliminant has degree 320");`);
lprint(`degree(std(H[1]));`);
lprint(`print("//  We compute the ideal of the eliminant and its first derivative");`);
lprint(`ideal K = H[1], diff(H[1],f);`);
lprint(`print("//  Since this is the unit ideal, the eliminant has no repeated factors");`);
lprint(`print("//  and thus the original ideal is reduced.");`);
lprint(`std(K);`);
lprint(`print("//  This is the time for the full calculation");`);
lprint(`timer-t;`);
lprint(`quit;`);

quit