Changeset 6ed15c in git for Singular/LIB/solve.lib

Timestamp:

May 2, 2005, 9:43:22 AM (19 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

6eff86e3417d554a62456e5d4f0634e215db4891

Parents:

d1a54b1d03b681db8f1032e2f169a09db7a0faf7

Message:

*lossen: new lib


git-svn-id: file:///usr/local/Singular/svn/trunk@7982 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/solve.lib (modified) (47 diffs)

Singular/LIB/solve.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: solve.lib,v 1.28 2005-04-25 16:57:36 Singular Exp $";
+version="$Id: solve.lib,v 1.29 2005-05-02 07:43:22 Singular Exp $";
 category="Symbolic-numerical solving";
 info="
@@ -46 +46 @@
          parameters. 
 RETURN:  list of (complex) roots of the polynomial f, depending on n. The 
-         result is of type
- string: if the basering is not complex,
- number: otherwise.
+         result is 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
@@ -516 +516 @@
 
 proc solve( ideal G, list # )
-"USAGE:   solve(G [, m, n, l] ); G = ideal,
-         m, n, l = integers (control parameters of the method)
+"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;
@@ -525 +525 @@
          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 )
+            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:  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).
- 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).
-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\".
+         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 (on the 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
 "
@@ -550 +557 @@
 // 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");}
+  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;
-  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;}}
+  // 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 (#[ii]=="oldring"){ oldr=1; }
+       if (#[ii]=="nodisplay"){ nodisp=1; }
+    }
+  }   
   if (outprec>prec){prec = outprec;}
-  string rinC = nameof(basering)+"C";
+
+  // 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);
@@ -578 +669 @@
   if (sb){if (dim(G)!=0){ERROR("ideal not zero-dimensional");}}
 
-// the trivial homog case
+// the trivial homogeneous case (unique solution: (0,...0))
   if (homog(G))
   {
@@ -587 +678 @@
       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);
+      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;");
+  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
-    {
+    else {
       def dphilb = basering;      
     }    
-    execute("ring lexhilb=("+charstr(rin)+"),("+
-    varstr(rin)+"),lp;");
+    execute("ring lexhilb=("+charstr(rin)+"),("+ varstr(rin)+"),lp;");
     option(redTail);
     ideal H = fglm(dphilb,G);
@@ -619 +728 @@
     H = simplify(H,2);
     if (lp){setring rin;}
-    else
-    {
+    else {
       execute("ring lplex=("+charstr(rin)+"),("+
       varstr(rin)+"),lp;");
@@ -631 +739 @@
 // only 1 variable
   def hr = basering;
-  if (nv==1)
-  {
-    if ((mu==0) and (charstr(basering)[1]=="0")) // special case 
-    {
+  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
-    {
+      if (oldr!=1) { 
+        execute("ring rinC =(complex,"+string(outprec)+
+                   "),("+varstr(basering)+"),lp;"); 
+        list SOL;
+        for (ii=1; ii<=size(L); ii++ ) { execute("SOL[ii]="+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)));
-      int jj = size(sp);
+      jj = size(sp);
       while(jj>0)
       {
@@ -659 +782 @@
       }
       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);
+      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);
+      }
     }    
   }    
@@ -695 +851 @@
   else
   {
-    changering(rinC,prec);
+    execute("ring rinC =(complex,"+string(outprec)+
+                 "),("+varstr(basering)+"),lp;"); 
   }
   list triC = imap(hr,sp);
@@ -724 +881 @@
 // final computations
   option(set,ovec);
-  if (outprec==prec)
-  {
-    keepring basering;
-    return(ret1);
-  }
-  changering(rinC,outprec);
-  keepring basering;
-  return(imap(internC,ret1));
+  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
@@ -737 +934 @@
     "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 d,t,s = 4,3,2 ; 
     int i;
-    ring A=0,(x(1..d)),dp;
+    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;}
+    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
+    // the multiplicity is
     vdim(std(I));
-    list l1=solve(I,6,0);
-    // the current ring is
-    AC;
+    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;
-    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];
+    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
 
