diff options
Diffstat (limited to 'include/2geom/orphan-code')
| -rw-r--r-- | include/2geom/orphan-code/arc-length.h | 58 | ||||
| -rw-r--r-- | include/2geom/orphan-code/chebyshev.h | 30 | ||||
| -rw-r--r-- | include/2geom/orphan-code/intersection-by-smashing.h | 78 | ||||
| -rw-r--r-- | include/2geom/orphan-code/linear-of.h | 269 | ||||
| -rw-r--r-- | include/2geom/orphan-code/linearN.h | 363 | ||||
| -rw-r--r-- | include/2geom/orphan-code/redblacktree.h | 121 | ||||
| -rw-r--r-- | include/2geom/orphan-code/rtree.h | 241 | ||||
| -rw-r--r-- | include/2geom/orphan-code/sbasis-of.h | 638 | ||||
| -rw-r--r-- | include/2geom/orphan-code/sbasisN.h | 1123 |
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
|