Merge from trunk

(bzr r9508.1.73)

7 files changed, 1414 insertions, 20 deletions

diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
index dcf0e8e1d..a44c1ddac 100644
--- a/src/2geom/numeric/fitting-model.h
+++ b/src/2geom/numeric/fitting-model.h
@@ -43,6 +43,7 @@
 #include <2geom/ellipse.h>
 #include <2geom/circle.h>
 #include <2geom/utils.h>
+#include <2geom/conicsec.h>
 
 
 namespace Geom { namespace NL {
@@ -83,6 +84,10 @@ namespace Geom { namespace NL {
  *   instance type:  the type of the objects produced by using
  *                   the fitting raw data solution
  */
+
+
+
+
 template< typename ParameterType, typename ValueType, typename InstanceType >
 class LinearFittingModel
 {
@@ -230,8 +235,32 @@ class LFMNormalizedPowerBasis
 // incomplete model, it can be inherited to make up different kinds of
 // instance type; the raw data is a vector of coefficients of the equation
 // of an ellipse curve
+//template< typename InstanceType >
+//class LFMEllipseEquation
+//    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+//{
+//  public:
+//    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+//    {
+//        coeff[0] = p[X] * p[Y];
+//        coeff[1] = p[Y] * p[Y];
+//        coeff[2] = p[X];
+//        coeff[3] = p[Y];
+//        coeff[4] = 1;
+//        fixed_term = p[X] * p[X];
+//    }
+//
+//    size_t size() const
+//    {
+//        return 5;
+//    }
+//};
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of a conic section
 template< typename InstanceType >
-class LFMEllipseEquation
+class LFMConicEquation
     : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
 {
   public:
@@ -251,10 +280,20 @@ class LFMEllipseEquation
     }
 };
 
+// this model generates Ellipse curves
+class LFMConicSection
+    : public LFMConicEquation<xAx>
+{
+  public:
+    void instance(xAx & c, ConstVectorView const& coeff) const
+    {
+        c.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+    }
+};
 
 // this model generates Ellipse curves
 class LFMEllipse
-    : public LFMEllipseEquation<Ellipse>
+    : public LFMConicEquation<Ellipse>
 {
   public:
     void instance(Ellipse & e, ConstVectorView const& coeff) const
@@ -431,13 +470,13 @@ class LFMBezier
 
 
 // this model generates Bezier curves
-template< unsigned int N >
-class LFMBezierCurve
-    : public LinearFittingModel< double, Point, BezierCurve<N> >
+template <unsigned degree>
+class LFMBezierCurveN
+    : public LinearFittingModel< double, Point, BezierCurveN<degree> >
 {
   public:
-    LFMBezierCurve( size_t _order )
-        : mob(_order)
+    LFMBezierCurveN()
+        : mob(degree+1)
     {
     }
 
@@ -451,13 +490,13 @@ class LFMBezierCurve
         return mob.size();
     }
 
-    void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const
+    void instance(BezierCurveN<degree> & bc, ConstMatrixView const& raw_data) const
     {
-        Bezier bx(size()-1);
-        Bezier by(size()-1);
+        Bezier bx(degree);
+        Bezier by(degree);
         mob.instance(bx, raw_data.column_const_view(X));
         mob.instance(by, raw_data.column_const_view(Y));
-        bc = BezierCurve<N>(bx, by);
+        bc = BezierCurveN<degree>(bx, by);
     }
 
   private:
diff --git a/src/2geom/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp
index 94a345fd5..98ff3b6ca 100644
--- a/src/2geom/numeric/matrix.cpp
+++ b/src/2geom/numeric/matrix.cpp
@@ -37,10 +37,8 @@
 #include <2geom/numeric/vector.h>
 
 
-
 namespace Geom { namespace NL {
 
-
 Vector operator*( detail::BaseMatrixImpl const& A,
                   detail::BaseVectorImpl const& v )
 {
@@ -101,6 +99,47 @@ Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
     return P;
 }
 
+
+double trace (detail::BaseMatrixImpl const& A)
+{
+    if (A.rows() != A.columns())
+    {
+        THROW_RANGEERROR ("NL::Matrix: computing trace: "
+                          "rows() != columns()");
+    }
+    double t = 0;
+    for (size_t i = 0; i < A.rows(); ++i)
+    {
+        t += A(i,i);
+    }
+    return t;
+}
+
+
+double det (detail::BaseMatrixImpl const& A)
+{
+    if (A.rows() != A.columns())
+    {
+        THROW_RANGEERROR ("NL::Matrix: computing determinant: "
+                          "rows() != columns()");
+    }
+
+    Matrix LU(A);
+    int s;
+    gsl_permutation * p = gsl_permutation_alloc(LU.rows());
+    gsl_linalg_LU_decomp (LU.get_gsl_matrix(), p, &s);
+
+    double t = 1;
+    for (size_t i = 0; i < LU.rows(); ++i)
+    {
+        t *= LU(i,i);
+    }
+
+    gsl_permutation_free(p);
+    return t;
+}
+
+
 } }  // end namespaces
 
 /*
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
index f2a934235..a130bd748 100644
--- a/src/2geom/numeric/matrix.h
+++ b/src/2geom/numeric/matrix.h
@@ -38,6 +38,7 @@
 #ifndef _NL_MATRIX_H_
 #define _NL_MATRIX_H_
 
+#include <2geom/exception.h>
 #include <2geom/numeric/vector.h>
 
 #include <cassert>
@@ -45,7 +46,6 @@
 #include <algorithm>  // for std::swap
 #include <sstream>
 #include <string>
-
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_linalg.h>
 
@@ -232,13 +232,13 @@ class MatrixImpl : public BaseMatrixImpl
 		gsl_matrix_set_identity(m_matrix);
 	}
 
-        using base_type::operator(); // VSC legacy support
-        const double & operator() (size_t i, size_t j) const
-        {
-            return base_type::operator ()(i, j);
-        }
+    using base_type::operator(); // VSC legacy support
+    const double & operator() (size_t i, size_t j) const
+    {
+        return base_type::operator ()(i, j);
+    }
 
-        double & operator() (size_t i, size_t j)
+    double & operator() (size_t i, size_t j)
 	{
             return *gsl_matrix_ptr(m_matrix, i, j);
 	}
@@ -310,6 +310,8 @@ using detail::operator==;
 using detail::operator<<;
 
 
+template <size_t N>
+class ConstBaseSymmetricMatrix;
 
 
 class Matrix: public detail::MatrixImpl
@@ -352,6 +354,18 @@ class Matrix: public detail::MatrixImpl
 		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 	}
 
+	template <size_t N>
+	explicit
+    Matrix(ConstBaseSymmetricMatrix<N> const& _smatrix)
+	{
+	    m_rows = N;
+	    m_columns = N;
+	    m_matrix = gsl_matrix_alloc(N, N);
+	    for (size_t i = 0; i < N; ++i)
+	        for (size_t j = 0; j < N ; ++j)
+	            (*gsl_matrix_ptr(m_matrix, i, j)) = _smatrix(i,j);
+	}
+
 	Matrix & operator=(Matrix const& _matrix)
 	{
 		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
@@ -366,6 +380,16 @@ class Matrix: public detail::MatrixImpl
 		return *this;
 	}
 
+	template <size_t N>
+	Matrix & operator=(ConstBaseSymmetricMatrix<N> const& _smatrix)
+	{
+	    assert (rows() == N && columns() ==  N);
+	    for (size_t i = 0; i < N; ++i)
+	        for (size_t j = 0; j < N ; ++j)
+	            (*this)(i,j) = _smatrix(i,j);
+	    return *this;
+	}
+
 	virtual ~Matrix()
 	{
 		gsl_matrix_free(m_matrix);
@@ -557,6 +581,10 @@ Matrix operator*( detail::BaseMatrixImpl const& A,
 
 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
 
+double trace (detail::BaseMatrixImpl const& A);
+
+double det (detail::BaseMatrixImpl const& A);
+
 } } // end namespaces
 
 #endif /*_NL_MATRIX_H_*/
diff --git a/src/2geom/numeric/symmetric-matrix-fs-operation.h b/src/2geom/numeric/symmetric-matrix-fs-operation.h
new file mode 100644
index 000000000..5222d2734
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs-operation.h
@@ -0,0 +1,107 @@
+/*
+ * SymmetricMatrix basic operation
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009  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 _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+#define _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+
+
+#include <2geom/numeric/symmetric-matrix-fs.h>
+#include <2geom/numeric/symmetric-matrix-fs-trace.h>
+
+
+
+
+namespace Geom { namespace NL {
+
+template <size_t N>
+inline
+SymmetricMatrix<N> adj(const ConstBaseSymmetricMatrix<N> & S)
+{
+    THROW_NOTIMPLEMENTED();
+    return SymmetricMatrix<N>();
+}
+
+template <>
+inline
+SymmetricMatrix<2> adj(const ConstBaseSymmetricMatrix<2> & S)
+{
+    SymmetricMatrix<2> result;
+    result.get<0,0>() = S.get<1,1>();
+    result.get<1,0>() = -S.get<1,0>();
+    result.get<1,1>() = S.get<0,0>();
+    return result;
+}
+
+template <>
+inline
+SymmetricMatrix<3> adj(const ConstBaseSymmetricMatrix<3> & S)
+{
+    SymmetricMatrix<3> result;
+
+    result.get<0,0>() = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+    result.get<1,0>() = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+    result.get<1,1>() = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+    result.get<2,0>() = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+    result.get<2,1>() = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+    result.get<2,2>() = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    return result;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> inverse(const ConstBaseSymmetricMatrix<N> & S)
+{
+    SymmetricMatrix<N> result = adj(S);
+    double d = det(S);
+    assert (d != 0);
+    result.scale (1/d);
+    return result;
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+#endif // _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+
+
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/symmetric-matrix-fs-trace.h b/src/2geom/numeric/symmetric-matrix-fs-trace.h
new file mode 100644
index 000000000..099c834a8
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs-trace.h
@@ -0,0 +1,433 @@
+/*
+ * SymmetricMatrix trace
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009  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 _NL_TRACE_H_
+#define _NL_TRACE_H_
+
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/symmetric-matrix-fs.h>
+
+
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+/*
+ *  helper routines
+ */
+
+inline
+int sgn_prod (int x, int y)
+{
+    if (x == 0 || y == 0)  return 0;
+    if (x == y) return 1;
+    return -1;
+}
+
+inline
+bool abs_less (double x, double y)
+{
+    return (std::fabs(x) < std::fabs(y));
+}
+
+
+/*
+ * trace K-th of symmetric matrix S of order N
+ */
+template <size_t K, size_t N>
+struct trace
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        THROW_NOTIMPLEMENTED();
+        return K;
+    }
+
+};
+
+template <size_t N>
+struct trace<1,N>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        double t = 0;
+        for (size_t i = 0; i < N; ++i)
+        {
+            t += S(i,i);
+        }
+        return t;
+    }
+};
+
+template <size_t N>
+struct trace<N,N>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        Matrix M(S);
+        return det(M);
+    }
+};
+
+/*
+ *  trace for symmetric matrix of order 2
+ */
+template <>
+struct trace<1,2>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        return (S.get<0,0>() + S.get<1,1>());
+    }
+};
+
+template <>
+struct trace<2,2>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>());
+    }
+};
+
+
+/*
+ *  trace for symmetric matrix of order 3
+ */
+template <>
+struct trace<1,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>());
+    }
+};
+
+template <>
+struct trace<2,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+        double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+        double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+        return (a00 + a11 + a22);
+    }
+};
+
+template <>
+struct trace<3,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+        d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>());
+        d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+        d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+        d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+        return d;
+    }
+};
+
+
+/*
+ *  sign of trace K-th
+ */
+template <size_t K, size_t N>
+struct trace_sgn
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        double d = trace<K, N>::evaluate(S);
+        return sgn(d);
+    }
+};
+
+
+/*
+ *  sign of trace for symmetric matrix of order 2
+ */
+template <>
+struct trace_sgn<2,2>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        double m00 = S.get<0,0>();
+        double m10 = S.get<1,0>();
+        double m11 = S.get<1,1>();
+
+        int sm00 = sgn (m00);
+        int sm10 = sgn (m10);
+        int sm11 = sgn (m11);
+
+        if (sm10 == 0)
+        {
+            return sgn_prod (sm00, sm11);
+        }
+        else
+        {
+            int sm00m11 = sgn_prod (sm00, sm11);
+            if (sm00m11 == 1)
+            {
+                int e00, e10, e11;
+                double f00 = std::frexp (m00, &e00);
+                double f10 = std::frexp (m10, &e10);
+                double f11 = std::frexp (m11, &e11);
+
+                int e0011 = e00 + e11;
+                int e1010 = e10 << 1;
+                int ed = e0011 - e1010;
+
+                if (ed > 1)
+                {
+                    return 1;
+                }
+                else if (ed < -1)
+                {
+                    return -1;
+                }
+                else
+                {
+                    double d = std::ldexp (f00 * f11, ed) - f10 * f10;
+                    //std::cout << "trace_sgn<2,2>: det = " << d << std::endl;
+                    double eps = std::ldexp (1, -50);
+                    if (std::fabs(d) < eps)  return 0;
+                    return sgn (d);
+                }
+            }
+            return -1;
+        }
+    }
+};
+
+
+/*
+ *  sign of trace for symmetric matrix of order 3
+ */
+template <>
+struct trace_sgn<2,3>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double eps = std::ldexp (1, -50);
+        double t[6];
+
+        t[0] = S.get<1,1>() * S.get<2,2>();
+        t[1] = - S.get<1,2>() * S.get<2,1>();
+        t[2] = S.get<0,0>() * S.get<2,2>();
+        t[3] = - S.get<0,2>() * S.get<2,0>();
+        t[4] = S.get<0,0>() * S.get<1,1>();
+        t[5] = - S.get<0,1>() * S.get<1,0>();
+
+
+        double* maxp = std::max_element (t, t+6, abs_less);
+        int em;
+        std::frexp(*maxp, &em);
+        double d = 0;
+        for (size_t i = 0; i < 6; ++i)
+        {
+            d += t[i];
+        }
+        double r = std::fabs (std::ldexp (d, -em));  // relative error
+        //std::cout << "trace_sgn<2,3>: d = " << d << std::endl;
+        //std::cout << "trace_sgn<2,3>: r = " << r << std::endl;
+        if (r < eps) return 0;
+        if (d > 0) return 1;
+        return -1;
+    }
+};
+
+template <>
+struct trace_sgn<3,3>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+
+        double eps = std::ldexp (1, -48);
+        double t[5];
+
+        t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+        t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>();
+        t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+        t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+        t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+
+        double* maxp = std::max_element (t, t+5, abs_less);
+        int em;
+        std::frexp(*maxp, &em);
+        double d = 0;
+        for (size_t i = 0; i < 5; ++i)
+        {
+            d += t[i];
+        }
+        //std::cout << "trace_sgn<3,3>: d = " << d << std::endl;
+        double r = std::fabs (std::ldexp (d, -em));  // relative error
+        //std::cout << "trace_sgn<3,3>: r = " << r << std::endl;
+
+        if (r < eps)  return 0;
+        if (d > 0) return 1;
+        return -1;
+    }
+}; // end struct trace_sgn<3,3>
+
+} // end namespace detail
+
+
+template <size_t K, size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double det (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<N, N>::evaluate(_matrix);
+}
+
+
+template <size_t K, size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<N, N>::evaluate(_matrix);
+}
+
+/*
+template <size_t N>
+inline
+size_t rank (const ConstBaseSymmetricMatrix<N> & S)
+{
+    THROW_NOTIMPLEMENTED();
+    return 0;
+}
+
+template <>
+inline
+size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S)
+{
+    if (S.is_zero())  return 0;
+    double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    if (d != 0)  return 2;
+    return 1;
+}
+
+template <>
+inline
+size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S)
+{
+    if (S.is_zero())  return 0;
+
+    double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+    double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+    double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>();
+
+    if (d != 0)  return 3;
+
+    if (a20 != 0 || a21 != 0 || a22 != 0)  return 2;
+
+    double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+    if (a00 != 0)  return 2;
+
+    double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+    if (a10 != 0)  return 2;
+
+    double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+    if (a11 != 0)  return 2;
+
+    return 1;
+}
+*/
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif // _NL_TRACE_H_
+
+
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/symmetric-matrix-fs.h b/src/2geom/numeric/symmetric-matrix-fs.h
new file mode 100644
index 000000000..c1de27afd
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs.h
@@ -0,0 +1,730 @@
+/*
+ * SymmetricMatrix, ConstSymmetricMatrixView, SymmetricMatrixView template
+ * classes implement fixed size symmetric matrix; "views" mimic the semantic
+ * of C++ references: any operation performed on a "view" is actually performed
+ * on the "viewed object"
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009  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 _NL_SYMMETRIC_MATRIX_FS_H_
+#define _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+#include <2geom/utils.h>
+#include <2geom/exception.h>
+
+#include <boost/static_assert.hpp>
+
+#include <cassert>
+#include <utility>     // for std::pair
+#include <algorithm>   // for std::swap, std::copy
+#include <sstream>
+#include <string>
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+template <size_t I, size_t J, bool B = (I < J)>
+struct index
+{
+    static const size_t K = index<J, I>::K;
+};
+
+template <size_t I, size_t J>
+struct index<I, J, false>
+{
+    static const size_t K = (((I+1) * I) >> 1) + J;
+};
+
+} // end namespace detail
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix;
+
+template <size_t N>
+class BaseSymmetricMatrix;
+
+template <size_t N>
+class SymmetricMatrix;
+
+template <size_t N>
+class ConstSymmetricMatrixView;
+
+template <size_t N>
+class SymmetricMatrixView;
+
+
+
+// declaration needed for friend clause
+template <size_t N>
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+                 ConstBaseSymmetricMatrix<N> const& _smatrix2);
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix
+{
+  public:
+    const static size_t DIM = N;
+    const static size_t DATA_SIZE = ((DIM+1) * DIM) / 2;
+
+  public:
+
+    ConstBaseSymmetricMatrix (VectorView const& _data)
+        : m_data(_data)
+    {
+    }
+
+    double operator() (size_t i, size_t j) const
+    {
+        return m_data[get_index(i,j)];
+    }
+
+    template <size_t I, size_t J>
+    double get() const
+    {
+        BOOST_STATIC_ASSERT ((I < N && J < N));
+        return m_data[detail::index<I, J>::K];
+    }
+
+
+    size_t rows() const
+    {
+        return DIM;
+    }
+
+    size_t columns() const
+    {
+        return DIM;
+    }
+
+    bool is_zero() const
+    {
+        return m_data.is_zero();
+    }
+
+    bool is_positive() const
+    {
+        return m_data.is_positive();
+    }
+
+    bool is_negative() const
+    {
+        return m_data.is_negative();
+    }
+
+    bool is_non_negative() const
+    {
+        return m_data.is_non_negative();
+    }
+
+    double min() const
+    {
+        return m_data.min();
+    }
+
+    double max() const
+    {
+        return m_data.max();
+    }
+
+    std::pair<size_t, size_t>
+    min_index() const
+    {
+        std::pair<size_t, size_t> indices(0,0);
+        double min_value = m_data[0];
+        for (size_t i = 1; i < DIM; ++i)
+        {
+            for (size_t j = 0; j <= i; ++j)
+            {
+                if (min_value > (*this)(i,j))
+                {
+                    min_value = (*this)(i,j);
+                    indices.first = i;
+                    indices.second = j;
+                }
+            }
+        }
+        return indices;
+    }
+
+    std::pair<size_t, size_t>
+    max_index() const
+    {
+        std::pair<size_t, size_t> indices(0,0);
+        double max_value = m_data[0];
+        for (size_t i = 1; i < DIM; ++i)
+        {
+            for (size_t j = 0; j <= i; ++j)
+            {
+                if (max_value < (*this)(i,j))
+                {
+                    max_value = (*this)(i,j);
+                    indices.first = i;
+                    indices.second = j;
+                }
+            }
+        }
+        return indices;
+    }
+
+    size_t min_on_row_index (size_t i) const
+    {
+        size_t idx = 0;
+        double min_value = (*this)(i,0);
+        for (size_t j = 1; j < DIM; ++j)
+        {
+            if (min_value > (*this)(i,j))
+            {
+                min_value = (*this)(i,j);
+                idx = j;
+            }
+        }
+        return idx;
+    }
+
+    size_t max_on_row_index (size_t i) const
+    {
+        size_t idx = 0;
+        double max_value = (*this)(i,0);
+        for (size_t j = 1; j < DIM; ++j)
+        {
+            if (max_value < (*this)(i,j))
+            {
+                max_value = (*this)(i,j);
+                idx = j;
+            }
+        }
+        return idx;
+    }
+
+    size_t min_on_column_index (size_t j) const
+    {
+        return min_on_row_index(j);
+    }
+
+    size_t max_on_column_index (size_t j) const
+    {
+        return max_on_row_index(j);
+    }
+
+    size_t min_on_diag_index () const
+    {
+        size_t idx = 0;
+        double min_value = (*this)(0,0);
+        for (size_t i = 1; i < DIM; ++i)
+        {
+            if (min_value > (*this)(i,i))
+            {
+                min_value = (*this)(i,i);
+                idx = i;
+            }
+        }
+        return idx;
+    }
+
+    size_t max_on_diag_index () const
+    {
+        size_t idx = 0;
+        double max_value = (*this)(0,0);
+        for (size_t i = 1; i < DIM; ++i)
+        {
+            if (max_value < (*this)(i,i))
+            {
+                max_value = (*this)(i,i);
+                idx = i;
+            }
+        }
+        return idx;
+    }
+
+    std::string str() const;
+
+    ConstSymmetricMatrixView<N-1> main_minor_const_view() const;
+
+    SymmetricMatrix<N> operator- () const;
+
+    Vector operator* (ConstVectorView _vector) const
+    {
+        assert (_vector.size() == DIM);
+        Vector result(DIM, 0.0);
+
+        for (size_t i = 0; i < DIM; ++i)
+        {
+            for (size_t j = 0; j < DIM; ++j)
+            {
+                result[i] += (*this)(i,j) * _vector[j];
+            }
+        }
+        return result;
+    }
+
+  protected:
+    static size_t get_index (size_t i, size_t j)
+    {
+        if (i < j) return get_index (j, i);
+        size_t k = (i+1) * i;
+        k >>= 1;
+        k += j;
+        return k;
+    }
+
+  protected:
+    ConstVectorView get_data() const
+    {
+      return m_data;
+    }
+
+    friend
+    bool operator==<N> (ConstBaseSymmetricMatrix const& _smatrix1,
+                        ConstBaseSymmetricMatrix const& _smatrix2);
+
+  protected:
+    VectorView m_data;
+
+}; //end ConstBaseSymmetricMatrix
+
+
+template <size_t N>
+class BaseSymmetricMatrix : public ConstBaseSymmetricMatrix<N>
+{
+  public:
+    typedef ConstBaseSymmetricMatrix<N> base_type;
+
+
+  public:
+
+    BaseSymmetricMatrix (VectorView const& _data)
+        : base_type(_data)
+    {
+    }
+
+    using base_type::operator();
+
+    double& operator() (size_t i, size_t j)
+    {
+        return m_data[base_type::get_index(i,j)];
+    }
+
+    template <size_t I, size_t J>
+    double& get()
+    {
+        BOOST_STATIC_ASSERT ((I < N && J < N));
+        return m_data[detail::index<I, J>::K];
+    }
+
+    void set_all (double x)
+    {
+        m_data.set_all(x);
+    }
+
+    SymmetricMatrixView<N-1> main_minor_view();
+
+    BaseSymmetricMatrix& transpose() const
+    {
+        return (*this);
+    }
+
+    BaseSymmetricMatrix& translate (double c)
+    {
+        m_data.translate(c);
+        return (*this);
+    }
+
+    BaseSymmetricMatrix& scale (double c)
+    {
+        m_data.scale(c);
+        return (*this);
+    }
+
+    BaseSymmetricMatrix& operator+= (base_type const& _smatrix)
+    {
+        m_data += (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+        return (*this);
+    }
+
+    BaseSymmetricMatrix& operator-= (base_type const& _smatrix)
+    {
+        m_data -= (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+        return (*this);
+    }
+
+    using base_type::DIM;
+    using base_type::DATA_SIZE;
+    using base_type::m_data;
+    using base_type::operator-;
+    using base_type::operator*;
+
+}; //end BaseSymmetricMatrix
+
+
+template <size_t N>
+class SymmetricMatrix : public BaseSymmetricMatrix<N>
+{
+  public:
+    typedef BaseSymmetricMatrix<N> base_type;
+    typedef typename  base_type::base_type base_base_type;
+
+    using base_type::DIM;
+    using base_type::DATA_SIZE;
+    using base_type::m_data;
+
+  public:
+    SymmetricMatrix ()
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+    }
+
+    explicit
+    SymmetricMatrix (ConstVectorView _data)
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+        assert (_data.size() == DATA_SIZE);
+        m_data = _data;
+    }
+
+    explicit
+    SymmetricMatrix (const double _data[DATA_SIZE])
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+        std::copy (_data, _data + DATA_SIZE, m_adata);
+    }
+
+    SymmetricMatrix (SymmetricMatrix const& _smatrix)
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+        m_data = _smatrix.m_data;
+    }
+
+    explicit
+    SymmetricMatrix (base_base_type const& _smatrix)
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+        m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+    }
+
+    explicit
+    SymmetricMatrix (ConstMatrixView const& _matrix)
+        : base_type (VectorView(m_adata, DATA_SIZE))
+    {
+        assert (_matrix.rows() == N && _matrix.columns() == N);
+        for (size_t i = 0; i < N; ++i)
+            for (size_t j = 0; j <= i ; ++j)
+                (*this)(i,j) = _matrix(i,j);
+    }
+
+    SymmetricMatrix& operator= (SymmetricMatrix const& _smatrix)
+    {
+        m_data = _smatrix.m_data;
+        return (*this);
+    }
+
+    SymmetricMatrix& operator= (base_base_type const& _smatrix)
+    {
+
+        m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+        return (*this);
+    }
+
+    SymmetricMatrix& operator= (ConstMatrixView const& _matrix)
+    {
+        assert (_matrix.rows() == N && _matrix.columns() == N);
+        for (size_t i = 0; i < N; ++i)
+            for (size_t j = 0; j <= i ; ++j)
+                (*this)(i,j) = _matrix(i,j);
+
+        return (*this);
+    }
+
+    // needed for accessing m_adata
+    friend class ConstSymmetricMatrixView<DIM>;
+    friend class SymmetricMatrixView<DIM>;
+  private:
+      double m_adata[DATA_SIZE];
+}; //end SymmetricMatrix
+
+
+template <size_t N>
+class ConstSymmetricMatrixView : public ConstBaseSymmetricMatrix<N>
+{
+  public:
+    typedef ConstBaseSymmetricMatrix<N> base_type;
+
+    using base_type::DIM;
+    using base_type::DATA_SIZE;
+    using base_type::m_data;
+
+
+  public:
+
+    explicit
+    ConstSymmetricMatrixView (ConstVectorView _data)
+        : base_type (const_vector_view_cast(_data))
+    {
+        assert (_data.size() == DATA_SIZE);
+    }
+
+    explicit
+    ConstSymmetricMatrixView (const double _data[DATA_SIZE])
+        : base_type (const_vector_view_cast (ConstVectorView (_data, DATA_SIZE)))
+    {
+    }
+
+    ConstSymmetricMatrixView (const ConstSymmetricMatrixView & _smatrix)
+        : base_type (_smatrix.m_data)
+    {
+    }
+
+    ConstSymmetricMatrixView (const base_type & _smatrix)
+        : base_type (static_cast<const ConstSymmetricMatrixView &>(_smatrix).m_data)
+    {
+    }
+
+}; //end SymmetricMatrix
+
+
+// declaration needed for friend clause
+template <size_t N>
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2);
+
+
+template <size_t N>
+class SymmetricMatrixView : public BaseSymmetricMatrix<N>
+{
+  public:
+    typedef BaseSymmetricMatrix<N> base_type;
+    typedef typename  base_type::base_type base_base_type;
+
+    using base_type::DIM;
+    using base_type::DATA_SIZE;
+    using base_type::m_data;
+
+  public:
+
+    explicit
+    SymmetricMatrixView (VectorView _data)
+        : base_type (_data)
+    {
+        assert (_data.size() == DATA_SIZE);
+    }
+
+    explicit
+    SymmetricMatrixView (double _data[DATA_SIZE])
+        : base_type (VectorView (_data, DATA_SIZE))
+    {
+    }
+
+    SymmetricMatrixView (const SymmetricMatrixView & _smatrix)
+        : base_type (_smatrix.m_data)
+    {
+    }
+
+    SymmetricMatrixView (SymmetricMatrix<DIM> & _smatrix)
+        : base_type (VectorView (_smatrix.m_adata, DATA_SIZE))
+    {
+    }
+
+    SymmetricMatrixView& operator= (const SymmetricMatrixView & _smatrix)
+    {
+        m_data = _smatrix.m_data;
+        return (*this);
+    }
+
+    SymmetricMatrixView& operator= (const base_base_type & _smatrix)
+    {
+        m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+        return (*this);
+    }
+
+    friend
+    void swap_view<N>(SymmetricMatrixView & m1, SymmetricMatrixView & m2);
+
+}; //end SymmetricMatrix
+
+
+
+
+/*
+ * class ConstBaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+std::string ConstBaseSymmetricMatrix<N>::str() const
+{
+    std::ostringstream oss;
+    oss << (*this);
+    return oss.str();
+}
+
+template <size_t N>
+inline
+ConstSymmetricMatrixView<N-1>
+ConstBaseSymmetricMatrix<N>::main_minor_const_view() const
+{
+    ConstVectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+    ConstSymmetricMatrixView<N-1> mm(data);
+    return mm;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> ConstBaseSymmetricMatrix<N>::operator- () const
+{
+    SymmetricMatrix<N> result;
+    for (size_t i = 0; i < DATA_SIZE; ++i)
+    {
+        result.m_data[i] = -m_data[i];
+    }
+    return result;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix friend free functions
+ */
+
+template <size_t N>
+inline
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+                 ConstBaseSymmetricMatrix<N> const& _smatrix2)
+{
+    return (_smatrix1.m_data == _smatrix2.m_data);
+}
+
+/*
+ * class ConstBaseSymmetricMatrix related free functions
+ */
+
+template< size_t N, class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os,
+            const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    os << "[[" << _matrix(0,0);
+    for (size_t j = 1; j < N; ++j)
+    {
+        os << ", " << _matrix(0,j);
+    }
+    os << "]";
+    for (size_t i = 1; i < N; ++i)
+    {
+        os << "\n [" << _matrix(i,0);
+        for (size_t j = 1; j < N; ++j)
+        {
+            os << ", " << _matrix(i,j);
+        }
+        os << "]";
+    }
+    os << "]";
+    return os;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix specialized methods
+ */
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<2>::get_index (size_t i, size_t j)
+{
+    return (i+j);
+}
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<3>::get_index (size_t i, size_t j)
+{
+    size_t k = i + j;
+    if (i == 2 || j == 2)  ++k;
+    return k;
+}
+
+
+/*
+ * class BaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+SymmetricMatrixView<N-1> BaseSymmetricMatrix<N>::main_minor_view()
+{
+    VectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+    SymmetricMatrixView<N-1> mm(data);
+    return mm;
+}
+
+
+/*
+ * class SymmetricMatrixView friend free functions
+ */
+
+template <size_t N>
+inline
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2)
+{
+    swap_view(m1.m_data, m2.m_data);
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif  // _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
index 46701645a..6ab898f29 100644
--- a/src/2geom/numeric/vector.h
+++ b/src/2geom/numeric/vector.h
@@ -555,6 +555,24 @@ void swap_view(VectorView & v1, VectorView & v2)
 	std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
 }
 
+inline
+const VectorView & const_vector_view_cast (const ConstVectorView & view)
+{
+    const detail::BaseVectorImpl & bvi
+        = static_cast<const detail::BaseVectorImpl &>(view);
+    const VectorView & vv = reinterpret_cast<const VectorView &>(bvi);
+    return vv;
+}
+
+inline
+VectorView & const_vector_view_cast (ConstVectorView & view)
+{
+    detail::BaseVectorImpl & bvi
+        = static_cast<detail::BaseVectorImpl &>(view);
+    VectorView & vv = reinterpret_cast<VectorView &>(bvi);
+    return vv;
+}
+
 
 } } // end namespaces