//last change: 13.02.2001 (Eric Westenberger)
///////////////////////////////////////////////////////////////////////////////
version="$Id: solve.lib,v 1.24 2001-12-11 09:42:15 pohl 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
laguerre_solve(p,[..]); find all roots of univariate polynomial p
mp_res_mat(i,[..]);     multipolynomial resultant matrix of ideal i
interpolate(p,v,d);     interpolate poly from evaluation points i and results j
fglm_solve(i,[..]);     find roots of 0-dim. ideal using FGLM and lex_solve
lex_solve(i,p,[..]);    find roots of reduced lexicographic standard basis
triangLf_solve(l,[..]); find roots using triangular sys. (factorizing Lazard)
triangM_solve(l,[..]);  find roots of given triangular system (Moeller)
triangL_solve(l,[..]);  find roots using triangular system (Lazard)
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, p] ); i = ideal, k, p = integers
 @format
         k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky,
         k=1: use resultant matrix of Macaulay which works only for
              homogeneous ideals,
         p>0: defines precision of the long floats for internal computation
              if the basering is not complex (in decimal digits),
              (default: k, p = 0, 30)
 @end format
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,
 @format
         the result is
         of type string: if the basering is not complex,
         of type number: otherwise.
 @end format
EXAMPLE: example ures_solve; shows an example
"
{
    int typ=0;// defaults
    int prec=30;

    if ( size(#) > 0 )
    {
        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(#) > 1 )
    {
        prec= #[2];
        if ( prec < 8 )
        {
            prec = 8;
        }
    }

    return(uressolve(gls,typ,prec,1));
    // the last nonzero parameter gives an extra run of
    // Laguerre's algorithm leading to better results
}
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,0,16);
    // 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,0,10);
    // 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(f [, l, m, n, s] ); f = polynomial,@*
         l, m, n, s = integers (control parameters of the method)
 @format
         l: precision of internal computation in decimal digits ( l >=8 )
            only if the basering is not complex or complex with smaller
            precision,
         m: precision of output in digits ( 4 <= m <= l ), if basering is not
            the ring of complex numbers;
            (if s != 0, m defines also the precision of the residual error
            i.e. | f(root) | should be <= 0.1^m ).
         n: control of multiplicity of roots or of splitting of f into
            squarefree factors
            n < 0, no split of f (good, if all roots are simple)
            n >= 0, try to split
            n = 0, return only different roots
            n > 0, find all roots (with multiplicity)
         s: s != 0, returns ERROR if  | f(root) | > 0.1^m
         ( default: l, m, n, s = 30, 10, 1, 0 )
 @end format
ASSUME:  f is an univariate polynomial
RETURN:  list of (complex) roots of the polynomial f, depending on n.
 @format
         The result is
         of type string: if the basering is not complex,
         of type number: otherwise.
 @end format
NOTE:    If printlevel >0: displays comments ( default = 0 ).
         If s != 0 and if the procrdure stops with ERROR, try a higher
         internal precision m.
EXAMPLE: example laguerre_solve; shows an example
"
{
  int OLD_COMPLEX=0;
  int iv=checkv(f);
  if(iv==0){ERROR("Wrong polynomial!");}
  poly v=var(iv);

  int numberprec=30;// set the control
  int solutionprec=10;
  int splitcontrol=1;
  int rootcheck=0;
  if(size(#)>0){numberprec=#[1];
    if(numberprec<8){numberprec=8;}}
  if(size(#)>1){solutionprec=#[2];
    if(solutionprec<4){solutionprec=4;}
    if(solutionprec>numberprec){solutionprec=numberprec;}}
  if(size(#)>2){splitcontrol=#[3];}
  if(size(#)>3){rootcheck=#[4];}
  int prot=printlevel-voice+2;
  int ringprec=0;

  poly p=divzero(f,iv);// divide out zeros as solution
  int iz=deg(f)-deg(p);
  if(prot!=0)
  {
    string pout;
    string nl=newline;
    pout="//BEGIN laguerre_solve";
    if(iz!=0){pout=pout+nl+"//zeros: divide out "+string(iz);}
    dbprint(prot,pout);
  }
  string ss,tt,oo;
  ss="";oo=ss;
  if(npars(basering)==1)
  {
    if(OLD_COMPLEX)
    {
      tt="1,"+string(par(1));
      if(tt==charstr(basering))
      {ss=tt;ringprec=system("getPrecDigits");}
    }
    else
    {
      tt=charstr(basering);
      if(size(tt)>7)
      {
        if(tt[1..7]=="complex")
        {
          ss=tt;
          ringprec=system("getPrecDigits");
        }
      }
    }
  }

  list roots,simple;
  if(deg(p)==0)// only zeros
  {
    roots=addzero(roots,ss,iz,splitcontrol);
    if(prot!=0){dbprint(prot,"//END laguerre_solve");}
    return(roots);
  }

  if(prot!=0)// more informations
  {
    pout="//control: complex ring with precision "+string(numberprec);
    if(size(ss)==0){pout=pout+nl+
        "//         basering not complex, hence solutiontype string";
    if(solutionprec<numberprec){pout=pout+nl+
        "//         with precision "+string(solutionprec);}}
    if(splitcontrol<0){pout=pout+nl+ "//       no spliting";}
    if(splitcontrol==0){pout=pout+nl+"//       output without multiple roots";}
    if(rootcheck){pout=pout+nl+
        "//         check roots with precision "+string(solutionprec);}
    dbprint(prot,pout);
  }

  def rn = basering;// set the complex groundfield
  if (ringprec<numberprec)
  {
    tt="ring lagc=(complex,"+string(numberprec)+","+string(numberprec)+
       "),"+string(var(iv))+",lp;";
    execute(tt);
    poly p=imap(rn,p);
    poly v=var(1);
  }
  int ima=0;
  if(size(ss)!=0){ima=checkim(p);}
  number prc=0.1;// set precision of the solution
  prc=prc^solutionprec;
  if(prot!=0)
  {
    if(ringprec<numberprec){pout="//working in:  "+tt;}
    if((size(ss)!=0)&&(ima!=0)){pout=pout+nl+
        "//         polynomial has complex coefficients";}
    dbprint(prot,pout);
  }

  int i1=1;
  int i2=1;
  ideal SPLIT=p;
  if(splitcontrol>=0)// spliting
  {
    if(prot!=0){dbprint(prot,"//split in working ring:");}
    SPLIT=splitsqrfree(p,v);
    i1=size(SPLIT);
    if((i1==1)&&(charstr(rn)=="0"))
    {
      if(prot!=0){dbprint(prot,"//split exact in basering:");}
      setring rn;
      if(v>1)
      {
        ideal SQQQQ=splitsqrfree(p,v);
        setring lagc;
        SPLIT=imap(rn,SQQQQ);
      }
      else
      {
        oo="ring exa=0,"+string(var(1))+",lp;";
        execute(oo);
        ideal SQQQQ=splitsqrfree(imap(rn,p),var(1));
        setring lagc;
        SPLIT=imap(exa,SQQQQ);
        kill exa;
      }
      i1=size(SPLIT);
    }
    if(prot!=0)
    {
      if(i1>1)
      {
        int i3=deg(SPLIT[1]);
        pout="//results of split(the squarefree factors):";
        if(i3>0){pout=pout+nl+
           "//  multiplicity "+string(i2)+", degree "+string(i3);}
        while(i2<i1)
        {
          i2++;
          i3=deg(SPLIT[i2]);
          if(i3>0){pout=pout+nl+
             "//  multiplicity "+string(i2)+", degree "+string(i3);}
        }
        dbprint(prot,pout);
        i2=1;
      }
      else
      {
        if(charstr(rn)=="0"){dbprint(prot,"// polynomial is squarefree");}
        else{dbprint(prot,"// split without result");}
      }
    }
  }

  p=SPLIT[1];// the first part
  if(deg(p)>0)
  {
    roots=laguerre(p,numberprec,1);// the ring is already complex, hence numberprec is dummy
    if((size(roots)==0)||(string(roots[1])=="0")){ERROR("laguerre: no roots found");}
    if(rootcheck){checkroots(p,v,roots,ima,prc);}
  }
  while(i2<i1)
  {
    i2++;
    p=SPLIT[i2];// the part with multiplicity i2
    if(deg(p)>0)
    {
      simple=laguerre(p,numberprec,1);
      if((size(simple)==0)||(string(simple[1])=="0")){ERROR("laguerre: no roots found");}
      if(rootcheck){checkroots(p,v,simple,ima,prc);}
      if(splitcontrol==0)// no multiple roots
      {
        roots=roots+simple;
      }
      else// multiple roots
      {
        roots=roots+makemult(simple,i2);
      }
    }
  }

  if((solutionprec<numberprec)&&(size(ss)==0))// shorter output
  {
    oo="ring lout=(complex,"+string(solutionprec)+",1),"
    +string(var(1))+",lp;";
    execute(oo);
    list roots=imap(lagc,roots);
    roots=transroots(roots);
    if(iz>0){roots=addzero(roots,ss,iz,splitcontrol);}
    if(prot!=0){dbprint(prot,"//END laguerre_solve");}
    return(roots);
  }
  if(size(ss)==0){roots=transroots(roots);}// transform to string
  else         // or map in basering
  {
    if(ringprec<numberprec)
    {
      setring rn;
      list roots=imap(lagc,roots);
    }
  }
  if(iz>0){roots=addzero(roots,ss,iz,splitcontrol);}
  if(prot!=0){dbprint(prot,"//END laguerre_solve");}
  return(roots);
}
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;
    // 10 digits precision
    laguerre_solve(f,10);

    // Now with complex coefficients,
    // internal precision is 30 digits (default)
    printlevel=2;
    ring rsc= (real,10,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]);
}
//////////////////////////////////////////////////////////////////////////////
//                subprocedures for laguerre_solve
/*
*  if p depends only on var(i)
*       returns i
*  otherwise    0
*/
static proc checkv(poly p)
{
  int n=nvars(basering);
  int i=0;
  int v;

  while (n>0)
  {
    if ((p-subst(p,var(n),0))!=0)
    {
      i++;
      if (i>1){return(0);}
      v=n;
    }
    n--;
  }
  return(v);
}
/*
*  if p hase only real coefficients
*       returns 0
*  otherwise    1
*/
static proc checkim(poly p)
{
  poly q=p;

  while(q!=0)
  {
    if(impart(leadcoef(q))!=0){return(1);}
    q=q-lead(q);
  }
  return(0);
}
/*
*  make multiplicity m
*/
static proc makemult(list si,int m)
{
  int k0=0;
  int k1=size(si);
  int k2,k3;
  number ro;
  list msi;

  for(k2=1;k2<=k1;k2++)
  {
    ro=si[k2];
    for(k3=m;k3>0;k3--){k0++;msi[k0]=ro;}
  }
  return(msi);
}
/*
*  returns 1 for n<1
*/
static proc cmp1(number n)
{
  number r=repart(n);
  number i=impart(n);
  number c=r*r+i*i;
  if(c>1){return(1);}
  else{return(0);}
}
/*
*  exact division of polys f/g
*  (should be internal)
*/
static proc exdiv(poly f,poly g,poly v)
{
  int d1=deg(f);
  int d2=deg(g);
  poly r0=f;
  poly rf=0;
  poly h;
  number n,m;

  m=leadcoef(g);
  while ((r0!=0)&&(d1>=d2))
  {
    n=leadcoef(r0)/m;
    h=n*v^(d1-d2);
    rf=rf+h;
    r0=r0-h*g;
    d1=deg(r0);
  }
  return(cleardenom(rf));
}
/*
*  p is uniform in x
*  perform a split of p into squarefree factors
*  such that the returned ideal 'split' consists of
*  the faktors, i.e.
*    p = n * product ( split[i]^i ) , n a number
*/
proc splitsqrfree(poly p, poly x)
{
  int i=1;
  int dd=deg(p);
  int j;
  ideal h,split;
  poly high;

  if(x<1){ERROR("wrog ringorder!");}
  h=p,diff(p,x);
  h=interred(h);
  if(deg(h[1])==0){return(p);}
  high=h[1];
  split[1]=exdiv(p,high,x);
  while(1)
  {
    h=split[i],high;
    h=interred(h);
    j=deg(h[1]);
    if(j==0){return(p);}
    if(deg(h[1])==deg(split[i]))
    {
      split=split,split[i];
      split[i]=1;
    }
    else
    {
      split[i]=exdiv(split[i],h[1],x);
      split=split,h[1];
      dd=dd-deg(split[i])*i;
    }
    j=j*(i+1);
    if(j==dd){break;}
    if(j>dd){return(p);}
    high=exdiv(high,h[1],x);
    if(deg(high)==0){return(p);};
    i++;
  }
  return(split);
}
/*
*  see checkroots
*/
static proc nerr(number n,number m)
{
  int r;
  number z=0;
  number nr=repart(n);
  number ni=impart(n);
  if(nr<z){nr=z-nr;}
  if(ni<z){ni=nr-ni;}
  else{ni=nr+ni;}
  if(ni<m){r=0;}
  else{r=1;}
  return(r);
}
/*
*  returns ERROR for nerr(p(r[i]))>=pr
*/
static proc checkroots(poly p,poly v,list r,int ima,number pr)
{
  int i=0;
  int j;
  number n,m;
  ideal li;

  while(i<size(r))
  {
    i++;
    n=r[i];
    j=cmp1(n);
    if(j!=0){li[1]=v/n-1;m=1;}
    else{li[1]=v-n;m=n;}
    if((ima==0)&&(impart(n)!=0))
    {
      i++;
      n=r[i];
      if(j!=0){li[1]=li[1]*(v/n-1);}
      else{li[1]=li[1]*(v-n);m=m*n;}
    }
    attrib(li,"isSB",1);
    n=leadcoef(reduce(p,li));n=n/m;
    if(n!=0)
    {if(nerr(n,pr)!=0){ERROR("Unsufficient accuracy!");}}
  }
}
/*
*  transforms thr result to string
*/
static proc transroots(list r)
{
  int i=size(r);
  while (i>0)
  {
    r[i]=string(r[i]);
    i--;
  }
  return(r);
}
/*
* returns a poly without zeroroots
*/
static proc divzero(poly f,int iv);
{
  poly p=f;
  poly q=p;
  poly r;
  while(p==q)
  {
    q=p/var(iv);
    r=q*var(iv);
    if(r==p){p=q;}
  }
  return(p);
}
/*
*  add zeros to solution
*/
static proc addzero(list zz,string ss,int iz,int a)
{
    int i=1;
    int j=size(zz);

    if(size(ss)==0){zz[j+1]="0";}
    else{zz[j+1]=number(0);}
    if(a==0){return(zz);}
    while(i<iz)
    {
      i++;
      if(size(ss)==0){zz[j+i]="0";}
      else{zz[j+i]=number(0);}
    }
    return(zz);
}
///////////////////////////////////////////////////////////////////////////////

proc mp_res_mat( ideal i, list # )
"USAGE:   mp_res_mat(i [, k] ); i ideal, k integer,
 @format
         k=0: sparse resultant matrix of Gelfand, Kapranov and Zelevinsky,
         k=1: resultant matrix of Macaulay (k=0 is default)
 @end format
ASSUME:
 @format
         the number of elements in the input system must be
         the number of variables in the basering plus one,
         if k=1 then i must be homogeneous,
 @end format
RETURN:  module representing the multipolynomial resultant matrix
EXAMPLE: example mp_res_mat; shows an example
"
{
    int typ=0;

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

    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 of numbers, d=integer
ASSUME:
@format
         ground field K are the rational numbers,
         p and v consists of elements of the ground field K
         with the following number of elements:
         size(p) = n and size(v)=N=(d+1)^n with n equals
         the number of variables,
         p is considered as point in K^n and with
         pi=(p[1]^i,..,p[n]^i), i=1,..,N the returned
         polynomial f satisfies f(pi) = v[i]
         (p[j] != -1,0,1);
@end format
RETURN:
@format
         the unique polynomial f of degree n*d with prescribed
         values v at certain points p1,..,pN derived from p;
@end format
NOTE:    mainly useful when n=1, i.e. f is satisfying
@format
         f(p^(i-1)) = v[i], i=1..d+1
@end format
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 );
}
///////////////////////////////////////////////////////////////////////////////

static proc psubst( int d, int dd, int n, list resl,
                    ideal fi, int elem, int nv, int prec )
{
    //   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
    // resl: 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 pruefen???
                    // 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(resl[k])+");";
            }
            ps = subst(ps,var(k),resl[k]);
        }

        if ( pdebug>=2 )
        { "// 5 substituted ps="+string(ps); }

        if ( ps != 0 )
        {
            lsr= laguerre_solve( ps, prec, prec/2, 0 );
        }
        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(resl[k])+");";
                            }
                            ps = subst(ps,var(k),resl[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),
                                        prec) )
                            {
                                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++ )
            {
                resl[nv-d]= nares[i];

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

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

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

proc fglm_solve( ideal fi, list # )
"USAGE:   fglm_solve(i [, p] ); i=ideal, p=integer,
 @format
         p>0: gives precision of complex numbers
         in decimal digits (default: p=30),
 @end format
ASSUME:  the ground field has char 0.
RETURN:
 @format
         a list of complex roots of type number, the procedure
         uses a standard basis of i to determine all complex
         roots of i.
 @end format
NOTE:
 @format
         The procedure creates a ring rC with the same number
         of variables but with complex coefficients (and precision p).
 @end format
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,10);
    rlist;
}

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

proc lex_solve( ideal fi, list # )
"USAGE:   lex_solve( i[,p] ); i=ideal, p=integer,
 @format
         p>0: gives precision of complex numbers
         in decimal digits (default: p=30),
 @end format
ASSUME:
 @format
         i is a reduced lexicographical Groebner bases of a
         zero-dimensional ideal (i), sorted by increasing leading terms.
 @end format
RETURN:  nothing
CREATE:
 @format
         The procedure creates a complec ring with the same variables but
         with complex coefficients (and precision p).
         In this ring a list @code{rlist} of numbers is created,
         in which the complex roots of i are stored.
 @format
EXAMPLE: example lex_solve; shows an example
"
{
    int prec=30;

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

    if ( !defined(pdebug) )
    {
        int pdebug;
        pdebug=0;
        export pdebug;
    }

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

    // change the ground field to complex numbers
    string nrings= "ring "+orings+"C=(complex,"+string(prec)
        +"),("+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;

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

    if ( !defined(rlist) )
    {
        list rlist;
        export rlist;
    }

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

    li= laguerre_solve(fi[1],prec,prec/2,0);
    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@++;
        resl[nv]= li[j];
        psubst( 1, 2, nv-1, resl, fi, idelem, nv, prec );
    }

    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),10);
    rlist;
}

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

proc triangLf_solve( ideal fi, list # )
"USAGE:   triangLf_solve(i [, p] ); i ideal, p integer,
p>0: gives precision of complex numbers in digits,@*
(default: p=30)
ASSUME:  the ground field has char 0;@*
         i zero-dimensional ideal
RETURN:  nothing
CREATE:  The procedure creates a ring rC with the same number of variables but
         with complex coefficients (and precision p).@*
         In rC a list @code{rlist} of numbers is created, in which the complex
         roots of i are stored.@*
         The proc 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.
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,10);
    rlist;
}

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

proc triangM_solve( ideal fi, list # )
"USAGE:   triangM_solve(i [, p ] ); i=ideal, p=integer,
p>0: gives precision of complex numbers in digits,@*
(default: p=30)
ASSUME:  the ground field has char 0;@*
         i zero-dimensional ideal
RETURN:  nothing
CREATE:  The procedure creates a ring rC with the same number of variables but
         with complex coefficients (and precision p).@*
         In rC a list @code{rlist} of numbers is created, in which the complex
         roots of i are stored.@*
         The proc 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.
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,10);
    rlist;
}

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

