///////////////////////////////////////////////////////////////////////////////
version="$Id: solve.lib,v 1.21.2.13 2002/10/21 14:38:44 lossen Exp $";
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 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
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)
@format
 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 )
@end format
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 
         result is of type
@format
 string: if the basering is not complex,
 number: otherwise.
@end format
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 poly 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 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 solve( ideal G, list # )
"USAGE:   solve(G [, m, n, l] ); G = ideal,
         m, n, l = integers (control parameters of the method)
 @format
         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
                    for (n != 0 and size(#) = 2) l = 60 )
 @end format
ASSUME:  the ideal is 0-dimensional;@*
         basering has characteristic 0 and is either complex or
         without parameters;
RETURN:  list of solutions of the ideal G, depending on n; one solution is a 
         list of complex numbers in the generated output ring (the new 
         basering).
@format
 The result is a list L
    n  = 0: a list of all different solutions (L[i]),
    n != 0: a list of two elements,
            L[i][1] contains all different solutions with the same multiplicity
            L[i][2] the multiplicity
 L is ordered w.r.t. multiplicity (the smallest first).
@end format
NOTE:    If the problem is not 0-dim. the procedure stops with ERROR, if the 
         ideal G is not a lex. standard basis, it is generated with internal 
         computation (Hilbert driven), if the input-ring (with char 0) has 
         the name "<A>", the lexicographical and complex output-ring has the 
         name "<A>C".
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 outprec = 8;
  int mu = 0;
  int prec = 30;
  if (size(#)>0){outprec = #[1];if (outprec<4){outprec = 4;}}
  if (size(#)>1){mu = #[2];}
  if (size(#)>2){prec = #[3];if (prec<8){prec = 8;}}
  else {if(mu!=0){prec = 60;}}
  if (outprec>prec){prec = outprec;}
  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 homog case
  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");}
    }
    changering(rinC,outprec);
    list ret0;
    if (mu==0){ret0[1] = zerolist(nv);}
    else{ret0[1] = list(zerolist(nv),list(vdim(G)));}
    option(set,ovec);
    keepring basering;
    return(ret0);
  }

// 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;      
    }    
    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
      changering(rinC,outprec);
      list sp;
      for (int ii=1; ii<=size(L); ii++ ) 
      { 
        execute("sp[ii]="+L[ii]); 
      }
      keepring basering;
      return(sp);
    }
    else
    {
      execute("ring internC=(complex,"+string(prec)+
                 "),("+varstr(basering)+"),lp;");

      ideal H = imap(hr,H);
      list sp = splittolist(splitsqrfree(H[1],var(1)));
      int jj = size(sp);
      while(jj>0)
      {
        sp[jj][1] = laguerre(sp[jj][1],prec,1);
        jj--;
      }
      setring hr;
      changering(rinC,outprec);
      list sp=imap(internC,sp);

      keepring basering;
      if(mu!=0){return(sp);}
      jj = size(sp);
      list ll=sp[jj][1];
      while(jj>1)
      {
        jj--;
        ll = sp[jj][1]+ll;
      }
      return(ll);
    }    
  }    
  

// 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
  {
    changering(rinC,prec);
  }
  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)
  {
    keepring basering;
    return(ret1);
  }
  changering(rinC,outprec);
  keepring basering;
  return(imap(internC,ret1));
}
example
{
    "EXAMPLE:";echo=2;
    // Find all roots of a multivariate ideal using triangular sets:
    int d=4;// with these 3 parameters you may construct
    int t=3;// very hard problems for 'solve'
    int s=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 mutiplicity is
    vdim(std(I));
    list l1=solve(I,6,0);
    // the current ring is
    AC;
    // you must start with char. 0
    setring A;
    list l2=solve(I,6,1);
    // the number of different solutions is
    size(l1);
    // this is equal to
    size(l2[1][1])+size(l2[2][1]);
    // the number of solutions with multiplicity is
    size(l2[1][1])*l2[1][2]+size(l2[2][1])*l2[2][2];
    // the solutions with multiplicity
    l2[2][2];
    // are
    l2[2][1];
}
//////////////////////////////////////////////////////////////////////////////
//                subprocedures for solve

/* ----------------------- support ----------------------- */
/*
* the complex ring with precision outprec
* has the well defined name: rinC
* 1. if such a ring exists with the precision outprec,
*    this will be the current ring
* 2. otherwise such a ring will be created
*/
static proc changering(string rinC, int outprec)
{
  string rinDC = "ring "+rinC+"=(complex,"+string(outprec)+
                 "),("+varstr(basering)+"),lp;";
  string h = "int ex=defined("+rinC+");";

  execute(h);
  if (ex)
  {
    h = "setring "+rinC+";";
    execute(h);
    if (system("getPrecDigits")==outprec)
    {"// name of current ring: "+rinC;}
    else
    {
      execute("kill "+rinC+";");
      execute(rinDC);
      execute("export "+rinC+";");
      "// name of new current ring: "+rinC;
    }
  }
  else
  {
    execute(rinDC);
    execute("export "+rinC+";");
    "// name of new current ring: "+rinC;
  }
  keepring basering;
}

/*
* 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));
    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
@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=0, p=30)
@end format
ASSUME:  i is a zerodimensional ideal with
         nvars(basering) = ncols(i) = number of vars
         actually occurring in i,
RETURN:  list of all (complex) roots of the polynomial system i = 0; the 
         result is of type
@format
   string: if the basering is not complex,
   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 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:  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 );
}
///////////////////////////////////////////////////////////////////////////////

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, 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
ASSUME:  the ground field has char 0.
RETURN:  a list of numbers, the complex roots of i; 
         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.
         It creates a ring rC with the same number of variables but with 
         complex coefficients (and precision p).
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:  i is a reduced lexicographical Groebner bases of a zero-dimensional 
         ideal, sorted by increasing leading terms.
RETURN:  nothing
CREATE:  The procedure creates a complec ring with the same variables but
         with complex coefficients (and precision p).
         In this ring a list rlist of numbers is created, in which the complex
         roots of i are stored.
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,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 is a 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 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 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 is a 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 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 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 resu,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,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@++;
            resu[nv]= li[j];
            psubst( 1, 2, nv-1, resu, 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 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: ***