Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

9 files changed, 2921 insertions, 0 deletions

diff --git a/include/2geom/orphan-code/arc-length.h b/include/2geom/orphan-code/arc-length.h
new file mode 100644
index 0000000..8029f04
--- /dev/null
+++ b/include/2geom/orphan-code/arc-length.h
@@ -0,0 +1,58 @@
+/**
+ * \file arc-length.h
+ * \brief Arc length computations for paths
+ *//*
+ * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *
+ * 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 __2GEOM_ARC_LENGTH_H
+#define __2GEOM_ARC_LENGTH_H
+
+#include <2geom/path.h>
+
+/* Routines in this group return a path that looks the same, but
+ * include extra knots for certain points of interest. */
+
+/* find_vector_extreme_points
+ * extreme points . dir.
+ */
+
+double arc_length_subdividing(Geom::Path const & p, double tol);
+double arc_length_integrating(Geom::Path const & p, double tol);
+
+#endif
+
+/*
+  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/include/2geom/orphan-code/chebyshev.h b/include/2geom/orphan-code/chebyshev.h
new file mode 100644
index 0000000..f729e1f
--- /dev/null
+++ b/include/2geom/orphan-code/chebyshev.h
@@ -0,0 +1,30 @@
+#ifndef _CHEBYSHEV
+#define _CHEBYSHEV
+
+#include <2geom/sbasis.h>
+#include <2geom/interval.h>
+
+/*** Conversion between Chebyshev approximation and SBasis.
+ * 
+ */
+
+namespace Geom{
+
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p=0);
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*), int order, Interval in, void* p=0);
+SBasis chebyshev(unsigned n);
+
+};
+
+/*
+  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 :
+
+#endif
diff --git a/include/2geom/orphan-code/intersection-by-smashing.h b/include/2geom/orphan-code/intersection-by-smashing.h
new file mode 100644
index 0000000..996ec99
--- /dev/null
+++ b/include/2geom/orphan-code/intersection-by-smashing.h
@@ -0,0 +1,78 @@
+/*
+ * Diffeomorphism-based intersector: given two curves
+ *  M(t)=(x(t),y(t)) and N(u)=(X(u),Y(u))
+ * and supposing M is a graph over the x-axis, we compute y(x) and the
+ * transformation:
+ *  X <- X
+ *  Y <- Y - y(X)
+ * smashes M on the x axis. The intersections are then given by the roots of:
+ *  Y(u) - y(X(u)) = 0
+ *//*
+ * Authors:
+ * 		J.-F. Barraud    <jfbarraud at gmail.com>
+ * Copyright 2010  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 SEEN_LIB2GEOM_INTERSECTION_BY_SMASHING_H
+#define SEEN_LIB2GEOM_INTERSECTION_BY_SMASHING_H
+
+#include <2geom/d2.h>
+#include <2geom/interval.h>
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-geometric.h>
+#include <cstdlib>
+#include <vector>
+#include <algorithm>
+
+
+namespace Geom{
+
+struct SmashIntersection {
+	Rect times;
+	Rect bbox;
+};
+
+std::vector<SmashIntersection> smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol);
+std::vector<SmashIntersection> monotonic_smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol);
+//std::vector<Intersection> monotonic_smash_intersect( Curve const &a, double a_from, double a_to,
+//			                                         Curve const &b, double b_from, double b_to, double tol);
+
+std::vector<Interval> monotonicSplit(D2<SBasis> const &p);
+
+} // end namespace Geom
+
+#endif // !SEEN_LIB2GEOM_INTERSECTION_BY_SMASHING_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/include/2geom/orphan-code/linear-of.h b/include/2geom/orphan-code/linear-of.h
new file mode 100644
index 0000000..9ba1fb2
--- /dev/null
+++ b/include/2geom/orphan-code/linear-of.h
@@ -0,0 +1,269 @@
+/**

+ * \file

+ * \brief  Linear fragment function class

+ *

+ *  Authors:

+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>

+ *   Michael Sloan <mgsloan@gmail.com>

+ *

+ * Copyright (C) 2006-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 SEEN_LINEAR_OF_H

+#define SEEN_LINEAR_OF_H

+#include <2geom/interval.h>

+#include <2geom/math-utils.h>

+

+namespace Geom{

+

+template <typename T>

+inline T lerp(double t, T a, T b) { return a*(1-t) + b*t; }

+

+template <typename T>

+class SBasisOf;

+

+template <typename T>

+class HatOf{

+public:

+    HatOf () {}

+    HatOf(T d) :d(d) {}

+    operator T() const { return d; }

+    T d;

+};

+

+template <typename T>

+class TriOf{

+public:

+    TriOf () {}

+    TriOf(double d) :d(d) {}

+    operator T() const { return d; }

+    T d;

+};

+

+

+//--------------------------------------------------------------------------

+#ifdef USE_SBASIS_OF

+template <typename T>

+class LinearOf;

+typedef Geom::LinearOf<double> Linear;

+#endif

+//--------------------------------------------------------------------------

+

+template <typename T>

+class LinearOf{

+public:

+    T a[2];

+    LinearOf() {}

+    LinearOf(T aa, T b) {a[0] = aa; a[1] = b;}

+    //LinearOf(double aa, double b) {a[0] = T(aa); a[1] = T(b);}

+    LinearOf(HatOf<T> h, TriOf<T> t) {

+        a[0] = T(h) - T(t)/2; 

+        a[1] = T(h) + T(t)/2;

+    }

+

+    LinearOf(HatOf<T> h) {

+        a[0] = T(h); 

+        a[1] = T(h);

+    }

+

+    unsigned input_dim(){return T::input_dim() + 1;}

+

+    T operator[](const int i) const {

+        assert(i >= 0);

+        assert(i < 2);

+        return a[i];

+    }

+    T& operator[](const int i) {

+        assert(i >= 0);

+        assert(i < 2);

+        return a[i];

+    }

+

+    //IMPL: FragmentConcept

+    typedef T output_type;

+    inline bool isZero() const { return a[0].isZero() && a[1].isZero(); }

+    inline bool isConstant() const { return a[0] == a[1]; }

+    inline bool isFinite() const { return std::isfinite(a[0]) && std::isfinite(a[1]); }

+

+    inline T at0() const { return a[0]; }

+    inline T at1() const { return a[1]; }

+

+    inline T valueAt(double t) const { return lerp(t, a[0], a[1]); }

+    inline T operator()(double t) const { return valueAt(t); }

+

+    //defined in sbasis.h

+    inline SBasisOf<T> toSBasis() const;

+

+//This is specific for T=double!!

+    inline OptInterval bounds_exact() const { return Interval(a[0], a[1]); }

+    inline OptInterval bounds_fast() const { return bounds_exact(); }

+    inline OptInterval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }

+

+    operator TriOf<T>() const {

+        return a[1] - a[0];

+    }

+    operator HatOf<T>() const {

+        return (a[1] + a[0])/2;

+    }

+};

+

+template <>

+unsigned LinearOf<double>::input_dim(){return 1;}

+template <>

+inline OptInterval LinearOf<double>::bounds_exact() const { return Interval(a[0], a[1]); }

+template <>

+inline OptInterval LinearOf<double>::bounds_fast() const { return bounds_exact(); }

+template <>

+inline OptInterval LinearOf<double>::bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }

+template <>

+inline bool LinearOf<double>::isZero() const { return a[0]==0 && a[1]==0; }

+

+template <typename T>

+inline LinearOf<T> reverse(LinearOf<T> const &a) { return LinearOf<T>(a[1], a[0]); }

+

+//IMPL: AddableConcept

+template <typename T>

+inline LinearOf<T> operator+(LinearOf<T> const & a, LinearOf<T> const & b) {

+    return LinearOf<T>(a[0] + b[0], a[1] + b[1]);

+}

+template <typename T>

+inline LinearOf<T> operator-(LinearOf<T> const & a, LinearOf<T> const & b) {

+    return LinearOf<T>(a[0] - b[0], a[1] - b[1]);

+}

+template <typename T>

+inline LinearOf<T>& operator+=(LinearOf<T> & a, LinearOf<T> const & b) {

+    a[0] += b[0]; a[1] += b[1];

+    return a;

+}

+template <typename T>

+inline LinearOf<T>& operator-=(LinearOf<T> & a, LinearOf<T> const & b) {

+    a[0] -= b[0]; a[1] -= b[1];

+    return a;

+}

+//IMPL: OffsetableConcept

+template <typename T>

+inline LinearOf<T> operator+(LinearOf<T> const & a, double b) {

+    return LinearOf<T>(a[0] + b, a[1] + b);

+}

+template <typename T>

+inline LinearOf<T> operator-(LinearOf<T> const & a, double b) {

+    return LinearOf<T>(a[0] - b, a[1] - b);

+}

+template <typename T>

+inline LinearOf<T>& operator+=(LinearOf<T> & a, double b) {

+    a[0] += b; a[1] += b;

+    return a;

+}

+template <typename T>

+inline LinearOf<T>& operator-=(LinearOf<T> & a, double b) {

+    a[0] -= b; a[1] -= b;

+    return a;

+}

+/*

+//We can in fact offset in coeff ring T...

+template <typename T>

+inline LinearOf<T> operator+(LinearOf<T> const & a, T b) {

+    return LinearOf<T>(a[0] + b, a[1] + b);

+}

+template <typename T>

+inline LinearOf<T> operator-(LinearOf<T> const & a, T b) {

+    return LinearOf<T>(a[0] - b, a[1] - b);

+}

+template <typename T>

+inline LinearOf<T>& operator+=(LinearOf<T> & a, T b) {

+    a[0] += b; a[1] += b;

+    return a;

+}

+template <typename T>

+inline LinearOf<T>& operator-=(LinearOf<T> & a, T b) {

+    a[0] -= b; a[1] -= b;

+    return a;

+}

+*/

+

+//IMPL: boost::EqualityComparableConcept

+template <typename T>

+inline bool operator==(LinearOf<T> const & a, LinearOf<T> const & b) {

+    return a[0] == b[0] && a[1] == b[1];

+}

+template <typename T>

+inline bool operator!=(LinearOf<T> const & a, LinearOf<T> const & b) {

+    return a[0] != b[0] || a[1] != b[1];

+}

+//IMPL: ScalableConcept

+template <typename T>

+inline LinearOf<T> operator-(LinearOf<T> const &a) {

+    return LinearOf<T>(-a[0], -a[1]);

+}

+template <typename T>

+inline LinearOf<T> operator*(LinearOf<T> const & a, double b) {

+    return LinearOf<T>(a[0]*b, a[1]*b);

+}

+template <typename T>

+inline LinearOf<T> operator/(LinearOf<T> const & a, double b) {

+    return LinearOf<T>(a[0]/b, a[1]/b);

+}

+template <typename T>

+inline LinearOf<T> operator*=(LinearOf<T> & a, double b) {

+    a[0] *= b; a[1] *= b;

+    return a;

+}

+template <typename T>

+inline LinearOf<T> operator/=(LinearOf<T> & a, double b) {

+    a[0] /= b; a[1] /= b;

+    return a;

+}

+/*

+//We can in fact rescale in coeff ring T... (but not divide!)

+template <typename T>

+inline LinearOf<T> operator*(LinearOf<T> const & a, T b) {

+    return LinearOf<T>(a[0]*b, a[1]*b);

+}

+template <typename T>

+inline LinearOf<T> operator/(LinearOf<T> const & a, T b) {

+    return LinearOf<T>(a[0]/b, a[1]/b);

+}

+template <typename T>

+inline LinearOf<T> operator*=(LinearOf<T> & a, T b) {

+    a[0] *= b; a[1] *= b;

+    return a;

+}

+*/

+

+};

+

+#endif //SEEN_LINEAR_OF_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/include/2geom/orphan-code/linearN.h b/include/2geom/orphan-code/linearN.h
new file mode 100644
index 0000000..bb27c30
--- /dev/null
+++ b/include/2geom/orphan-code/linearN.h
@@ -0,0 +1,363 @@
+/**

+ * @file

+ * @brief  LinearN fragment function class

+ *//*

+ *  Authors:

+ *   JF Barraud <jf.barraud@gmail.com>

+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>

+ *   Michael Sloan <mgsloan@gmail.com>

+ *

+ * Copyright (C) 2006-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 SEEN_LINEARN_H

+#define SEEN_LINEARN_H

+#include <2geom/interval.h>

+#include <2geom/math-utils.h>

+#include <2geom/linear.h> //for conversion purpose ( + lerp() )

+

+#include <iostream>

+

+namespace Geom{

+

+//TODO: define this only once!! (see linear.h)

+inline double lerpppp(double t, double a, double b) { return a*(1-t) + b*t; }

+

+template<unsigned n>

+class SBasisN;

+

+template<unsigned n>

+class LinearN{

+public:

+    double a[1<<n];// 1<<n is 2^n

+    LinearN() {

+        for (unsigned i=0; i < (1<<n); i++){

+            a[i] = 0.;

+        }

+    }

+    LinearN(double aa[]) {

+        for (unsigned i=0; i < (1<<n); i++){

+            a[i] = aa[i];

+        }

+    }

+    LinearN(double c) {

+        for (unsigned i=0; i<(1<<n); i++){

+            a[i] = c;

+        }

+    }

+    LinearN(LinearN<n-1> const &aa, LinearN<n-1> const &b, unsigned var=0) {

+//         for (unsigned i=0; i<(1<<n-1); i++){

+//             a[i] = aa[i];

+//             a[i+(1<<(n-1))] = b[i];

+//         }

+        unsigned mask = (1<<var)-1;

+        for (unsigned i=0; i < (1<<(n-1)); i++){

+            unsigned low_i = i & mask, high_i = i & ~mask;

+            unsigned idx0 = (high_i<<1)|low_i;

+            unsigned idx1 = (high_i<<1)|(1<<var)|low_i;

+            a[idx0] = aa[i];

+            a[idx1] = b[i];

+        }

+

+    }

+    double operator[](const int i) const {

+        assert( i >= 0 );

+        assert( i < (1<<n) );

+        return a[i];

+    }

+    double& operator[](const int i) {

+        assert(i >= 0);

+        assert(i < (1<<n) );

+        return a[i];

+    }

+

+    //IMPL: FragmentConcept

+    typedef double output_type;

+    unsigned input_dim() const {return n;}

+    inline bool isZero() const { 

+        for (unsigned i=0; i < (1<<n); i++){

+            if (a[i] != 0) return false;

+        }

+        return true; }

+    inline bool isConstant() const { 

+        for (unsigned i=1; i < (1<<n); i++){

+            if (a[i] != a[0]) return false;

+        }

+        return true; }

+    inline bool isConstant(unsigned var) const { 

+        unsigned mask = (1<<var)-1;

+        for (unsigned i=0; i < (1<<(n-1)); i++){

+            unsigned low_i = i & mask, high_i = i & ~mask;

+            unsigned idx0 = (high_i<<1)|low_i;

+            unsigned idx1 = (high_i<<1)|(1<<var)|low_i;

+            if (a[idx0] != a[idx1]) return false;

+        }

+        return true; 

+    }

+    inline bool isFinite() const {

+        for (unsigned i=0; i < (1<<n); i++){

+            if ( !std::isfinite(a[i]) ) return false;

+        }

+        return true; }

+    //value if k-th variable is set to 0.

+    inline LinearN<n-1> at0(unsigned k=0) const {

+        LinearN<n-1> res;

+        unsigned mask = (1<<k)-1;

+        for (unsigned i=0; i < (1<<(n-1)); i++){

+            unsigned low_i = i & mask, high_i = i & ~mask;

+            unsigned idx = (high_i<<1)|low_i;

+            res[i] = a[idx];

+        }

+        return res;

+    }

+    //value if k-th variable is set to 1.

+    inline LinearN<n-1> at1(unsigned k=0) const {

+        LinearN<n-1> res;

+        for (unsigned i=0; i < (1<<(n-1)); i++){

+            unsigned mask = (1<<k)-1;

+            unsigned low_i = i & mask, high_i = i & ~mask;

+            unsigned idx = (high_i<<1)|(1<<k)|low_i;

+            res[i] = a[idx];

+        }

+        return res;

+    }

+    inline double atCorner(unsigned k) const {

+        assert( k < (1<<n) );

+        return a[k];

+    }

+    inline double atCorner(double t[]) const { 

+        unsigned k=0;

+        for(unsigned i=0; i<n; i++){

+            if (t[i] == 1.) k = k | (1<<i);

+            else assert( t[i] == 0. );

+        }

+        return atCorner(k);

+    }

+    inline LinearN<n-1> partialEval(double t, unsigned var=0 ) const {

+        LinearN<n-1> res;

+        res = at0(var)*(1-t) + at1(var)*t;

+        return res;

+    }

+    

+    //fixed and flags are used for recursion.

+    inline double valueAt(double t[], unsigned fixed=0, unsigned flags=0 ) const {

+        if (fixed == n) {

+            return a[flags]; 

+        }else{

+            double a0 = valueAt(t, fixed+1, flags);

+            double a1 = valueAt(t, fixed+1, flags|(1<<fixed));

+            return lerpppp( t[fixed], a0, a1 );

+        }

+    }

+    inline double operator()(double t[]) const { return valueAt(t); }

+

+    //defined in sbasisN.h

+    inline SBasisN<n> toSBasisN() const;

+

+    inline OptInterval bounds_exact() const {

+        double min=a[0], max=a[0];

+        for (unsigned i=1; i < (1<<n); i++){

+            if (a[i] < min) min = a[i];

+            if (a[i] > max) max = a[i];

+        }

+        return Interval(min, max); 

+    }

+    inline OptInterval bounds_fast() const { return bounds_exact(); }

+    //inline OptInterval bounds_local(double u, double v) const { return Interval(valueAt(u), valueAt(v)); }

+};

+

+//LinearN<0> are doubles. Specialize them out.

+template<>

+class LinearN<0>{

+public:

+    double d;

+    LinearN () {}

+    LinearN(double d) :d(d) {}

+    operator double() const { return d; }

+    double operator[](const int i) const {assert (i==0); return d;}

+    double& operator[](const int i) {assert (i==0); return d;}

+    typedef double output_type;

+    unsigned input_dim() const {return 0;}

+    inline bool isZero() const { return d==0; }

+    inline bool isConstant() const { return true; }

+    inline bool isFinite() const { return std::isfinite(d); }

+};

+

+//LinearN<1> are usual Linear. Allow conversion.

+Linear toLinear(LinearN<1> f){

+    return Linear(f[0],f[1]);

+}

+

+

+

+//inline Linear reverse(Linear const &a) { return Linear(a[1], a[0]); }

+

+//IMPL: AddableConcept

+template<unsigned n>

+inline LinearN<n> operator+(LinearN<n> const & a, LinearN<n> const & b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] + b[i];

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n> operator-(LinearN<n> const & a, LinearN<n> const & b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] - b[i];

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n>& operator+=(LinearN<n> & a, LinearN<n> const & b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] += b[i];

+    }

+    return a;

+}

+template<unsigned n>

+inline LinearN<n>& operator-=(LinearN<n> & a, LinearN<n> const & b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] -= b[i];

+    }

+    return a;

+}

+//IMPL: OffsetableConcept

+template<unsigned n>

+inline LinearN<n> operator+(LinearN<n> const & a, double b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] + b;

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n> operator-(LinearN<n> const & a, double b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] - b;

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n>& operator+=(LinearN<n> & a, double b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] += b;

+    }

+    return a;

+}

+template<unsigned n>

+inline LinearN<n>& operator-=(LinearN<n> & a, double b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] -= b;

+    }

+    return a;

+}

+//IMPL: boost::EqualityComparableConcept

+template<unsigned n>

+inline bool operator==(LinearN<n> const & a, LinearN<n> const & b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        if (a[i] != b[i]) return false;

+    }

+    return true;

+}

+template<unsigned n>

+inline bool operator!=(LinearN<n> const & a, LinearN<n> const & b) {

+    return !(a==b);

+}

+//IMPL: ScalableConcept

+template<unsigned n>

+inline LinearN<n> operator-(LinearN<n> const &a) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = -a[i];

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n> operator*(LinearN<n> const & a, double b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] * b;

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n> operator/(LinearN<n> const & a, double b) {

+    LinearN<n> res;

+    for (unsigned i=0; i < (1<<n); i++){

+        res[i] = a[i] / b;

+    }

+    return res;

+}

+template<unsigned n>

+inline LinearN<n> operator*=(LinearN<n> & a, double b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] *= b;

+    }

+    return a;

+}

+template<unsigned n>

+inline LinearN<n> operator/=(LinearN<n> & a, double b) {

+    for (unsigned i=0; i < (1<<n); i++){

+        a[i] /= b;

+    }

+    return a;

+}

+

+template<unsigned n>

+void setToVariable(LinearN<n> &x, unsigned k){;

+    x = LinearN<n>(0.);

+    unsigned mask = 1<<k;

+    for (unsigned i=0; i < (1<<n); i++){

+        if ( i & mask ) x[i] = 1;

+    }

+}

+

+template<unsigned n>

+inline std::ostream &operator<< (std::ostream &out_file, const LinearN<n> &bo) {

+    out_file << "{";

+    for (unsigned i=0; i < (1<<n); i++){

+        out_file << bo[i]<<(i == (1<<n)-1 ? "}" : ",");

+    }

+    return out_file;

+}

+

+

+}

+#endif //SEEN_LINEAR_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/include/2geom/orphan-code/redblacktree.h b/include/2geom/orphan-code/redblacktree.h
new file mode 100644
index 0000000..9d79342
--- /dev/null
+++ b/include/2geom/orphan-code/redblacktree.h
@@ -0,0 +1,121 @@
+/**
+ * \file
+ * \brief 
+ * Implementation of Red-Black Tree as described in 
+ * Intorduction to Algorithms. Cormen et al. Mc Grow Hill. 1990. pp 263-280
+ * 
+ * The intention is to implement interval trees mentioned in the same book, after the red-black.
+ * Interval are heavily based on red-black trees (most operations are the same). So, we begin first 
+ * with implementing red-black!
+ *
+ * Authors:
+ *      ? <?@?.?>
+ * 
+ * Copyright 2009-2009 Evangelos Katsikaros
+ *
+ * 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 SEEN_LIB2GEOM_REDBLACKTREE_H
+#define SEEN_LIB2GEOM_REDBLACKTREE_H
+
+#include <vector>
+//#include <cassert>
+#include <limits>
+#include <cfloat>
+
+#include <2geom/d2.h>
+#include <2geom/interval.h>
+
+namespace Geom{
+
+class RedBlack{    
+public:
+    Interval interval;    // Key of the redblack tree will be the min of the interval
+    RedBlack *left, *right, *parent;
+    bool isRed;
+    Coord subtree_max;    // subtree_max = max( x->left->subtree_max, x->right->subtree_max, x->high )
+    int data;
+
+    RedBlack(): left(0), right(0), parent(0), isRed(false), subtree_max(0.0), data(0) {
+        Interval interval(0.0, 0.0);
+    }
+/*
+    RedBlack(Coord min, Coord max): left(0), right(0), parent(0), isRed(false), subtree_max(0.0), data(0) {
+        Interval interval( min, max );
+    }
+*/
+    inline Coord key(){ return interval.min(); };
+    inline Coord high(){ return interval.max(); };
+};
+
+
+class RedBlackTree{
+public:
+    RedBlack* root;
+
+    RedBlackTree(): root(0) {}
+
+    void insert(Rect const &r, int shape, int dimension);
+    void insert(Coord dimension_min, Coord dimension_max, int shape);
+
+    void erase(Rect const &r);
+    void erase(int shape);
+
+    RedBlack* search(Rect const &r, int dimension);
+    RedBlack* search(Interval i);
+    RedBlack* search(Coord a, Coord b);
+
+    void print_tree();
+private:
+    void inorder_tree_walk(RedBlack* x);
+    RedBlack* tree_minimum(RedBlack* x);
+    RedBlack* tree_successor(RedBlack* x);
+
+    void left_rotate(RedBlack* x);
+    void right_rotate(RedBlack* x);
+    void tree_insert(RedBlack* x);
+
+    void update_max(RedBlack* x);
+
+    RedBlack* erase(RedBlack* x); // TODO why rerutn pointer? to collect garbage ???
+    void erase_fixup(RedBlack* x);
+
+};
+
+} // end namespace Geom
+
+#endif // !SEEN_LIB2GEOM_REDBLACKTREE_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/include/2geom/orphan-code/rtree.h b/include/2geom/orphan-code/rtree.h
new file mode 100644
index 0000000..3ffae8e
--- /dev/null
+++ b/include/2geom/orphan-code/rtree.h
@@ -0,0 +1,241 @@
+/**
+ * \file
+ * \brief 
+ * Implementation of Red-Black Tree as described in 
+ * Intorduction to Algorithms. Cormen et al. Mc Grow Hill. 1990. pp 263-280
+ * 
+ * The intention is to implement interval trees mentioned in the same book, after the red-black.
+ * Interval are heavily based on red-black trees (most operations are the same). So, we begin first 
+ * with implementing red-black!
+ *
+ * Authors:
+ *      ? <?@?.?>
+ * 
+ * Copyright 2009-2009 Evangelos Katsikaros
+ *
+ * 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 SEEN_LIB2GEOM_RTREE_H
+#define SEEN_LIB2GEOM_RTREE_H
+
+#include <vector>
+#include <utility>
+
+
+#include <2geom/d2.h>
+#include <2geom/interval.h>
+
+namespace Geom{
+
+// used only in pick_next( )
+enum enum_add_to_group { 
+    ADD_TO_GROUP_A = 0, 
+    ADD_TO_GROUP_B
+};
+
+
+enum enum_split_strategy { 
+    QUADRATIC_SPIT = 0,
+    LINEAR_COST,
+    TOTAL_STRATEGIES // this one must be the last item
+};
+
+
+
+template <typename T>
+class pedantic_vector:public std::vector<T> {
+public:
+    pedantic_vector(size_t s=0) : std::vector<T>(s) {}
+    T& operator[](unsigned i) {
+        //assert(i >= 0);
+        assert(i < std::vector<T>::size());
+        return std::vector<T>::operator[](i);
+    }
+    T const& operator[](unsigned i) const {
+        //assert(i >= 0);
+        assert(i < std::vector<T>::size());
+        return std::vector<T>::operator[](i);
+    }
+/*
+    erase( std::vector<T>::iterator it ) {
+        //assert(i >= 0);
+        assert( it < std::vector<T>::size());
+        return std::vector<T>::erase(it);
+    }
+*/
+};
+
+class RTreeNode;
+
+class RTreeRecord_Leaf{
+public:
+    Rect bounding_box;
+    int data;
+
+    RTreeRecord_Leaf(): bounding_box(), data(0)
+    {}
+
+    RTreeRecord_Leaf(Rect bb, int d): bounding_box(bb), data(d)
+    {}
+};
+
+class RTreeRecord_NonLeaf{
+public:
+    Rect bounding_box;
+    RTreeNode* data;    
+
+    RTreeRecord_NonLeaf(): bounding_box(), data(0)
+    {}
+
+    RTreeRecord_NonLeaf(Rect bb, RTreeNode* d): bounding_box(bb), data(d)
+    {}
+};
+
+/*
+R-Tree has 2 kinds of nodes
+* Leaves which store:
+  - the actual data
+  - the bounding box of the data
+
+* Non-Leaves which store:
+  - a child node (data)
+  - the bounding box of the child node
+
+This causes some code duplication in rtree.cpp. There are 2 cases:
+- we care whether we touch a leaf/non-leaf node, since we write data in the node, so we want to 
+  write the correct thing (int or RTreeNode*)
+- we do NOT care whether we touch a leaf/non-leaf node, because we only read/write the bounding 
+  boxes which is the same in both cases.
+
+TODO:
+A better design would eliminate the duplication in the 2nd case, but we can't avoid the 1st probably.
+*/
+class RTreeNode{    
+public:
+    // first: bounding box
+    // second: "data" (leaf-node) or node (NON leaf-node)
+    //pedantic_vector< RTreeRecord_Leaf > children_leaves; // if this is empty, then node is leaf-node
+    //pedantic_vector< RTreeRecord_NonLeaf > children_nodes;  // if this is empty, then node is NON-leaf node
+
+    std::vector< RTreeRecord_Leaf > children_leaves; // if this is empty, then node is leaf-node
+    std::vector< RTreeRecord_NonLeaf > children_nodes;  // if this is empty, then node is NON-leaf node
+
+    RTreeNode(): children_leaves(0), children_nodes(0)
+    {}
+
+};
+
+
+class RTree{
+public:
+    RTreeNode* root;
+
+    // min/max records per node
+    unsigned min_records;
+    unsigned max_records; // allow +1 (used during insert)
+
+    enum_split_strategy split_strategy;
+
+
+    RTree( unsigned n, unsigned m, enum_split_strategy split_s ): 
+        root(0), min_records( n ), max_records( m ), split_strategy( split_s ),
+        tree_height(0)
+    {}
+
+    void insert( Rect const &r, unsigned shape);
+    void search( const Rect &search_area, std::vector< int >* result, const RTreeNode* subtree ) const;
+    //int erase( const RTreeRecord_Leaf & search );
+    int erase( const Rect &search_area, const int shape_to_delete );
+
+//  update
+
+    void print_tree(RTreeNode* subtree_root, int depth ) const;
+
+private:
+    unsigned tree_height; // 0 is the root level
+
+    void insert(  //Rect const &r, 
+                //int shape, 
+                const RTreeRecord_Leaf &leaf_record,
+                const bool &insert_high = false, 
+                const unsigned &stop_height = 0,
+                const RTreeRecord_NonLeaf &nonleaf_record = RTreeRecord_NonLeaf()
+                );
+    // I1
+    RTreeNode* choose_node( const Rect &r, const bool &insert_high = false, const unsigned &stop_height=0 ) const;
+    double find_waste_area( const Rect &a, const Rect &b ) const;
+    double find_enlargement( const Rect &a, const Rect &b ) const;
+
+    // I2
+    std::pair<RTreeNode*, RTreeNode*> split_node( RTreeNode *s );
+        // QUADRATIC_SPIT
+    std::pair<RTreeNode*, RTreeNode*> quadratic_split( RTreeNode* s );
+    std::pair<unsigned, unsigned> pick_seeds( RTreeNode* s ) const;
+    std::pair<unsigned, enum_add_to_group>  pick_next( RTreeNode* group_a, RTreeNode* group_b, RTreeNode* s, std::vector<bool> &assigned_v );
+        // others...
+
+    // I3
+    bool adjust_tree(   RTreeNode* position, 
+                        std::pair<RTreeNode*, RTreeNode*>  &node_division, 
+                        bool split_performed
+                        );
+    std::pair< RTreeNode*, bool > find_parent( RTreeNode* subtree_root, Rect search_area, RTreeNode* wanted ) const;
+
+    void recalculate_bounding_box( RTreeNode* parent, RTreeNode* child, unsigned &child_in_parent );
+
+    void copy_group_a_to_existing_node( RTreeNode *position, RTreeNode* group_a );
+    RTreeRecord_NonLeaf create_nonleaf_record_from_rtreenode( Rect &new_entry_bounding, RTreeNode *rtreenode );
+    RTreeRecord_Leaf create_leaf_record_from_rtreenode( Rect &new_entry_bounding, RTreeNode *rtreenode );
+
+    // erase
+//    RTreeNode* find_leaf( RTreeNode* subtree, const RTreeRecord_Leaf &search ) const;
+    RTreeNode* find_leaf( RTreeNode* subtree, const Rect &search_area, const int shape_to_delete ) const;
+
+    bool condense_tree( RTreeNode* position
+    //                  std::pair<RTreeNode*, RTreeNode*>  &node_division, // modified: it holds the last split group
+    //                  bool initial_split_performed, 
+    //                    const unsigned min_nodes
+    //                    const unsigned max_nodes 
+                        );
+    int remove_record_from_parent( RTreeNode* parent, RTreeNode* child );
+    void sanity_check(RTreeNode* subtree_root, int depth, bool used_during_insert = false ) const;
+
+};
+
+} // end namespace Geom
+
+#endif // !SEEN_LIB2GEOM_RTREE_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/include/2geom/orphan-code/sbasis-of.h b/include/2geom/orphan-code/sbasis-of.h
new file mode 100644
index 0000000..e5b76d6
--- /dev/null
+++ b/include/2geom/orphan-code/sbasis-of.h
@@ -0,0 +1,638 @@
+/**

+ * \file

+ * \brief Defines S-power basis function class 

+ * with coefficient in arbitrary ring

+ *

+ *  Authors:

+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>

+ *   Michael Sloan <mgsloan@gmail.com>

+ *

+ * Copyright (C) 2006-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 SEEN_SBASIS_OF_H

+#define SEEN_SBASIS_OF_H

+#include <vector>

+#include <cassert>

+#include <iostream>

+

+#include <2geom/interval.h>

+#include <2geom/utils.h>

+#include <2geom/exception.h>

+

+#include <2geom/orphan-code/linear-of.h>

+

+namespace Geom {

+

+template<typename T>

+class SBasisOf;

+

+#ifdef USE_SBASIS_OF

+typedef Geom::SBasisOf<double> SBasis;

+#endif

+

+/*** An empty SBasisOf<T> is identically 0. */

+template<typename T>

+class SBasisOf : public std::vector<LinearOf<T> >{

+public:

+    SBasisOf() {}

+    explicit SBasisOf(T a) {

+        this->push_back(LinearOf<T>(a,a));

+    }

+    SBasisOf(SBasisOf<T> const & a) :

+        std::vector<LinearOf<T> >(a)

+    {}

+    SBasisOf(LinearOf<T> const & bo) {

+        this->push_back(bo);

+    }

+    SBasisOf(LinearOf<T>* bo) {

+        this->push_back(*bo);

+    }

+    //static unsigned input_dim(){return T::input_dim()+1;}

+

+    //IMPL: FragmentConcept

+    typedef T output_type;

+    inline bool isZero() const {

+        if(this->empty()) return true;

+        for(unsigned i = 0; i < this->size(); i++) {

+            if(!(*this)[i].isZero()) return false;

+        }

+        return true;

+    }

+    inline bool isConstant() const {

+        if (this->empty()) return true;

+        for (unsigned i = 0; i < this->size(); i++) {

+            if(!(*this)[i].isConstant()) return false;

+        }

+        return true;

+    }

+

+    //TODO: code this...

+    bool isFinite() const;

+

+    inline T at0() const { 

+        if(this->empty()) return T(0); else return (*this)[0][0];

+    }

+    inline T at1() const{

+        if(this->empty()) return T(0); else return (*this)[0][1];

+    }

+

+    T valueAt(double t) const {

+        double s = t*(1-t);

+        T p0 = T(0.), p1 = T(0.);

+        for(unsigned k = this->size(); k > 0; k--) {

+            const LinearOf<T> &lin = (*this)[k-1];

+            p0 = p0*s + lin[0];

+            p1 = p1*s + lin[1];

+        }

+        return p0*(1-t) + p1*t;

+    }

+

+    T operator()(double t) const {

+        return valueAt(t);

+    }

+

+    /**

+    *  The size of the returned vector equals n+1.

+    */

+    std::vector<T> valueAndDerivatives(double t, unsigned n) const{

+        std::vector<T> ret(n+1);

+        ret[0] = valueAt(t);

+        SBasisOf<T> tmp = *this;

+        for(unsigned i = 0; i < n; i++) {

+            tmp.derive();

+            ret[i+1] = tmp.valueAt(t);

+        }

+        return ret;

+    }

+

+    //The following lines only makes sens if T=double!

+    SBasisOf<T> toSBasis() const { return SBasisOf<T>(*this); }

+    double tailError(unsigned tail) const{

+        Interval bs = *bounds_fast(*this, tail);

+        return std::max(fabs(bs.min()),fabs(bs.max()));

+    }

+

+// compute f(g)

+    SBasisOf<T> operator()(SBasisOf<T> const & g) const;

+

+    LinearOf<T> operator[](unsigned i) const {

+        assert(i < this->size());

+        return std::vector<LinearOf<T> >::operator[](i);

+    }

+

+//MUTATOR PRISON

+    LinearOf<T>& operator[](unsigned i) { return this->at(i); }

+

+    //remove extra zeros

+    void normalize() {

+        while(!this->empty() && this->back().isZero())

+            this->pop_back();

+    }

+

+    void truncate(unsigned k) { if(k < this->size()) this->resize(k); }

+private:

+    void derive(); // in place version

+    unsigned dim;

+};

+

+//template<>

+//inline unsigned SBasisOf<double>::input_dim() { return 1; }

+

+//--------------------------------------------------------------------------

+#ifdef USE_SBASIS_OF

+

+//implemented in sbasis-roots.cpp

+OptInterval bounds_exact(SBasis const &a);

+OptInterval bounds_fast(SBasis const &a, int order = 0);

+OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0);

+

+std::vector<double> roots(SBasis const & s);

+std::vector<std::vector<double> > multi_roots(SBasis const &f,

+                                 std::vector<double> const &levels,

+                                 double htol=1e-7,

+                                 double vtol=1e-7,

+                                 double a=0,

+                                 double b=1);

+#endif

+//--------------------------------------------------------------------------

+

+

+//TODO: figure out how to stick this in linear, while not adding an sbasis dep

+template<typename T>

+inline SBasisOf<T> LinearOf<T>::toSBasis() const { return SBasisOf<T>(*this); }

+

+template<typename T>

+inline SBasisOf<T> reverse(SBasisOf<T> const &a) {

+    SBasisOf<T> result;

+    result.reserve(a.size());

+    for(unsigned k = 0; k < a.size(); k++)

+       result.push_back(reverse(a[k]));

+    return result;

+}

+

+//IMPL: ScalableConcept

+template<typename T>

+inline SBasisOf<T> operator-(const SBasisOf<T>& p) {

+    if(p.isZero()) return SBasisOf<T>();

+    SBasisOf<T> result;

+    result.reserve(p.size());

+        

+    for(unsigned i = 0; i < p.size(); i++) {

+        result.push_back(-p[i]);

+    }

+    return result;

+}

+

+template<typename T>

+SBasisOf<T> operator*(SBasisOf<T> const &a, double k){

+    SBasisOf<T> c;

+    //TODO: what does this mean for vectors of vectors??

+    //c.reserve(a.size());

+    for(unsigned i = 0; i < a.size(); i++)

+        c.push_back(a[i] * k);

+    return c;

+}

+

+template<typename T>

+inline SBasisOf<T> operator*(double k, SBasisOf<T> const &a) { return a*k; }

+template<typename T>

+inline SBasisOf<T> operator/(SBasisOf<T> const &a, double k) { return a*(1./k); }

+template<typename T>

+SBasisOf<T>& operator*=(SBasisOf<T>& a, double b){

+    if (a.isZero()) return a;

+    if (b == 0)

+        a.clear();

+    else

+        for(unsigned i = 0; i < a.size(); i++)

+            a[i] *= b;

+    return a;

+}

+

+template<typename T>

+inline SBasisOf<T>& operator/=(SBasisOf<T>& a, double b) { return (a*=(1./b)); }

+

+/*

+//We can also multiply by element of ring coeff T: 

+template<typename T>

+SBasisOf<T> operator*(SBasisOf<T> const &a, T k){

+    SBasisOf<T> c;

+    //TODO: what does this mean for vectors of vectors??

+    //c.reserve(a.size());

+    for(unsigned i = 0; i < a.size(); i++)

+        c.push_back(a[i] * k);

+    return c;

+}

+

+template<typename T>

+inline SBasisOf<T> operator*(T k, SBasisOf<T> const &a) { return a*k; }

+template<typename T>

+SBasisOf<T>& operator*=(SBasisOf<T>& a, T b){

+    if (a.isZero()) return a;

+    if (b == 0)

+        a.clear();

+    else

+        for(unsigned i = 0; i < a.size(); i++)

+            a[i] *= b;

+    return a;

+}

+*/

+

+//IMPL: AddableConcept

+template<typename T>

+inline SBasisOf<T> operator+(const SBasisOf<T>& a, const SBasisOf<T>& b){

+    SBasisOf<T> result;

+    const unsigned out_size = std::max(a.size(), b.size());

+    const unsigned min_size = std::min(a.size(), b.size());

+    //TODO: what does this mean for vector<vector>;

+    //result.reserve(out_size);

+

+    for(unsigned i = 0; i < min_size; i++) {

+        result.push_back(a[i] + b[i]);

+    }

+    for(unsigned i = min_size; i < a.size(); i++)

+        result.push_back(a[i]);

+    for(unsigned i = min_size; i < b.size(); i++)

+        result.push_back(b[i]);

+

+    assert(result.size() == out_size);

+    return result;

+}

+

+template<typename T>

+SBasisOf<T> operator-(const SBasisOf<T>& a, const SBasisOf<T>& b){

+    SBasisOf<T> result;

+    const unsigned out_size = std::max(a.size(), b.size());

+    const unsigned min_size = std::min(a.size(), b.size());

+    //TODO: what does this mean for vector<vector>;

+    //result.reserve(out_size);

+

+    for(unsigned i = 0; i < min_size; i++) {

+        result.push_back(a[i] - b[i]);

+    }

+    for(unsigned i = min_size; i < a.size(); i++)

+        result.push_back(a[i]);

+    for(unsigned i = min_size; i < b.size(); i++)

+        result.push_back(-b[i]);

+

+    assert(result.size() == out_size);

+    return result;

+}

+

+template<typename T>

+SBasisOf<T>& operator+=(SBasisOf<T>& a, const SBasisOf<T>& b){

+    const unsigned out_size = std::max(a.size(), b.size());

+    const unsigned min_size = std::min(a.size(), b.size());

+    //TODO: what does this mean for vectors of vectors

+    //a.reserve(out_size);

+    for(unsigned i = 0; i < min_size; i++)

+        a[i] += b[i];

+    for(unsigned i = min_size; i < b.size(); i++)

+        a.push_back(b[i]);

+

+    assert(a.size() == out_size);

+    return a;

+}

+

+template<typename T>

+SBasisOf<T>& operator-=(SBasisOf<T>& a, const SBasisOf<T>& b){

+    const unsigned out_size = std::max(a.size(), b.size());

+    const unsigned min_size = std::min(a.size(), b.size());

+    //TODO: what does this mean for vectors of vectors

+    //a.reserve(out_size);

+    for(unsigned i = 0; i < min_size; i++)

+        a[i] -= b[i];

+    for(unsigned i = min_size; i < b.size(); i++)

+        a.push_back(-b[i]);

+

+    assert(a.size() == out_size);

+    return a;

+}

+

+//TODO: remove?

+template<typename T>

+inline SBasisOf<T> operator+(const SBasisOf<T> & a, LinearOf<T> const & b) {

+    if(b.isZero()) return a;

+    if(a.isZero()) return b;

+    SBasisOf<T> result(a);

+    result[0] += b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T> operator-(const SBasisOf<T> & a, LinearOf<T> const & b) {

+    if(b.isZero()) return a;

+    SBasisOf<T> result(a);

+    result[0] -= b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T>& operator+=(SBasisOf<T>& a, const LinearOf<T>& b) {

+    if(a.isZero())

+        a.push_back(b);

+    else

+        a[0] += b;

+    return a;

+}

+template<typename T>

+inline SBasisOf<T>& operator-=(SBasisOf<T>& a, const LinearOf<T>& b) {

+    if(a.isZero())

+        a.push_back(-b);

+    else

+        a[0] -= b;

+    return a;

+}

+

+//IMPL: OffsetableConcept

+/*

+template<typename T>

+inline SBasisOf<T> operator+(const SBasisOf<T> & a, double b) {

+    if(a.isZero()) return LinearOf<T>(b, b);

+    SBasisOf<T> result(a);

+    result[0] += b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T> operator-(const SBasisOf<T> & a, double b) {

+    if(a.isZero()) return LinearOf<T>(-b, -b);

+    SBasisOf<T> result(a);

+    result[0] -= b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T>& operator+=(SBasisOf<T>& a, double b) {

+    if(a.isZero())

+        a.push_back(LinearOf<T>(b,b));

+    else

+        a[0] += b;

+    return a;

+}

+template<typename T>

+inline SBasisOf<T>& operator-=(SBasisOf<T>& a, double b) {

+    if(a.isZero())

+        a.push_back(LinearOf<T>(-b,-b));

+    else

+        a[0] -= b;

+    return a;

+}

+*/

+//We can also offset by elements of coeff ring T

+template<typename T>

+inline SBasisOf<T> operator+(const SBasisOf<T> & a, T b) {

+    if(a.isZero()) return LinearOf<T>(b, b);

+    SBasisOf<T> result(a);

+    result[0] += b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T> operator-(const SBasisOf<T> & a, T b) {

+    if(a.isZero()) return LinearOf<T>(-b, -b);

+    SBasisOf<T> result(a);

+    result[0] -= b;

+    return result;

+}

+template<typename T>

+inline SBasisOf<T>& operator+=(SBasisOf<T>& a, T b) {

+    if(a.isZero())

+        a.push_back(LinearOf<T>(b,b));

+    else

+        a[0] += b;

+    return a;

+}

+template<typename T>

+inline SBasisOf<T>& operator-=(SBasisOf<T>& a, T b) {

+    if(a.isZero())

+        a.push_back(LinearOf<T>(-b,-b));

+    else

+        a[0] -= b;

+    return a;

+}

+

+

+template<typename T>

+SBasisOf<T> shift(SBasisOf<T> const &a, int sh){

+    SBasisOf<T> c = a;

+    if(sh > 0) {

+        c.insert(c.begin(), sh, LinearOf<T>(0,0));

+    } else {

+        //TODO: truncate

+    }

+    return c;

+}

+

+template<typename T>

+SBasisOf<T> shift(LinearOf<T> const &a, int sh) {

+    SBasisOf<T> c;

+    if(sh > 0) {

+        c.insert(c.begin(), sh, LinearOf<T>(0,0));

+        c.push_back(a);

+    }

+    return c;

+}

+

+template<typename T>

+inline SBasisOf<T> truncate(SBasisOf<T> const &a, unsigned terms) {

+    SBasisOf<T> c;

+    c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));

+    return c;

+}

+

+template<typename T>

+SBasisOf<T> multiply_add(SBasisOf<T> const &a, SBasisOf<T> const &b, SBasisOf<T> c) {

+    if(a.isZero() || b.isZero())

+        return c;

+    c.resize(a.size() + b.size(), LinearOf<T>(T(0.),T(0.)));

+    for(unsigned j = 0; j < b.size(); j++) {

+        for(unsigned i = j; i < a.size()+j; i++) {

+            T tri = (b[j][1]-b[j][0])*(a[i-j][1]-a[i-j][0]);

+            c[i+1/*shift*/] += LinearOf<T>(-tri);

+        }

+    }

+    for(unsigned j = 0; j < b.size(); j++) {

+        for(unsigned i = j; i < a.size()+j; i++) {

+            for(unsigned dim = 0; dim < 2; dim++)

+                c[i][dim] += b[j][dim]*a[i-j][dim];

+        }

+    }

+    c.normalize();

+    //assert(!(0 == c.back()[0] && 0 == c.back()[1]));

+    return c;

+}

+

+template<typename T>

+SBasisOf<T> multiply(SBasisOf<T> const &a, SBasisOf<T> const &b) {

+    SBasisOf<T> c;

+    if(a.isZero() || b.isZero())

+        return c;

+    return multiply_add(a, b, c);

+}

+

+template<typename T>

+SBasisOf<T> integral(SBasisOf<T> const &c){

+    SBasisOf<T> a;

+    T aTri = T(0.);

+    for(int k = c.size()-1; k >= 0; k--) {

+        aTri = (HatOf<T>(c[k]).d + (k+1)*aTri/2)/(2*k+1);

+        a[k][0] -= aTri/2;

+        a[k][1] += aTri/2;

+    }

+    a.normalize();

+    return a;

+}

+

+template<typename T>

+SBasisOf<T> derivative(SBasisOf<T> const &a){

+    SBasisOf<T> c;

+    c.resize(a.size(), LinearOf<T>());

+    if(a.isZero())

+        return c;

+

+    for(unsigned k = 0; k < a.size()-1; k++) {

+        T d = (2*k+1)*(a[k][1] - a[k][0]);

+        

+        c[k][0] = d + (k+1)*a[k+1][0];

+        c[k][1] = d - (k+1)*a[k+1][1];

+    }

+    int k = a.size()-1;

+    T d = (2*k+1)*(a[k][1] - a[k][0]);

+    //TODO: do a real test to know if d==0!

+    if(d == T(0.0))

+        c.pop_back();

+    else {

+        c[k][0] = d;

+        c[k][1] = d;

+    }

+

+    return c;

+}

+

+template<typename T>

+void SBasisOf<T>::derive() { // in place version

+    if(isZero()) return;

+    for(unsigned k = 0; k < this->size()-1; k++) {

+        T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);

+        

+        (*this)[k][0] = d + (k+1)*(*this)[k+1][0];

+        (*this)[k][1] = d - (k+1)*(*this)[k+1][1];

+    }

+    int k = this->size()-1;

+    T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);

+    if(d == 0)//TODO: give this a meaning for general coeff ring.

+        this->pop_back();

+    else {

+        (*this)[k][0] = d;

+        (*this)[k][1] = d;

+    }

+}

+

+

+template<typename T>

+inline SBasisOf<T> operator*(SBasisOf<T> const & a, SBasisOf<T> const & b) {

+    return multiply(a, b);

+}

+

+template<typename T>

+inline SBasisOf<T>& operator*=(SBasisOf<T>& a, SBasisOf<T> const & b) {

+    a = multiply(a, b);

+    return a;

+}

+

+// a(b(t))

+//TODO: compose all compatibles types!

+template<typename T>

+SBasisOf<T> compose(SBasisOf<T> const &a, SBasisOf<T> const &b){

+    SBasisOf<double> s = multiply((SBasisOf<T>(LinearOf<T>(1,1))-b), b);

+    SBasisOf<T> r;

+

+    for(int i = a.size()-1; i >= 0; i--) {

+        r = multiply_add(r, s, SBasisOf<T>(LinearOf<T>(HatOf<T>(a[i][0]))) - b*a[i][0] + b*a[i][1]);

+    }

+    return r;

+}

+

+template<typename T>

+SBasisOf<T> compose(SBasisOf<T> const &a, SBasisOf<T> const &b, unsigned k){

+    SBasisOf<T> s = multiply((SBasisOf<T>(LinearOf<T>(1,1))-b), b);

+    SBasisOf<T> r;

+

+    for(int i = a.size()-1; i >= 0; i--) {

+        r = multiply_add(r, s, SBasisOf<T>(LinearOf<T>(HatOf<T>(a[i][0]))) - b*a[i][0] + b*a[i][1]);

+    }

+    r.truncate(k);

+    return r;

+}

+template<typename T>

+SBasisOf<T> compose(LinearOf<T> const &a, SBasisOf<T> const &b){

+    return compose(SBasisOf<T>(a),b);

+}

+template<typename T>

+SBasisOf<T> compose(SBasisOf<T> const &a, LinearOf<T> const &b){

+    return compose(a,SBasisOf<T>(b));

+}

+template<typename T>//TODO: don't be so lazy!!

+SBasisOf<T> compose(LinearOf<T> const &a, LinearOf<T> const &b){

+    return compose(SBasisOf<T>(a),SBasisOf<T>(b));

+}

+

+

+

+template<typename T>

+inline SBasisOf<T> portion(const SBasisOf<T> &t, double from, double to) { return compose(t, LinearOf<T>(from, to)); }

+

+// compute f(g)

+template<typename T>

+inline SBasisOf<T>

+SBasisOf<T>::operator()(SBasisOf<T> const & g) const {

+    return compose(*this, g);

+}

+ 

+template<typename T>

+inline std::ostream &operator<< (std::ostream &out_file, const LinearOf<T> &bo) {

+    out_file << "{" << bo[0] << ", " << bo[1] << "}";

+    return out_file;

+}

+

+template<typename T>

+inline std::ostream &operator<< (std::ostream &out_file, const SBasisOf<T> & p) {

+    for(unsigned i = 0; i < p.size(); i++) {

+        out_file << p[i] << "s^" << i << " + ";

+    }

+    return out_file;

+}

+

+};

+#endif

+

+

+/*

+  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/include/2geom/orphan-code/sbasisN.h b/include/2geom/orphan-code/sbasisN.h
new file mode 100644
index 0000000..0b5a48f
--- /dev/null
+++ b/include/2geom/orphan-code/sbasisN.h
@@ -0,0 +1,1123 @@
+/**

+ * \file

+ * \brief Multi-dimensional symmetric power basis function.

+ * A SBasisN<n> is a polynomial f of n variables (t0,...,tn-1), 

+ * written in a particular form. Let si = ti(1-t_i). f is written as

+ *

+ *  f = sum_p  s^p a_{p}(t0,...,t_{n-1})

+ *

+ * where p=(p0,...,p_{n-1}) is a multi index (called MultiDegree<n> in the code)

+ *  s^p = prod_i si^pi, and a_p is a LinearN<n>.

+ * Recall a LinearN<n> is sum over all choices xi = ti or (1-ti) of terms of form

+ *   a * x0*...*x_{n-1}

+ *

+ * Caution: degrees are expressed as degrees of s=t*(1-t). The real degree

+ * (with respect to t) of the polynomial is twice that + 0 or 1 depending 

+ * whether the relevant LinearN<n> coeff is constant or not.

+ *//*

+ *

+ *  Authors:

+ *   JF Barraud <jf.barraud@gmail.com>

+ *   Nathan Hurst <njh@mail.csse.monash.edu.au>

+ *   Michael Sloan <mgsloan@gmail.com>

+ *

+ * Copyright (C) 2006-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 SEEN_SBASISN_H

+#define SEEN_SBASISN_H

+#include <vector>

+#include <cassert>

+#include <iostream>

+

+#include <2geom/orphan-code/linearN.h>

+#include <2geom/linear.h>//for conversion purpose

+#include <2geom/sbasis.h>//for conversion purpose

+#include <2geom/interval.h>

+#include <2geom/utils.h>

+#include <2geom/exception.h>

+

+

+namespace Geom{

+

+/** MultiDegree:

+ *  \brief multi-degree (p0,...,p^{n-1}) of a s0^p0...s_{n-1}p^{n-1} monomial.

+ *

+ *  a "Multi_deg" is a sequence p={p0,...,p_n-1} representing the monomial

+ *  s^p = s_0^{p_0}*...*s_{n-1}^{p_{n-1}}.

+ *  Caution: the degrees are expressed with respect to si! ( in SBasis code

+ *  below, si = ti*(1-ti) ).

+ */

+

+template<unsigned n>

+class MultiDegree{

+public:

+    unsigned p[n];

+    MultiDegree(){

+        for (unsigned i = 0; i <n; i++) {

+            p[i]=0;

+        }

+    }

+    MultiDegree( unsigned const other_p[] ){

+        p = other_p;

+    }

+    MultiDegree(unsigned const idx, unsigned const sizes[]){

+        unsigned q = idx;

+        for (unsigned i = n-1; i >0; i--) {

+            div_t d = std::div(int(q), int(sizes[i]));

+            p[i] = d.rem;

+            q = d.quot;

+        }

+        p[0] = q;

+    }

+    unsigned operator[](const unsigned i) const {

+        assert(i < n); return p[i];

+    }

+    unsigned& operator[](const unsigned i) {

+        assert(i < n); return p[i]; 

+    }

+

+    unsigned asIdx(unsigned const sizes[]) const{

+        unsigned ret = p[0];

+        bool in_range = (p[0]<sizes[0]);

+        for (unsigned i = 1; i < n; i++) {

+            in_range = in_range && (p[i]<sizes[i]);

+            ret = ret*sizes[i] + p[i];

+        }

+        if (in_range) return ret;

+        //TODO: find a better warning than returning out of range idx!

+        ret =1;

+        for (unsigned i = 0; i < n; i++) {

+            ret *= sizes[i];

+        }

+        return ret;

+    }

+    bool stepUp(unsigned const sizes[], unsigned frozen_mask = 0){

+        unsigned i = 0;

+        while ( i < n && ( (1<<i) & frozen_mask ) ) i++;

+        while ( i <n && p[i] == sizes[i]-1 ) {

+            i++;

+            while (i<n && ( (1<<i) & frozen_mask ) ) i++;

+        }

+        if (i<n){

+            p[i]+=1; 

+            for (unsigned j = 0; j <  i; j++) {

+                if ( !( (1<<j) & frozen_mask ) ) p[j] = 0;

+            }

+            return true; 

+        }else{

+            return false;

+        }

+    }

+    bool stepDown(unsigned const sizes[], unsigned frozen_mask = 0){

+        int i = n-1;

+        while (i>=0 && ( (1<<i) & frozen_mask ) ) i--;

+        while ( i >= 0 && p[i] == 0 ) {

+            i--;

+            while (i>=0 && ( (1<<i) & frozen_mask ) ) i--;

+        }

+        if ( i >= 0 ){

+            p[i]-=1; 

+            for (unsigned j = i+1; j <  n; j++) {

+                if ( !( (1<<j) & frozen_mask ) ) p[j] = sizes[j]-1;

+            }

+            return true; 

+        }else{

+            return false;

+        }

+    }

+};

+

+/**

+ * Returns the maximal degree appearing in the two arguments for each variables.

+ */

+template <unsigned n>

+MultiDegree<n> max(MultiDegree<n> const &p, MultiDegree<n> const &q){

+    MultiDegree<n> ret;

+    for (unsigned i = 0; i <n; i++) {

+        ret.p[i] = (p[i]>q[i] ? p[i] : q[i]);

+    }

+    return ret;

+}

+

+template <unsigned n>

+MultiDegree<n> operator + (MultiDegree<n> const &p, MultiDegree<n> const &q){

+    MultiDegree<n> ret;

+    for (unsigned i = 0; i <n; i++) {

+        ret.p[i] = p[i] + q[i];

+    }

+    return ret;

+}

+template <unsigned n>

+MultiDegree<n> operator += (MultiDegree<n> const &p, MultiDegree<n> const &q){

+    for (unsigned i = 0; i <n; i++) {

+        p[i] += q[i];

+    }

+    return p;

+}

+

+/**

+ *  \brief MultiDegree comparison.

+ * A MultiDegree \param p is smaller than another \param q

+ * if all it's smaller for all variables.

+ *

+ * In particular, p<=q and q<=p can both be false! 

+ */

+template<unsigned n>

+bool operator<=(MultiDegree<n> const &p, MultiDegree<n> const &q){

+    for (unsigned i = 0; i <n; i++) {

+        if (p[i]>q[i]) return false;

+    }

+    return true;

+}

+

+

+/**

+ *  \brief Polynomials of n variables, written in SBasis form.

+ * An SBasisN<n> f is internaly represented as a vector of LinearN<n>.

+ * It should be thought of as an n-dimensional vector: the coef of s0^p0...s_{n-1}p^{n-1}

+ * is soterd in f[p0,...,p_{n-1}]. The sizes of each dimension is stored in "sizes". 

+ * Note: an empty SBasis is identically 0.

+ */

+template<unsigned n>

+class SBasisN : public std::vector<LinearN<n> >{

+public:

+    unsigned sizes[n];

+    SBasisN() {

+        for (unsigned i = 0; i < n; i++) {

+            sizes[i] = 0;

+        }

+    }

+    explicit SBasisN(double a) {

+        for (unsigned i = 0; i < n; i++) {

+            sizes[i] = 1;

+        }

+        this->push_back(LinearN<n>(a));

+    }

+    SBasisN(SBasisN<n> const & a) : std::vector<LinearN<n> >(a){

+        //TODO: efficient array copy??

+        for (unsigned i = 0; i < n; i++) {

+            sizes[i] = a.sizes[i];

+        }

+    }

+    SBasisN(LinearN<n> const & bo) {

+        for (unsigned i = 0; i < n; i++) {

+            sizes[i] = 1;

+        }

+        this->push_back(bo);

+    }

+    SBasisN(LinearN<n>* bo) {

+        for (unsigned i = 0; i < n; i++) {

+            sizes[i] = 1;

+        }

+        this->push_back(*bo); 

+    }

+

+//----------------------------------------------

+//-- Degree/Sizing facilities ------------------

+//----------------------------------------------

+/**

+ * Internal recursive function used to compute partial degrees.

+ */

+    bool find_non_empty_level(unsigned var, MultiDegree<n> &fixed_degrees)const{

+        if (this->size()==0){

+            for (unsigned i = 0; i < n; i++) {

+                fixed_degrees[i] = 0;//FIXME this should be -infty!!

+            }

+            return false;

+        }

+        if ( !((*this)[fixed_degrees.asIdx(sizes)].isZero()) ) return true;

+

+        unsigned frozen = (1<<var);

+        if ( fixed_degrees.stepDown(sizes, frozen) ){

+            if ( find_non_empty_level(var, fixed_degrees) ) return true;

+        }

+        if ( fixed_degrees[var] > 0 ){

+            fixed_degrees[var] -= 1;

+            for (unsigned i = 0; i < n; i++) {

+                if (i!=var) fixed_degrees[i] = sizes[i]-1;

+            }

+            if (find_non_empty_level(var, fixed_degrees)) return true;

+        }

+        return false;//FIXME: this should return -infty in all variables!

+    }

+

+/**

+ * Returns the degree of an SBasisN<n> with respect to a given variable form its sizes.

+ * All terms are taken into account, even eventual trailing zeros.

+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.

+ */

+    unsigned quick_degree(unsigned var) const{

+        return ( sizes[var] > 0 ? sizes[var]-1 : 0 );//this should be -infty.

+     }

+/**

+ * Computes the multi degree of the SBasis from it's sizes.

+ * All terms are taken into account, even eventual trailing zeros.

+ * Note: degrees are expressed with respect to s = t*(1-t), not t itself.

+ */

+    MultiDegree<n> quick_multi_degree() const{

+        MultiDegree<n> ret;

+        if (this->size()==0) return ret;//should be -infty for all vars.

+        for (unsigned i = 0; i < n; i++) {

+            assert( sizes[i]>0 );

+            ret.p[i] = sizes[i]-1;

+       }

+        return ret;

+    }

+/**

+ * Returns the degree of an SBasisN<n> with respect to a given variable.

+ * Trailing zeros are not taken into account.

+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.

+ */

+    unsigned degree(unsigned var)const{

+        MultiDegree<n> degrees;

+        for(unsigned i = 0; i < n; i++) {

+            degrees[i] = sizes[i]-1;

+        }

+        if ( find_non_empty_level(var, degrees) ) return degrees[var];

+        else return 0;//should be -infty.

+     }

+/**

+ * Returns the *real* degree of an SBasisN<n> with respect to a given variable.

+ * Trailing zeros are not taken into account.

+ * Note: degree is expressed with respect to t itself, not s = t*(1-t).

+ * In particular: real_t_degree() = 2*degree() + 0 or 1.

+ */

+    unsigned real_t_degree(unsigned var)const{

+        unsigned deg = 0;

+        bool even = true;

+        bool notfound = true;

+        unsigned frozen = (1<<var);

+        MultiDegree<n> degrees;

+        for(unsigned i = 0; i < n; i++) {

+            degrees[i] = sizes[i]-1;

+        }

+        while( notfound ){

+            if ( find_non_empty_level(var, degrees) && degrees[var]>= deg ){

+                deg = degrees[var];

+                even = (*this)[degrees.asIdx(sizes)].isConstant(var);

+            }

+            notfound =  even && degrees.stepDown(sizes, frozen);

+        }

+        return 2*deg + ( even ? 0 : 1 );

+     }

+/**

+ * Returns the *real* degrees of an SBasisN<n>.

+ * Trailing zeros are not taken into account.

+ * Note: degree is expressed with respect to t itself, not s = t*(1-t).

+ * In particular: real_t_degree() = 2*degree() + 0 or 1.

+ */

+    MultiDegree<n> real_t_degrees()const{

+        MultiDegree<n>res;

+        for(unsigned i = 0; i < n; i++) {

+            res[i] = real_t_degree(i);

+        }

+        return res;

+     }

+/**

+ * Computes the multi degree of the SBasis.

+ * Trailing zeros are not taken into account.

+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.

+ */

+    MultiDegree<n> multi_degree() const{

+        MultiDegree<n> ret;

+        if (this->size()==0) return ret;//should be -infty for all vars.

+        for (unsigned i = 0; i < n; i++) {

+            ret[i] = this->degree(i);

+        }

+        return ret;

+    }

+/**

+ * Returns the highest degree over all variables.

+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.

+ */

+    unsigned max_degree() const {

+        if (this->size()==0) return 0;//should be -infty!

+        unsigned d=0;

+        for (unsigned i = 0; i < n; i++) {

+            assert( sizes[i]>0 );

+            if (d < sizes[i]-1) d = sizes[i]-1;

+        }

+        return d;

+    }

+

+/**

+ * Resize an SBasisN<n> to match new sizes.

+ *

+ * Caution: if a new size is smaller, the corresponding coefficients are discarded.

+ */

+    void multi_resize(unsigned new_sizes[], LinearN<n> def_value = LinearN<n>(0.)){

+        SBasisN<n> result;

+        bool nothing_todo = true;

+        unsigned tot_size = 1;

+        for(unsigned i = 0; i < n; i++) {

+            nothing_todo = nothing_todo && (sizes[i] == new_sizes[i]);

+            result.sizes[i] = new_sizes[i];

+            tot_size *= new_sizes[i];

+        }

+        if (nothing_todo) return;

+        result.resize(tot_size, def_value);

+        for(unsigned i = 0; i < tot_size; i++) {

+            MultiDegree<n> d( i, result.sizes );

+            unsigned j = d.asIdx(sizes);

+            if ( j < this->size() ){

+                result[i] = (*this)[j];

+            }

+        }

+        *this = result;

+    }

+

+    //remove extra zeros

+    void normalize() {

+        MultiDegree<n> max_p = multi_degree();

+        unsigned new_sizes[n];

+        for (unsigned i=0; i<n; i++){

+            new_sizes[i] = max_p[i]+1;

+        }

+        multi_resize(new_sizes);

+    }

+

+//-----------------------------

+//-- Misc. --------------------

+//-----------------------------

+

+/**

+ * Returns the number of variables this function takes as input: n.

+ */

+    unsigned input_dim(){return n;};

+

+    //IMPL: FragmentConcept

+    typedef double output_type;

+

+    inline bool isZero() const {

+        if(this->size()==0) return true;

+        for(unsigned i = 0; i < this->size(); i++) {

+            if(!(*this)[i].isZero()) return false;

+        }

+        return true;

+    }

+    inline bool isConstant() const {

+        if (this->size()==0) return true;

+        if(!(*this)[0].isConstant()) return false;

+        for (unsigned i = 1; i < this->size(); i++) {

+            if(!(*this)[i].isZero()) return false;

+        }

+        return true;

+    }

+

+    bool isFinite() const{

+        for (unsigned i = 0; i < this->size(); i++) {

+            if(!(*this)[i].isFinite()) return false;

+        }

+        return true;

+    }

+

+

+//------------------------------------------

+//-- Evaluation methods --------------------

+//------------------------------------------

+/**

+ * Returns the value of the SBasis at a given corner of [0,1]^n.

+ * \param k describes the corner: if i-th bit is 0, ti=0, otherwise ti=1. 

+ */

+    inline double atCorner(unsigned k) const { 

+        if(this->size()==0) return 0.;

+        return (*this)[0].atCorner(k);

+    }

+/**

+ * Returns the value of the SBasis at a given corner of [0,1]^n.

+ * \param t[n] describes the corner: the values should be 0's and 1's. 

+ */

+    inline double atCorner(double t[]) const { 

+        if(this->size()==0) return 0.;

+        return (*this)[0].atCorner(t);

+    }

+/**

+ * Returns a "slice" of the array.

+ * Returns an SBasis containing all the coeff of (s-)degree \param deg in variable \param var

+ */

+    //TODO: move by bigger blocks (?) but they are broken into pieces in general...

+    SBasisN<n> slice(unsigned const var, unsigned const deg) const{

+        if (deg >= sizes[var] ) return SBasisN<n>();

+        SBasisN<n> res;

+        unsigned tot_size = 1;

+        for (unsigned i = 0; i < n; i++) {

+            res.sizes[i] = (i==var ? 1 : sizes[i]);

+            tot_size *= res.sizes[i];

+        }

+        res.resize( tot_size, LinearN<n>(0.));

+        for (unsigned i = 0; i < tot_size; i++) {

+            MultiDegree<n> d(i,res.sizes);

+            d.p[var] = deg;

+            res[i] = (*this)[d.asIdx(sizes)];

+        }

+        return res;

+    }

+/**

+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to 0.

+ */

+    inline SBasisN<n-1> at0(unsigned var=0, unsigned deg=0) const {

+        SBasisN<n> sl = slice(var,deg);

+        SBasisN<n-1> res;

+        res.reserve(sl.size());

+        for (unsigned i = 0; i < n-1; i++) {

+            res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];

+        }

+        for (unsigned i = 0; i < sl.size(); i++) {

+            res.push_back( sl[i].at0(var) );

+        }

+        return res;

+    }

+/**

+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to 1.

+ */

+    inline SBasisN<n-1> at1(unsigned var=0, unsigned deg=0) const {

+        SBasisN<n> sl = slice(var,deg);

+        SBasisN<n-1> res;

+        res.reserve(sl.size());

+        for (unsigned i = 0; i < n-1; i++) {

+            res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];

+        }

+        for (unsigned i = 0; i < sl.size(); i++) {

+            res.push_back( sl[i].at1(var) );

+        }

+        return res;

+    }

+/**

+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to \param t.

+ */

+    inline SBasisN<n-1> partialEval(double t, unsigned var=0 ) const {

+        SBasisN<n> sl;

+        double s = t*(1-t);

+        double si = 1;

+        for (unsigned i = 0; i <sizes[var]; i++) {

+            sl = sl + slice(var, i)*si;

+            si *= s;

+        }

+        SBasisN<n-1> res;

+        res.resize(sl.size(), LinearN<n-1>(0.));

+        for (unsigned i = 0; i < n-1; i++) {

+            res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];

+        }

+        for (unsigned i = 0; i < sl.size(); i++) {

+            res[i] = sl[i].partialEval(t,var);

+        }

+        return res;

+    }

+

+/**

+ * \brief Internal recursive function.

+ * Replace each variable by it's value in the 's=t*(1-t)' factor 

+ * but not in the LinearN<n> coeffs. Then sum up all coefficients.

+ * \param t[n]: values of the variables.

+ */

+    LinearN<n> sumCoefs( double t[], unsigned const k, unsigned const idx) const{

+        LinearN<n> a;

+        if (k == n){

+            a = (*this)[idx];

+            return (*this)[idx];

+        }

+        double s = t[k]*(1-t[k]);

+        double si=1;

+        for (unsigned i=0; i<sizes[k]; i++){

+            a += sumCoefs(t,k+1,idx*sizes[k]+i)*si;;

+            si *= s;

+        }

+        return a;

+    }

+/**

+ * Evaluate at given n-dimensional point.

+ * \param t[n]: values of the variables.

+ */

+    double valueAt(double t[]) const {

+        LinearN<n> a = sumCoefs(t,0,0);

+        return a.valueAt(t);

+    }

+

+    double operator()(double t[]) const {

+        return valueAt(t);

+    }

+

+    //double valueAndDerivative(double t, double &der) const;

+    //std::vector<double> valueAndDerivatives(double t, unsigned n) const;

+    //SBasisN toSBasisN() const { return SBasisN(*this); }

+    //double tailError(unsigned tail) const;

+

+

+//--------------------------------------------------

+//-- Coeff. manipulation ---------------------------

+//--------------------------------------------------

+

+/**

+ * Accessing the SBasisN<n> coefficients.

+ */

+    LinearN<n> operator[](unsigned i) const {

+        assert(i < this->size());

+        return std::vector<LinearN<n> >::operator[](i);

+    }

+    LinearN<n> operator[](MultiDegree<n> const &p) const {

+        unsigned i = p.asIdx(sizes);

+        assert(i < this->size());        

+        return std::vector<LinearN<n> >::operator[](i);

+    }

+

+//MUTATOR PRISON

+    LinearN<n>& operator[](unsigned i) { return this->at(i); }

+//    LinearN<n>& operator[](MultiDegree const &p) { 

+//        unsigned i = p.asIdx(sizes);

+//        return this->at(i); 

+//    }

+

+    void appendCoef(const SBasisN<n-1> &a, const SBasisN<n-1> &b, unsigned var=0){

+        unsigned new_sizes[n];

+        MultiDegree<n-1> deg_a = a.multi_degree(), deg_b = b.multi_degree();

+        MultiDegree<n-1> dcoef = max( deg_a, deg_b );

+        for (unsigned i=0; i<n; i++){

+            if ( i == var ){

+                new_sizes[var] = sizes[var] + 1;

+            }else{

+                unsigned coef_size = dcoef[(i<var?i:i-1)] + 1;

+                new_sizes[i] = ( sizes[i]>coef_size ? sizes[i] : coef_size );

+            }

+        }

+        multi_resize(new_sizes);

+        

+        MultiDegree<n> d;

+        d[var] = sizes[var]-1;

+        unsigned frozen_mask = (1<<var);

+        do{

+            for (unsigned i=0; i<n-1; i++){

+                dcoef.p[i] = d.p[ ( i<var ? i : i+1) ];

+            }

+            LinearN<n-1> a_d,b_d;

+            unsigned ia = dcoef.asIdx(a.sizes);

+            if ( ia < a.size() ) a_d = a[ia];

+            unsigned ib = dcoef.asIdx(b.sizes);

+            if ( ib < b.size() ) b_d = b[ib];

+            (*this)[d.asIdx(sizes)] = LinearN<n>(a_d,b_d);

+        }while (d.stepUp(sizes,frozen_mask));

+    }

+

+//private:

+    //void derive(); // in place version

+};

+

+//SBasisN<0> is a double. Specialize it out.

+template<>

+class SBasisN<0>{

+public:

+    double d;

+    SBasisN () {}

+    SBasisN(double d) :d(d) {}

+    operator double() const { return d; }

+};

+

+

+//SBasisN<1> are usual SBasis. Allow conversion.

+SBasis toSBasis(SBasisN<1> f){

+    SBasis res(f.size(), Linear());

+    for (unsigned i = 0; i < f.size(); i++) {

+        res[i] = toLinear(f[i]);

+    }

+    return res;

+}

+

+//TODO: figure out how to stick this in linear, while not adding an sbasis dep

+template<unsigned n>

+inline SBasisN<n> LinearN<n>::toSBasisN() const { return SBasisN<n>(*this); }

+

+

+

+

+//implemented in sbasis-roots.cpp

+//OptInterval bounds_exact(SBasisN const &a);

+//OptInterval bounds_fast(SBasisN const &a, int order = 0);

+//OptInterval bounds_local(SBasisN const &a, const OptInterval &t, int order = 0);

+

+/** Returns a function which reverses the domain of a.

+ \param a sbasis function

+

+useful for reversing a parameteric curve.

+*/

+//template<unsigned n>

+//inline SBasisN<n> reverse(SBasisN<n> const &a);

+

+//IMPL: ScalableConcept

+template<unsigned n>

+inline SBasisN<n> operator-(const SBasisN<n>& p) {

+    if(p.isZero()) return SBasisN<n>();

+    SBasisN<n> result;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = p.sizes[i];

+    }

+    result.reserve(p.size());

+    for(unsigned i = 0; i < p.size(); i++) {

+        result.push_back(-p[i]);

+    }

+    return result;

+}

+template<unsigned n>

+SBasisN<n> operator*(SBasisN<n> const &a, double c){

+    if(a.isZero()) return SBasisN<n>();

+    SBasisN<n> result;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = a.sizes[i];

+    }

+    result.reserve(a.size());

+    for(unsigned i = 0; i < a.size(); i++) {

+        result.push_back(a[i] * c);

+    }

+    return result;

+}

+template<unsigned n>

+inline SBasisN<n> operator*(double k, SBasisN<n> const &a) { return a*k; }

+template<unsigned n>

+inline SBasisN<n> operator/(SBasisN<n> const &a, double k) { return a*(1./k); }

+template<unsigned n>

+SBasisN<n>& operator*=(SBasisN<n>& a, double c){

+    for(unsigned i = 0; i < a.size(); i++) a[i] *= c;

+    return a;

+}

+template<unsigned n>

+inline SBasisN<n>& operator/=(SBasisN<n>& a, double b) { return (a*=(1./b)); }

+

+//IMPL: AddableConcept

+template<unsigned n>

+SBasisN<n> operator + (const SBasisN<n>& a, const SBasisN<n>& b){

+    if( a.isZero() ) return b;

+    if( b.isZero() ) return a;

+    SBasisN<n> result;

+    MultiDegree<n> deg = max(a.quick_multi_degree(),b.quick_multi_degree());

+    unsigned max_size = 1;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = deg[i]+1;

+        max_size *= result.sizes[i];

+    }

+    result.resize( max_size, LinearN<n>(0.) );

+    for(unsigned i = 0; i < result.size(); i++) {

+        MultiDegree<n> p(i,result.sizes);

+        unsigned ia = p.asIdx(a.sizes);

+        unsigned ib = p.asIdx(b.sizes);

+        if (ia<a.size()) {

+            result[i] += a[ia];

+        }

+        if (ib<b.size()) {

+            result[i] += b[ib];

+        }

+    }

+    return result;

+}

+template<unsigned n>

+SBasisN<n> operator-(const SBasisN<n>& a, const SBasisN<n>& b){return a+(-b);}

+template<unsigned n>

+SBasisN<n>& operator+=(SBasisN<n>& a, const SBasisN<n>& b){

+    if(b.isZero()) return a;

+    a = a + b;

+    return a;

+}

+template<unsigned n>

+SBasisN<n>& operator-=(SBasisN<n>& a, const SBasisN<n>& b){

+    a += -b;

+    return a;

+}

+

+//TODO: remove?

+template<unsigned n>

+inline SBasisN<n> operator+(const SBasisN<n> & a, LinearN<n> const & b) {

+    if(b.isZero()) return a;

+    if(a.isZero()) return b;

+    SBasisN<n> result(a);

+    result[0] += b;

+    return result;

+}

+template<unsigned n>

+

+inline SBasisN<n> operator-(const SBasisN<n> & a, LinearN<n> const & b) {

+    if(b.isZero()) return a;

+    if(a.isZero()) return -b;

+    SBasisN<n> result(a);

+    result[0] -= b;

+    return result;

+}

+template<unsigned n>

+inline SBasisN<n>& operator+=(SBasisN<n>& a, const LinearN<n>& b) {

+    if(a.size()==0)

+        a.push_back(b);

+    else

+        a[0] += b;

+    return a;

+}

+template<unsigned n>

+inline SBasisN<n>& operator-=(SBasisN<n>& a, const LinearN<n>& b) {

+    if(a.size()==0)

+        a.push_back(-b);

+    else

+        a[0] -= b;

+    return a;

+}

+

+//IMPL: OffsetableConcept

+template<unsigned n>

+inline SBasisN<n> operator+(const SBasisN<n> & a, double b) {

+    if(a.isZero()) return LinearN<n>(b);

+    SBasisN<n> result(a);

+    result[0] += b;

+    return result;

+}

+template<unsigned n>

+inline SBasisN<n> operator-(const SBasisN<n> & a, double b) {

+    if(a.isZero()) return LinearN<n>(-b);

+    SBasisN<n> result(a);

+    result[0] -= b;

+    return result;

+}

+template<unsigned n>

+inline SBasisN<n>& operator+=(SBasisN<n>& a, double b) {

+    if(a.size()==0)

+        a.push_back(LinearN<n>(b));

+    else

+        a[0] += b;

+    return a;

+}

+template<unsigned n>

+inline SBasisN<n>& operator-=(SBasisN<n>& a, double b) {

+    if(a.size()==0)

+        a.push_back(LinearN<n>(-b));

+    else

+        a[0] -= b;

+    return a;

+}

+

+template<unsigned n>

+SBasisN<n> shift(SBasisN<n> const &a, MultiDegree<n> sh){

+    SBasisN<n> result;

+    MultiDegree<n> deg = a.quick_multi_degree() + sh;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = deg[i]+1;

+    }

+    unsigned max_size = deg.asIdx(result.sizes);

+    result.resize( max_size, LinearN<n>(0.) );

+    for(unsigned i = 0; i < a.size(); i++) {

+        MultiDegree<n> p(i,a.sizes);

+        p+=sh;

+        result[p.asIdx(result.sizes)]=a[i];

+    }

+    return result;

+}

+template<unsigned n>

+SBasisN<n> shift(LinearN<n> const &a, MultiDegree<n> sh){

+    SBasisN<n> result;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = sh[i]+1;

+    }

+    unsigned max_size = sh.asIdx(result.sizes);

+    result.resize( max_size, LinearN<n>(0.) );

+    result[max_size-1]=a;

+    return result;

+}

+//shift only in one variable

+template<unsigned n>

+SBasisN<n> shift(LinearN<n> const &a, unsigned sh, unsigned var){

+    assert( var < n );

+    SBasisN<n> result;

+    for(unsigned i = 0; i < n; i++) {

+        result.sizes[i] = 1;

+    }

+    result.sizes[var] = sh+1;

+    result.resize( sh+1, LinearN<n>(0.) );

+    result[sh]=a;

+    return result;

+}

+

+//truncate only in first variable

+template<unsigned n>

+inline SBasisN<n> truncate(SBasisN<n> const &a, unsigned first_size) {

+    if ( first_size <= a.sizes[0] ) return a;

+    SBasisN<n> c;

+    for (unsigned i = 0; i < n; i++) {

+        c.sizes[i] = a.sizes[i];

+    }

+    c.sizes[0] = first_size;

+    unsigned tot_size = 1;

+    for(unsigned i = 0; i < n; i++) {

+        tot_size*=c.sizes[i];

+    }

+    c.insert(c.begin(), a.begin(), a.begin() + tot_size);

+    return c;

+}

+

+template<unsigned n>

+SBasisN<n> multiply(SBasisN<n> const &a, SBasisN<n> const &b){

+    SBasisN<n> c;

+    MultiDegree<n> d;

+    MultiDegree<n> t_deg = a.real_t_degrees() + b.real_t_degrees();

+    for(unsigned i = 0; i < n; i++) {

+        d[i] = ( t_deg[i]%2 == 0 ? t_deg[i]/2 : (t_deg[i]-1)/2 ) ;

+    }

+    unsigned new_sizes[n], tot_size = 1;

+    for(unsigned i = 0; i < n; i++) {

+        //c.sizes[i] = d[i] + 1+1;//product of linears might give 1 more s in each dir!!

+        new_sizes[i] = d[i] + 1;

+        tot_size*=new_sizes[i];

+    }

+    c.resize( tot_size, LinearN<n>(0.) );

+    for(unsigned i = 0; i < n; i++) {

+        c.sizes[i] = new_sizes[i];

+    }

+

+    for(unsigned ai = 0; ai < a.size(); ai++) {

+        for(unsigned bj = 0; bj < b.size(); bj++) {

+            MultiDegree<n> di( ai, a.sizes ); 

+            MultiDegree<n> dj( bj, b.sizes ); 

+            //compute a[ai]*b[bj]:

+            for(unsigned p = 0; p < (1<<n); p++) {

+                for(unsigned q = 0; q < (1<<n); q++) {

+

+                    //compute a[ai][p]*b[bj][q]:

+                    unsigned m = p^q;//m has ones for factors s, 0 for (t-s) or ((1-t)-s).

+                    for(unsigned r = 0; r < (1<<n); r++) {

+                        if (!(r&m)) {// a 1 in r means take t (or (1-t)), otherwise take -s.

+                            int sign = 1;

+                            MultiDegree<n> dr; 

+                            unsigned t0 = 0, t1 = 0;

+                            for (unsigned var = 0; var < n; var++) {

+                                //if var is in mask m, no choice, take s

+                                if ( m & (1<<var) ){

+                                    dr.p[var] = 1;

+                                }//if var is in mask r, take t or (1-t)

+                                else  if ( r & (1<<var) ){

+                                    dr.p[var] = 0;

+                                    if ( p&(1<<var) ) {

+                                        t0 = t0 | (1<<var);

+                                    }else{

+                                        t1 = t1 | (1<<var);

+                                    }

+                                }//ohterwise take -s

+                                else{

+                                    dr.p[var] = 1;

+                                    sign *= -1;

+                                }

+                            }

+                            unsigned idx = (di+dj+dr).asIdx(c.sizes);

+                            if (idx < c.size()){

+                                for(unsigned s = 0; s < (1<<n); s++) {

+                                    if ( (t0 & ~s) || (t1 & s) ){

+                                        c[idx][s] += 0;

+                                    }else{

+                                        c[idx][s] += sign * a[ai][p] * b[bj][q];

+                                    }

+                                }

+                            }

+                        }

+                    }//r loop: all choices have been expanded in the product a[ai][p]*b[bj][q]

+                }//q loop

+            }//p loop: all products a[ai][p]*b[bj][q] have been computed.

+        }//bj loop

+    }//ai loop: all a[ai]b[bj] have been computed.

+

+    //TODO: normalize c, or even better, compute with the right size from scratch

+    return c;

+}

+

+

+template<unsigned n>

+inline SBasisN<n> operator*(SBasisN<n> const & a, SBasisN<n> const & b) {

+    return multiply(a, b);

+}

+

+template<unsigned n>

+inline SBasisN<n>& operator*=(SBasisN<n>& a, SBasisN<n> const & b) {

+    a = multiply(a, b);

+    return a;

+}

+

+template<unsigned m,unsigned n>

+SBasisN<m> compose(LinearN<n> const &f, std::vector<SBasisN<m> > const &t, unsigned fixed=0, unsigned flags=0 ){

+    assert (t.size() == n );

+    if (fixed == n) {

+        return SBasisN<m>(1.) * f[flags]; 

+    }else{

+        SBasisN<m> a0 = compose(f, t, fixed+1, flags);

+        SBasisN<m> a1 = compose(f, t, fixed+1, flags|(1<<fixed));

+        return (-t[fixed]+1) * a0 + t[fixed] * a1;

+    }

+}

+

+template<unsigned m,unsigned n>

+SBasisN<m> compose(SBasisN<n> const &f, std::vector<SBasisN<m> > const &t, unsigned const k=0, unsigned const idx = 0){

+    assert (t.size() == n );

+    if (k == n){

+        return compose( f[idx], t);

+    }

+    SBasisN<m> a;

+    SBasisN<m> s = multiply( t[k], (-t[k]+1.) );

+    SBasisN<m> si= SBasisN<m>(1.);

+    for (unsigned i=0; i<f.sizes[k]; i++){

+        a += compose(f, t,k+1,idx*f.sizes[k]+i)*si;;

+        si *= s;

+    }

+    return a;

+}

+

+template <unsigned n>

+inline std::ostream &operator<< (std::ostream &out_file, const MultiDegree<n> & d) {

+    out_file << "s^{";

+    for(unsigned i = 0; i < n; i++) {

+        out_file << d[i] << (i == n-1 ? "}" : ",");        

+    }

+    return out_file;

+}

+template <unsigned n>

+inline std::ostream &operator<< (std::ostream &out_file, const SBasisN<n> & p) {

+    for(unsigned i = 0; i < p.size(); i++) {

+        MultiDegree<n> d(i, p.sizes);

+        out_file << d << " " << p[i] << " + ";

+    }

+    return out_file;

+}

+

+

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+

+#if 0

+

+

+// This performs a multiply and accumulate operation in about the same time as multiply.  return a*b + c

+template<unsigned n>

+SBasisN<n> multiply_add(SBasisN<n> const &a, SBasisN<n> const &b, SBasisN<n> c);

+

+template<unsigned n>

+SBasisN<n> integral(SBasisN<n> const &c);

+template<unsigned n>

+SBasisN<n> derivative(SBasisN<n> const &a);

+

+template<unsigned n>

+SBasisN<n> sqrt(SBasisN<n> const &a, int k);

+

+// return a kth order approx to 1/a)

+template<unsigned n>

+SBasisN<n> reciprocal(LinearN<n> const &a, int k);

+template<unsigned n>

+SBasisN<n> divide(SBasisN<n> const &a, SBasisN<n> const &b, int k);

+

+

+/** Returns the degree of the first non zero coefficient.

+ \param a sbasis function

+ \param tol largest abs val considered 0

+ \returns first non zero coefficient

+*/

+template<unsigned n>

+inline unsigned 

+valuation(SBasisN<n> const &a, double tol=0){

+    unsigned val=0;

+    while( val<a.size() &&

+           fabs(a[val][0])<tol &&

+           fabs(a[val][1])<tol ) 

+        val++;

+    return val;

+}

+

+// a(b(t))

+template<unsigned n>

+SBasisN<n> compose(SBasisN<n> const &a, SBasisN<n> const &b);

+template<unsigned n>

+SBasisN<n> compose(SBasisN<n> const &a, SBasisN<n> const &b, unsigned k);

+template<unsigned n>

+SBasisN<n> inverse(SBasisN<n> a, int k);

+//compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...

+//TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.

+template<unsigned n>

+SBasisN<n> compose_inverse(SBasisN<n> const &f, SBasisN<n> const &g, unsigned order=2, double tol=1e-3);

+

+/** Returns the sbasis on domain [0,1] that was t on [from, to]

+ \param a sbasis function

+ \param from,to interval

+ \returns sbasis

+

+*/

+template<unsigned n>

+inline SBasisN<n> portion(const SBasisN<n> &t, double from, double to) { return compose(t, LinearN<n>(from, to)); }

+

+// compute f(g)

+template<unsigned n>

+inline SBasisN<n>

+SBasisN<n>::operator()(SBasisN<n> const & g) const {

+    return compose(*this, g);

+}

+ 

+template<unsigned n>

+inline std::ostream &operator<< (std::ostream &out_file, const LinearN<n> &bo) {

+    out_file << "{" << bo[0] << ", " << bo[1] << "}";

+    return out_file;

+}

+

+template<unsigned n>

+inline std::ostream &operator<< (std::ostream &out_file, const SBasisN<n> & p) {

+    for(unsigned i = 0; i < p.size(); i++) {

+        out_file << p[i] << "s^" << i << " + ";

+    }

+    return out_file;

+}

+

+// These are deprecated, use sbasis-math.h versions if possible

+template<unsigned n>

+SBasisN<n> sin(LinearN<n> bo, int k);

+template<unsigned n>

+SBasisN<n> cos(LinearN<n> bo, int k);

+

+template<unsigned n>

+std::vector<double> roots(SBasisN<n> const & s);

+template<unsigned n>

+std::vector<std::vector<double> > multi_roots(SBasisN<n> const &f,

+                                 std::vector<double> const &levels,

+                                 double htol=1e-7,

+                                 double vtol=1e-7,

+                                 double a=0,

+                                 double b=1);

+    

+#endif

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+//--------------------------------------------------

+}

+

+/*

+  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 :

+#endif