//search.singular
//
//	This SINGULAR script was used to study feedback laws for m=2, 
//	p=4, n=8 real linear systems.

//	It generates 250 systems, and creates a MAPLE file (which must
//	be edited, see below) that determines the number of real
//	feedback laws for each system.

//	First, it generates 10 pairs [D(s):N(S)] with integral
//	polynomial entries representing systems, and for each, generates
//	25 characteristic polynomials whose roots are integers.
//
//	For each of the 250, it creates the ideal defining all possible
//	feedback laws, and computes an elimination Gr\"obner basis.
//	
//	The eliminant is written to the MAPLE file, which determines the
//	number of real roots.   The nature of the elimination ideal is
//	that the number of real solutions is determined solely by the
//	number of real roots of the eliminant.
//
//	The MAPLE file creates a 10 x 25 matrix which records the
//	numbers of real roots of the eliminant, hence the number of real
//	feedback laws.
//
//	To edit the maple file, remove all occurrences of `>' and `.'.
//	Then the    SINGULAR output: x14+2345x13+...       must be 
//	changed to      MAPLE input: x^14 + 23*x^13 +...
//
//	This may be accomplished using a global change x |--> *x^,
//	However, the linear terms will look like either x^+ or x^-,
//      and this must be corrected.   Also, in the unlikely event that
//      some term has coefficient 1, there will be a stray *x, and the 
//      `*' must be removed.
//

int t=timer;

ring S = 0,(a,b,c,d,e,f,g,h,s), (dp(7),dp(2));

ideal I;
ideal G;
poly Det;
poly Eliminant;
matrix F[6][6];
matrix Co[2][6];
matrix Li[2][6];
matrix Qu[2][6];
matrix Cu[2][6];
matrix Y[25][8];
matrix X[25][8];
int j;
int ind;


print(" ");
print("interface(quiet=true): ");
print("with(linalg): ");
print("readlib( realroot ):");
print("NumberReal:= matrix(10,25,0): ");
print("ii:=1:");
print("jj:=0:");
print(" ");

//
//	These matrices Co, Li, Qu, Cu are the matrices of coefficients
//	for the Consant, Linear, Quadratic, and Cubic terms in the
//	matrix  F(s):= [D(s): N(s)], which represents the system.
//
//	We used several different distributions during the course of our
//	search.   The ranges below are an example of one.
//	For example, random(70,2,6), means choose a 2X6 matrix whose 
//	entries range from -70 to +70.
//

for (ind = 1; ind<=10; ind=ind+1)
{
 Co = random(70,2,6);
 Li = random(35,2,6);
 Qu = random(14,2,6);
 Cu = random(7,2,6);

for (int i = 1; i<=6; i=i+1)
{
F[1,i] = Co[1,i]+ s*Li[1,i]+s*s*Qu[1,i]+s*s*s*Cu[1,i];
F[2,i] = Co[2,i]+ s*Li[2,i]+s*s*Qu[2,i]+s*s*s*Cu[2,i];
}
F[1,1] = F[1,1]+ s*s*s*s;
F[2,2] = F[2,2]+ s*s*s*s;

//
//	The entries below are the coordinates of the potential feedback
//	laws 
//

F[3,1]=a;
F[4,1]=b;
F[5,1]=c;
F[6,1]=d;
F[3,2]=e;
F[4,2]=f;
F[5,2]=g;
F[6,2]=h;

//
//	This is the Identity matrix, embedded in F(s).
//

F[3,3]=1;
F[4,4]=1;
F[5,5]=1;
F[6,6]=1;

Det=det(F);

//
//	The poles are chosen in the following manner (several diffferent
//	methods were used for different parts of the search).   That
//	displayed below is one example.
//
//	Here, we choose a 25X8 matrix Y filled with numbers chosen from 
//      {-1, 0, 1}


Y=random(1,25,8);

//
//	Then the poles are numbers that begin with 0, and the successive
//	differences are the corresponding entry of Y, plus 2.  That is,
//	they range from 1 to 3.
//



for (j = 1; j<=25; j=j+1 )
 {
  X[j,1]=0;
  for ( i = 1; i<=7; i=i+1 )
   {
    X[j,i+1]= X[j,i]+2+Y[j,i];
   }
 }

//	This creates the ideal


for (j = 1; j<=25; j=j+1 )
{
    I =  subst(Det,s,X[j,1]), subst(Det,s,X[j,2]),subst(Det,s,X[j,3]),
	subst(Det,s,X[j,4]),subst(Det,s,X[j,5]),subst(Det,s,X[j,6]),
	 subst(Det,s,X[j,7]),subst(Det,s,X[j,8]);


G = std(I);
Eliminant = G[1];

//
//	This is the heart of teh MAPLE code for determiing the number of
//	real feedback laws
/m


print("jj:=jj+1: ");
print(" ");
print("eliminant :=");
print(Eliminant);
print(":");
print("NumberReal[ii,jj] := nops(realroot( elimiant, 1/10)):");
print(" ");
print("############################################################");
print("#                                                          # ");
print("#          Next One                                        #");
print("#                                                          # ");
print("############################################################");
print(" ");

}
print("ii:= ii+1:");
print("jj:=0: ");

}
print("evalm(NunberReal); ");
print("time(); ");
print("quit; ");
timer-t;
quit;