proc triangL_solve( ideal fi, list # )
"USAGE:   triangL_solve(i [, p] ); i=ideal, p=integer,@*
p>0: gives precision of complex numbers in digits,@*
(default: p=30)
ASSUME:  the ground field has char 0;@*
         i zero-dimensional ideal
RETURN:  nothing
CREATE:  The procedure creates a ring rC with the same number of variables but
         with complex coefficients (and precision p).@*
         In rC a list @code{rlist} of numbers is created, in which the complex
         roots of i are stored.@*
         The proc 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.
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,10);
    rlist;
}


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

proc triang_solve( list lfi, int prec, list # )
"USAGE:   triang_solve(l,p [, d] ); l=list, p,d=integers,@*
l a list of finitely many triangular systems, such that
the union of their varieties equals the variety of the
initial ideal.@*
p>0: gives precision of complex numbers in digits,@*
d>0: gives precision (1<d<p) for near-zero-determination,@*
(default: d=1/2*p)

ASSUME:  the ground field has char 0;@*
         l was computed using Algorithm of Lazard or Algorithm of Moeller
         (see triang.lib).
RETURN:  nothing
CREATE:  The procedure creates a ring rC with the same number of variables but
         with complex coefficients (and precision p).@*
         In rC a list @code{rlist} of numbers is created, in which the complex
         roots of i are stored.@*
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"));
    }

    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,myCompDigits/2,0);
        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, myCompDigits );
        }
    }

    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)),10);
    rlist;
}

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

proc pcheck( ideal fi, list sols, list # )
"USAGE:   pcheck(i,l [, d] ); i=ideal, l=list, d=integer,@*
d>0: precision in digits for near-zero determination
ASSUME: the ground field has char 0;@*
        l is a list of numbers
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),10);
    rlist;
    ideal s=imap(r,s);
    pcheck(s,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: ***