#2e6p3.maple
#
# Frank Sottile
# 3 November 1998
# Berkeley, CA
#
#	This computes instances of Shapiro's Conjecture for the case of
#    3-planes in C^7 satisfying (J_2)^6 = 16.
#
#Tested 1001 instances 
#    choose([1,2,3,4,5,6,7,8,9,10,11,12,13,14],4):
#
#     Task			Time(sec.)     Size
#
#    Singular input file:        91           398734
#    Compute eliminants:      58979            94823
#    Verify all real roots:    5846   
#
#Tested 1001 instances 
#    choose([-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,15,16,17],4):
#
#     Task			Time(sec.)     Size
#
#    Singular input file:         18         3800975
#    Compute eliminants:       11601          878125
#    Verify all real roots:     4810
#
#realroots =    10   Vals =    [-3,-2,2,3],[-3,-1,1,3] [-2,-1,1,2]
#	   Dealt with these by a linear change of coordinates
#
#		[ 0 , 1 , c ,  d+h ,  e+h, f+h  ,  0   ],
#		[ 0 , 0 , 0 ,  0  ,  1  ,  g  ,  h+g   ]]):
#
interface(quiet=true):
with(Ore_algebra):
with(linalg):
with(Groebner):
with(combinat):
#trials:=choose([1,2,3,4,5,6,7,8,9,10,11,12,13,14],4):
#trials:=choose([-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,15,16,17],4):
#nops(trials);   quit;
trials:=[[-3,-2,2,3],[-3,-1,1,3],[-2,-1,1,2]]:
ST_Nreal:="lprint(`realroots = `,NREAL,`Vals = `,":
ST_end:=");":

SmMat := s -> matrix([
		[ 1 , s ,s^2, s^3 , s^4 ,  s^5 ,   s^6],
		[ 0 , 1 ,2*s,3*s^2,4*s^3, 5*s^4, 6*s^5],
		[ 0 , 0 , 1 , 3*s ,6*s^2,10*s^3,15*s^4],
		[ 1 , a , b ,  0  ,  0  ,  0   ,  0   ],
		[ 0 , 1 , c ,  d  ,  e  ,  f   ,  0   ],
		[ 0 , 0 , 0 ,  0  ,  1  ,  g   ,  h   ]]):


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

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

#for i from 1 to nops(GB) do indets(GB[i]);  od;

lprint(`option(redSB);`);
lprint(`int t=timer;`);
lprint(`ring R= 0, (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), (dp(7),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<16 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;