-/* ----------------------- 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;
-}
 
 /*
@@ -1335 +1495 @@
     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 )
+                    ideal fi, int elem, int nv, int prec, int rn@, list rlist)
 {
     //   nv: number of ring variables         (fixed value)
@@ -1349 +1511 @@
     //    d: actual variable
 
+    list LL;
+    int pdebug;
     int olddd=dd;
 
-    if ( pdebug>=1 )
-    {"// 0 step "+string(dd)+" of "+string(elem);}
+    dbprint(printlevel-voice+2, "// 0 step "+string(dd)+" of "+string(elem) );
 
     if ( dd <= elem )
@@ -1362 +1525 @@
         int thedd;
 
-        if ( pdebug>=1 )
-        {"// 1 dd = "+string(dd);}
+        dbprint( printlevel-voice+1,"// 1 dd = "+string(dd) );
 
         thedd=0;
@@ -1372 +1534 @@
             if ( n-1 > 0 )
             {
-                if ( pdebug>=2 )
-                {
+                dbprint( printlevel-voice,
                     "// 2 ps=fi["+string(dd+1)+"]"+" size="
                         +string(size(coeffs(ps,var(n-1))))
-                        +"  leadexp(ps)="+string(leadexp(ps));
-                }
+                        +"  leadexp(ps)="+string(leadexp(ps)) );
+
                 if ( size(coeffs(ps,var(n-1))) == 1 )
                 {
@@ -1392 +1553 @@
             else
             {
-                if ( pdebug>=2 )
-                {
+                dbprint( printlevel-voice,
                     "// 2 ps=fi["+string(dd+1)+"]"+"  leadexp(ps)="
-                        +string(leadexp(ps));
-                }
+                        +string(leadexp(ps)) );
                 dd++;
             }
@@ -1403 +1562 @@
         ps= fi[thedd];
 
-        if ( pdebug>=1 )
-        {
+        dbprint( printlevel-voice+1,
             "// 3    fi["+string(thedd-1)+"]"+"  leadexp(fi[thedd-1])="
-                +string(leadexp(fi[thedd-1]));
+                +string(leadexp(fi[thedd-1])) );
+        dbprint( printlevel-voice+1,
             "// 3 ps=fi["+string(thedd)+"]"+"  leadexp(ps)="
-                +string(leadexp(ps));
-        }
+                +string(leadexp(ps)) );
 
         for ( k= nv; k > nv-d; k-- )
         {
-            if ( pdebug>=2 )
-            {
+            dbprint( printlevel-voice,
                 "// 4 subst(fi["+string(thedd)+"],"
-                    +string(var(k))+","+string(resl[k])+");";
-            }
+                    +string(var(k))+","+string(resl[k])+");" );
             ps = subst(ps,var(k),resl[k]);
         }
 
-        if ( pdebug>=2 )
-        { "// 5 substituted ps="+string(ps); }
+        dbprint( printlevel-voice, "// 5 substituted ps="+string(ps) ); 
 
         if ( ps != 0 )
@@ -1430 +1585 @@
         else
         {
-            if ( pdebug>=1 )
-            { "// 30 ps == 0, thats not cool..."; }
-            lsr=@ln; // lsr=number(0);
+            dbprint( printlevel-voice+1,"// 30 ps == 0, thats not cool...");
+            lsr=list(number(0)); 
         }
 
-        if ( pdebug>=1 )
-        { "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]"; }
+        dbprint( printlevel-voice+1, 
+         "// 6 laguerre_solve found roots: lsr["+string(size(lsr))+"]" ); 
 
         if ( size(lsr) > 1 )
         {
-            if ( pdebug>=1 )
-            {
+            dbprint( printlevel-voice+1,  
                 "// 10 checking roots found before, range "
-                    +string(dd-olddd)+" -- "+string(dd);
-                "// 10 thedd = "+string(thedd);
-            }
+                    +string(dd-olddd)+" -- "+string(dd) );
+            dbprint( printlevel-voice+1,  
+                "// 10 thedd = "+string(thedd) );
 
             int i,j,l;
@@ -1525 +1678 @@
                 "Numerical problem: No root found...";
                 "Output may be incorrect!";
-                nares=@ln;
+                nares=list(number(0));
             }
 
@@ -1541 +1694 @@
                         rn@++;  
                     }
-                    psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec );
+                    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;
+                   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;
                 }
             }
@@ -1562 +1717 @@
             if ( dd < elem )
             {
-                psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec );
+                LL= psubst( d+1, dd+1, n-1, resl, fi, elem, nv, prec, 
+                            rn@, rlist );
+                rlist = LL[1];
+                rn@ = LL[2];
             }
             else
@@ -1574 +1732 @@
         }
     }
+    return(list(rlist,rn@));
 }
 
@@ -1581 +1740 @@
 "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). 
+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.
-         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
 "
@@ -1598 +1757 @@
     }
 
-    lex_solve(stdfglm(fi),prec);
-    keepring basering;
+    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
@@ -1606 +1771 @@
     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;
+    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
+    def R = fglm_solve(s,10);
+    setring R; rlist;
 }
 
@@ -1618 +1783 @@
 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.
+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")
@@ -1633 +1797 @@
     }
 
-    if ( !defined(pdebug) )
-    {
-        int pdebug;
-        pdebug=0;
-        export pdebug;
-    }
-
-    string orings= nameof(basering);
+    if ( !defined(pdebug) ) { int pdebug; }
     def oring= basering;
 
     // change the ground field to complex numbers
-    string nrings= "ring "+orings+"C=(complex,"+string(prec)
+    string nrings= "ring RC =(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);
@@ -1674 +1819 @@
     }
 
-    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));}
+    dbprint(printlevel-voice+2,"// laguerre found roots: "+string(size(li)));
+    int rn@;
 
     for ( j= 1; j <= lis; j++ )
     {
-        if ( pdebug>=1 )
-        {"// root "+string(j);}
+        dbprint(printlevel-voice+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;
+        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
@@ -1707 +1849 @@
     ring r = 0,(x,y),lp;
     // compute the intersection points of two curves
-    ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
-    lex_solve(stdfglm(s),10);
-    rlist;
+    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
+    def R = lex_solve(stdfglm(s),10);
+    setring R; rlist;
 }
 
@@ -1718 +1860 @@
          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.
+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 proc 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
 "
@@ -1736 +1877 @@
     }
 
-    triang_solve(triangLfak(stdfglm(fi)),prec);
-    keepring basering;
+    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
@@ -1744 +1891 @@
     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;
+    ideal s = x2 + y2 - 10, x2 + xy + 2y2 - 16;
+    def R = triangLf_solve(s,10);
+    setring R; rlist;
 }
 
@@ -1756 +1903 @@
 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
+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 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.
@@ -1774 +1919 @@
     }
 
-    triang_solve(triangM(stdfglm(fi)),prec);
-    keepring basering;
+    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
@@ -1782 +1933 @@
     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;
+    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
+    def R = triangM_solve(s,10);
+    setring R; rlist;
 }
 
@@ -1793 +1944 @@
          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 
+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 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.
@@ -1811 +1960 @@
     }
 
-    triang_solve(triangL(stdfglm(fi)),prec);
-    keepring basering;
+    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
@@ -1819 +1974 @@
     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;
+    ideal s =  x2 + y2 - 10, x2 + xy + 2y2 - 16;
+    def R = triangL_solve(s,10);
+    setring R; rlist;
 }
 
@@ -1828 +1983 @@
 
 proc triang_solve( list lfi, int prec, list # )
-"USAGE:   triang_solve(l,p [, d] ); l=list, p,d=integers,@*
+"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.@*
@@ -1837 +1992 @@
          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.@*
+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
 "
 {
-    if ( !defined(pdebug) )
-    {
-        int pdebug;
-        export pdebug;
-    }
-    pdebug=0;
-
-    string orings= nameof(basering);
     def oring= basering;
+    list LL;
 
     // change the ground field to complex numbers
-    string nrings= "ring "+orings+"C=(real,"+string(prec)
+    string nrings= "ring RC =(complex,"+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
@@ -1876 +2013 @@
     {
         int myCompDigits;
-        export  myCompDigits;
     }
     if ( size(#)>=1  && typeof(#[1])=="int" )
@@ -1887 +2023 @@
     }
 
-    if ( pdebug>=1 )
-    {"// myCompDigits="+string(myCompDigits);}
+    dbprint( printlevel-voice+2,"// myCompDigits="+string(myCompDigits) );
 
     int idelem;
     int nv= nvars(basering);
-    int i,j,k,lis;
+    int i,j,lis;
     list resu,li;
 
@@ -1901 +2036 @@
     }
 
-    if ( !defined(rn@) )
-    {
-        int rn@;
-        export rn@;
-    }
-    rn@=0;
+    int rn@=0;
 
     // map the list
     list lfi= imap(oring,lfi);
-
     int slfi= size(lfi);
+
     ideal fi;
-
     for ( i= 1; i <= slfi; i++ )
     {
@@ -1925 +2054 @@
         lis= size(li);
 
-        if ( pdebug>=1 )
-        {"// laguerre found roots: "+string(size(li));}
+        dbprint( printlevel-voice+2,"// laguerre found roots: "+string(lis) );
 
         for ( j= 1; j <= lis; j++ )
         {
-            if ( pdebug>=1 )
-            {"// root "+string(j);}
+            dbprint( printlevel-voice+2,"// root "+string(j) );
             rn@++;
             resu[nv]= li[j];
-            psubst( 1, 2, nv-1, resu, fi, idelem, nv, myCompDigits );
+            LL = psubst( 1, 2, nv-1, resu, fi, idelem, nv, myCompDigits, 
+                         rn@, rlist );
+            rlist = LL[1];
+            rn@ = LL[2];
         }
     }
 
-    if ( pdebug>=0 )
-    {"// list of roots: "+nameof(rlist);}
-    // keep the ring and exit
-    keepring basering;
+    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
@@ -1949 +2083 @@
     // compute the intersection points of two curves
     ideal s=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
-    triang_solve(triangLfak(stdfglm(s)),10);
-    rlist;
+    def R=triang_solve(triangLfak(stdfglm(s)),10);
+    setring R; rlist;
 }