Index: Singular/LIB/rootisolation.lib
===================================================================
--- Singular/LIB/rootisolation.lib	(revision ab7016d343470b5a509daea472404e794c96199f)
+++ Singular/LIB/rootisolation.lib	(revision b7e650aa28103f9685a9dd9313ce231b375f1ac5)
@@ -97,5 +97,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,x,lp;
+    ring R = 0,x,dp;
     interval I = bounds(1);
     I;
@@ -199,5 +199,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
 
     box B = list(bounds(0,1),
@@ -265,5 +265,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,x(1..5),lp;
+    ring R = 0,x(1..5),dp;
 
     ivmat A = ivmatInit(4, 5); A;
@@ -280,5 +280,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,x,lp;
+    ring R = 0,x,dp;
 
     ivmat A = ivmatInit(2,3);
@@ -319,5 +319,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,x,lp;
+    ring R = 0,x,dp;
     ivmat A = ivmatInit(2,2);             A;
     A = ivmatSet(A, 1, 2, bounds(1, 2));  A;
@@ -339,5 +339,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,x,lp;
+    ring R = 0,x,dp;
     ivmat A = diagMatrix(2, bounds(1, 2)); A;
 }
@@ -352,5 +352,6 @@
 example
 {
-    ring R = 0,(x,y),lp;
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y),dp;
     unitMatrix(2);
     unitMatrix(3);
@@ -463,5 +464,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
 
     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;
@@ -521,5 +522,5 @@
 
 proc evalJacobianAtBox(ideal I, box B)
-"USAGE: @code{evalJacobianAtBox(I, B)}; @code{I ideal B box}
+"USAGE: @code{evalJacobianAtBox(I, B)}; @code{I ideal, B box}
 RETURN: Jacobian matrix of @code{I} where polynomials are evaluated at @code{B}
 PURPOSE: evaluates each polynomial of the Jacobian matrix of @code{I} using
@@ -546,5 +547,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
     ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
 
@@ -560,9 +561,11 @@
         0 if test is inconclusive;
         box is intersection of Newton step and supplied box if applicable
-NOTE:   rounding is performed on fractions obtained by intersecting to prevent
-        the size of denominators and numerators from increasing dramatically
+NOTE:   rounding is performed on fractions obtained by matrix inversion to
+        prevent the size of denominators and numerators from increasing
+        dramatically
 EXAMPLE: example testPolyBox; tests the above for intersection of ellipses."
 {
     int N = nvars(basering);
+    int isreal = find(charstr(basering), "Float(");
     int i;
 
@@ -614,41 +617,46 @@
         }
 
-        // in this case, fB may have horrible fractions, so try to simplify
-        //B = Bint;
-
-        // attempt simplification of fractions
-        // add check if we work over reals
-        list bb;
-        for (i = 1; i <= N; i++)
-        {
-            lo = B[i][1];
-            up = B[i][2];
-
-            // modify numerators of B to tighten box
-            if (lo < Bint[i][1])
+        if (isreal)
+        {
+            // fraction simplification over real basering is pointless
+            B = Bint;
+        }
+        else
+        {
+            // attempt simplification of fractions
+            // add check if we work over reals
+            list bb;
+            for (i = 1; i <= N; i++)
             {
-                n = denominator(lo);
-                // floor
-                lo = round(Bint[i][1]*n - 1/2)/n;
+                lo = B[i][1];
+                up = B[i][2];
+
+                // modify numerators of B to tighten box
+                if (lo < Bint[i][1])
+                {
+                    n = denominator(lo);
+                    // floor
+                    lo = round(Bint[i][1]*n - 1/2)/n;
+                }
+                if (up > Bint[i][2])
+                {
+                    n = denominator(up);
+                    // ceil
+                    up = round(Bint[i][2]*n + 1/2)/n;
+                }
+
+                // make sure box does not grow
+                if (lo >= B[i][1] && up <= B[i][2])
+                {
+                    bb[i] = bounds(lo, up);
+                }
+                else
+                {
+                    bb[i] = Bint[i];
+                }
             }
-            if (up > Bint[i][2])
-            {
-                n = denominator(up);
-                // ceil
-                up = round(Bint[i][2]*n + 1/2)/n;
-            }
-
-            // make sure box does not grow
-            if (lo >= B[i][1] && up <= B[i][2])
-            {
-                bb[i] = bounds(lo, up);
-            }
-            else
-            {
-                bb[i] = Bint[i];
-            }
-        }
-
-        B = bb;
+
+            B = bb;
+        }
 
         if (laststep) { return(1, B); }
