update to latest 2geom

diff --git a/src/2geom/angle.h b/src/2geom/angle.h
index 05364fb222b57187767315b83223c0c41d3a8c80..4d548ab49b5750a1a5c6a72e514efd5d852c83ab 100644 (file)
--- a/src/2geom/angle.h
+++ b/src/2geom/angle.h
@@ -1,50 +1,50 @@
-/**
- *  \file angle.h
- *  \brief Various trigoniometric helper functions
- *
- *  Authors:
- *   Johan Engelen <goejendaagh@zonnet.nl>
- *
- * Copyright (C) 2007 authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- *
- */
-
-#ifndef LIB2GEOM_SEEN_ANGLE_H
-#define LIB2GEOM_SEEN_ANGLE_H
-
-namespace Geom {
-
-#ifndef M_PI
-# define M_PI 3.14159265358979323846
-#endif
-
-inline double deg_to_rad(double deg) { return deg*M_PI/180.0;}
-
-inline double rad_to_deg(double rad) { return rad*180.0/M_PI;}
-
-}
-
-#endif
+/**\r
+ *  \file angle.h\r
+ *  \brief Various trigoniometric helper functions\r
+ *\r
+ *  Authors:\r
+ *   Johan Engelen <goejendaagh@zonnet.nl>\r
+ *\r
+ * Copyright (C) 2007 authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ *\r
+ */\r
+ \r
+#ifndef LIB2GEOM_SEEN_ANGLE_H\r
+#define LIB2GEOM_SEEN_ANGLE_H\r
+\r
+namespace Geom {\r
+\r
+#ifndef M_PI\r
+# define M_PI 3.14159265358979323846\r
+#endif\r
+\r
+inline double deg_to_rad(double deg) { return deg*M_PI/180.0;}\r
+\r
+inline double rad_to_deg(double rad) { return rad*180.0/M_PI;}\r
+\r
+}\r
+\r
+#endif\r

diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp
index 6d7c772c067a1dc693ffb730f04eb8517bb6200f..f90bb6d420543a5c9f17b77ec663d60226b446a2 100644 (file)
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
@@ -108,9 +108,9 @@ void Matrix::setIdentity() {
 //TODO: use eps
 
 bool Matrix::isIdentity(Coord const eps) const {
-    return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) &&
-           are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
-           are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 }
 
 /** Answers the question "Does this matrix perform a translation, and \em{only} a translation?"
@@ -118,9 +118,9 @@ bool Matrix::isIdentity(Coord const eps) const {
  \return A bool representing yes/no.
  */
 bool Matrix::isTranslation(Coord const eps) const {
-    return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) &&
-           are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
-           (!are_near(_c[4], 0.0) || !are_near(_c[5], 0.0));
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           (!are_near(_c[4], 0.0, eps) || !are_near(_c[5], 0.0, eps));
 }
 
 /** Answers the question "Does this matrix perform a scale, and \em{only} a Scale?"
@@ -128,9 +128,9 @@ bool Matrix::isTranslation(Coord const eps) const {
  \return A bool representing yes/no.
  */
 bool Matrix::isScale(Coord const eps) const {
-    return !are_near(_c[0], 1.0) || !are_near(_c[3], 1.0) &&  //NOTE: these are the diags, and the next line opposite diags
-           are_near(_c[1], 0.0) && are_near(_c[2], 0.0) && 
-           are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
+    return !are_near(_c[0], 1.0, eps) || !are_near(_c[3], 1.0, eps) &&  //NOTE: these are the diags, and the next line opposite diags
+           are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && 
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 }
 
 /** Answers the question "Does this matrix perform a uniform scale, and \em{only} a uniform scale?"
@@ -138,9 +138,9 @@ bool Matrix::isScale(Coord const eps) const {
  \return A bool representing yes/no.
  */
 bool Matrix::isUniformScale(Coord const eps) const {
-    return !are_near(_c[0], 1.0) && are_near(_c[0], _c[3]) &&
-           are_near(_c[1], 0.0) && are_near(_c[2], 0.0) &&  
-           are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
+    return !are_near(_c[0], 1.0, eps) && are_near(_c[0], _c[3], eps) &&
+           are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) &&  
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 }
 
 /** Answers the question "Does this matrix perform a rotation, and \em{only} a rotation?"
@@ -148,13 +148,17 @@ bool Matrix::isUniformScale(Coord const eps) const {
  \return A bool representing yes/no.
  */
 bool Matrix::isRotation(Coord const eps) const {
-    return !are_near(_c[0], _c[3]) && are_near(_c[1], -_c[2]) &&
-           are_near(_c[4], 0.0) && are_near(_c[5], 0.0) &&
-           are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0);
+    return !are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps) &&
+           are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps);
 }
 
 bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
