#2e7p2.maple
#
# Frank Sottile
# 3 November 1998
# Berkeley, CA
#
#	This computes instances of Conjecture for the case of 2-planes in
#    C^9 satisfying (13)^7 = 36.
#
#	This files generates a Singular input file, which Singular
#  uses to compute an eliminant for the system, and then create a MAPLE
#  input file, which checks to see if the eliminant has 36 real roots.
#  If the eliminant does not, then it prints an error message and
#  the values of the parameters for that instance.
#
#	We first checked the 252 instances of values in 
#   choose([1,2,3,4,5,6,7,8,9,10],5)
#                                 Time (sec.)   Size
#   Singular input file:              10        3987347        
#   Computing eliminants:	   29202        1108890
#   Checking for real roots:       85087
#
#	Next the 252 instances of values in 
#   choose([1,2,3,4,5,6,7,8,9,10],5)
#                                 Time (sec.)    Size
#   Singular input file:              4.33     1124095
#   Computing eliminants:	   8059         894705
#   Checking for real roots:      30738
# 
#  Ambiguity, (only 21 roots found) but computing an eliminant in x worked
#
#    [-4, -2, -1, 1, 4],  [-4, -1, 1, 2, 4]
#
interface(quiet=true):
with(Ore_algebra):
with(linalg):
with(Groebner):
with(combinat):
#trials:=choose([-6,-5,-4,-3,-2,-1,1,2,3,4],5):
#nops(trials);   quit;
trials:=[ [-4, -2, -1, 1, 4], [-4, -1, 1, 2, 4]]:
ST_Nreal:="lprint(`realroots = `,NREAL,`Vals = `,":
ST_end:=");":

#SmMat := s -> matrix([
#	[ 1 , s ,s^2, s^3 , s^4 ,  s^5 ,   s^6,   s^7,    s^8],
#	[ 0 , 1 ,2*s,3*s^2,4*s^3, 5*s^4, 6*s^5, 7*s^6,  8*s^7],
#	[ 0 , 0 , 1 , 3*s ,6*s^2,10*s^3,15*s^4,21*s^5, 28*s^6],
#	[ 0 , 0 , 0 ,  1  , 4*s ,10*s^2,20*s^3,35*s^4, 56*s^5],
#	[ 0 , 0 , 0 ,  0  ,  1  ,  5*s ,15*s^2,35*s^3, 70*s^4],
#	[ 0 , 0 , 0 ,  0  ,  0  ,   1  ,  6*s ,21*s^2, 56*s^3],
#	[ 1 , a , b ,  c  ,  d  ,  e   ,  0   ,  0   ,   0   ],
#	[ 0 , 0 , 0 ,  1  ,  f  ,  g   ,  h   ,   x  ,   y   ]]):
##

#A:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,4,5,6,7,8]))/s):
#B:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,4,5,6,7,9]))/s):
#AA:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,4,5,6,8,9]))/s/s):
#BB:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,4,5,7,8,9]))/s^4):
#CC:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,4,6,7,8,9]))/s^4):
#DD:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,3,5,6,7,8,9]))/s^5):
#EE:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,2,4,5,6,7,8,9]))/s^6):
#FF:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[1,3,4,5,6,7,8,9]))/s^7):
#GG:=simplify(det(submatrix(SmMat(s),[1,2,3,4,5,6,7,8],[2,3,4,5,6,7,8,9]))/s^8):


#GB:=gbasis([A,B,AA,BB,CC,DD,EE,FF,GG],
#	wdeg([1,2,3,4,5,1,2,3,4,5],[a,b,c,d,e,f,g,h,x,y])):

#for i from 1 to nops(GB) do indets(GB[i]); lprint(GB[i]); od;
#time();
#quit;
GB:=
-560*s^5*a*h+840*s^7*a*f+1176*s^4*b*h+210*s^9-70*x*s^5-175*s^8*f-896*s^8*a+1176*s^7*b-980*s^5*d+15*y*s^4-1470*s^6*b*f+980*s^5*c*f-840*x*s^3*b+896*x*s^2*c-420*x*s*d+384*x*s^4*a-84*y*s^3*a+189*y*s^2*b-210*y*s*c+105*y*d+105*s^6*h-1176*s^3*c*h+490*s^2*d*h,
20*x*s^2*c-105*s^7*a*f+84*s^4*b*h+70*s^2*d*h-35*s^5*a*h-21*s*e*h+84*s^6*a*g+210*s^6*b*f-189*s^5*b*g-175*s^5*c*f+210*s^4*c*g-105*s^3*d*g+105*s^3*e*f-105*s^3*c*h+6*x*s^4*a-15*x*s^3*b-15*x*s*d+6*x*e-15*s^9+6*s^6*h+20*s^8*f-15*s^7*g+70*s^8*a-105*s^7*b+175*s^5*d-210*s^4*e-x*s^5,
-420*s^4*b*x+448*s^3*c*x-210*s^2*d*x+192*s^5*a*x-35*s^6*x-1470*s^8*a*f+1568*s^5*b*h+980*s^3*d*h-700*s^6*a*h-196*s^2*e*h+1344*s^7*a*g+2940*s^7*b*f-2940*s^6*b*g-2450*s^6*c*f+3136*s^5*c*g-1470*s^4*d*g+1470*s^4*e*f-1764*s^4*c*h-189*s^10+126*s^7*h+280*s^9*f-245*s^8*g+896*s^9*a-1372*s^8*b+2450*s^6*d-3136*s^5*e-35*y*s^2*c+21*y*e-21*y*s^4*a+42*y*s^3*b+4*y*s^5:

lprint(`timer = 1;`);
lprint(`option(redSB);`);
lprint(`int t=timer;`);
lprint(`ring R= 0, (y,x,h,g,f,e,d,c,b,a), dp;`);
lprint(`ideal I;`);
lprint(`ideal H;`);
lprint(`ring S= 0, (a,b,c,d,e,f,g,h,x,y), (dp(9),dp(1));`);
lprint(`ideal G;`);


lprint(`print("interface(quiet=true):");`);
lprint(`print("readlib(realroot):");`);

for ntimes from 1 to nops(trials) do

 Vals:=trials[ntimes]:
 lprint(`setring R;`);
 lprint(`I = `);

 equations:=[]:
 for P in GB do
  for ii in Vals do
   equations:=[equations[],subs(s=ii,P)]:
 od:od:

 for EQS in equations do
 lprint(EQS,`,`);  od:
 lprint(`0;`);
 lprint(`H = std(I);`);
 lprint(`setring S;`);
 lprint(`G  = fglm(R,H);`);
 lprint(`print("POLY:=");`);
 lprint(`short=0;`);
 lprint(`G[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLY, 1/100)): ");`);
 lprint(`print(" if NREAL<36 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol));
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);
 lprint(`print(" ");`);
 lprint(`print(" ");`);
 lprint(`print(" ");`);
 lprint(`print("############################################################");`);
 lprint(`print("#                                                          #");`);
 lprint(`print("#          Next One                                        #");`);
 lprint(`print("#                                                          #");`);
 lprint(`print("############################################################");`);
 lprint(`print(" ");`);
od:

lprint(`print("time(); ");`);
lprint(`print("quit; ");`);
lprint(`timer-t;`);
lprint(`quit;`);

time();

quit;
quit;