@@ -661,5 +669,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
     ideal I = 2x2-xy+2y2-2, 2x2-3xy+3y2-2;
 
@@ -690,5 +698,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
     interval I1 = bounds(0, 1); I1;
     interval I2 = bounds(0, 1); I2;
@@ -704,10 +712,10 @@
 
 proc rootIsolationNoPreprocessing(ideal I, def start, number eps)
-"USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)}; @code{I ideal, B box/list of boxes,
-        eps number};
+"USAGE:  @code{rootIsolationNoPreprocessing(I, B, eps)}; @code{I ideal,
+        B box/list of boxes, eps number};
 ASSUME: @code{I} is a zero-dimensional radical ideal
-RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which may contain an
-        element of V(@code{I}), i.e. a root and @code{L2} contains boxes which contain
-        exactly one element of V(@code{I})
+RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
+        may contain an element of V(@code{I}), i.e. a root and @code{L2}
+        contains boxes which contain exactly one element of V(@code{I})
 PURPOSE: Given input box(es) @code{start} we try to find all roots of @code{I}
         lying in @code{start} by computing boxes that contain exactly one root.
@@ -798,5 +806,5 @@
     "EXAMPLE:"; echo = 2;
 
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
 
@@ -848,5 +856,5 @@
 {
     "EXAMPLE:"; echo = 2;
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
 
     ideal I = x+1,y-2;
@@ -861,17 +869,22 @@
 //im Moment geht das nur mit eingegebener eliminationsordnung
 proc rootIsolation(ideal I, box start, number eps)
-"USAGE:  @code{rootIsolation(I, start, eps)}; @code{I ideal, start box, eps number}
+"USAGE:  @code{rootIsolation(I, start, eps)}; @code{I ideal, start box,
+        eps number}
 ASSUME: @code{I} is a zero-dimensional radical ideal
-        and @code{basering} is defined with an elimination ordering.
-RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which may contain an
-        element of V(@code{I}), i.e. a root and @code{L2} contains boxes which contain
-        exactly one element of V(@code{I})
-PURPOSE: same as @code{rootIsolationNoPreprocessing}, but speeds up computation by
-        preprocessing starting box
+RETURN: @code{(L1, L2)}, where @code{L1} contains boxes smaller than eps which
+        may contain an element of V(@code{I}), i.e. a root and @code{L2}
+        contains boxes which contain exactly one element of V(@code{I})
+PURPOSE: same as @code{rootIsolationNoPreprocessing}, but speeds up computation
+        by preprocessing starting box
 THEORY: As every root of @code{I} is a root of the polynomials @code{I[i]}, we
         use Groebner elimination to find univariate polynomials for every
-        variable which have these roots as well. Applying root isolation to
-        these univariate polynomials then provides smaller starting boxes which
-        speed up computations in the multivariate case.
+        variable which have these roots as well. Using that @code{I} is
+        zero-dimensional these Groebner bases may be quickly computed using
+        FGLM. Applying root isolation to these univariate polynomials then
+        provides smaller starting boxes which speed up computations in the
+        multivariate case.
+NOTE:   This algorithm and some procedures used therein perform Groebner basis
+        computations in @code{basering}. It is thus advised to define @code{I}
+        w.r.t. a fast monomial ordering.
 EXAMPLE: example rootIsolation; for intersection of two ellipses"
 {
@@ -907,52 +920,64 @@
     }
 
-    // construct single variable ring from basering
+    // compute reduced GB in (hopefully) fast ordering
+    option(redSB);
+    ideal fastGB = groebner(I);
+
     ring rSource = basering;
     list rList = ringlist(rSource);
-    // delete all but first variable
+    // first construct univariate ring
+    int numords = size(rList[3]);
+    // remove all but first variables
     rList[2] = list(rList[2][1]);
