Update to 2Geom rev. 1113

(bzr r3622)

1 files changed, 147 insertions, 147 deletions

diff --git a/src/2geom/poly-laguerre-solve.cpp b/src/2geom/poly-laguerre-solve.cpp
index 83f049286..fb3c2075b 100644
--- a/src/2geom/poly-laguerre-solve.cpp
+++ b/src/2geom/poly-laguerre-solve.cpp
@@ -1,147 +1,147 @@
-#include "poly-laguerre-solve.h"
-#include <iterator>
-
-typedef std::complex<double> cdouble;
-
-cdouble laguerre_internal_complex(Poly const & p,
-                                  double x0,
-                                  double tol,
-                                  bool & quad_root) {
-    cdouble a = 2*tol;
-    cdouble xk = x0;
-    double n = p.degree();
-    quad_root = false;
-    const unsigned shuffle_rate = 10;
-    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
-    unsigned shuffle_counter = 0;
-    while(std::norm(a) > (tol*tol)) {
-        //std::cout << "xk = " << xk << std::endl;
-        cdouble b = p.back();
-        cdouble d = 0, f = 0;
-        double err = abs(b);
-        double abx = abs(xk);
-        for(int j = p.size()-2; j >= 0; j--) {
-            f = xk*f + d;
-            d = xk*d + b;
-            b = xk*b + p[j];
-            err = abs(b) + abx*err;
-        }
-        
-        err *= 1e-7; // magic epsilon for convergence, should be computed from tol
-        
-        cdouble px = b;
-        if(abs(b) < err)
-            return xk;
-        //if(std::norm(px) < tol*tol)
-        //    return xk;
-        cdouble G = d / px;
-        cdouble H = G*G - f / px;
-        
-        //std::cout << "G = " << G << "H = " << H;
-        cdouble radicand = (n - 1)*(n*H-G*G);
-        //assert(radicand.real() > 0);
-        if(radicand.real() < 0)
-            quad_root = true;
-        //std::cout << "radicand = " << radicand << std::endl;
-        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
-            a = - sqrt(radicand);
-        else
-            a = sqrt(radicand);
-        //std::cout << "a = " << a << std::endl;
-        a = n / (a + G);
-        //std::cout << "a = " << a << std::endl;
-        if(shuffle_counter % shuffle_rate == 0)
-            ;//a *= shuffle[shuffle_counter / shuffle_rate];
-        xk -= a;
-        shuffle_counter++;
-        if(shuffle_counter >= 90)
-            break;
-    }
-    //std::cout << "xk = " << xk << std::endl;
-    return xk;
-}
-
-double laguerre_internal(Poly const & p,
-                         Poly const & pp,
-                         Poly const & ppp,
-                         double x0,
-                         double tol,
-                         bool & quad_root) {
-    double a = 2*tol;
-    double xk = x0;
-    double n = p.degree();
-    quad_root = false;
-    while(a*a > (tol*tol)) {
-        //std::cout << "xk = " << xk << std::endl;
-        double px = p(xk);
-        if(px*px < tol*tol)
-            return xk;
-        double G = pp(xk) / px;
-        double H = G*G - ppp(xk) / px;
-        
-        //std::cout << "G = " << G << "H = " << H;
-        double radicand = (n - 1)*(n*H-G*G);
-        assert(radicand > 0);
-        //std::cout << "radicand = " << radicand << std::endl;
-        if(G < 0) // here we try to maximise the denominator avoiding cancellation
-            a = - sqrt(radicand);
-        else
-            a = sqrt(radicand);
-        //std::cout << "a = " << a << std::endl;
-        a = n / (a + G);
-        //std::cout << "a = " << a << std::endl;
-        xk -= a;
-    }
-    //std::cout << "xk = " << xk << std::endl;
-    return xk;
-}
-
-
-std::vector<cdouble > 
-laguerre(Poly p, const double tol) {
-    std::vector<cdouble > solutions;
-    //std::cout << "p = " << p << " = ";
-    while(p.size() > 1)
-    {
-        double x0 = 0;
-        bool quad_root = false;
-        cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
-        //if(abs(sol) > 1) break;
-        Poly dvs;
-        if(quad_root) {
-            dvs.push_back((sol*conj(sol)).real());
-            dvs.push_back(-(sol + conj(sol)).real());
-            dvs.push_back(1.0);
-            //std::cout << "(" <<  dvs << ")";
-	    //solutions.push_back(sol);
-	    //solutions.push_back(conj(sol));
-        } else {
-            //std::cout << sol << std::endl;
-            dvs.push_back(-sol.real());
-            dvs.push_back(1.0);
-	    solutions.push_back(sol);
-            //std::cout << "(" <<  dvs << ")";
-        }
-        Poly r;
-        p = divide(p, dvs, r);
-        //std::cout << r << std::endl;
-    }
-    return solutions;
-}
-
-std::vector<double> 
-laguerre_real_interval(Poly const & ply, 
-		       const double lo, const double hi,
-		       const double tol) {
-}
-
-/*
-  Local Variables:
-  mode:c++
-  c-file-style:"stroustrup"
-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
-  indent-tabs-mode:nil
-  fill-column:99
-  End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+#include "poly-laguerre-solve.h"

+#include <iterator>

+

+typedef std::complex<double> cdouble;

+

+cdouble laguerre_internal_complex(Poly const & p,

+                                  double x0,

+                                  double tol,

+                                  bool & quad_root) {

+    cdouble a = 2*tol;

+    cdouble xk = x0;

+    double n = p.degree();

+    quad_root = false;

+    const unsigned shuffle_rate = 10;

+    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};

+    unsigned shuffle_counter = 0;

+    while(std::norm(a) > (tol*tol)) {

+        //std::cout << "xk = " << xk << std::endl;

+        cdouble b = p.back();

+        cdouble d = 0, f = 0;

+        double err = abs(b);

+        double abx = abs(xk);

+        for(int j = p.size()-2; j >= 0; j--) {

+            f = xk*f + d;

+            d = xk*d + b;

+            b = xk*b + p[j];

+            err = abs(b) + abx*err;

+        }

+        

+        err *= 1e-7; // magic epsilon for convergence, should be computed from tol

+        

+        cdouble px = b;

+        if(abs(b) < err)

+            return xk;

+        //if(std::norm(px) < tol*tol)

+        //    return xk;

+        cdouble G = d / px;

+        cdouble H = G*G - f / px;

+        

+        //std::cout << "G = " << G << "H = " << H;

+        cdouble radicand = (n - 1)*(n*H-G*G);

+        //assert(radicand.real() > 0);

+        if(radicand.real() < 0)

+            quad_root = true;

+        //std::cout << "radicand = " << radicand << std::endl;

+        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation

+            a = - sqrt(radicand);

+        else

+            a = sqrt(radicand);

+        //std::cout << "a = " << a << std::endl;

+        a = n / (a + G);

+        //std::cout << "a = " << a << std::endl;

+        if(shuffle_counter % shuffle_rate == 0)

+            ;//a *= shuffle[shuffle_counter / shuffle_rate];

+        xk -= a;

+        shuffle_counter++;

+        if(shuffle_counter >= 90)

+            break;

+    }

+    //std::cout << "xk = " << xk << std::endl;

+    return xk;

+}

+

+double laguerre_internal(Poly const & p,

+                         Poly const & pp,

+                         Poly const & ppp,

+                         double x0,

+                         double tol,

+                         bool & quad_root) {

+    double a = 2*tol;

+    double xk = x0;

+    double n = p.degree();

+    quad_root = false;

+    while(a*a > (tol*tol)) {

+        //std::cout << "xk = " << xk << std::endl;

+        double px = p(xk);

+        if(px*px < tol*tol)

+            return xk;

+        double G = pp(xk) / px;

+        double H = G*G - ppp(xk) / px;

+        

+        //std::cout << "G = " << G << "H = " << H;

+        double radicand = (n - 1)*(n*H-G*G);

+        assert(radicand > 0);

+        //std::cout << "radicand = " << radicand << std::endl;

+        if(G < 0) // here we try to maximise the denominator avoiding cancellation

+            a = - sqrt(radicand);

+        else

+            a = sqrt(radicand);

+        //std::cout << "a = " << a << std::endl;

+        a = n / (a + G);

+        //std::cout << "a = " << a << std::endl;

+        xk -= a;

+    }

+    //std::cout << "xk = " << xk << std::endl;

+    return xk;

+}

+

+

+std::vector<cdouble > 

+laguerre(Poly p, const double tol) {

+    std::vector<cdouble > solutions;

+    //std::cout << "p = " << p << " = ";

+    while(p.size() > 1)

+    {

+        double x0 = 0;

+        bool quad_root = false;

+        cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);

+        //if(abs(sol) > 1) break;

+        Poly dvs;

+        if(quad_root) {

+            dvs.push_back((sol*conj(sol)).real());

+            dvs.push_back(-(sol + conj(sol)).real());

+            dvs.push_back(1.0);

+            //std::cout << "(" <<  dvs << ")";

+	    //solutions.push_back(sol);

+	    //solutions.push_back(conj(sol));

+        } else {

+            //std::cout << sol << std::endl;

+            dvs.push_back(-sol.real());

+            dvs.push_back(1.0);

+	    solutions.push_back(sol);

+            //std::cout << "(" <<  dvs << ")";

+        }

+        Poly r;

+        p = divide(p, dvs, r);

+        //std::cout << r << std::endl;

+    }

+    return solutions;

+}

+

+std::vector<double> 

+laguerre_real_interval(Poly const & ply, 

+		       const double lo, const double hi,

+		       const double tol) {

+}

+

+/*

+  Local Variables:

+  mode:c++

+  c-file-style:"stroustrup"

+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

+  indent-tabs-mode:nil

+  fill-column:99

+  End:

+*/

+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :