source: git/Singular/LIB/solve.lib @ 7708934

///////////////////////////////////////////////////////////////////////////////
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: ***