-    // adjust weight vector
-    rList[3][1][2] = intvec(1);
-
-    // compute reduced groebner bases
-    option(redSB);
-
-    ideal gbUnivarPolys, gbTemp, Ipermutation;
-    intvec perm;
-
-    // permute every variable to the end once and compute Groebner basis to
-    // find univariate polynomial in that variable
-    for (i = 1; i <= N; i++)
-    {
-        perm = 1..N;
-        perm[i] = N;
-        perm[N] = i;
-
-        Ipermutation = fetch(rSource, I, perm);
-        gbTemp = groebner(Ipermutation);
-
-        // first poly in GB now only depends on last variable
-        gbUnivarPolys[i] = gbTemp[1];
-    }
-
-    // do computations in univariate ring
+    // change ordering accordingly (keep last block)
+    rList[3] = list( list(rList[3][1][1], intvec(1)), rList[3][numords] );
+
+    // construct and switch to univariate ring
     ring rUnivar = ring(rList);
-    setring rUnivar;
-
-    // gbUnivarPoly only contains var(N),
-    // make sure variables are properly mapped
-    ideal gbUnivarPolys = fetch(rSource, gbUnivarPolys, intvec(0:(N-1),1));
+
+    // some necessary variables
+    ideal gbUnivar;
+    // maps var(N) in rElimOrd(i) to var(1) in rUnivar
+    intvec varmap = (0:(N-1)),1;
     number eps = fetch(rSource, eps);
     list univarResult, startBoxesPerDim;
 
+    // now prepare lp-ring
+    setring rSource;
+    rList = ringlist(rSource);
+    // set ordering to lp
+    rList[3] = list( list( "lp", 1:N ), list( "C", 0 ) );
+    // save variable order for later
+    list varstrs = rList[2];
+
     for (i = 1; i <= N; i++)
     {
-        univarResult = rootIsolationNoPreprocessing(ideal(gbUnivarPolys[i]),
+        // permute variables in ringlist,
+        // create and switch to ring with elimination ordering
+        rList[2][i] = varstrs[N];
+        rList[2][N] = varstrs[i];
+        ring rElimOrd = ring(rList);
+        // get GB in elimination ordering, note that GB[1] only depends on
+        // var(i) of rSource, which is var(N) in rElimOrd
+        ideal GB = fglm(rSource, fastGB);
+        GB = GB[1];
+
+        setring rUnivar;
+        gbUnivar = fetch(rElimOrd, GB, varmap);
+        // clean up ring and its elements
+        kill rElimOrd;
+
+        // get boxes containing roots of univariate polynomial
+        univarResult = rootIsolationNoPreprocessing(gbUnivar,
             box(list(start[i])), eps);
         // maybe result[1] is not empty, so take both
         startBoxesPerDim[i] = univarResult[1] + univarResult[2];
         // debug:
-        print(string("Sieved variable ", varstr(rSource, i), " to ",
+        print(string("Sieved variable ", varstrs[i], " to ",
             size(startBoxesPerDim[i]), " intervals."));
+
+        // reset ring variable order
+        setring rSource;
+        rList[2] = varstrs;
     }
 
@@ -973,4 +998,5 @@
     list startBoxes, sbTemp;
 
+    // prepare a list of lists
     for (i = 1; i <= numBoxes; i++)
     {
@@ -997,6 +1023,5 @@
     }
 
-    // switch back for final computation and so that boxes have proper size
-    setring rSource;
+    // since we're back in rSource box(...) will return box of proper size
     for (i = 1; i <= size(sbTemp); i++)
     {
@@ -1010,5 +1035,5 @@
     "EXAMPLE:"; echo = 2;
 
-    ring R = 0,(x,y),lp;
+    ring R = 0,(x,y),dp;
     ideal I = 2x2-xy+2y2-2,2x2-3xy+3y2-2;  // V(I) has four elements
 
@@ -1016,7 +1041,5 @@
     box B = list(i, i);
 
-    list result = rootIsolation(I, B, 1/512);
-    size(result[1]);
-    size(result[2]);
+    list result = rootIsolation(I, B, 0);
 
     result;
