#32.maple
#
# Frank Sottile
# 10 October 1998
# Berkeley, CA
#
#   This maple script is intended to test many instances of 
#   Shapiro's conjecture for totally positive
#   in the case m=3, p=2.
#
#   Instances tested   
#
# _seed:=1: 	  1500 	  6484.460  seconds  =  1:48:04
# _seed:=123:     6000   39493.720  seconds  = 10:58:13
# _seed:=121212:  1500    6394.460  seconds  =
# _seed:=2934:	  3000   15387.479  seconds  =
#_seed:=818408:
#                                  
#

interface(quiet=true):
with(linalg):
with(Groebner):
die:=rand(1..14):
readlib(realroot):
_seed:=818408:

Pos:=proc(a,A,B,x,y,z,X,Y,Z,W)
matrix([
[1,X+x+A+a,(x+A+a)*Y+(A+a)*y+a*B,((A+a)*y+a*B)*Z+a*B*z,a*B*z*W],
[0,1,Y+y+B,(y+B)*Z+B*z,B*z*W],
[0,0,1,Z+z,z*W],
[0,0,0,1,W],
[0,0,0,0,1]])
end:				

Coords:=matrix([
	[1,a,b,0,0],
	[0,0,1,c,d]]);

M.0:=Pos(0,0,0,0,0,0,0,0,0,0):

for inst from 1 to 3000 do

for j from 1 to 4 do
Params:=[]:
SEED:=_seed:

for i from 1 to 10 do Params:=[Params[],die()]: od:
M.j:=multiply(Pos(Params[]),M.(j-1));
od:

for j from 1 to 4 do
F.j:=det(stackmatrix(submatrix(M.j,[1,2,3],[1,2,3,4,5]),Coords)):
od:
N:=nops(realroot(gbasis([F1,F2,F3,F4],wdeg([1,6,6,6],[a,b,c,d]))[1],1/10)):

if (N<5) then lprint(`realroots = `,N,` inst=`,inst,` seed=`,SEED); fi:
od:

time();

print(_seed);

quit;