///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Symbolic-numerical solving";
info="
LIBRARY: solve.lib     Complex Solving of Polynomial Systems
AUTHOR:  Moritz Wenk,  email: wenk@mathematik.uni-kl.de
         Wilfred Pohl, email: pohl@mathematik.uni-kl.de

PROCEDURES:
laguerre_solve(p,[..]); find all roots of univariate polynomial p
solve(i,[..]);          all roots of 0-dim. ideal i using triangular sets
ures_solve(i,[..]);     find all roots of 0-dimensional ideal i with resultants
mp_res_mat(i,[..]);     multipolynomial resultant matrix of ideal i
interpolate(p,v,d);     interpolate polynomial 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
simplexOut(l);          prints solution of simplex in nice format
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
";

LIB "triang.lib";    // needed for triang_solve

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

proc laguerre_solve( poly f, list # )
"USAGE:   laguerre_solve(f [, m, l, n, s] ); f = polynomial,@*
         m, l, n, s = integers (control parameters of the method)@*
 m: precision of output in digits ( 4 <= m), if basering is not ring of
      complex numbers;
 l: precision of internal computation in decimal digits ( l >=8 )
      only if the basering is not complex or complex with smaller precision;@*
 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 (when computing in the
      current ring)
 ( default: m, l, n, s = 8, 30, 1, 0 )
ASSUME:  f is a univariate polynomial;@*
         basering has characteristic 0 and is either complex or without
         parameters.
RETURN:  list of (complex) roots of the polynomial f, depending on n. The
         entries of the result are of type@*
          string: if the basering is not complex,@*
          number: otherwise.
NOTE:    If printlevel >0: displays comments ( default = 0 ).
         If s != 0 and if the procedure stops with ERROR, try a higher
         internal precision m.
EXAMPLE: example laguerre_solve; shows an example
"
{
  if (char(basering)!=0){ERROR("characteristic of basering not 0");}
  if ((charstr(basering)[1]=="0") and (npars(basering)!=0))
    {ERROR("basering has parameters");}
  int OLD_COMPLEX=0;
  int iv=checkv(f);  // check for variable appearing in f
  if(iv==0){ERROR("Wrong polynomial!");}
  poly v=var(iv);    // f univariate in v

  int solutionprec=8;// set the control
  int numberprec=30;
  int splitcontrol=1;
  int rootcheck=0;
  if(size(#)>0){solutionprec=#[1];if(solutionprec<4){solutionprec=4;}}
  if(size(#)>1){numberprec=#[2];if(numberprec<8){numberprec=8;}}
  if(solutionprec>numberprec){numberprec=solutionprec;}
  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); // multiplicity of zero solution
  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(string(tt[1..7])=="complex")
        {
          ss=tt;
          ringprec=system("getPrecDigits");
        }
      }
    }
  }

  list roots,simple;
  if(deg(p)==0) // only zero was root
  {
    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 ground field
  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)// splitting
  {
    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 polynomial should be near to zero
    // remember that l contains a list of strings
    // in the case of a different ring
    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 has 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 univariant 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
*/
static proc splitsqrfree(poly p, poly x)
{
  int dd=deg(p);
  if(dd==1){return(p);}
  int i=1;
  int j;
  ideal h,split;
  poly high;

  h=interred(ideal(p,diff(p,x)));
  if(deg(h[1])==0){return(p);}
  high=h[1];
  split[1]=exdiv(p,high,x);
  while(1)
  {
    h=interred(ideal(split[i],high));
    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 polynomial 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 solve( ideal G, list # )
"USAGE:   solve(G [, m, n [, l]] [,\"oldring\"] [,\"nodisplay\"] ); G = ideal,
         m, n, l = integers (control parameters of the method), outR ring,@*
         m: precision of output in digits ( 4 <= m) and of the generated ring
            of complex numbers;
         n: control of multiplicity
            n = 0, return all different roots
            n != 0, find all roots (with multiplicity)
         l: precision of internal computation in decimal digits ( l >=8 )
            only if the basering is not complex or complex with smaller
            precision, @*
         [default: (m,n,l) = (8,0,30), or if only (m,n) are set explicitly
          with n!=0, then (m,n,l) = (m,n,60) ]
ASSUME:  the ideal is 0-dimensional;@*
         basering has characteristic 0 and is either complex or without
         parameters;
RETURN:  (1) If called without the additional parameter @code{\"oldring\"}: @*
         ring @code{R} with the same number of variables but with complex
         coefficients (and precision m). @code{R} comes with a list
         @code{SOL} of numbers, in which complex roots of G are stored: @*
         * If n  = 0, @code{SOL} is the list of all different solutions, each
         of them being represented by a list of numbers. @*
         * If n != 0, @code{SOL} is a list of two list: SOL[i][1] is the list
         of all different solutions with the multiplicity SOL[i][2].@*
         SOL is ordered w.r.t. multiplicity (the smallest first). @*
         (2) If called with the additional parameter @code{\"oldring\"}, the
         procedure looks for an appropriate ring (at top level) in which
         the solutions can be stored (interactive). @*
         The user may then select an appropriate ring and choose a name for
         the output list in this ring. The list is exported directly to the
         selected ring and the return value is a string \"result exported to\"
         + name of the selected ring.
NOTE:    If the problem is not 0-dim. the procedure stops with ERROR. If the
         ideal G is not a lexicographic Groebner basis, the lexicographic
         Groebner basis is computed internally (Hilbert driven).  @*
         The computed solutions are displayed, unless @code{solve} is called
         with the additional parameter @code{\"nodisplay\"}.
EXAMPLE: example solve; shows an example
"
{
// test if basering admissible
  if (char(basering)!=0){ERROR("characteristic of basering not 0");}
  if ((charstr(basering)[1]=="0") and (npars(basering)!=0))
  { ERROR("basering has parameters"); }

// some global settings and control
  int oldr, nodisp, ii, jj;
  list LL;
  int outprec = 8;
  int mu = 0;
  int prec = 30;
  // check additional parameters...
  if (size(#)>0)
  {
    int sofar=1;
    if (typeof(#[1])=="int")
    {
      outprec = #[1];
      if (outprec<4){outprec = 4;}
      if (size(#)>1)
      {
        if (typeof(#[2])=="int")
        {
          mu = #[2];
          if (size(#)>2)
          {
            if (typeof(#[3])=="int")
            {
              prec = #[3];
              if (prec<8){prec = 8;}
            }
            else
            {
              if(mu!=0){prec = 60;}
              if (#[3]=="oldring"){ oldr=1; }
              if (#[3]=="nodisplay"){ nodisp=1; }
            }
            sofar=3;
          }
        }
        else
        {
          if (#[2]=="oldring"){ oldr=1; }
          if (#[2]=="nodisplay"){ nodisp=1; }
        }
        sofar=2;
      }
    }
    else
    {
      if (#[1]=="oldring"){ oldr=1; }
      if (#[1]=="nodisplay"){ nodisp=1; }
    }
    for (ii=sofar+1;ii<=size(#);ii++)
    { // check for additional strings
      if (typeof(#[ii])=="string")
      {
        if (#[ii]=="oldring"){ oldr=1; }
        if (#[ii]=="nodisplay"){ nodisp=1; }
      }
    }
  }
  if (outprec>prec){prec = outprec;}
  // if interaktive version is chosen -- choice of basering (Top::`outR`)
  // and name for list of solutions (outL):
  if (oldr==1)
  {
    list Out;
    LL=names(Top);
    for (ii=1;ii<=size(LL);ii++)
    {
      if (typeof(`LL[ii]`)=="ring")
      {
        if (find(charstr(`LL[ii]`),"complex,"+string(outprec)))
        {
          jj++;
          Out[jj]=LL[ii];
        }
      }
    }
    if (size(Out)>0)
    {
      print("// *** You may select between the following rings for storing "+
            "the list of");
      print("// *** complex solutions:");
      Out;
      print("// *** Enter the number of the chosen ring");
      print("// ***  (0: none of them => new ring created and returned)");
      string chosen;
      while (chosen=="") { chosen=read(""); }
      execute("def tchosen = "+chosen);
      if (typeof(tchosen)=="int")
      {
        if ((tchosen>0) and (tchosen<=size(Out)))
        {
          string outR = Out[tchosen];
          print("// *** You have chosen the ring "+ outR +". In this ring"
                +" the following objects");
          print("//*** are defined:");
          listvar(Top::`outR`);
          print("// *** Enter a name for the list of solutions (different "+
                "from existing names):");
          string outL;
          while (outL==""){ outL=read(""); }
        }
      }
    }
    else
    {
      print("No appropriate ring for storing the list of solutions found " +
             "=> new ring created and returned");
    }
    if (not(defined(outR))) { oldr=0; }
  }

//  string rinC = nameof(basering)+"C";
  string sord = ordstr(basering);
  int nv = nvars(basering);
  def rin = basering;
  intvec ovec = option(get);
  option(redSB);
  option(returnSB);
  int sb = attrib(G,"isSB");
  int lp = 0;
  if (size(sord)==size("C,lp()"+string(nv)))
  {
    lp = find(sord,"lp");
  }

// ERROR
  if (sb){if (dim(G)!=0){ERROR("ideal not zero-dimensional");}}

// the trivial homogeneous case (unique solution: (0,...0))
  if (homog(G))
  {
    if (sb==0)
    {
      execute("ring dphom=("+charstr(rin)+"),("+
      varstr(rin)+"),dp;");
      ideal G = std(imap(rin,G));
      if (dim(G)!=0){ERROR("ideal not zero-dimensional");}
      int vdG=vdim(G);
    }
    if (oldr!=1)
    {
      execute("ring rinC =(complex,"+string(outprec)+
                 "),("+varstr(basering)+"),lp;");
      list SOL;
      if (mu==0){SOL[1] = zerolist(nv);}
      else{SOL[1] = list(zerolist(nv),list(vdG));}
      export SOL;
      if (nodisp==0) { print(SOL); }
      option(set,ovec);
      dbprint( printlevel-voice+3,"
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
      return(rinC);
    }
    else
    {
      setring (Top::`outR`);
      list SOL;
      if (mu==0){SOL[1] = zerolist(nv);}
      else{SOL[1] = list(zerolist(nv),list(vdG));}
      execute("def "+outL + "=SOL;");
      execute("export "+outL+";");
      if (nodisp==0) { print(SOL); }
      option(set,ovec);
      kill SOL;
      return("result exported to "+outR+" as list "+outL);
    }
  }

// look for reduced standard basis in lex
  if (sb*lp==0)
  {
    if (sb==0)
    {
      execute("ring dphilb=("+charstr(rin)+"),("+ varstr(rin)+"),dp;");
      ideal G = imap(rin,G);
      G = std(G);
      if (dim(G)!=0){ERROR("ideal not zero-dimensional");}
    }
    else
    {
      def dphilb = basering;
      G=interred(G);
      attrib(G,"isSB",1);
    }
    execute("ring lexhilb=("+charstr(rin)+"),("+ varstr(rin)+"),lp;");
    option(redTail);
    ideal H = fglm(dphilb,G);
    kill dphilb;
    H = simplify(H,2);
    if (lp){setring rin;}
    else
    {
      execute("ring lplex=("+charstr(rin)+"),("+varstr(rin)+"),lp;");
    }
    ideal H = imap(lexhilb,H);
    kill lexhilb;
  }
  else{ideal H = interred(G);}

// only 1 variable
  def hr = basering;
  if (nv==1)
  {
    if ((mu==0) and (charstr(basering)[1]=="0"))
    { // special case
      list L = laguerre_solve(H[1],prec,prec,mu,0); // list of strings
      if (oldr!=1)
      {
        execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
        list SOL;
        for (ii=1; ii<=size(L); ii++ ) { execute("SOL[ii]=number("+L[ii]+");"); }
        export SOL;
        if (nodisp==0) { print(SOL); }
        option(set,ovec);
        dbprint( printlevel-voice+3,"
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
        return(rinC);
      }
      else
      {
        setring (Top::`outR`);
        list SOL;
        for (ii=1; ii<=size(L); ii++ ) { execute("SOL[ii]="+L[ii]+";"); }
        execute("def "+outL + "=SOL;");
        execute("export "+outL+";");
        if (nodisp==0) { print(SOL); }
        option(set,ovec);
        kill SOL;
        return("result exported to "+outR+" as list "+outL);
      }
    }
    else
    {
      execute("ring internC=(complex,"+string(prec)+"),("+varstr(basering)+"),lp;");
      ideal H = imap(hr,H);
      list sp = splittolist(splitsqrfree(H[1],var(1)));
      jj = size(sp);
      while(jj>0)
      {
        sp[jj][1] = laguerre(sp[jj][1],prec,1);
        jj--;
      }
      setring hr;
      if (oldr!=1)
      {
        execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
        list SOL;
        list sp=imap(internC,sp);
        if(mu!=0){ SOL=sp; }
        else
        {
          jj = size(sp);
          SOL=sp[jj][1];
          while(jj>1)
          {
            jj--;
            SOL = sp[jj][1]+SOL;
          }
        }
        export SOL;
        if (nodisp==0) { print(SOL); }
        option(set,ovec);
        dbprint( printlevel-voice+3,"
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
        return(rinC);
      }
      else
      {
        setring (Top::`outR`);
        list SOL;
        list sp=imap(internC,sp);
        if(mu!=0){ SOL=sp; }
        else
        {
          jj = size(sp);
          SOL=sp[jj][1];
          while(jj>1)
          {
            jj--;
            SOL = sp[jj][1]+SOL;
          }
        }
        kill sp;
        execute("def "+outL + "=SOL;");
        execute("export "+outL+";");
        if (nodisp==0) { print(SOL); }
        option(set,ovec);
        kill SOL;
        return("result exported to "+outR+" as list "+outL);
      }
    }
  }

// the triangular sets (not univariate case)
  attrib(H,"isSB",1);
  if (mu==0)
  {
    list sp = triangMH(H); // faster, but destroy multiplicity
  }
  else
  {
    list sp = triangM(H);
  }

//   create the complex ring and map the result
  if (outprec<prec)
  {
    execute("ring internC=(complex,"+string(prec)+"),("+varstr(hr)+"),lp;");
  }
  else
  {
    execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
  }
  list triC = imap(hr,sp);

// solve the tridiagonal systems
  int js = size(triC);
  list ret1;
  if (mu==0)
  {
    ret1 = trisolve(list(),triC[1],prec);
    while (js>1)
    {
      ret1 = trisolve(list(),triC[js],prec)+ret1;
      js--;
    }
  }
  else
  {
    ret1 = mutrisolve(list(),triC[1],prec);
    while (js>1)
    {
      ret1 = addlist(mutrisolve(list(),triC[js],prec),ret1,1);
      js--;
    }
    ret1 = finalclear(ret1);
  }

// final computations
  option(set,ovec);
  if (outprec==prec)
  { // we are in ring rinC
    if (oldr!=1)
    {
      list SOL=ret1;
      export SOL;
      if (nodisp==0) { print(SOL); }
      dbprint( printlevel-voice+3,"
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
      return(rinC);
    }
    else
    {
      setring (Top::`outR`);
      list SOL=imap(rinC,ret1);
      execute("def "+outL + "=SOL;");
      execute("export "+outL+";");
      if (nodisp==0) { print(SOL); }
      kill SOL;
      return("result exported to "+outR+" as list "+outL);
    }
  }
  else
  {
    if (oldr!=1)
    {
      execute("ring rinC =(complex,"+string(outprec)+"),("+varstr(basering)+"),lp;");
      list SOL=imap(internC,ret1);
      export SOL;
      if (nodisp==0) { print(SOL); }
      dbprint( printlevel-voice+3,"
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
      return(rinC);
    }
    else
    {
      setring (Top::`outR`);
      list SOL=imap(internC,ret1);
      execute("def "+outL + "=SOL;");
      execute("export "+outL+";");
      if (nodisp==0) { print(SOL); }
      kill SOL;
      return("result exported to "+outR+" as list "+outL);
    }
  }
}
example
{
    "EXAMPLE:";echo=2;
    // Find all roots of a multivariate ideal using triangular sets:
    int d,t,s = 4,3,2 ;
    int i;
    ring A=0,x(1..d),dp;
    poly p=-1;
    for (i=d; i>0; i--) { p=p+x(i)^s; }
    ideal I = x(d)^t-x(d)^s+p;
    for (i=d-1; i>0; i--) { I=x(i)^t-x(i)^s+p,I; }
    I;
    // the multiplicity is
    vdim(std(I));
    def AC=solve(I,6,0,"nodisplay");  // solutions should not be displayed
    // list of solutions is stored in AC as the list SOL (default name)
    setring AC;
    size(SOL);               // number of different solutions
    SOL[5];                  // the 5th solution
    // you must start with char. 0
    setring A;
    def AC1=solve(I,6,1,"nodisplay");
    setring AC1;
    size(SOL);               // number of different multiplicities
    SOL[1][1][1];            // a solution with
    SOL[1][2];               // multiplicity 1
    SOL[2][1][1];            // a solution with
    SOL[2][2];               // multiplicity 12
    // the number of different solutions is equal to
    size(SOL[1][1])+size(SOL[2][1]);
    // the number of complex solutions (counted with multiplicities) is
    size(SOL[1][1])*SOL[1][2]+size(SOL[2][1])*SOL[2][2];
}
//////////////////////////////////////////////////////////////////////////////
//                subprocedures for solve


/*
* return one zero-solution
*/
static proc zerolist(int nv)
{
  list ret;
  int i;
  number o=0;

  for (i=nv;i>0;i--){ret[i] = o;}
  return(ret);
}

/* ----------------------- check solution ----------------------- */
static proc multsol(list ff, int c)
{
  int i,j;

  i = 0;
  j = size(ff);
  while (j>0)
  {
    if(c){i = i+ff[j][2]*size(ff[j][1]);}
    else{i = i+size(ff[j][1]);}
    j--;
  }
  return(i);
}

/*
* the inputideal A => zero ?
*/
static proc checksol(ideal A, list lr)
{
  int d = nvars(basering);
  list ro;
  ideal re,h;
  int i,j,k;

  for (i=size(lr);i>0;i--)
  {
    ro = lr[i];
    for (j=d;j>0;j--)
    {
      re[j] = var(j)-ro[j];
    }
    attrib(re,"isSB",1);
    k = size(reduce(A,re,1));
    if (k){return(i);}
  }
  return(0);
}

/*
*  compare 2 solutions: returns 0 for equal
*/
static proc cmpn(list a,list b)
{
  int ii;

  for(ii=size(a);ii>0;ii--){if(a[ii]!=b[ii]) break;}
  return(ii);
}

/*
*  delete equal solutions in the list
*/
static proc delequal(list r, int w)
{
  list h;
  int i,j,k,c;

  if (w)
  {
    k = size(r);
    h = r[k][1];
    k--;
    while (k>0)
    {
      h = r[k][1]+h;
      k--;
    }
  }
  else{h = r;}
  k=size(h);
  i=1;
  while(i<k)
  {
    j=k;
    while(j>i)
    {
      c=cmpn(h[i],h[j]);
      if(c==0)
      {
        h=delete(h,j);
        k--;
      }
      j--;
    }
    i++;
  }
  return(h);
}

/* ----------------------- substitution ----------------------- */
/*
* instead of subst(T,var(v),n), much faster
*   need option(redSB) !
*/
static proc linreduce(ideal T, int v, number n)
{
  ideal re = var(v)-n;
  attrib (re,"isSB",1);
  return (reduce(T,re));
}

/* ----------------------- triangular solution ----------------------- */
/*
* solution of one tridiagonal system T
*   with precision prec
*   T[1] is univariant in var(1)
*   list o is empty for the first call
*/
static proc trisolve(list o, ideal T, int prec)
{
  list lroots,ll;
  ideal S;
  int i,d;

  d = size(T);
  S = interred(ideal(T[1],diff(T[1],var(d))));
  if (deg(S[1]))
  {
    T[1] = exdiv(T[1],S[1],var(d));
  }
  ll = laguerre(T[1],prec,1);
  for (i=size(ll);i>0;i--){ll[i] = list(ll[i])+o;}
  if (d==1){return(ll);}
  for (i=size(ll);i>0;i--)
  {
    S = linreduce(ideal(T[2..d]),d,ll[i][1]);
    lroots = trisolve(ll[i],S,prec)+lroots;
  }
  return(lroots);
}

/* ------------------- triangular solution (mult) ------------------- */
/*
*  recompute equal solutions w.r.t. multiplicity
*/
static proc finalclear(list b)
{
  list a = b;
  list r;
  int i,l,ju,j,k,ku,mu,c;

// a[i] only
  i = 1;
  while (i<=size(a))
  {
    ju = size(a[i][1]);
    j = 1;
    while (j<=ju)
    {
      mu = 1;
      k = j+1;
      while (k<=ju)
      {
        c = cmpn(a[i][1][j],a[i][1][k]);
        if (c==0)
        {
          a[i][1] = delete(a[i][1],k);
          ju--;
          mu++;
        }
        else{k++;}
      }
      if (mu>1)
      {
        r[1] = a[i];
        r[1][1] = list(a[i][1][j]);
        a[i][1] = delete(a[i][1],j);
        a = addlist(r,a,mu);
        ju--;
      }
      else{j++;}
    }
    if (ju==0){a = delete(a,i);}
    else{i++;}
  }
// a[i], a[l]
  i = 1;
  while (i<size(a))
  {
    ju = size(a[i][1]);
    l = i+1;
    while (l<=size(a))
    {
      ku = size(a[l][1]);
      j = 1;
      while (j<=ju)
      {
        mu = 0;
        k = 1;
        while (k<=ku)
        {
          c = cmpn(a[i][1][j],a[l][1][k]);
          if (c==0)
          {
            mu = a[i][2]+a[l][2];
            r[1] = a[l];
            r[1][1] = list(a[l][1][k]);
            r[1][2] = mu;
            a[l][1] = delete(a[l][1],k);
            a = addlist(r,a,1);
            ku--;
            break;
          }
          else{k++;}
        }
        if (mu)
        {
          a[i][1] = delete(a[i][1],j);
          ju--;
        }
        else{j++;}
      }
      if (ku){l++;}
      else
      {
        a = delete(a,l);
      }
    }
    if (ju){i++;}
    else
    {
      a = delete(a,i);
    }
  }
  return(a);
}

/*
* convert to list
*/
static proc splittolist(ideal sp)
{
  int j = size(sp);
  list spl = list(list(sp[j],j));

  j--;
  while (j>0)
  {
    if (deg(sp[j]))
    {
      spl = list(list(sp[j],j))+spl;
    }
    j--;
  }
  return(spl);
}

/*
*  multiply the multiplicity
*/
static proc multlist(list a, int m)
{
  int i;
  for (i=size(a);i>0;i--){a[i][2] = a[i][2]*m;}
  return(a);
}

/*
*  a+b w.r.t. to multiplicity as ordering
*    (programming like spolys)
*/
static proc addlist(list a, list b, int m)
{
  int i,j,k,l,s;
  list r = list();

  if (m>1){a = multlist(a,m);}
  k = size(a);
  l = size(b);
  i = 1;
  j = 1;
  while ((i<=k)&&(j<=l))
  {
    s = a[i][2]-b[j][2];
    if (s>=0)
    {
      r = r+list(b[j]);
      j++;
      if (s==0)
      {
        s = size(r);
        r[s][1] = r[s][1]+a[i][1];
        i++;
      }
    }
    else
    {
      r = r+list(a[i]);
      i++;
    }
  }
  if (i>k)
  {
    if (j<=l){r = r+list(b[j..l]);}
  }
  else{r = r+list(a[i..k]);}
  return(r);
}

/*
* solution of one tridiagonal system T with multiplicity
*   with precision prec
*   T[1] is univariant in var(1)
*   list o is empty for the first call
*/
static proc mutrisolve(list o, ideal T, int prec)
{
  list lroots,ll,sp;
  ideal S,h;
  int i,d,m,z;

  d = size(T);
  sp = splittolist(splitsqrfree(T[1],var(d)));
  if (d==1){return(l_mutrisolve(sp,o,prec));}
  z = size(sp);
  while (z>0)
  {
    m = sp[z][2];
    ll = laguerre(sp[z][1],prec,1);
    i = size(ll);
    while(i>0)
    {
      h = linreduce(ideal(T[2..d]),d,ll[i]);
      if (size(lroots))
      {
        lroots = addlist(mutrisolve(list(ll[i])+o,h,prec),lroots,m);
      }
      else
      {
        lroots = mutrisolve(list(ll[i])+o,h,prec);
        if (m>1){lroots=multlist(lroots,m);}
      }
      i--;
    }
    z--;
  }
  return(lroots);
}

/*
*  the last call, we are ready
*/
static proc l_mutrisolve(list sp, list o, int prec)
{
  list lroots,ll;
  int z,m,i;

  z = size(sp);
  while (z>0)
  {
    m = sp[z][2];
    ll = laguerre(sp[z][1],prec,1);
    for (i=size(ll);i>0;i--){ll[i] = list(ll[i])+o;}
    if (size(lroots))
    {
      lroots = addlist(list(list(ll,m)),lroots,1);
    }
    else
    {
      lroots = list(list(ll,m));
    }
    z--;
  }
  return(lroots);
}
///////////////////////////////////////////////////////////////////////////////

proc ures_solve( ideal gls, list # )
"USAGE:   ures_solve(i [, k, p] ); i = ideal, k, p = integers
   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=0, p=30)
ASSUME:  i is a zerodimensional ideal given by a quadratic system, that is,@*
         nvars(basering) = ncols(i) = number of vars actually occurring in i,
RETURN:  If the ground field is the field of complex numbers: list of numbers
         (the complex roots of the polynomial system i=0). @*
         Otherwise: ring @code{R} with the same number of variables but with
         complex coefficients (and precision p). @code{R} comes with a list
         @code{SOL} of numbers, in which complex roots of the polynomial
         system i are stored: @*
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;
    }
  }

  list LL=uressolve(gls,typ,prec,1);
  int sizeLL=size(LL);
  if (sizeLL==0)
  {
    dbprint(printlevel-voice+3,"No solution found!");
    return(list());
  }
  if (typeof(LL[1][1])=="string")
  {
    int ii,jj;
    int nv=size(LL[1]);
    execute("ring rinC =(complex,"+string(prec)+",I),("
                           +varstr(basering)+"),lp;");
    list SOL,SOLnew;
    for (ii=1; ii<=sizeLL; ii++)
    {
      SOLnew=list();
      for (jj=1; jj<=nv; jj++)
      {
        execute("SOLnew["+string(jj)+"]="+LL[ii][jj]+";");
      }
      SOL[ii]=SOLnew;
    }
    kill SOLnew;
    export SOL;
    dbprint( printlevel-voice+3,"
// 'ures_solve' created a ring, in which a list SOL of numbers (the complex
// solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; SOL; ");
    return(rinC);
  }
  else
  {
    return(LL);
  }
}
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;
    def R=ures_solve(gls,0,16);
    setring R; SOL;
}
///////////////////////////////////////////////////////////////////////////////

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:  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.
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:  Ground field K is the field of rational numbers, p and v are lists
         of elements of the ground field K with p[j] != -1,0,1, size(p) = n
         (= number of vars) and size(v)=N=(d+1)^n.
RETURN:  poly f, the unique polynomial f of degree n*d with prescribed values
         v[i] at the points p(i)=(p[1]^(i-1),..,p[n]^(i-1)), i=1,..,N.
NOTE:    mainly useful when n=1, i.e. f is satisfying f(p^(i-1)) = v[i],
         i=1..d+1.
SEE ALSO: vandermonde.
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 );
}

///////////////////////////////////////////////////////////////////////////////
// changed for Singular 3
// Return value is now a list: (rlist, rn@)
static proc psubst( int d, int dd, int n, list resl,
                    ideal fi, int elem, int nv, int prec, int rn@, list rlist)
{
    //   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

    list LL;
    int pdebug;
    int olddd=dd;

    dbprint(printlevel-voice+2, "// 0 step "+string(dd)+" of "+string(elem) );

    if ( dd <= elem )
    {
        int loop = 1;
        int k;
        list lsr,lh;
        poly ps;
        int thedd;

        dbprint( printlevel-voice+1,"// 1 dd = "+string(dd) );

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

            if ( n-1 > 0 )
            {
                dbprint( printlevel-voice,
                    "// 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 Polynom immer als letztes in der Liste?!?
                    // leadexp(ps)
                }
                else
                {
                    loop=0;
                }
            }
            else
            {
                dbprint( printlevel-voice,
                    "// 2 ps=fi["+string(dd+1)+"]"+"  leadexp(ps)="
                        +string(leadexp(ps)) );
                dd++;
            }
        }
        thedd=dd;
        ps= fi[thedd];

        dbprint( printlevel-voice+1,
            "// 3    fi["+string(thedd-1)+"]"+"  leadexp(fi[thedd-1])="
                +string(leadexp(fi[thedd-1])) );
        dbprint( printlevel-voice+1,
            "// 3 ps=fi["+string(thedd)+"]"+"  leadexp(ps)="
                +string(leadexp(ps)) );

        for ( k= nv; k > nv-d; k-- )
        {
            dbprint( printlevel-voice,
                "// 4 subst(fi["+string(thedd)+"],"
                    +string(var(k))+","+string(resl[k])+");" );
            ps = subst(ps,var(k),resl[k]);
        }

        dbprint( printlevel-voice, "// 5 substituted ps="+string(ps) );

        if ( ps != 0 )
        {
            lsr= laguerre_solve( ps, prec, prec, 0 );
        }
        else
        {
            dbprint( printlevel-voice+1,"// 30 ps == 0, thats not cool...");
            lsr=list(number(0));
        }

        dbprint( printlevel-voice+1,
         "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]" );

        if ( size(lsr) > 1 )
        {
            dbprint( printlevel-voice+1,
                "// 10 checking roots found before, range "
                    +string(dd-olddd)+" -- "+string(dd) );
            dbprint( printlevel-voice+1,
                "// 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=list(number(0));
            }

            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@++;
                    }
                    LL = psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec,
                                 rn@, rlist );
                    rlist = LL[1];
                    rn@ = LL[2];
                }
                else
                {
                   if ( i > 1 ) {  rn@++; }  //bug found by O.Labs
                   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 )
            {
                LL= psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec,
                            rn@, rlist );
                rlist = LL[1];
                rn@ = LL[2];
            }
            else
            {
                if ( pdebug>=1 )
                { "// 30_2 <"+string(rn@)+"> "+string(size(resl))+" <-----";}
                if ( pdebug>=2 )
                { resl; }
                rlist[rn@]=resl;
            }
        }
    }
    return(list(rlist,rn@));
}

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

proc fglm_solve( ideal fi, list # )
"USAGE:   fglm_solve(i [, p] ); i ideal, p integer
ASSUME:  the ground field has char 0.
RETURN:  ring @code{R} with the same number of variables but with complex
         coefficients (and precision p). @code{R} comes with a list
         @code{rlist} of numbers, in which the complex roots of i are stored.@*
         p>0: gives precision of complex numbers in decimal digits [default:
         p=30].
NOTE:    The procedure uses a standard basis of i to determine all complex
         roots of i.
EXAMPLE: example fglm_solve; shows an example
"
{
    int prec=30;

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

    def R = lex_solve(stdfglm(fi),prec);
    dbprint( printlevel-voice+3,"
// 'fglm_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");
    return(R);
}
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;
    def R = fglm_solve(s,10);
    setring R; rlist;
}

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

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

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

    if ( !defined(pdebug) ) { int pdebug; }
    def oring= basering;

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

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

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

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

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

    dbprint(printlevel-voice+2,"// laguerre found roots: "+string(size(li)));
    int rn@;

    for ( j= 1; j <= lis; j++ )
    {
        dbprint(printlevel-voice+1,"// root "+string(j) );
        rn@++;
        resl[nv]= li[j];
        LL = psubst( 1, 2, nv-1, resl, fi, idelem, nv, prec, rn@, rlist );
        rlist=LL[1];
        rn@=LL[2];
    }

    dbprint( printlevel-voice+3,"
// 'lex_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");

    return(RC);
}
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;
    def R = lex_solve(stdfglm(s),10);
    setring R; 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 is a zero-dimensional ideal
RETURN:  ring @code{R} with the same number of variables but with complex
         coefficients (and precision p). @code{R} comes with a list
         @code{rlist} of numbers, in which the complex roots of i are stored.
NOTE:    The procedure uses a triangular system (Lazard's Algorithm with
         factorization) computed from a standard basis to determine
         recursively all complex roots of the input ideal i with Laguerre's
         algorithm.
EXAMPLE: example triangLf_solve; shows an example
"
{
    int prec=30;

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

    def R=triang_solve(triangLfak(stdfglm(fi)),prec);
    dbprint( printlevel-voice+3,"
// 'triangLf_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");
    return(R);
}
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;
    def R = triangLf_solve(s,10);
    setring R; 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:  ring @code{R} with the same number of variables but with complex
         coefficients (and precision p). @code{R} comes with a list
         @code{rlist} of numbers, in which the complex roots of i are stored.
NOTE:    The procedure uses a triangular system (Moellers Algorithm) computed
         from a standard basis  of input ideal i to determine recursively all
         complex roots with Laguerre's algorithm.
EXAMPLE: example triangM_solve; shows an example
"
{
    int prec=30;

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

    def R = triang_solve(triangM(stdfglm(fi)),prec);
    dbprint( printlevel-voice+3,"
// 'triangM_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");
    return(R);
}
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;
    def R = triangM_solve(s,10);
    setring R; 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 is a zero-dimensional ideal.
RETURN:  ring @code{R} with the same number of variables, but with complex
         coefficients (and precision p). @code{R} comes with a list
         @code{rlist} of numbers, in which the complex roots of i are stored.
NOTE:    The procedure uses a triangular system (Lazard's Algorithm) computed
         from a standard basis of input ideal i to determine recursively all
         complex roots with Laguerre's algorithm.
EXAMPLE: example triangL_solve; shows an example
"
{
    int prec=30;

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

    def R=triang_solve(triangL(stdfglm(fi)),prec);
    dbprint( printlevel-voice+3,"
// 'triangL_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");
    return(R);
}
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;
    def R = triangL_solve(s,10);
    setring R; rlist;
}


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

proc triang_solve( list lfi, int prec, list # )
"USAGE:   triang_solve(l,p [,d] ); l=list, p,d=integers@*
         l is 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 the algorithm of Lazard or the algorithm of
         Moeller (see triang.lib).
RETURN:  ring @code{R} with the same number of variables, but with complex
         coefficients (and precision p). @code{R} comes with a list
         @code{rlist} of numbers, in which the complex roots of l are stored.@*
EXAMPLE: example triang_solve; shows an example
"
{
    def oring= basering;
    list LL;

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

    // list with entry 0 (number)
    number nn=0;

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

    dbprint( printlevel-voice+2,"// myCompDigits="+string(myCompDigits) );

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

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

    int 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,0);
        lis= size(li);

        dbprint( printlevel-voice+2,"// laguerre found roots: "+string(lis) );

        for ( j= 1; j <= lis; j++ )
        {
            dbprint( printlevel-voice+2,"// root "+string(j) );
            rn@++;
            resu[nv]= li[j];
            LL = psubst( 1, 2, nv-1, resu, fi, idelem, nv, myCompDigits,
                         rn@, rlist );
            rlist = LL[1];
            rn@ = LL[2];
        }
    }

    dbprint( printlevel-voice+3,"
// 'triang_solve' created a ring, in which a list rlist of numbers (the
// complex solutions) is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):
        setring R; rlist; ");

    return(RC);
}
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;
    def R=triang_solve(triangLfak(stdfglm(s)),10);
    setring R; rlist;
}

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

proc simplexOut(list l)
"USAGE:   simplexOut(l); l list
ASSUME:  l is the output of simplex.
RETURN:  Nothing. The procedure prints the computed solution of simplex
         (as strings) in a nice format.
SEE ALSO: simplex
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 )
      {
        "x"+string(iposv[i-1])+" = "+string(m[i,1]);
      }
//        else
//        {
//               "Y"; iposv[i-1]-N+1;
//        }
    }
  }
}
example
{
  "EXAMPLE:";echo=2;
  ring r = (real,10),(x),lp;

  // consider the max. problem:
  //
  //    maximize  x(1) + x(2) + 3*x(3) - 0.5*x(4)
  //
  //  with constraints:   x(1) +          2*x(3)          <= 740
  //                             2*x(2)          - 7*x(4) <=   0
  //                               x(2) -   x(3) + 2*x(4) >=   0.5
  //                      x(1) +   x(2) +   x(3) +   x(4)  =   9
  //
  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;

  int n = 4;  // number of constraints
  int m = 4;  // number of variables
  int m1= 2;  // number of <= constraints
  int m2= 1;  // number of >= constraints
  int m3= 1;  // number of == constraints

  list sol=simplex(sm, n, m, m1, m2, m3);
  simplexOut(sol);
}


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