#2e6p2.maple
#
# Frank Sottile
# 11 October 1998
# Berkeley, CA
#
#   This computes many instances of Shapiro's Conjecture for the case of
# 2-planes in C^8 satisfying (J_2)^6 = 15. Of the 1821 cases checked 26 were
# anamolous; these are dealt with in the file at the bottom by computing a
# different eliminant.
#
#########################################################################
#
#Tested 1001 instances:   on Schubert.laptop
#    choose([1,2,3,4,5,6,7,8,9,10,11,12,13,14],4):
#     Task                      Time(sec.)     Size
#    Singular input file:
#    Compute eliminants. 
#    test for real roots:      14761.840  
#
# There were 23 anamolous cases: 
#  [1,2,3,6],  [1,2,4,8],   [1,2,5,10],  [1,2,6,12],  [1,2,7,14], [1,3,4,12],
#  [2,3,4,6],  [2,3,6,9],   [2,3,8,12],  [2,4,5,10],  [2,4,6,12], [2,4,7,14],
#  [3,4,6,8],  [3,4,9,12],  [3,5,6,10],  [3,6,7,14],  [4,5,8,10], [4,6,8,12],
#  [4,7,8,14], [5,6,10,12], [5,7,10,14], [6,7,12,14], [6,8,9,12]
#
#  which were dealt with by computing a different eliminant
#
#Tested 1820 instances:   On Blinn
#   choose([-14,-12,-11,-7,-5,-3,-2,2,3,5,17,19,23,29,31,37],4):
#     Task                      Time(sec.)      Size
#    Singular input file:          77.530      5 154 987
#    Compute eliminants:         4612          2 168 101
#    Verify all real roots:    170534.589
#
#  There were 3 anamolous cases: 
#     [-5, -3, 3, 5]   [-5, -2, 2, 5]    [-3, -2, 2, 3]
#
#
#
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([-14,-12,-11,-7,-5,-3,-2,2,3,5,17,19,23,29,31,37],4):
nops(trials);   quit;

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],
		[ 0 , 1 ,2*s,3*s^2,4*s^3, 5*s^4, 6*s^5, 7*s^6],
		[ 0 , 0 , 1 , 3*s ,6*s^2,10*s^3,15*s^4,21*s^5],
		[ 0 , 0 , 0 ,  1  , 4*s ,10*s^2,20*s^3,35*s^4],
		[ 0 , 0 , 0 ,  0  ,  1  ,  5*s ,15*s^2,35*s^3],
		[ 1 , a , b ,  c  ,  d  ,  0   ,  0   ,   0  ],
		[ 0 , 0 , 0 ,  1  ,  e  ,  f   ,  g   ,   h  ]]):


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

GB:=gbasis([A,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(`timer = 1;`);
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<15 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;
###############################################################
#2e6p2.anamolous
#
#	This checks the anamolous cases of Shapiro's
#	Conjecture for the case of 2-planes in C^8 satisfying
#	(J_2)^6 = 15.
#
#
interface(quiet=true):
with(Ore_algebra):
with(linalg):
with(Groebner):
with(combinat):
trials:=[[1,2,3,6],[1,2,4,8],[1,2,5,10],[1,2,6,12],[1,2,7,14],[1,3,4,12],[2,3,4,6],[2,3,6,9],[2,3,8,12],[2,4,5,10],[2,4,6,12],[2,4,7,14],[3,4,6,8],[3,4,9,12],[3,5,6,10],[3,6,7,14],[4,5,8,10],[4,6,8,12],[4,7,8,14],[5,6,10,12],[5,7,10,14],[6,7,12,14],[6,8,9,12]]:
#                                   14761.840


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],
		[ 0 , 1 ,2*s,3*s^2,4*s^3, 5*s^4, 6*s^5, 7*s^6],
		[ 0 , 0 , 1 , 3*s ,6*s^2,10*s^3,15*s^4,21*s^5],
		[ 0 , 0 , 0 ,  1  , 4*s ,10*s^2,20*s^3,35*s^4],
		[ 0 , 0 , 0 ,  0  ,  1  ,  5*s ,15*s^2,35*s^3],
		[ 1 , a , b ,  c  ,  d  ,  0   ,  0   ,   0  ],
		[ 0 , 0 , 0 ,  1  ,  e  ,  f   ,  g   ,   h  ]]):


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

GB:=gbasis([A,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(`timer = 1;`);
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 Sh= 0, (a,b,c,d,e,f,g,h), (dp(7),dp(1));`);
lprint(`ideal Gh;`);
lprint(`ring Sg= 0, (h,a,b,c,d,e,f,g), (dp(7),dp(1));`);
lprint(`ideal Gg;`);
lprint(`ring Sf= 0, (g,h,a,b,c,d,e,f), (dp(7),dp(1));`);
lprint(`ideal Gf;`);
lprint(`ring Se= 0, (f,g,h,a,b,c,d,e), (dp(7),dp(1));`);
lprint(`ideal Ge;`);
lprint(`ring Sd= 0, (e,f,g,h,a,b,c,d), (dp(7),dp(1));`);
lprint(`ideal Gd;`);
lprint(`ring Sc= 0, (d,e,f,g,h,a,b,c), (dp(7),dp(1));`);
lprint(`ideal Gc;`);
lprint(`ring Sb= 0, (c,d,e,f,g,h,a,b), (dp(7),dp(1));`);
lprint(`ideal Gb;`);
lprint(`ring Sa= 0, (b,c,d,e,f,g,h,a), (dp(7),dp(1));`);
lprint(`ideal Ga;`);


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 Sg;`);
 lprint(`Gg  = fglm(R,H);`);
 lprint(`print("POLYg:=");`);
 lprint(`short=0;`);
 lprint(`Gg[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYg, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`g`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);


 lprint(`setring Sf;`);
 lprint(`Gf  = fglm(R,H);`);
 lprint(`print("POLYf:=");`);
 lprint(`short=0;`);
 lprint(`Gf[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYf, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`f`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);


 lprint(`setring Se;`);
 lprint(`Ge  = fglm(R,H);`);
 lprint(`print("POLYe:=");`);
 lprint(`short=0;`);
 lprint(`Ge[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYe, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`e`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);

 lprint(`setring Sd;`);
 lprint(`Gd  = fglm(R,H);`);
 lprint(`print("POLYd:=");`);
 lprint(`short=0;`);
 lprint(`Gd[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYd, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`d`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);


 lprint(`setring Sc;`);
 lprint(`Gc  = fglm(R,H);`);
 lprint(`print("POLYc:=");`);
 lprint(`short=0;`);
 lprint(`Gc[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYc, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`c`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);


 lprint(`setring Sb;`);
 lprint(`Gb  = fglm(R,H);`);
 lprint(`print("POLYb:=");`);
 lprint(`short=0;`);
 lprint(`Gb[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYb, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`b`);
 lprint(convert(ST_end,symbol),` ");`);
 lprint(`print("fi: ");`);


 lprint(`setring Sa;`);
 lprint(`Ga  = fglm(R,H);`);
 lprint(`print("POLYa:=");`);
 lprint(`short=0;`);
 lprint(`Ga[1];`);
 lprint(`print(":");`);
 lprint(`print("NREAL:= nops(realroot(POLYa, 1/100)): ");`);
 lprint(`print(" if NREAL<15 then");`);
 lprint(`print("`,convert(ST_Nreal,symbol));
 lprint(convert(Vals,symbol),`a`);
 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;