-    return are_near(_c[0], _c[3]) && are_near(_c[1], 0) && are_near(_c[2], 0);
+    return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps);
+}
+
+bool Matrix::isSingular(Coord const eps) const {
+    return are_near(det(), 0.0, eps);
 }
 
 bool Matrix::flips() const {

diff --git a/src/2geom/matrix.h b/src/2geom/matrix.h
index c39c997166a51ef4a1a6a2c71f4b6abddecfe061..c9f244d6243ac7a8e1a38cc56e9f913ecb577aae 100644 (file)
--- a/src/2geom/matrix.h
+++ b/src/2geom/matrix.h
@@ -91,6 +91,7 @@ class Matrix {
     bool isScale(double eps = EPSILON) const;
     bool isUniformScale(double eps = EPSILON) const;
     bool onlyScaleAndTranslation(double eps = EPSILON) const;
+    bool isSingular(double eps = EPSILON) const;
 
     bool flips() const;
 

diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp
index fdc1cefe55a585d2205a63fcf925b1a74a587058..87d238f14294c2692924ec7b9d1c97d53e19807d 100644 (file)
--- a/src/2geom/poly-dk-solve.cpp
+++ b/src/2geom/poly-dk-solve.cpp
@@ -1,64 +1,64 @@
-#include "poly-dk-solve.h"
-#include <iterator>
-
-/*** implementation of the Durand-Kerner method.  seems buggy*/
-
-std::complex<double> evalu(Poly const & p, std::complex<double> x) {
-    std::complex<double> result = 0;
-    std::complex<double> xx = 1;
-
-    for(unsigned i = 0; i < p.size(); i++) {
-        result += p[i]*xx;
-        xx *= x;
-    }
-    return result;
-}
-
-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
-    std::vector<std::complex<double> > roots;
-    const int N = ply.degree();
-
-    std::complex<double> b(0.4, 0.9);
-    std::complex<double> p = 1;
-    for(int i = 0; i < N; i++) {
-        roots.push_back(p);
-        p *= b;
-    }
-    assert(roots.size() == ply.degree());
-
-    double error = 0;
-    int i;
-    for( i = 0; i < 30; i++) {
-        error = 0;
-        for(int r_i = 0; r_i < N; r_i++) {
-            std::complex<double> denom = 1;
-            std::complex<double> R = roots[r_i];
-            for(int d_i = 0; d_i < N; d_i++) {
-                if(r_i != d_i)
-                    denom *= R-roots[d_i];
-            }
-            assert(norm(denom) != 0);
-            std::complex<double> dr = evalu(ply, R)/denom;
-            error += norm(dr);
-            roots[r_i] = R - dr;
-        }
-        /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
-          std::cout << std::endl;*/
-        if(error < tol)
-            break;
-    }
-    //std::cout << error << ", " << i<< std::endl;
-    return roots;
-}
-
-
-/*
-  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-dk-solve.h"\r
+#include <iterator>\r
+\r
+/*** implementation of the Durand-Kerner method.  seems buggy*/\r
+\r
+std::complex<double> evalu(Poly const & p, std::complex<double> x) {\r
+    std::complex<double> result = 0;\r
+    std::complex<double> xx = 1;\r
+\r
+    for(unsigned i = 0; i < p.size(); i++) {\r
+        result += p[i]*xx;\r
+        xx *= x;\r
+    }\r
+    return result;\r
+}\r
+\r
+std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {\r
+    std::vector<std::complex<double> > roots;\r
+    const int N = ply.degree();\r
+\r
+    std::complex<double> b(0.4, 0.9);\r
+    std::complex<double> p = 1;\r
+    for(int i = 0; i < N; i++) {\r
+        roots.push_back(p);\r
+        p *= b;\r
+    }\r
+    assert(roots.size() == ply.degree());\r
+\r
+    double error = 0;\r
+    int i;\r
+    for( i = 0; i < 30; i++) {\r
+        error = 0;\r
+        for(int r_i = 0; r_i < N; r_i++) {\r
+            std::complex<double> denom = 1;\r
+            std::complex<double> R = roots[r_i];\r
+            for(int d_i = 0; d_i < N; d_i++) {\r
+                if(r_i != d_i)\r
+                    denom *= R-roots[d_i];\r
+            }\r
+            assert(norm(denom) != 0);\r
+            std::complex<double> dr = evalu(ply, R)/denom;\r
+            error += norm(dr);\r
+            roots[r_i] = R - dr;\r
+        }\r
+        /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));\r
+          std::cout << std::endl;*/\r
+        if(error < tol)\r
+            break;\r
+    }\r
+    //std::cout << error << ", " << i<< std::endl;\r
+    return roots;\r
+}\r
+\r
+\r
+/*\r
+  Local Variables:\r
+  mode:c++\r
+  c-file-style:"stroustrup"\r
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+  indent-tabs-mode:nil\r
+  fill-column:99\r
+  End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r

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

diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index b9ef71b5c0cd3cd165a066bc93599dbdee536a71..c0e1c399713c3a6252b1816d10fab04a8c7ef715 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
@@ -1,105 +1,105 @@
-#include "sweep.h"
-
-#include <algorithm>
-
-namespace Geom {
-
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {
-    std::vector<Event> events; events.reserve(rs.size()*2);
-    std::vector<std::vector<unsigned> > pairs(rs.size());
-
-    for(unsigned i = 0; i < rs.size(); i++) {
-        events.push_back(Event(rs[i].left(), i, false));
-        events.push_back(Event(rs[i].right(), i, true));
-    }
-    std::sort(events.begin(), events.end());
-
-    std::vector<unsigned> open;
-    for(unsigned i = 0; i < events.size(); i++) {
-        unsigned ix = events[i].ix;
-        if(events[i].closing) {
-            std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);
-            //if(iter != open.end())
-            open.erase(iter);
-        } else {
-            for(unsigned j = 0; j < open.size(); j++) {
-                unsigned jx = open[j];
-                if(rs[jx][Y].intersects(rs[ix][Y])) {
-                    pairs[jx].push_back(ix);
-                }
-            }
-            open.push_back(ix);
-        }
-    }
-    return pairs;
-}
-
-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {
-    std::vector<std::vector<unsigned> > pairs(a.size());
-    if(a.empty() || b.empty()) return pairs;
-    std::vector<Event> events[2];
-    events[0].reserve(a.size()*2);
-    events[1].reserve(b.size()*2);
-
-    for(unsigned n = 0; n < 2; n++) {
-        unsigned sz = n ? b.size() : a.size();
-        events[n].reserve(sz*2);
-        for(unsigned i = 0; i < sz; i++) {
-            events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));
-            events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));
-        }
-        std::sort(events[n].begin(), events[n].end());
-    }
-
-    std::vector<unsigned> open[2];
-    bool n = events[1].front() < events[0].front();
-    for(unsigned i[] = {0,0}; i[n] < events[n].size();) {
-        unsigned ix = events[n][i[n]].ix;
-        bool closing = events[n][i[n]].closing;
-        //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";
-        if(closing) {
-            open[n].erase(std::find(open[n].begin(), open[n].end(), ix));
-        } else {
-            if(n) {
-                //n = 1
-                //opening a B, add to all open a
-                for(unsigned j = 0; j < open[0].size(); j++) {
-                    unsigned jx = open[0][j];
-                    if(a[jx][Y].intersects(b[ix][Y])) {
-                        pairs[jx].push_back(ix);
-                    }
-                }
-            } else {
-                //n = 0
-                //opening an A, add all open b
-                for(unsigned j = 0; j < open[1].size(); j++) {
-                    unsigned jx = open[1][j];
-                    if(b[jx][Y].intersects(a[ix][Y])) {
-                        pairs[ix].push_back(jx);
-                    }
-                }
-            }
-            open[n].push_back(ix);
-        }
-        i[n]++;
-       if(i[n]>=events[n].size()) {break;}
-        n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;
-    }
-    return pairs;
-}
-
-//Fake cull, until the switch to the real sweep is made.
-std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {
-    std::vector<std::vector<unsigned> > ret;
-
-    std::vector<unsigned> all;
-    for(unsigned j = 0; j < b; j++)
-        all.push_back(j);
-
-    for(unsigned i = 0; i < a; i++)
-        ret.push_back(all);
-
-    return ret;
-}
-
-}
+#include "sweep.h"\r
+\r
+#include <algorithm>\r
+\r
+namespace Geom {\r
+\r
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {\r
+    std::vector<Event> events; events.reserve(rs.size()*2);\r
+    std::vector<std::vector<unsigned> > pairs(rs.size());\r
+\r
+    for(unsigned i = 0; i < rs.size(); i++) {\r
+        events.push_back(Event(rs[i].left(), i, false));\r
+        events.push_back(Event(rs[i].right(), i, true));\r
+    }\r
+    std::sort(events.begin(), events.end());\r
+\r
+    std::vector<unsigned> open;\r
+    for(unsigned i = 0; i < events.size(); i++) {\r
+        unsigned ix = events[i].ix;\r
+        if(events[i].closing) {\r
+            std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);\r
+            //if(iter != open.end())\r
+            open.erase(iter);\r
+        } else {\r
+            for(unsigned j = 0; j < open.size(); j++) {\r
+                unsigned jx = open[j];\r
+                if(rs[jx][Y].intersects(rs[ix][Y])) {\r
+                    pairs[jx].push_back(ix);\r
+                }\r
+            }\r
+            open.push_back(ix);\r
+        }\r
+    }\r
+    return pairs;\r
+}\r
+\r
+std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {\r
+    std::vector<std::vector<unsigned> > pairs(a.size());\r
+    if(a.empty() || b.empty()) return pairs;\r
+    std::vector<Event> events[2];\r
+    events[0].reserve(a.size()*2);\r
+    events[1].reserve(b.size()*2);\r
+\r
+    for(unsigned n = 0; n < 2; n++) {\r
+        unsigned sz = n ? b.size() : a.size();\r
+        events[n].reserve(sz*2);\r
+        for(unsigned i = 0; i < sz; i++) {\r
+            events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));\r
+            events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));\r
+        }\r
+        std::sort(events[n].begin(), events[n].end());\r
+    }\r
+\r
+    std::vector<unsigned> open[2];\r
+    bool n = events[1].front() < events[0].front();\r
+    for(unsigned i[] = {0,0}; i[n] < events[n].size();) {\r
+        unsigned ix = events[n][i[n]].ix;\r
+        bool closing = events[n][i[n]].closing;\r
+        //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";\r
+        if(closing) {\r
+            open[n].erase(std::find(open[n].begin(), open[n].end(), ix));\r
+        } else {\r
+            if(n) {\r
+                //n = 1\r
+                //opening a B, add to all open a\r
+                for(unsigned j = 0; j < open[0].size(); j++) {\r
+                    unsigned jx = open[0][j];\r
+                    if(a[jx][Y].intersects(b[ix][Y])) {\r
+                        pairs[jx].push_back(ix);\r
+                    }\r
+                }\r
+            } else {\r
+                //n = 0\r
+                //opening an A, add all open b\r
+                for(unsigned j = 0; j < open[1].size(); j++) {\r
+                    unsigned jx = open[1][j];\r
+                    if(b[jx][Y].intersects(a[ix][Y])) {\r
+                        pairs[ix].push_back(jx);\r
+                    }\r
+                }\r
+            }\r
+            open[n].push_back(ix);\r
+        }\r
+        i[n]++;\r
+       if(i[n]>=events[n].size()) {break;}\r
+        n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;\r
+    }\r
+    return pairs;\r
+}\r
+\r
+//Fake cull, until the switch to the real sweep is made.\r
+std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {\r
+    std::vector<std::vector<unsigned> > ret;\r
+\r
+    std::vector<unsigned> all;\r
+    for(unsigned j = 0; j < b; j++)\r
+        all.push_back(j);\r
+\r
+    for(unsigned i = 0; i < a; i++)\r
+        ret.push_back(all);\r
+\r
+    return ret;\r
+}\r
+\r
+}\r