///////////////////////////////////////////////////////////////////////////////
version="$Id: solve.lib,v 1.18 2000-12-22 14:48:14 greuel Exp $";
category="Symbolic-numerical solving";
info="
LIBRARY: solve.lib     Complex Solving of Polynomial Systems
AUTHOR:  Moritz Wenk,  email: wenk@mathematik.uni-kl.de

PROCEDURES:
ures_solve(i,..);       find all roots of 0-dimensional ideal i with resultants
mp_res_mat(i,..);       multipolynomial resultant matrix of ideal i
laguerre_solve(p,..);   find all roots of univariate polynom p
interpolate(i,j,d);     interpolate poly from evaluation points i and results j
fglm_solve(i,p,...);    find roots of 0-dim. ideal using FGLM and lex_solve
triangL_solve(l,p,...); find roots using triangular system (Lazard) 
triangLf_solve(l,p,..); find roots using triangular sys. (factorizing Lazard)
triangM_solve(l,p,...); find roots of given triangular system (Moeller)
lex_solve(i,p,...);     find roots of reduced lexicographic standard basis
triang_solve(l,p,...);  find roots of given triangular system
pcheck(i,l,...);        checks if elements (numbers) of l are roots of ideal i
";

LIB "triang.lib";    // needed for triang*_solve

///////////////////////////////////////////////////////////////////////////////

proc ures_solve( ideal gls, list # )
"USAGE:   ures_solve(i[,k,l,m]); i ideal, k,l,m integers
         k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky
         k=1: use resultant matrix of Macaulay (k=0 is default)
         l>0: defines precision of fractional part if ground field is Q
         m=0,1,2: number of iterations for approximation of roots (default=2)
ASSUME:  i is a zerodimensional ideal with
         nvars(basering) = ncols(i) = number of vars actually occuring in i
RETURN:  list of all (complex) roots of the polynomial system i = 0,
         of type number if the ground field is the complex numbers,
         of type string if the ground field is the rational or real numbers
EXAMPLE: example ures_solve; shows an example
"
{
    int typ=0;
    int polish=2;
    int prec=30;

    if ( size(#) >= 1 )
    {
	typ= #[1];
	if ( typ < 0 || typ > 1 )
	{
	    ERROR("Valid values for second parameter k are:
                   0: use sparse Resultant (default)
                   1: use Macaulay Resultant");
	}
    }
    if ( size(#) >= 2 )
    {
	prec= #[2];
	if ( prec == 0 ) { prec = 30; }
	if ( prec < 0 )
	{
	    ERROR("Third parameter l must be positive!");
	}    }
    if ( size(#) >= 3 )
    {
	polish= #[3];
	if ( polish < 0 || polish > 2 )
	{
	    ERROR("Valid values for fourth parameter m are:
                   0,1,2: number of iterations for approximation of roots");
	}
    }

    if ( size(#) > 3 )
    {
	ERROR("only three parameters allowed!");
    }

    return(uressolve(gls,typ,prec,polish));

}
example
{
    "EXAMPLE:";echo=2;
    // compute the intersection points of two curves
    ring rsq = 0,(x,y),lp;
    ideal gls=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    ures_solve(gls);
    // result is a list (x,y)-coordinates as strings

    // now with complex coefficient field, precision is 20 digits
    ring rsc= (real,20,I),(x,y),lp;
    ideal i = (2+3*I)*x2 + (0.35+I*45.0e-2)*y2 - 8, x2 + xy + (42.7)*y2;
    list l= ures_solve(i);
    // result is a list of (x,y)-coordinates of complex numbers
    l;
    // check the result
    subst(subst(i[1],x,l[1][1]),y,l[1][2]);
}
///////////////////////////////////////////////////////////////////////////////

proc laguerre_solve( poly f, list # )
"USAGE:   laguerre_solve( p[,l,m]); f poly, l,m integers
         l>0: defines precision of fractional part if ground field is Q
         m=0,1,2: number of iterations for approximation of roots (default=2)
ASSUME:  p is an univariate polynom
RETURN:  list of all (complex) roots of the polynomial p;
         of type number if the ground field is the complex numbers,
         of type string if the ground field is the rational or real numbers
EXAMPLE: example laguerre_solve; shows an example
"
{
    int polish=2;
    int prec=30;

    if ( size(#) >= 1 )
    {
	prec= #[1];
	if ( prec == 0 ) { prec = 30; }
	if ( prec < 0 )
	{
	    ERROR("Fisrt parameter must be positive!");
	}
    }
    if ( size(#) >= 2 )
    {
	polish= #[2];
	if ( polish < 0 || polish > 2 )
	{
	    ERROR("Valid values for third parameter are:
                   0,1,2: number of iterations for approximation of roots");
	}
    }
    if ( size(#) > 2 )
    {
	ERROR("only two parameters allowed!");
    }

    return(laguerre(f,prec,polish));

}
example
{
    "EXAMPLE:";echo=2;
    // Find all roots of an univariate polynomial using Laguerre's method:
    ring rs1= 0,(x,y),lp;
    poly f = 15x5 + x3 + x2 - 10;
    laguerre_solve(f);

    // Now with 10 digits precision:
    laguerre_solve(f,10);

    // Now with complex coefficients, precision is 20 digits:
    ring rsc= (real,20,I),x,lp;
    poly f = (15.4+I*5)*x^5 + (25.0e-2+I*2)*x^3 + x2 - 10*I;
    list l = laguerre_solve(f);
    l;
    // check result, value of substituted poly should be near to zero
    subst(f,x,l[1]);
    subst(f,x,l[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc mp_res_mat( ideal i, list # )
"USAGE:   mp_res_mat(i[,k]); i ideal, k integer
         k=0: sparse resultant matrix of Gelfand, Kapranov and Zelevinsky
         k=1: resultant matrix of Macaulay (k=0 is default)
ASSUME:  nvars(basering) = ncols(i)-1 = number of vars actually occuring in i,
         for k=1 i must be homogeneous
RETURN:  module representing the multipolynomial resultant matrix
EXAMPLE: example mp_res_mat; shows an example
"
{
    int typ=2;

    if ( size(#) == 1 )
    {
	typ= #[1];
	if ( typ < 0 || typ > 1 )
	{
	    ERROR("Valid values for third parameter are:
                   0: sparse resultant (default)
                   1: Macaulay resultant");
	}
    }
    if ( size(#) > 1 )
    {
	ERROR("only two parameters allowed!");
    }

    return(mpresmat(i,typ));

}
example
{
    "EXAMPLE:";echo=2;
    // compute resultant matrix in ring with parameters (sparse resultant matrix)
    ring rsq= (0,u0,u1,u2),(x1,x2),lp;
    ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
    module m = mp_res_mat(i);
    print(m);
    // computing sparse resultant
    det(m);

    // compute resultant matrix (Macaulay resultant matrix)
    ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
    ideal h=  homog(imap(rsq,i),x0);
    h;
  
    module m = mp_res_mat(h,1);
    print(m);
    // computing Macaulay resultant (should be the same as above!)
    det(m);

    // compute numerical sparse resultant matrix
    setring rsq;
    ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
    module mn = mp_res_mat(ir);
    print(mn);
    // computing sparse resultant
    det(mn);
}
///////////////////////////////////////////////////////////////////////////////

proc interpolate( ideal p, ideal w, int d )
"USAGE:   interpolate(p,v,d); p,v ideals, d int
ASSUME:  ground field K is the rational numbers,
         p and v consist of numbers of the ground filed K, p must have
         n elements, v must have N=(d+1)^n elements where n=nvars(basering)
         and d=deg(f) (f is the unknown polynomial in K[x1,...,xn])
COMPUTE: polynomial f with values given by v at points p1,..,pN derived from p;
         more precisely: consider p as point in K^n and v as N elements in K,
         let p1,..,pN be the points in K^n obtained by evaluating all monomials
         of degree 0,1,...,N at p in lexicographical order,
         then the procedure computes the polynomial f satisfying f(pi) = v[i]
RETURN:  polynomial f of degree d
NOTE:    mainly useful for n=1, with f satisfying f(p^(i-1)) = v[i], i=1..d+1
EXAMPLE: example interpolate; shows an example
"
{
    return(vandermonde(p,w,d));
}
example
{
    "EXAMPLE:";  echo=2;
    ring r1 = 0,(x),lp;
    // determine f with deg(f) = 4 and
    // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
    ideal v=16,0,11376,1046880,85949136;
    interpolate( 3, v, 4 );

    ring r2 = 0,(x,y),dp;
    // determine f with deg(f) = 3 and
    // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15)
    // valuation point (2,3)
    ideal p = 2,3;
    ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
    poly ip= interpolate( p,v,3 );
    ip;
    // compute poly at point 2,3, result must be 2
    subst(subst(ip,x,2),y,3);
    // compute poly at point 2^15,3^15, result must be 16
    subst(subst(ip,x,2^15),y,3^15);
}
///////////////////////////////////////////////////////////////////////////////

static proc psubst( int d, int dd, int n, list res,
		    ideal fi, int elem, int nv )
{
    //   nv: number of ring variables         (fixed value)
    // elem: number of elements in ideal fi   (fixed value)
    //   fi: input ideal                      (fixed value)
    //   rl: output list of roots
    //  res: actual list of roots 
    //    n: 
    //   dd: actual element of fi
    //    d: actual variable  

    int olddd=dd;

    if ( pdebug>=1 ) 
    {"// 0 step "+string(dd)+" of "+string(elem);}
  
    if ( dd <= elem ) 
    {
	int loop = 1;
	int k;
	list lsr,lh;
	poly ps;
	int thedd;

	if ( pdebug>=1 ) 
	{"// 1 dd = "+string(dd);}

	thedd=0;
	while ( (dd+1 <= elem) && loop ) 
	{
	    ps= fi[dd+1];

	    if ( n-1 > 0 ) 
	    {
		if ( pdebug>=2 ) 
		{
		    "// 2 ps=fi["+string(dd+1)+"]"+" size="
			+string(size(coeffs(ps,var(n-1))))
			+"  leadexp(ps)="+string(leadexp(ps));
		}
		if ( size(coeffs(ps,var(n-1))) == 1 ) 
		{
		    dd++;
		    // hier Leading-Exponent puefen???
		    // oder ist das Poly immer als letztes in der Liste?!?
		    // leadexp(ps)
		}
		else
		{
		    loop=0;
		}
	    }
	    else
	    {
		if ( pdebug>=2 ) 
		{
		    "// 2 ps=fi["+string(dd+1)+"]"+"  leadexp(ps)="
			+string(leadexp(ps));
		}
		dd++;
	    }
	} 
	thedd=dd;
	ps= fi[thedd];

	if ( pdebug>=1 ) 
	{
	    "// 3    fi["+string(thedd-1)+"]"+"  leadexp(fi[thedd-1])="
		+string(leadexp(fi[thedd-1]));
	    "// 3 ps=fi["+string(thedd)+"]"+"  leadexp(ps)="
		+string(leadexp(ps));
	}

	for ( k= nv; k > nv-d; k-- )
	{
	    if ( pdebug>=2 ) 
	    { 
		"// 4 subst(fi["+string(thedd)+"],"
		    +string(var(k))+","+string(res[k])+");"; 
	    }
	    ps = subst(ps,var(k),res[k]);
	}
	      
	if ( pdebug>=2 ) 
	{ "// 5 substituted ps="+string(ps); }

	if ( ps != 0 ) 
	{
	    lsr= laguerre_solve( ps );
	}
	else
	{
	    if ( pdebug>=1 ) 
	    { "// 30 ps == 0, thats not cool..."; }
	    lsr=@ln; // lsr=number(0);
	}

	if ( pdebug>=1 ) 
	{ "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]"; }
	  
	if ( size(lsr) > 1 ) 
	{
	    if ( pdebug>=1 ) 
	    {
		"// 10 checking roots found before, range "
		    +string(dd-olddd)+" -- "+string(dd);
		"// 10 thedd = "+string(thedd);
	    }
	      
	    int i,j,l;
	    int ls=size(lsr);
	    int lss;
	    poly pss;
	    list nares;
	    int rroot;
	    int nares_size;

	      
	    for ( i = 1; i <= ls; i++ ) // lsr[1..ls]
	    {
		rroot=1;

		if ( pdebug>=2 ) 
		{"// 13 root lsr["+string(i)+"] = "+string(lsr[i]);}
		for ( l = 0; l <= dd-olddd; l++ ) 
		{
		    if ( l+olddd != thedd )
		    {
			if ( pdebug>=2 ) 
			{"// 11 checking ideal element "+string(l+olddd);}
			ps=fi[l+olddd];
			if ( pdebug>=3 ) 
			{"// 14 ps=fi["+string(l+olddd)+"]";}
			for ( k= nv; k > nv-d; k-- )
			{
			    if ( pdebug>=3 ) 
			    {
				"// 11 subst(fi["+string(olddd+l)+"],"
				    +string(var(k))+","+string(res[k])+");";
			    }
			    ps = subst(ps,var(k),res[k]);
			  
			}
			
			pss=subst(ps,var(k),lsr[i]); // k=nv-d
			if ( pdebug>=3 ) 
			{ "// 15 0 == "+string(pss); }
			if ( pss != 0 ) 
			{
			    if ( system("complexNearZero",
					leadcoef(pss),
					myCompDigits) ) 
			    {
				if ( pdebug>=2 ) 
				{ "// 16 root "+string(i)+" is a real root"; }
			    }
			    else
			    {
				if ( pdebug>=2 ) 
				{ "// 17 0 == "+string(pss); }
				rroot=0;
			    }
			}
		      
		    }
		}
		
		if ( rroot == 1 ) // add root to list ?
		{
		    if ( size(nares) > 0 )
		    {
			nares=nares[1..size(nares)],lsr[i];
		    }
		    else
		    {
			nares=lsr[i];
		    }
		    if ( pdebug>=2 ) 
		    { "// 18 added root to list nares"; }
		}
	    }

	    nares_size=size(nares);
	    if ( nares_size == 0 ) 
	    {
		"Numerical problem: No root found...";
		"Output may be incorrect!";
		nares=@ln;
	    }
	      
	    if ( pdebug>=1 ) 
	    { "// 20 found <"+string(size(nares))+"> roots"; }

	    for ( i= 1; i <= nares_size; i++ )
	    {
		res[nv-d]= nares[i];

		if ( dd < elem ) 
		{
		    if ( i > 1 ) 
		    { 
			rn@++; 
		    }
		    psubst( d+1, dd+1, n-1, res, fi, elem, nv );
		}
		else
		{
		    if ( pdebug>=1 ) 
		    {"// 30_1 <"+string(rn@)+"> "+string(size(res))+" <-----";}
		    if ( pdebug>=2 ) 
		    { res; }
		    rlist[rn@]=res;
		}
	    }
	}
	else
	{
	    if ( pdebug>=2 ) 
	    { "// 21 found root to be: "+string(lsr[1]); }
	    res[nv-d]= lsr[1];

	    if ( dd < elem ) 
	    {
		psubst( d+1, dd+1, n-1, res, fi, elem, nv );
	    }
	    else
	    {
		if ( pdebug>=1 ) 
		{ "// 30_2 <"+string(rn@)+"> "+string(size(res))+" <-----";}
		if ( pdebug>=2 ) 
		{ res; }
		rlist[rn@]=res;
	    }
	}
    }
}

///////////////////////////////////////////////////////////////////////////////

proc fglm_solve( ideal fi, list # )
"USAGE:   fglm_solve( i [, d ] ); i ideal, p integer.
         p gives precision of complex numbers in digits,
         uses a standard basis computed from fi to determine
         recursively all complex roots of fi
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example fglm_solve; shows an example
"
{
    int prec=30;

    if ( size(#)>=1  && typeof(#[1])=="int")
    {
	prec=#[1];
    }

    lex_solve(stdfglm(fi),prec);
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    fglm_solve(s);
    rlist;
}

///////////////////////////////////////////////////////////////////////////////

proc lex_solve( ideal fi, int prec, list # )
"USAGE:   lex_solve( i, d [, l] ); i ideal, p,d,l integers.
         i a reduced lexicographical Groebner bases of a zero-dimensional 
         ideal (i), sorted by increasing leading terms. 
         p gives precision of complex numbers in digits,
         d gives precision (1<d<p) for near-zero-determination (default:1/2*p),
         l gives debug level (-1: no output,..., 3: lots of output, 0: default)
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example lex_solve; shows an example
"
{
    if ( !defined(pdebug) ) 
    {
	int pdebug;
	export pdebug;
    }
    pdebug=0;

    string orings= nameof(basering);
    def oring= basering;

    // change the ground field to complex numbers  
    string nrings= "ring "+orings+"C=(real,"+string(prec)
	+",I),("+varstr(basering)+"),lp;";
    execute(nrings);

    if ( pdebug>=0 )
    { "// name of new ring: "+string(nameof(basering));}

    // map fi from old to new ring
    ideal fi= imap(oring,fi);

    // list with entry 0 (number)
    number nn=0;
    if ( !defined(@ln) ) 
    {
	list @ln;
	export @ln;
    }
    @ln=nn;

    // set number of digits for zero-comparison of roots
    if ( !defined(myCompDigits) )
    { 
	int myCompDigits;
	export  myCompDigits;
    }
    if ( size(#)>=1  && typeof(#[1])=="int")
    {
	myCompDigits=#[1];
    }
    else
    {
	myCompDigits=(system("getPrecDigits")/2);
    }

    if ( pdebug>=1 )
    {"// myCompDigits="+string(myCompDigits);}

    int idelem= size(fi);
    int nv= nvars(basering);
    int i,j,k,lis;
    list res,li;

    if ( !defined(rlist) )
    {
	list rlist;
	export rlist;
    }
    
    if ( !defined(rn@) )
    {
	int rn@;
	export rn@;
    }
    rn@=0;

    li= laguerre_solve(fi[1]);
    lis= size(li);

    if ( pdebug>=1 )
    {"// laguerre found roots: "+string(size(li));}
    
    for ( j= 1; j <= lis; j++ )
    {
	if ( pdebug>=1 )
	{"// root "+string(j);}
	rn@++;
	res[nv]= li[j];
	psubst( 1, 2, nv-1, res, fi, idelem, nv );
    }
    
    if ( pdebug>=0 )
    {"// list of roots: "+nameof(rlist);}

    // keep the ring and exit
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    lex_solve(stdfglm(s),30);
    rlist;
}

///////////////////////////////////////////////////////////////////////////////

proc triangLf_solve( ideal fi, list # )
"USAGE:   triangLf_solve( i [, d ] ); i ideal, p integer.
         i zero-dimensional ideal
         p gives precision of complex numbers in digits,
         uses a triangular system (Lazard's Algorithm with factorization) 
         computed from a standard basis to determine recursively all complex 
         roots with Laguerre's algorithm of input ideal i
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example triangLf_solve; shows an example
"
{
    int prec=30;

    if ( size(#)>=1  && typeof(#[1])=="int")
    {
	prec=#[1];
    }

    triang_solve(triangLfak(stdfglm(fi)),prec);
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    triangLf_solve(s);
    rlist;
}

///////////////////////////////////////////////////////////////////////////////

proc triangM_solve( ideal fi, list # )
"USAGE:   triangM_solve( i [, d ] ); i ideal, p integer.
         i zero-dimensional ideal
         p gives precision of complex numbers in digits,
         uses a triangular system (Moellers Algorithm) computed from a standard
         basis to determine recursively all complex roots with Laguerre's 
         algorithm of input ideal i
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example triangM_solve; shows an example
"
{
    int prec=30;

    if ( size(#)>=1  && typeof(#[1])=="int")
    {
	prec=#[1];
    }

    triang_solve(triangM(stdfglm(fi)),prec);
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    triangM_solve(s);
    rlist;
}

///////////////////////////////////////////////////////////////////////////////

proc triangL_solve( ideal fi, list # )
"USAGE:   triangL_solve( i [, d ] ); i ideal, p integer.
         i zero-dimensional ideal
         p gives precision of complex numbers in digits,
         uses a triangular system (Lazard's Algorithm) 
         computed from a standard basis to determine recursively all complex 
         roots with Laguerre's algorithm of input ideal i
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example triangL_solve; shows an example
"
{
    int prec=30;

    if ( size(#)>=1  && typeof(#[1])=="int")
    {
	prec=#[1];
    }

    triang_solve(triangL(stdfglm(fi)),prec);
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    triangL_solve(s);
    rlist;
}


///////////////////////////////////////////////////////////////////////////////

proc triang_solve( list lfi, int prec, list # )
"USAGE:   triang_solve( i, d [, l] ); l list, p,d,l integers.
         l a list of finitely many triangular systems, such that
         the union of their varieties equals the variety of the
         initial ideal. 
         p gives precision of complex numbers in digits,
         d gives precision (1<d<p) for near-zero-determination (default:1/2*p),
         l gives debug level (-1: no output,..., 3: lots of output, 0: default)
ASSUME:  l was computed using Algorithm of Lazard or 
         Algorithm of Moeller (see triang.lib).
RETURN:  a list of complex roots of type number.
NOTE:    changes the given ring to the ring ring with complex coefficient
         field with precision p in digits.
EXAMPLE: example triang_solve; shows an example
"
{
    if ( !defined(pdebug) ) 
    {
	int pdebug;
	export pdebug;
    }
    pdebug=0;

    string orings= nameof(basering);
    def oring= basering;

    // change the ground field to complex numbers  
    string nrings= "ring "+orings+"C=(real,"+string(prec)
	+",I),("+varstr(basering)+"),lp;";
    execute(nrings);

    if ( pdebug>=0 )
    { "// name of new ring: "+string(nameof(basering));}

    // list with entry 0 (number)
    number nn=0;
    if ( !defined(@ln) ) 
    {
	list @ln;
	export @ln;
    }
    @ln=nn;

    // set number of digits for zero-comparison of roots
    if ( !defined(myCompDigits) )
    { 
	int myCompDigits;
	export  myCompDigits;
    }
    if ( size(#)>=1  && typeof(#[1])=="int" )
    {
	myCompDigits=#[1];
    }
    else
    {
	myCompDigits=(system("getPrecDigits")/2);
    }

    if ( pdebug>=1 )
    {"// myCompDigits="+string(myCompDigits);}

    int idelem;
    int nv= nvars(basering);
    int i,j,k,lis;
    list res,li;

    if ( !defined(rlist) )
    {
	list rlist;
	export rlist;
    }
    
    if ( !defined(rn@) )
    {
	int rn@;
	export rn@;
    }
    rn@=0;

    // map the list
    list lfi= imap(oring,lfi);

    int slfi= size(lfi);
    ideal fi;

    for ( i= 1; i <= slfi; i++ ) 
    {
	// map fi from old to new ring
	fi= lfi[i]; //imap(oring,lfi[i]);
	
	idelem= size(fi);

	// solve fi[1]
	li= laguerre_solve(fi[1],myCompDigits,2);
	lis= size(li);

	if ( pdebug>=1 )
	{"// laguerre found roots: "+string(size(li));}
	
	for ( j= 1; j <= lis; j++ )
	{
	    if ( pdebug>=1 )
	    {"// root "+string(j);}
	    rn@++;
	    res[nv]= li[j];
	    psubst( 1, 2, nv-1, res, fi, idelem, nv );
	}
    }

    if ( pdebug>=0 )
    {"// list of roots: "+nameof(rlist);}
	
    // keep the ring and exit
    keepring basering;
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    triang_solve(triangLfak(stdfglm(s)),30);
    rlist;
}

///////////////////////////////////////////////////////////////////////////////

proc pcheck( ideal fi, list sols, list # )
"USAGE:   pcheck( i, l [, n] ); i ideal, l list, n integer
          i system of equations, l list of numers assumend
          to be roots of i.
RETURN:  1 iff all elements of l are roots of i, else 0
EXAMPLE: example pcheck; shows an example
"
{
    int i,j,k,nnrfound;
    int ss= size(sols);
    int nv= nvars(basering);
    poly ps;
    number nn;
    int cprec=0;
    
    if ( size(#)>=1  && typeof(#[1])=="int" )
    {
	cprec=#[1];
    }
    if ( cprec <= 0 ) 
    {
	cprec=system("getPrecDigits")/2;
    } 

    nnrfound=1;
    for ( j= 1; j <= size(fi); j++ )
    {
	for ( i= 1; i <= ss; i++ )
	{
	    ps= fi[j];
	    for ( k= 1; k <= nv; k++ )
	    {
		ps = subst(ps,var(k),sols[i][k]);
	    }
	    //ps;
	    nn= leadcoef(ps);
	    if ( nn != 0 ) 
	    {
		if ( !system("complexNearZero",nn,cprec) )
		{
		    " no root: ideal["+string(j)+"]( root "
			+string(i)+"): 0 != "+string(ps);
		    nnrfound=0; 
		}
	    }
	}
    }
    return(nnrfound);
}
example
{
    "EXAMPLE:";echo=2;
    ring r = 0,(x,y),lp;
    // compute the intersection points of two curves
    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
    lex_solve(stdfglm(s),30);
    rlist;
    pcheck(rlist);
}

///////////////////////////////////////////////////////////////////////////////

proc simplexOut(list l)
"USAGE:   simpelxOut(l); l list
ASSUME:
RETURN:
EXAMPLE: example simplexOut; shows an example
"
{
  int i,j;
  matrix m= l[1];
  intvec iposv= l[3];
  int icase= l[2];

  int cols= ncols(m);
  int rows= nrows(m);

  int N= l[6];

  if ( -1 == icase )  // objective function is unbound
  {
    "objective function is unbound";
    return;
  }
  if ( 1 == icase )  // no solution satisfies the given constraints
  {
    "no solution satisfies the given constraints";
    return;
  }
  if ( -2 == icase )  // other error
  {
    "an error occurred during simplex computation!";
    return;
  }

  for ( i = 1; i <= rows; i++ )
  {
    if (i == 1) 
    {
      "z = "+string(m[1][1]);
    }
    else 
    {
      if ( iposv[i-1] <= N+1 ) 
      {
	"x"+string(i-1)+" = "+string(m[1][iposv[i-1]]);
      }
//        else 
//        {
//  	     "Y"; iposv[i-1]-N+1;
//        }
    }
  } 
}
example
{
    "EXAMPLE:";echo=2;
    ring r = (real,10),(x),lp;

    // suppose we have the 
    
    matrix sm[5][5]=(  0, 1, 1, 3,-0.5,
		     740,-1, 0,-2, 0,
		       0, 0,-2, 0, 7,
		     0.5, 0,-1, 1,-2,
		       9,-1,-1,-1,-1);

    // simplex input:
    // 1: matrix
    // 2: number of variables
    // 3: total number of constraints (#4+#5+#6)
    // 4: number of <= constraints
    // 5: number of >= constraints
    // 6: number of == constraints

    list sol= simplex(sm, 4, 4, 2, 1, 1);
    simplexOut(sol);
}

///////////////////////////////////////////////////////////////////////////////

// local Variables: ***
// c-set-style: bsd ***
// End: ***