CRLF fix

(bzr r4247)

7 files changed, 1511 insertions, 1511 deletions

diff --git a/src/2geom/ord.h b/src/2geom/ord.h
index d0e348aec..735de6a3a 100644
--- a/src/2geom/ord.h
+++ b/src/2geom/ord.h
@@ -4,11 +4,11 @@
 
 namespace {
 
-enum Cmp {

-  LESS_THAN=-1,

-  GREATER_THAN=1,

-  EQUAL_TO=0

-};

+enum Cmp {
+  LESS_THAN=-1,
+  GREATER_THAN=1,
+  EQUAL_TO=0
+};
 
 inline Cmp operator-(Cmp x) {
   switch(x) {
@@ -20,16 +20,16 @@ inline Cmp operator-(Cmp x) {
     return EQUAL_TO;
   }
 }
-

-template <typename T1, typename T2>

-inline Cmp cmp(T1 const &a, T2 const &b) {

-  if ( a < b ) {

-    return LESS_THAN;

-  } else if ( b < a ) {

-    return GREATER_THAN;

-  } else {

-    return EQUAL_TO;

-  }

+
+template <typename T1, typename T2>
+inline Cmp cmp(T1 const &a, T2 const &b) {
+  if ( a < b ) {
+    return LESS_THAN;
+  } else if ( b < a ) {
+    return GREATER_THAN;
+  } else {
+    return EQUAL_TO;
+  }
 }
 
 }
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index c70ecd42c..171599768 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -1,690 +1,690 @@
-/*

- * piecewise.h - Piecewise function class

- *

- * Copyright 2007 Michael Sloan <mgsloan@gmail.com>

- *

- * 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, output 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_GEOM_PW_SB_H

-#define SEEN_GEOM_PW_SB_H

-

-#include "sbasis.h"

-#include <vector>

-#include <map>

-

-#include "concepts.h"

-#include "isnan.h"

-#include <boost/concept_check.hpp>

-

-namespace Geom {

-

-template <typename T>

-class Piecewise {

-  BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);

-

-  public:

-    std::vector<double> cuts;

-    std::vector<T> segs;

-    //segs[i] stretches from cuts[i] to cuts[i+1].

-

-    Piecewise() {}

-

-    explicit Piecewise(const T &s) {

-        push_cut(0.);

-        push_seg(s);

-        push_cut(1.);

-    }

-

-    typedef typename T::output_type output_type;

-

-    explicit Piecewise(const output_type & v) {

-        push_cut(0.);

-        push_seg(T(v));

-        push_cut(1.);

-    }

-

-    inline T operator[](unsigned i) const { return segs[i]; }

-    inline T &operator[](unsigned i) { return segs[i]; }

-    inline output_type operator()(double t) const { return valueAt(t); }

-    inline output_type valueAt(double t) const {

-        unsigned n = segN(t);

-        return segs[n](segT(t, n));

-    }

-    //TODO: maybe it is not a good idea to have this?

-    Piecewise<T> operator()(SBasis f);

-    Piecewise<T> operator()(Piecewise<SBasis>f);

-

-    inline unsigned size() const { return segs.size(); }

-    inline bool empty() const { return segs.empty(); }

-

-    /**Convenience/implementation hiding function to add segment/cut pairs.

-     * Asserts that basic size and order invariants are correct

-     */

-    inline void push(const T &s, double to) {

-        assert(cuts.size() - segs.size() == 1);

-        push_seg(s);

-        push_cut(to);

-    }

-    //Convenience/implementation hiding function to add cuts.

-    inline void push_cut(double c) {

-        assert_invariants(cuts.empty() || c > cuts.back()); 

-        cuts.push_back(c);

-    }

-    //Convenience/implementation hiding function to add segments.

-    inline void push_seg(const T &s) { segs.push_back(s); }

-

-    /**Returns the segment index which corresponds to a 'global' piecewise time.

-     * Also takes optional low/high parameters to expedite the search for the segment.

-     */

-    inline unsigned segN(double t, int low = 0, int high = -1) const {

-        high = (high == -1) ? size() : high;

-        if(t < cuts[0]) return 0;

-        if(t >= cuts[size()]) return size() - 1;

-        while(low < high) {

-            int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences

-            double mv = cuts[mid];

-            if(mv < t) {

-                if(t < cuts[mid + 1]) return mid; else low = mid + 1;

-            } else if(t < mv) {

-                if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;

-            } else {

-                return mid;

-            }

-        }

-        return low;

-    }

-

-    /**Returns the time within a segment, given the 'global' piecewise time.

-     * Also takes an optional index parameter which may be used for efficiency or to find the time on a

-     * segment outside its range.  If it is left to its default, -1, it will call segN to find the index.

-     */

-    inline double segT(double t, int i = -1) const {

-        if(i == -1) i = segN(t);

-        assert(i >= 0);

-        return (t - cuts[i]) / (cuts[i+1] - cuts[i]);

-    }

-

-    inline double mapToDomain(double t, unsigned i) const {

-        return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]

-    }    

-

-    //Offsets the piecewise domain

-    inline void offsetDomain(double o) {

-        assert(is_finite(o));

-        if(o != 0)

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

-                cuts[i] += o;

-    }

-

-    //Scales the domain of the function by a value. 0 will result in an empty Piecewise.

-    inline void scaleDomain(double s) {

-        assert(s > 0);

-        if(s == 0) {

-            cuts.clear(); segs.clear();

-            return;

-        }

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

-            cuts[i] *= s;

-    }

-

-    //Retrieves the domain in interval form

-    inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }

-

-    //Transforms the domain into another interval

-    inline void setDomain(Interval dom) {

-        if(empty()) return;

-        if(dom.isEmpty()) {

-            cuts.clear(); segs.clear();

-            return;

-        }

-        double cf = cuts.front();

-        double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);

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

-            cuts[i] = (cuts[i] - cf) * s + o;

-    }

-

-    //Concatenates this Piecewise function with another, offseting time of the other to match the end.

-    inline void concat(const Piecewise<T> &other) {

-        if(other.empty()) return;

-

-        if(empty()) { 

-            cuts = other.cuts; segs = other.segs;

-            return;

-        }

-

-        segs.insert(segs.end(), other.segs.begin(), other.segs.end());

-        double t = cuts.back() - other.cuts.front();        

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

-            push_cut(other.cuts[i + 1] + t);

-    }

-

-    //Like concat, but ensures continuity.

-    inline void continuousConcat(const Piecewise<T> &other) {

-        boost::function_requires<AddableConcept<typename T::output_type> >();

-        if(other.empty()) return;

-        typename T::output_type y = segs.back().at1() - other.segs.front().at0();

-

-        if(empty()) {

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

-                push_seg(other[i] + y);

-            cuts = other.cuts;

-            return;

-        }

-

-        double t = cuts.back() - other.cuts.front();

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

-            push(other[i] + y, other.cuts[i + 1] + t);

-    }

-

-    //returns true if the Piecewise<T> meets some basic invariants.

-    inline bool invariants() const {

-        // segs between cuts

-        if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))

-            return false;

-        // cuts in order

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

-            if(cuts[i] >= cuts[i+1])

-                return false;

-        return true;

-    }

-

-};

-

-template<typename T>

-inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {

-    boost::function_requires<FragmentConcept<T> >();

-    

-    if(f.empty()) return typename FragmentConcept<T>::BoundsType();

-    typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));

-    for(unsigned i = 1; i < f.size(); i++)

-        ret.unionWith(bounds_fast(f[i])); 

-    return ret;

-}

-

-template<typename T>

-inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {

-    boost::function_requires<FragmentConcept<T> >();

-    

-    if(f.empty()) return typename FragmentConcept<T>::BoundsType();

-    typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));

-    for(unsigned i = 1; i < f.size(); i++)

-        ret.unionWith(bounds_exact(f[i])); 

-    return ret;

-}

-

-template<typename T>

-inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {

-    boost::function_requires<FragmentConcept<T> >();

-    

-    if(f.empty()) return typename FragmentConcept<T>::BoundsType();

-    if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));

-    

-    unsigned fi = f.segN(m.min()), ti = f.segN(m.max());

-    double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);

-    

-    if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));

-    

-    typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));

-    for(unsigned i = fi + 1; i < ti; i++)

-        ret.unionWith(bounds_exact(f[i]));

-    if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));

-    

-    return ret;

-}

-

-//returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.

-template<typename T>

-T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {

-    assert(i < a.size());

-    double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);

-    return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );

-}

-

-/**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);

- * Further subdivides the Piecewise<T> such that there is a cut at every value in c.

- * Precondition: c sorted lower to higher.

- * 

- * //Given Piecewise<T> a and b:

- * Piecewise<T> ac = a.partition(b.cuts);

- * Piecewise<T> bc = b.partition(a.cuts);

- * //ac.cuts should be equivalent to bc.cuts

- */

-template<typename T>

-Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {

-    assert(pw.invariants());

-    if(c.empty()) return Piecewise<T>(pw);

-

-    Piecewise<T> ret = Piecewise<T>();

-    ret.cuts.reserve(c.size() + pw.cuts.size());

-    ret.segs.reserve(c.size() + pw.cuts.size() - 1);

-

-    if(pw.empty()) {

-        ret.cuts = c;

-        for(unsigned i = 0; i < c.size() - 1; i++)

-            ret.push_seg(T());

-        return ret;

-    }

-

-    unsigned si = 0, ci = 0;     //Segment index, Cut index

-

-    //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment

-    while(c[ci] < pw.cuts.front() && ci < c.size()) {

-        bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());

-        ret.push_cut(c[ci]);

-        ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );

-        ci++;

-    }

-

-    ret.push_cut(pw.cuts[0]);

-    double prev = pw.cuts[0];    //previous cut

-    //Loop which handles cuts within the Piecewise<T> domain

-    //Should have the cuts = segs + 1 invariant

-    while(si < pw.size() && ci <= c.size()) {

-        if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest

-            ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());

-            ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());

-            return ret;

-        }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) {  //no more cuts within this segment, finalize

-            if(prev > pw.cuts[si]) {      //segment already has cuts, so portion is required

-                ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));

-            } else {                     //plain copy is fine

-                ret.push_seg(pw[si]);

-            }

-            ret.push_cut(pw.cuts[si + 1]);

-            prev = pw.cuts[si + 1];

-            si++;

-        } else if(c[ci] == pw.cuts[si]){                  //coincident

-            //Already finalized the seg with the code immediately above

-            ci++;

-        } else {                                         //plain old subdivision

-            ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);

-            prev = c[ci];

-            ci++;

-        }

-    }

-    

-    //input cuts extend further than this Piecewise<T>, extend the last segment.

-    while(ci < c.size()) {

-        if(c[ci] > prev) {

-            ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);

-            prev = c[ci];

-        }

-        ci++;

-    }

-    return ret;

-}

-

-/**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);

- * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].

- */

-template<typename T>

-Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {

-    if(pw.empty() || from == to) return Piecewise<T>();

-

-    Piecewise<T> ret;

-

-    double temp = from;

-    from = std::min(from, to);

-    to = std::max(temp, to);

-    

-    unsigned i = pw.segN(from);

-    ret.push_cut(from);

-    if(i == pw.size() - 1 || to < pw.cuts[i + 1]) {    //to/from inhabit the same segment

-        ret.push(elem_portion(pw, i, from, to), to);

-        return ret;

-    }

-    ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));

-    i++;

-    unsigned fi = pw.segN(to, i);

-

-    ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi);  //copy segs

-    ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1);  //and their cuts

-

-    ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));

-    if(to != ret.cuts.back()) ret.push_cut(to);

-    ret.invariants();

-    return ret;

-}

-

-template<typename T>

-Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {

-    if(f.empty()) return f;

-    Piecewise<T> ret;

-    ret.push_cut(f.cuts[0]);

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

-        if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {

-            ret.push(f[i], f.cuts[i+1]);

-        }

-    }

-    return ret;

-}

-

-template<typename T>

-Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {

-    if(f.empty()) return f;

-    Piecewise<T> ret;

-    ret.push_cut(f.cuts[0]);

-    double last = f.cuts[0]; // last cut included

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

-        if (f.cuts[i+1]-f.cuts[i] >= tol) {

-            ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);

-            last = f.cuts[i+1];

-        }

-    }

-    return ret;

-}

-

-template<typename T>

-std::vector<double> roots(const Piecewise<T> &pw) {

-    std::vector<double> ret;

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

-        std::vector<double> sr = roots(pw[i]);

-        for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);

-

-    }

-    return ret;

-}

-

-//IMPL: OffsetableConcept

-template<typename T>

-Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {

-    boost::function_requires<OffsetableConcept<T> >();

-//TODO:empty

-    Piecewise<T> ret = Piecewise<T>();

-    ret.cuts = a.cuts;

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

-        ret.push_seg(a[i] + b);

-    return ret;

-}

-template<typename T>

-Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {

-    boost::function_requires<OffsetableConcept<T> >();

-//TODO: empty

-    Piecewise<T> ret = Piecewise<T>();

-    ret.cuts = a.cuts;

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

-        ret.push_seg(a[i] - b);

-    return ret;

-}

-template<typename T>

-Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {

-    boost::function_requires<OffsetableConcept<T> >();

-

-    if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }

-

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

-        a[i] += b;

-    return a;

-}

-template<typename T>

-Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {

-    boost::function_requires<OffsetableConcept<T> >();

-

-    if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }

-

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

-        a[i] -= b;

-    return a;

-}

-

-//IMPL: ScalableConcept

-template<typename T>

-Piecewise<T> operator-(Piecewise<T> const &a) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    Piecewise<T> ret;

-    ret.cuts = a.cuts;

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

-        ret.push_seg(- a[i]);

-    return ret;

-}

-template<typename T>

-Piecewise<T> operator*(Piecewise<T> const &a, double b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    if(a.empty()) return Piecewise<T>();

-

-    Piecewise<T> ret;

-    ret.cuts = a.cuts;

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

-        ret.push_seg(a[i] * b);

-    return ret;

-}

-template<typename T>

-Piecewise<T> operator/(Piecewise<T> const &a, double b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    //FIXME: b == 0?

-    if(a.empty()) return Piecewise<T>();

-

-    Piecewise<T> ret;

-    ret.cuts = a.cuts;

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

-        ret.push_seg(a[i] / b);

-    return ret;

-}

-template<typename T>

-Piecewise<T> operator*=(Piecewise<T>& a, double b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    if(a.empty()) return Piecewise<T>();

-

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

-        a[i] *= b;

-    return a;

-}

-template<typename T>

-Piecewise<T> operator/=(Piecewise<T>& a, double b) {

-    boost::function_requires<ScalableConcept<T> >();

-

-    //FIXME: b == 0?

-    if(a.empty()) return Piecewise<T>();

-

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

-        a[i] /= b;

-    return a;

-}

-

-//IMPL: AddableConcept

-template<typename T>

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

-    boost::function_requires<AddableConcept<T> >();

-

-    Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);

-    Piecewise<T> ret = Piecewise<T>();

-    assert(pa.size() == pb.size());

-    ret.cuts = pa.cuts;

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

-        ret.push_seg(pa[i] + pb[i]);

-    return ret;

-}

-template<typename T>

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

-    boost::function_requires<AddableConcept<T> >();

-

-    Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);

-    Piecewise<T> ret = Piecewise<T>();

-    assert(pa.size() == pb.size());

-    ret.cuts = pa.cuts;

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

-        ret.push_seg(pa[i] - pb[i]);

-    return ret;

-}

-template<typename T>

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

-    a = a+b;

-    return a;

-}

-template<typename T>

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

-    a = a-b;

-    return a;

-}

-

-template<typename T1,typename T2>

-Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {

-    //function_requires<MultiplicableConcept<T1> >();

-    //function_requires<MultiplicableConcept<T2> >();

-

-    Piecewise<T1> pa = partition(a, b.cuts);

-    Piecewise<T2> pb = partition(b, a.cuts);

-    Piecewise<T2> ret = Piecewise<T2>();

-    assert(pa.size() == pb.size());

-    ret.cuts = pa.cuts;

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

-        ret.push_seg(pa[i] * pb[i]);

-    return ret;

-}

-

-template<typename T>

-inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {

-    a = a * b;

-    return a;

-}

-

-Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);

-//TODO: replace divide(a,b,k) by divide(a,b,tol,k)?

-//TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.

-//Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero). 

-Piecewise<SBasis> 

-divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);

-Piecewise<SBasis> 

-divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);

-Piecewise<SBasis> 

-divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);

-Piecewise<SBasis> 

-divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);

-

-//Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.

-std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);

-int compose_findSegIdx(std::map<double,unsigned>::iterator  const &cut,

-                       std::map<double,unsigned>::iterator  const &next,

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

-                       SBasis const &g);

-

-//TODO: add concept check

-template<typename T>

-Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){

-    Piecewise<T> result;

-    if (f.empty()) return result;

-    if (g.isZero()) return Piecewise<T>(f(0));

-    if (f.size()==1){

-        double t0 = f.cuts[0], width = f.cuts[1] - t0;

-        return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));

-    }

-    

-    //first check bounds...

-    Interval bs = bounds_fast(g);

-    if (f.cuts.front() > bs.max()  || bs.min() > f.cuts.back()){

-        int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;

-        double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;

-        return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));  

-    }

-    

-    std::vector<double> levels;//we can forget first and last cuts...

-    levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);

-    //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.

-    std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);

-    

-    //-- Compose each piece of g with the relevant seg of f.

-    result.cuts.push_back(0.);

-    std::map<double,unsigned>::iterator cut=cuts_pb.begin();

-    std::map<double,unsigned>::iterator next=cut; next++;

-    while(next!=cuts_pb.end()){

-        //assert(std::abs(int((*cut).second-(*next).second))<1);

-        //TODO: find a way to recover from this error? the root finder missed some root;

-        //  the levels/variations of f might be too close/fast...

-        int idx = compose_findSegIdx(cut,next,levels,g);

-        double t0=(*cut).first;

-        double t1=(*next).first;

-        

-        SBasis sub_g=compose(g, Linear(t0,t1));

-        sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),

-                             (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);

-        result.push(compose(f[idx],sub_g),t1);

-        cut++;

-        next++;

-    }

-    return(result);

-} 

-

-//TODO: add concept check for following composition functions

-template<typename T>

-Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){

-  Piecewise<T> result;

-  for(unsigned i = 0; i < g.segs.size(); i++){

-      Piecewise<T> fgi=compose(f, g.segs[i]);

-      fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));

-      result.concat(fgi);

-  }

-  return result;

-}

-

-template <typename T>

-Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}

-template <typename T>

-Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}

-

-template<typename T>

-Piecewise<T> integral(Piecewise<T> const &a) {

-    Piecewise<T> result;

-    result.segs.resize(a.segs.size());

-    result.cuts = a.cuts;

-    typename T::output_type c = a.segs[0].at0();

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

-        result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);

-        result.segs[i]+= c-result.segs[i].at0();

-        c = result.segs[i].at1();

-    }

-    return result;

-}

-

-template<typename T>

-Piecewise<T> derivative(Piecewise<T> const &a) {

-    Piecewise<T> result;

-    result.segs.resize(a.segs.size());

-    result.cuts = a.cuts;

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

-        result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);

-    }

-    return result;

-}

-

-std::vector<double> roots(Piecewise<SBasis> const &f);

-

-}

-

-#endif //SEEN_GEOM_PW_SB_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:encoding=utf-8:textwidth=99 :

+/*
+ * piecewise.h - Piecewise function class
+ *
+ * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
+ *
+ * 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, output 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_GEOM_PW_SB_H
+#define SEEN_GEOM_PW_SB_H
+
+#include "sbasis.h"
+#include <vector>
+#include <map>
+
+#include "concepts.h"
+#include "isnan.h"
+#include <boost/concept_check.hpp>
+
+namespace Geom {
+
+template <typename T>
+class Piecewise {
+  BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);
+
+  public:
+    std::vector<double> cuts;
+    std::vector<T> segs;
+    //segs[i] stretches from cuts[i] to cuts[i+1].
+
+    Piecewise() {}
+
+    explicit Piecewise(const T &s) {
+        push_cut(0.);
+        push_seg(s);
+        push_cut(1.);
+    }
+
+    typedef typename T::output_type output_type;
+
+    explicit Piecewise(const output_type & v) {
+        push_cut(0.);
+        push_seg(T(v));
+        push_cut(1.);
+    }
+
+    inline T operator[](unsigned i) const { return segs[i]; }
+    inline T &operator[](unsigned i) { return segs[i]; }
+    inline output_type operator()(double t) const { return valueAt(t); }
+    inline output_type valueAt(double t) const {
+        unsigned n = segN(t);
+        return segs[n](segT(t, n));
+    }
+    //TODO: maybe it is not a good idea to have this?
+    Piecewise<T> operator()(SBasis f);
+    Piecewise<T> operator()(Piecewise<SBasis>f);
+
+    inline unsigned size() const { return segs.size(); }
+    inline bool empty() const { return segs.empty(); }
+
+    /**Convenience/implementation hiding function to add segment/cut pairs.
+     * Asserts that basic size and order invariants are correct
+     */
+    inline void push(const T &s, double to) {
+        assert(cuts.size() - segs.size() == 1);
+        push_seg(s);
+        push_cut(to);
+    }
+    //Convenience/implementation hiding function to add cuts.
+    inline void push_cut(double c) {
+        assert_invariants(cuts.empty() || c > cuts.back());
+        cuts.push_back(c);
+    }
+    //Convenience/implementation hiding function to add segments.
+    inline void push_seg(const T &s) { segs.push_back(s); }
+
+    /**Returns the segment index which corresponds to a 'global' piecewise time.
+     * Also takes optional low/high parameters to expedite the search for the segment.
+     */
+    inline unsigned segN(double t, int low = 0, int high = -1) const {
+        high = (high == -1) ? size() : high;
+        if(t < cuts[0]) return 0;
+        if(t >= cuts[size()]) return size() - 1;
+        while(low < high) {
+            int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences
+            double mv = cuts[mid];
+            if(mv < t) {
+                if(t < cuts[mid + 1]) return mid; else low = mid + 1;
+            } else if(t < mv) {
+                if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;
+            } else {
+                return mid;
+            }
+        }
+        return low;
+    }
+
+    /**Returns the time within a segment, given the 'global' piecewise time.
+     * Also takes an optional index parameter which may be used for efficiency or to find the time on a
+     * segment outside its range.  If it is left to its default, -1, it will call segN to find the index.
+     */
+    inline double segT(double t, int i = -1) const {
+        if(i == -1) i = segN(t);
+        assert(i >= 0);
+        return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
+    }
+
+    inline double mapToDomain(double t, unsigned i) const {
+        return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]
+    }
+
+    //Offsets the piecewise domain
+    inline void offsetDomain(double o) {
+        assert(is_finite(o));
+        if(o != 0)
+            for(unsigned i = 0; i <= size(); i++)
+                cuts[i] += o;
+    }
+
+    //Scales the domain of the function by a value. 0 will result in an empty Piecewise.
+    inline void scaleDomain(double s) {
+        assert(s > 0);
+        if(s == 0) {
+            cuts.clear(); segs.clear();
+            return;
+        }
+        for(unsigned i = 0; i <= size(); i++)
+            cuts[i] *= s;
+    }
+
+    //Retrieves the domain in interval form
+    inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
+
+    //Transforms the domain into another interval
+    inline void setDomain(Interval dom) {
+        if(empty()) return;
+        if(dom.isEmpty()) {
+            cuts.clear(); segs.clear();
+            return;
+        }
+        double cf = cuts.front();
+        double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
+        for(unsigned i = 0; i <= size(); i++)
+            cuts[i] = (cuts[i] - cf) * s + o;
+    }
+
+    //Concatenates this Piecewise function with another, offseting time of the other to match the end.
+    inline void concat(const Piecewise<T> &other) {
+        if(other.empty()) return;
+
+        if(empty()) {
+            cuts = other.cuts; segs = other.segs;
+            return;
+        }
+
+        segs.insert(segs.end(), other.segs.begin(), other.segs.end());
+        double t = cuts.back() - other.cuts.front();
+        for(unsigned i = 0; i < other.size(); i++)
+            push_cut(other.cuts[i + 1] + t);
+    }
+
+    //Like concat, but ensures continuity.
+    inline void continuousConcat(const Piecewise<T> &other) {
+        boost::function_requires<AddableConcept<typename T::output_type> >();
+        if(other.empty()) return;
+        typename T::output_type y = segs.back().at1() - other.segs.front().at0();
+
+        if(empty()) {
+            for(unsigned i = 0; i < other.size(); i++)
+                push_seg(other[i] + y);
+            cuts = other.cuts;
+            return;
+        }
+
+        double t = cuts.back() - other.cuts.front();
+        for(unsigned i = 0; i < other.size(); i++)
+            push(other[i] + y, other.cuts[i + 1] + t);
+    }
+
+    //returns true if the Piecewise<T> meets some basic invariants.
+    inline bool invariants() const {
+        // segs between cuts
+        if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
+            return false;
+        // cuts in order
+        for(unsigned i = 0; i < segs.size(); i++)
+            if(cuts[i] >= cuts[i+1])
+                return false;
+        return true;
+    }
+
+};
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {
+    boost::function_requires<FragmentConcept<T> >();
+
+    if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+    typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
+    for(unsigned i = 1; i < f.size(); i++)
+        ret.unionWith(bounds_fast(f[i]));
+    return ret;
+}
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {
+    boost::function_requires<FragmentConcept<T> >();
+
+    if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+    typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
+    for(unsigned i = 1; i < f.size(); i++)
+        ret.unionWith(bounds_exact(f[i]));
+    return ret;
+}
+
+template<typename T>
+inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {
+    boost::function_requires<FragmentConcept<T> >();
+
+    if(f.empty()) return typename FragmentConcept<T>::BoundsType();
+    if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
+
+    unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
+    double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
+
+    if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));
+
+    typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
+    for(unsigned i = fi + 1; i < ti; i++)
+        ret.unionWith(bounds_exact(f[i]));
+    if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
+
+    return ret;
+}
+
+//returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.
+template<typename T>
+T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {
+    assert(i < a.size());
+    double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
+    return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
+}
+
+/**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);
+ * Further subdivides the Piecewise<T> such that there is a cut at every value in c.
+ * Precondition: c sorted lower to higher.
+ *
+ * //Given Piecewise<T> a and b:
+ * Piecewise<T> ac = a.partition(b.cuts);
+ * Piecewise<T> bc = b.partition(a.cuts);
+ * //ac.cuts should be equivalent to bc.cuts
+ */
+template<typename T>
+Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {
+    assert(pw.invariants());
+    if(c.empty()) return Piecewise<T>(pw);
+
+    Piecewise<T> ret = Piecewise<T>();
+    ret.cuts.reserve(c.size() + pw.cuts.size());
+    ret.segs.reserve(c.size() + pw.cuts.size() - 1);
+
+    if(pw.empty()) {
+        ret.cuts = c;
+        for(unsigned i = 0; i < c.size() - 1; i++)
+            ret.push_seg(T());
+        return ret;
+    }
+
+    unsigned si = 0, ci = 0;     //Segment index, Cut index
+
+    //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment
+    while(c[ci] < pw.cuts.front() && ci < c.size()) {
+        bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
+        ret.push_cut(c[ci]);
+        ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
+        ci++;
+    }
+
+    ret.push_cut(pw.cuts[0]);
+    double prev = pw.cuts[0];    //previous cut
+    //Loop which handles cuts within the Piecewise<T> domain
+    //Should have the cuts = segs + 1 invariant
+    while(si < pw.size() && ci <= c.size()) {
+        if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest
+            ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
+            ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
+            return ret;
+        }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) {  //no more cuts within this segment, finalize
+            if(prev > pw.cuts[si]) {      //segment already has cuts, so portion is required
+                ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));
+            } else {                     //plain copy is fine
+                ret.push_seg(pw[si]);
+            }
+            ret.push_cut(pw.cuts[si + 1]);
+            prev = pw.cuts[si + 1];
+            si++;
+        } else if(c[ci] == pw.cuts[si]){                  //coincident
+            //Already finalized the seg with the code immediately above
+            ci++;
+        } else {                                         //plain old subdivision
+            ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);
+            prev = c[ci];
+            ci++;
+        }
+    }
+
+    //input cuts extend further than this Piecewise<T>, extend the last segment.
+    while(ci < c.size()) {
+        if(c[ci] > prev) {
+            ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
+            prev = c[ci];
+        }
+        ci++;
+    }
+    return ret;
+}
+
+/**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);
+ * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].
+ */
+template<typename T>
+Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {
+    if(pw.empty() || from == to) return Piecewise<T>();
+
+    Piecewise<T> ret;
+
+    double temp = from;
+    from = std::min(from, to);
+    to = std::max(temp, to);
+
+    unsigned i = pw.segN(from);
+    ret.push_cut(from);
+    if(i == pw.size() - 1 || to < pw.cuts[i + 1]) {    //to/from inhabit the same segment
+        ret.push(elem_portion(pw, i, from, to), to);
+        return ret;
+    }
+    ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));
+    i++;
+    unsigned fi = pw.segN(to, i);
+
+    ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi);  //copy segs
+    ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1);  //and their cuts
+
+    ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));
+    if(to != ret.cuts.back()) ret.push_cut(to);
+    ret.invariants();
+    return ret;
+}
+
+template<typename T>
+Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {
+    if(f.empty()) return f;
+    Piecewise<T> ret;
+    ret.push_cut(f.cuts[0]);
+    for(unsigned i=0; i<f.size(); i++){
+        if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
+            ret.push(f[i], f.cuts[i+1]);
+        }
+    }
+    return ret;
+}
+
+template<typename T>
+Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {
+    if(f.empty()) return f;
+    Piecewise<T> ret;
+    ret.push_cut(f.cuts[0]);
+    double last = f.cuts[0]; // last cut included
+    for(unsigned i=0; i<f.size(); i++){
+        if (f.cuts[i+1]-f.cuts[i] >= tol) {
+            ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
+            last = f.cuts[i+1];
+        }
+    }
+    return ret;
+}
+
+template<typename T>
+std::vector<double> roots(const Piecewise<T> &pw) {
+    std::vector<double> ret;
+    for(unsigned i = 0; i < pw.size(); i++) {
+        std::vector<double> sr = roots(pw[i]);
+        for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
+
+    }
+    return ret;
+}
+
+//IMPL: OffsetableConcept
+template<typename T>
+Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {
+    boost::function_requires<OffsetableConcept<T> >();
+//TODO:empty
+    Piecewise<T> ret = Piecewise<T>();
+    ret.cuts = a.cuts;
+    for(unsigned i = 0; i < a.size();i++)
+        ret.push_seg(a[i] + b);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {
+    boost::function_requires<OffsetableConcept<T> >();
+//TODO: empty
+    Piecewise<T> ret = Piecewise<T>();
+    ret.cuts = a.cuts;
+    for(unsigned i = 0; i < a.size();i++)
+        ret.push_seg(a[i] - b);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {
+    boost::function_requires<OffsetableConcept<T> >();
+
+    if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
+
+    for(unsigned i = 0; i < a.size();i++)
+        a[i] += b;
+    return a;
+}
+template<typename T>
+Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {
+    boost::function_requires<OffsetableConcept<T> >();
+
+    if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
+
+    for(unsigned i = 0;i < a.size();i++)
+        a[i] -= b;
+    return a;
+}
+
+//IMPL: ScalableConcept
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    Piecewise<T> ret;
+    ret.cuts = a.cuts;
+    for(unsigned i = 0; i < a.size();i++)
+        ret.push_seg(- a[i]);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator*(Piecewise<T> const &a, double b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    if(a.empty()) return Piecewise<T>();
+
+    Piecewise<T> ret;
+    ret.cuts = a.cuts;
+    for(unsigned i = 0; i < a.size();i++)
+        ret.push_seg(a[i] * b);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator/(Piecewise<T> const &a, double b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    //FIXME: b == 0?
+    if(a.empty()) return Piecewise<T>();
+
+    Piecewise<T> ret;
+    ret.cuts = a.cuts;
+    for(unsigned i = 0; i < a.size();i++)
+        ret.push_seg(a[i] / b);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator*=(Piecewise<T>& a, double b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    if(a.empty()) return Piecewise<T>();
+
+    for(unsigned i = 0; i < a.size();i++)
+        a[i] *= b;
+    return a;
+}
+template<typename T>
+Piecewise<T> operator/=(Piecewise<T>& a, double b) {
+    boost::function_requires<ScalableConcept<T> >();
+
+    //FIXME: b == 0?
+    if(a.empty()) return Piecewise<T>();
+
+    for(unsigned i = 0; i < a.size();i++)
+        a[i] /= b;
+    return a;
+}
+
+//IMPL: AddableConcept
+template<typename T>
+Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
+    Piecewise<T> ret = Piecewise<T>();
+    assert(pa.size() == pb.size());
+    ret.cuts = pa.cuts;
+    for (unsigned i = 0; i < pa.size(); i++)
+        ret.push_seg(pa[i] + pb[i]);
+    return ret;
+}
+template<typename T>
+Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {
+    boost::function_requires<AddableConcept<T> >();
+
+    Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
+    Piecewise<T> ret = Piecewise<T>();
+    assert(pa.size() == pb.size());
+    ret.cuts = pa.cuts;
+    for (unsigned i = 0; i < pa.size(); i++)
+        ret.push_seg(pa[i] - pb[i]);
+    return ret;
+}
+template<typename T>
+inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T> const &b) {
+    a = a+b;
+    return a;
+}
+template<typename T>
+inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T> const &b) {
+    a = a-b;
+    return a;
+}
+
+template<typename T1,typename T2>
+Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {
+    //function_requires<MultiplicableConcept<T1> >();
+    //function_requires<MultiplicableConcept<T2> >();
+
+    Piecewise<T1> pa = partition(a, b.cuts);
+    Piecewise<T2> pb = partition(b, a.cuts);
+    Piecewise<T2> ret = Piecewise<T2>();
+    assert(pa.size() == pb.size());
+    ret.cuts = pa.cuts;
+    for (unsigned i = 0; i < pa.size(); i++)
+        ret.push_seg(pa[i] * pb[i]);
+    return ret;
+}
+
+template<typename T>
+inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {
+    a = a * b;
+    return a;
+}
+
+Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);
+//TODO: replace divide(a,b,k) by divide(a,b,tol,k)?
+//TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.
+//Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero).
+Piecewise<SBasis>
+divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
+Piecewise<SBasis>
+divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
+
+//Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.
+std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);
+int compose_findSegIdx(std::map<double,unsigned>::iterator  const &cut,
+                       std::map<double,unsigned>::iterator  const &next,
+                       std::vector<double>  const &levels,
+                       SBasis const &g);
+
+//TODO: add concept check
+template<typename T>
+Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
+    Piecewise<T> result;
+    if (f.empty()) return result;
+    if (g.isZero()) return Piecewise<T>(f(0));
+    if (f.size()==1){
+        double t0 = f.cuts[0], width = f.cuts[1] - t0;
+        return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
+    }
+
+    //first check bounds...
+    Interval bs = bounds_fast(g);
+    if (f.cuts.front() > bs.max()  || bs.min() > f.cuts.back()){
+        int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
+        double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
+        return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
+    }
+
+    std::vector<double> levels;//we can forget first and last cuts...
+    levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
+    //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.
+    std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
+
+    //-- Compose each piece of g with the relevant seg of f.
+    result.cuts.push_back(0.);
+    std::map<double,unsigned>::iterator cut=cuts_pb.begin();
+    std::map<double,unsigned>::iterator next=cut; next++;
+    while(next!=cuts_pb.end()){
+        //assert(std::abs(int((*cut).second-(*next).second))<1);
+        //TODO: find a way to recover from this error? the root finder missed some root;
+        //  the levels/variations of f might be too close/fast...
+        int idx = compose_findSegIdx(cut,next,levels,g);
+        double t0=(*cut).first;
+        double t1=(*next).first;
+
+        SBasis sub_g=compose(g, Linear(t0,t1));
+        sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
+                             (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
+        result.push(compose(f[idx],sub_g),t1);
+        cut++;
+        next++;
+    }
+    return(result);
+}
+
+//TODO: add concept check for following composition functions
+template<typename T>
+Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){
+  Piecewise<T> result;
+  for(unsigned i = 0; i < g.segs.size(); i++){
+      Piecewise<T> fgi=compose(f, g.segs[i]);
+      fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
+      result.concat(fgi);
+  }
+  return result;
+}
+
+template <typename T>
+Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
+template <typename T>
+Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
+
+template<typename T>
+Piecewise<T> integral(Piecewise<T> const &a) {
+    Piecewise<T> result;
+    result.segs.resize(a.segs.size());
+    result.cuts = a.cuts;
+    typename T::output_type c = a.segs[0].at0();
+    for(unsigned i = 0; i < a.segs.size(); i++){
+        result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
+        result.segs[i]+= c-result.segs[i].at0();
+        c = result.segs[i].at1();
+    }
+    return result;
+}
+
+template<typename T>
+Piecewise<T> derivative(Piecewise<T> const &a) {
+    Piecewise<T> result;
+    result.segs.resize(a.segs.size());
+    result.cuts = a.cuts;
+    for(unsigned i = 0; i < a.segs.size(); i++){
+        result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
+    }
+    return result;
+}
+
+std::vector<double> roots(Piecewise<SBasis> const &f);
+
+}
+
+#endif //SEEN_GEOM_PW_SB_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:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/poly-dk-solve.cpp b/src/2geom/poly-dk-solve.cpp
index 136687ae8..fdc1cefe5 100644
--- a/src/2geom/poly-dk-solve.cpp
+++ b/src/2geom/poly-dk-solve.cpp
@@ -1,64 +1,64 @@
-#include "poly-dk-solve.h"

-#include <iterator>

-

-/*** implementation of the Durand-Kerner method.  seems buggy*/

-

-std::complex<double> evalu(Poly const & p, std::complex<double> x) {

-    std::complex<double> result = 0;

-    std::complex<double> xx = 1;

-    

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

-        result += p[i]*xx;

-        xx *= x;

-    }

-    return result;

-}

-

-std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {

-    std::vector<std::complex<double> > roots;

-    const int N = ply.degree();

-    

-    std::complex<double> b(0.4, 0.9);

-    std::complex<double> p = 1;

-    for(int i = 0; i < N; i++) {

-        roots.push_back(p);

-        p *= b;

-    }

-    assert(roots.size() == ply.degree());

-    

-    double error = 0;

-    int i;

-    for( i = 0; i < 30; i++) {

-        error = 0;

-        for(int r_i = 0; r_i < N; r_i++) {

-            std::complex<double> denom = 1;

-            std::complex<double> R = roots[r_i];

-            for(int d_i = 0; d_i < N; d_i++) {

-                if(r_i != d_i)

-                    denom *= R-roots[d_i];

-            }

-            assert(norm(denom) != 0);

-	    std::complex<double> dr = evalu(ply, R)/denom;

-            error += norm(dr);

-            roots[r_i] = R - dr;

-        }

-        /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));

-          std::cout << std::endl;*/

-        if(error < tol)

-            break;

-    }

-    //std::cout << error << ", " << i<< std::endl;

-    return roots;

-}

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

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

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

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

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

-#include <iterator>

-

-typedef std::complex<double> cdouble;

-

-cdouble laguerre_internal_complex(Poly const & p,

-                                  double x0,

-                                  double tol,

-                                  bool & quad_root) {

-    cdouble a = 2*tol;

-    cdouble xk = x0;

-    double n = p.degree();

-    quad_root = false;

-    const unsigned shuffle_rate = 10;

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

-    unsigned shuffle_counter = 0;

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

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

-        cdouble b = p.back();

-        cdouble d = 0, f = 0;

-        double err = abs(b);

-        double abx = abs(xk);

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

-            f = xk*f + d;

-            d = xk*d + b;

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

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

-        }

-        

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

-        

-        cdouble px = b;

-        if(abs(b) < err)

-            return xk;

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

-        //    return xk;

-        cdouble G = d / px;

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

-        

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

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

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

-        if(radicand.real() < 0)

-            quad_root = true;

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

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

-            a = - sqrt(radicand);

-        else

-            a = sqrt(radicand);

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

-        a = n / (a + G);

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

-        if(shuffle_counter % shuffle_rate == 0)

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

-        xk -= a;

-        shuffle_counter++;

-        if(shuffle_counter >= 90)

-            break;

-    }

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

-    return xk;

-}

-

-double laguerre_internal(Poly const & p,

-                         Poly const & pp,

-                         Poly const & ppp,

-                         double x0,

-                         double tol,

-                         bool & quad_root) {

-    double a = 2*tol;

-    double xk = x0;

-    double n = p.degree();

-    quad_root = false;

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

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

-        double px = p(xk);

-        if(px*px < tol*tol)

-            return xk;

-        double G = pp(xk) / px;

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

-        

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

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

-        assert(radicand > 0);

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

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

-            a = - sqrt(radicand);

-        else

-            a = sqrt(radicand);

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

-        a = n / (a + G);

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

-        xk -= a;

-    }

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

-    return xk;

-}

-

-

-std::vector<cdouble > 

-laguerre(Poly p, const double tol) {

-    std::vector<cdouble > solutions;

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

-    while(p.size() > 1)

-    {

-        double x0 = 0;

-        bool quad_root = false;

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

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

-        Poly dvs;

-        if(quad_root) {

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

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

-            dvs.push_back(1.0);

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

-	    //solutions.push_back(sol);

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

-        } else {

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

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

-            dvs.push_back(1.0);

-	    solutions.push_back(sol);

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

-        }

-        Poly r;

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

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

-    }

-    return solutions;

-}

-

-std::vector<double> 

-laguerre_real_interval(Poly const & ply, 

-		       const double lo, const double hi,

-		       const double tol) {

-}

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

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

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

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

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

+#include <vector>
 #include <cassert>
 
 class Quad{
@@ -20,7 +20,7 @@ public:
     double by0, by1;
 
     QuadTree() : root(0), scale(1) {}
-    
+
     Quad* search(double x0, double y0, double x1, double y1);
     void insert(double x0, double y0, double x1, double y1, int shape);
     void erase(Quad *q, int shape);
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index 5bf0d2876..7157bc885 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -1,490 +1,490 @@
-/*

- *  sbasis.cpp - S-power basis function class + supporting classes

- *

- *  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.

- */

-

-#include <cmath>

-

-#include "sbasis.h"

-#include "isnan.h"

-

-namespace Geom{

-

-/*** At some point we should work on tighter bounds for the error.  It is clear that the error is

- * bounded by the L1 norm over the tail of the series, but this is very loose, leading to far too

- * many cubic beziers.  I've changed this to be \sum _i=tail ^\infty |hat a_i| 2^-i but I have no

- * evidence that this is correct.

- */

-

-/*

-double SBasis::tail_error(unsigned tail) const {

-    double err = 0, s = 1./(1<<(2*tail)); // rough

-    for(unsigned i = tail; i < size(); i++) {

-        err += (fabs((*this)[i][0]) + fabs((*this)[i][1]))*s;

-        s /= 4;

-    }

-    return err;

-}

+/*
+ *  sbasis.cpp - S-power basis function class + supporting classes
+ *
+ *  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.
+ */
+
+#include <cmath>
+
+#include "sbasis.h"
+#include "isnan.h"
+
+namespace Geom{
+
+/*** At some point we should work on tighter bounds for the error.  It is clear that the error is
+ * bounded by the L1 norm over the tail of the series, but this is very loose, leading to far too
+ * many cubic beziers.  I've changed this to be \sum _i=tail ^\infty |hat a_i| 2^-i but I have no
+ * evidence that this is correct.
+ */
+
+/*
+double SBasis::tail_error(unsigned tail) const {
+    double err = 0, s = 1./(1<<(2*tail)); // rough
+    for(unsigned i = tail; i < size(); i++) {
+        err += (fabs((*this)[i][0]) + fabs((*this)[i][1]))*s;
+        s /= 4;
+    }
+    return err;
+}
 */
-

-double SBasis::tailError(unsigned tail) const {

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

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

-}

-

-bool SBasis::isFinite() const {

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

-        if(!(*this)[i].isFinite())

-            return false;

-    }

-    return true;

-}

-

-SBasis operator+(const SBasis& a, const SBasis& b) {

-    SBasis result;

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

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

-    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;

-}

-

-SBasis operator-(const SBasis& a, const SBasis& b) {

-    SBasis result;

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

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

-    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;

-}

-

-SBasis& operator+=(SBasis& a, const SBasis& b) {

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

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

-    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;

-}

-

-SBasis& operator-=(SBasis& a, const SBasis& b) {

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

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

-    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;

-}

-

-SBasis operator*(SBasis const &a, double k) {

-    SBasis c;

-    c.reserve(a.size());

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

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

-    return c;

-}

-

-SBasis& operator*=(SBasis& 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;

-}

-

-SBasis shift(SBasis const &a, int sh) {

-    SBasis c = a;

-    if(sh > 0) {

-        c.insert(c.begin(), sh, Linear(0,0));

-    } else {

-        //TODO: truncate

-    }

-    return c;

-}

-

-SBasis shift(Linear const &a, int sh) {

-    SBasis c;

-    if(sh > 0) {

-        c.insert(c.begin(), sh, Linear(0,0));

-        c.push_back(a);

-    }

-    return c;

-}

-

-SBasis multiply(SBasis const &a, SBasis const &b) {

-    // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}

-    

-    // shift(1, a.Tri*b.Tri)

-    SBasis c;

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

-        return c;

-    c.resize(a.size() + b.size(), Linear(0,0));

-    c[0] = Linear(0,0);

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

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

-            double tri = Tri(b[j])*Tri(a[i-j]);

-            c[i+1/*shift*/] += Linear(Hat(-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;

-}

-

-SBasis integral(SBasis const &c) {

-    SBasis a;

-    a.resize(c.size() + 1, Linear(0,0));

-    a[0] = Linear(0,0);

-    

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

-        double ahat = -Tri(c[k-1])/(2*k);

-        a[k] = Hat(ahat);

-    }

-    double aTri = 0;

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

-        aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);

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

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

-    }

-    a.normalize();

-    return a;

-}

-

-SBasis derivative(SBasis const &a) {

-    SBasis c;

-    c.resize(a.size(), Linear(0,0));

-    

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

-        double d = (2*k+1)*Tri(a[k]);

-        

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

-            c[k][dim] = d;

-            if(k+1 < a.size()) {

-                if(dim)

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

-                else

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

-            }

-        }

-    }

-    

-    return c;

-}

-

-//TODO: convert int k to unsigned k, and remove cast

-SBasis sqrt(SBasis const &a, int k) {

-    SBasis c;

-    if(a.isZero() || k == 0)

-        return c;

-    c.resize(k, Linear(0,0));

-    c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));

-    SBasis r = a - multiply(c, c); // remainder

-    

-    for(unsigned i = 1; i <= (unsigned)k and i<r.size(); i++) {

-        Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));

-        SBasis cisi = shift(ci, i);

-        r -= multiply(shift((c*2 + cisi), i), SBasis(ci));

-        r.truncate(k+1);

-        c += cisi;

-        if(r.tailError(i) == 0) // if exact

-            break;

-    }

-    

-    return c;

-}

-

-// return a kth order approx to 1/a)

-SBasis reciprocal(Linear const &a, int k) {

-    SBasis c;

-    assert(!a.isZero());

-    c.resize(k, Linear(0,0));

-    double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);

-    double r_s0k = 1;

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

-        c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);

-        r_s0k *= r_s0;

-    }

-    return c;

-}

-

-SBasis divide(SBasis const &a, SBasis const &b, int k) {

-    SBasis c;

-    assert(!a.isZero());

-    SBasis r = a; // remainder

-    

-    k++;

-    r.resize(k, Linear(0,0));

-    c.resize(k, Linear(0,0));

-

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

-        Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0

-        c[i] += ci;

-        r -= shift(multiply(ci,b), i);

-        r.truncate(k+1);

-        if(r.tailError(i) == 0) // if exact

-            break;

-    }

-    

-    return c;

-}

-

-// a(b)

-// return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k

-SBasis compose(SBasis const &a, SBasis const &b) {

-    SBasis s = multiply((SBasis(Linear(1,1))-b), b);

-    SBasis r;

-    

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

-        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);

-    }

-    return r;

-}

-

-// a(b)

-// return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k

-SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {

-    SBasis s = multiply((SBasis(Linear(1,1))-b), b);

-    SBasis r;

-    

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

-        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);

-    }

-    r.truncate(k);

-    return r;

-}

-

-/*

-Inversion algorithm. The notation is certainly very misleading. The

-pseudocode should say:

-

-c(v) := 0

-r(u) := r_0(u) := u

-for i:=0 to k do

-  c_i(v) := H_0(r_i(u)/(t_1)^i; u)

-  c(v) := c(v) + c_i(v)*t^i

-  r(u) := r(u) ? c_i(u)*(t(u))^i

-endfor

-*/

-

-//#define DEBUG_INVERSION 1 

-

-SBasis inverse(SBasis a, int k) {

-    assert(a.size() > 0);

-// the function should have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.

-    double a0 = a[0][0];

-    if(a0 != 0) {

-        a -= a0;

-    }

-    double a1 = a[0][1];

-    assert(a1 != 0); // not invertable.

-    

-    if(a1 != 1) {

-        a /= a1;

-    }

-    SBasis c;                           // c(v) := 0

-    if(a.size() >= 2 && k == 2) {

-        c.push_back(Linear(0,1));

-        Linear t1(1+a[1][0], 1-a[1][1]);    // t_1

-        c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));

-    } else if(a.size() >= 2) {                      // non linear

-        SBasis r = Linear(0,1);             // r(u) := r_0(u) := u

-        Linear t1(1./(1+a[1][0]), 1./(1-a[1][1]));    // 1./t_1

-        Linear one(1,1);

-        Linear t1i = one;                   // t_1^0

-        SBasis one_minus_a = SBasis(one) - a;

-        SBasis t = multiply(one_minus_a, a); // t(u)

-        SBasis ti(one);                     // t(u)^0

-#ifdef DEBUG_INVERSION

-        std::cout << "a=" << a << std::endl;

-        std::cout << "1-a=" << one_minus_a << std::endl;

-        std::cout << "t1=" << t1 << std::endl;

-        //assert(t1 == t[1]);

-#endif

-    

-        c.resize(k+1, Linear(0,0));

-        for(unsigned i = 0; i < (unsigned)k; i++) {   // for i:=0 to k do

-#ifdef DEBUG_INVERSION

-            std::cout << "-------" << i << ": ---------" <<std::endl;

-            std::cout << "r=" << r << std::endl

-                      << "c=" << c << std::endl

-                      << "ti=" << ti << std::endl

-                      << std::endl;

-#endif

-            if(r.size() <= i)                // ensure enough space in the remainder, probably not needed

-                r.resize(i+1, Linear(0,0));

-            Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)

-#ifdef DEBUG_INVERSION

-            std::cout << "t1i=" << t1i << std::endl;

-            std::cout << "ci=" << ci << std::endl;

-#endif

-            for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1

-                t1i[dim] *= t1[dim];

-            c[i] = ci; // c(v) := c(v) + c_i(v)*t^i

-            // change from v to u parameterisation

-            SBasis civ = one_minus_a*ci[0] + a*ci[1]; 

-            // r(u) := r(u) - c_i(u)*(t(u))^i

-            // We can truncate this to the number of final terms, as no following terms can

-            // contribute to the result.

-            r -= multiply(civ,ti);

-            r.truncate(k);

-            if(r.tailError(i) == 0)

-                break; // yay!

-            ti = multiply(ti,t);

-        }

-#ifdef DEBUG_INVERSION

-        std::cout << "##########################" << std::endl;

-#endif

-    } else

-        c = Linear(0,1); // linear

-    c -= a0; // invert the offset

-    c /= a1; // invert the slope

-    return c;

-}

-

-SBasis sin(Linear b, int k) {

-    SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));

-    Tri tr(s[0]);

-    double t2 = Tri(b);

-    s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));

-    

-    t2 *= t2;

-    for(int i = 0; i < k; i++) {

-        Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],

-                  -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);

-        bo -= s[i]*(t2/(i+1));

-        

-        

-        s.push_back(bo/double(i+2));

-    }

-    

-    return s;

-}

-

-SBasis cos(Linear bo, int k) {

-    return sin(Linear(bo[0] + M_PI/2,

-                      bo[1] + M_PI/2),

-               k);

-}

-

-//compute fog^-1. ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)

-//TODO: compute order according to tol?

-//TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!

-SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){

-    SBasis result; //result

-    SBasis r=f; //remainder

-    SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;

-    Pk.truncate(order);

-    Qk.truncate(order);

-    Pk.resize(order,Linear(0.));

-    Qk.resize(order,Linear(0.));

-    r.resize(order,Linear(0.));

-

-    int vs= valuation(sg,zero);

-    

-    for (unsigned k=0; k<order; k+=vs){

-        double p10 = Pk.at(k)[0];// we have to solve the linear system:

-        double p01 = Pk.at(k)[1];//

-        double q10 = Qk.at(k)[0];//   p10*a + q10*b = r10

-        double q01 = Qk.at(k)[1];// &                    

-        double r10 =  r.at(k)[0];//   p01*a + q01*b = r01

-        double r01 =  r.at(k)[1];//

-        double a,b;

-        double det = p10*q01-p01*q10;

-

-        //TODO: handle det~0!! 

-        if (fabs(det)<zero){

-            det = zero;

-            a=b=0;

-        }else{

-            a=( q01*r10-q10*r01)/det;

-            b=(-p01*r10+p10*r01)/det;

-        }

-        result.push_back(Linear(a,b));

-        r=r-Pk*a-Qk*b;

-        

-        Pk=Pk*sg;

-        Qk=Qk*sg;

-        Pk.truncate(order);

-        Qk.truncate(order);

-        r.truncate(order);

-    }

-    result.normalize();

-    return result;

-}

-

-}

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

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

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

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

+
+double SBasis::tailError(unsigned tail) const {
+  Interval bs = bounds_fast(*this, tail);
+  return std::max(fabs(bs.min()),fabs(bs.max()));
+}
+
+bool SBasis::isFinite() const {
+    for(unsigned i = 0; i < size(); i++) {
+        if(!(*this)[i].isFinite())
+            return false;
+    }
+    return true;
+}
+
+SBasis operator+(const SBasis& a, const SBasis& b) {
+    SBasis result;
+    const unsigned out_size = std::max(a.size(), b.size());
+    const unsigned min_size = std::min(a.size(), b.size());
+    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;
+}
+
+SBasis operator-(const SBasis& a, const SBasis& b) {
+    SBasis result;
+    const unsigned out_size = std::max(a.size(), b.size());
+    const unsigned min_size = std::min(a.size(), b.size());
+    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;
+}
+
+SBasis& operator+=(SBasis& a, const SBasis& b) {
+    const unsigned out_size = std::max(a.size(), b.size());
+    const unsigned min_size = std::min(a.size(), b.size());
+    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;
+}
+
+SBasis& operator-=(SBasis& a, const SBasis& b) {
+    const unsigned out_size = std::max(a.size(), b.size());
+    const unsigned min_size = std::min(a.size(), b.size());
+    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;
+}
+
+SBasis operator*(SBasis const &a, double k) {
+    SBasis c;
+    c.reserve(a.size());
+    for(unsigned i = 0; i < a.size(); i++)
+        c.push_back(a[i] * k);
+    return c;
+}
+
+SBasis& operator*=(SBasis& 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;
+}
+
+SBasis shift(SBasis const &a, int sh) {
+    SBasis c = a;
+    if(sh > 0) {
+        c.insert(c.begin(), sh, Linear(0,0));
+    } else {
+        //TODO: truncate
+    }
+    return c;
+}
+
+SBasis shift(Linear const &a, int sh) {
+    SBasis c;
+    if(sh > 0) {
+        c.insert(c.begin(), sh, Linear(0,0));
+        c.push_back(a);
+    }
+    return c;
+}
+
+SBasis multiply(SBasis const &a, SBasis const &b) {
+    // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}
+
+    // shift(1, a.Tri*b.Tri)
+    SBasis c;
+    if(a.isZero() || b.isZero())
+        return c;
+    c.resize(a.size() + b.size(), Linear(0,0));
+    c[0] = Linear(0,0);
+    for(unsigned j = 0; j < b.size(); j++) {
+        for(unsigned i = j; i < a.size()+j; i++) {
+            double tri = Tri(b[j])*Tri(a[i-j]);
+            c[i+1/*shift*/] += Linear(Hat(-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;
+}
+
+SBasis integral(SBasis const &c) {
+    SBasis a;
+    a.resize(c.size() + 1, Linear(0,0));
+    a[0] = Linear(0,0);
+
+    for(unsigned k = 1; k < c.size() + 1; k++) {
+        double ahat = -Tri(c[k-1])/(2*k);
+        a[k] = Hat(ahat);
+    }
+    double aTri = 0;
+    for(int k = c.size()-1; k >= 0; k--) {
+        aTri = (Hat(c[k]).d + (k+1)*aTri/2)/(2*k+1);
+        a[k][0] -= aTri/2;
+        a[k][1] += aTri/2;
+    }
+    a.normalize();
+    return a;
+}
+
+SBasis derivative(SBasis const &a) {
+    SBasis c;
+    c.resize(a.size(), Linear(0,0));
+
+    for(unsigned k = 0; k < a.size(); k++) {
+        double d = (2*k+1)*Tri(a[k]);
+
+        for(unsigned dim = 0; dim < 2; dim++) {
+            c[k][dim] = d;
+            if(k+1 < a.size()) {
+                if(dim)
+                    c[k][dim] = d - (k+1)*a[k+1][dim];
+                else
+                    c[k][dim] = d + (k+1)*a[k+1][dim];
+            }
+        }
+    }
+
+    return c;
+}
+
+//TODO: convert int k to unsigned k, and remove cast
+SBasis sqrt(SBasis const &a, int k) {
+    SBasis c;
+    if(a.isZero() || k == 0)
+        return c;
+    c.resize(k, Linear(0,0));
+    c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));
+    SBasis r = a - multiply(c, c); // remainder
+
+    for(unsigned i = 1; i <= (unsigned)k and i<r.size(); i++) {
+        Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));
+        SBasis cisi = shift(ci, i);
+        r -= multiply(shift((c*2 + cisi), i), SBasis(ci));
+        r.truncate(k+1);
+        c += cisi;
+        if(r.tailError(i) == 0) // if exact
+            break;
+    }
+
+    return c;
+}
+
+// return a kth order approx to 1/a)
+SBasis reciprocal(Linear const &a, int k) {
+    SBasis c;
+    assert(!a.isZero());
+    c.resize(k, Linear(0,0));
+    double r_s0 = (Tri(a)*Tri(a))/(-a[0]*a[1]);
+    double r_s0k = 1;
+    for(unsigned i = 0; i < (unsigned)k; i++) {
+        c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);
+        r_s0k *= r_s0;
+    }
+    return c;
+}
+
+SBasis divide(SBasis const &a, SBasis const &b, int k) {
+    SBasis c;
+    assert(!a.isZero());
+    SBasis r = a; // remainder
+
+    k++;
+    r.resize(k, Linear(0,0));
+    c.resize(k, Linear(0,0));
+
+    for(unsigned i = 0; i < (unsigned)k; i++) {
+        Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0
+        c[i] += ci;
+        r -= shift(multiply(ci,b), i);
+        r.truncate(k+1);
+        if(r.tailError(i) == 0) // if exact
+            break;
+    }
+
+    return c;
+}
+
+// a(b)
+// return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
+SBasis compose(SBasis const &a, SBasis const &b) {
+    SBasis s = multiply((SBasis(Linear(1,1))-b), b);
+    SBasis r;
+
+    for(int i = a.size()-1; i >= 0; i--) {
+        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+    }
+    return r;
+}
+
+// a(b)
+// return a0 + s(a1 + s(a2 +...  where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
+SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
+    SBasis s = multiply((SBasis(Linear(1,1))-b), b);
+    SBasis r;
+
+    for(int i = a.size()-1; i >= 0; i--) {
+        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+    }
+    r.truncate(k);
+    return r;
+}
+
+/*
+Inversion algorithm. The notation is certainly very misleading. The
+pseudocode should say:
+
+c(v) := 0
+r(u) := r_0(u) := u
+for i:=0 to k do
+  c_i(v) := H_0(r_i(u)/(t_1)^i; u)
+  c(v) := c(v) + c_i(v)*t^i
+  r(u) := r(u) ? c_i(u)*(t(u))^i
+endfor
+*/
+
+//#define DEBUG_INVERSION 1
+
+SBasis inverse(SBasis a, int k) {
+    assert(a.size() > 0);
+// the function should have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.
+    double a0 = a[0][0];
+    if(a0 != 0) {
+        a -= a0;
+    }
+    double a1 = a[0][1];
+    assert(a1 != 0); // not invertable.
+
+    if(a1 != 1) {
+        a /= a1;
+    }
+    SBasis c;                           // c(v) := 0
+    if(a.size() >= 2 && k == 2) {
+        c.push_back(Linear(0,1));
+        Linear t1(1+a[1][0], 1-a[1][1]);    // t_1
+        c.push_back(Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]));
+    } else if(a.size() >= 2) {                      // non linear
+        SBasis r = Linear(0,1);             // r(u) := r_0(u) := u
+        Linear t1(1./(1+a[1][0]), 1./(1-a[1][1]));    // 1./t_1
+        Linear one(1,1);
+        Linear t1i = one;                   // t_1^0
+        SBasis one_minus_a = SBasis(one) - a;
+        SBasis t = multiply(one_minus_a, a); // t(u)
+        SBasis ti(one);                     // t(u)^0
+#ifdef DEBUG_INVERSION
+        std::cout << "a=" << a << std::endl;
+        std::cout << "1-a=" << one_minus_a << std::endl;
+        std::cout << "t1=" << t1 << std::endl;
+        //assert(t1 == t[1]);
+#endif
+
+        c.resize(k+1, Linear(0,0));
+        for(unsigned i = 0; i < (unsigned)k; i++) {   // for i:=0 to k do
+#ifdef DEBUG_INVERSION
+            std::cout << "-------" << i << ": ---------" <<std::endl;
+            std::cout << "r=" << r << std::endl
+                      << "c=" << c << std::endl
+                      << "ti=" << ti << std::endl
+                      << std::endl;
+#endif
+            if(r.size() <= i)                // ensure enough space in the remainder, probably not needed
+                r.resize(i+1, Linear(0,0));
+            Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)
+#ifdef DEBUG_INVERSION
+            std::cout << "t1i=" << t1i << std::endl;
+            std::cout << "ci=" << ci << std::endl;
+#endif
+            for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1
+                t1i[dim] *= t1[dim];
+            c[i] = ci; // c(v) := c(v) + c_i(v)*t^i
+            // change from v to u parameterisation
+            SBasis civ = one_minus_a*ci[0] + a*ci[1];
+            // r(u) := r(u) - c_i(u)*(t(u))^i
+            // We can truncate this to the number of final terms, as no following terms can
+            // contribute to the result.
+            r -= multiply(civ,ti);
+            r.truncate(k);
+            if(r.tailError(i) == 0)
+                break; // yay!
+            ti = multiply(ti,t);
+        }
+#ifdef DEBUG_INVERSION
+        std::cout << "##########################" << std::endl;
+#endif
+    } else
+        c = Linear(0,1); // linear
+    c -= a0; // invert the offset
+    c /= a1; // invert the slope
+    return c;
+}
+
+SBasis sin(Linear b, int k) {
+    SBasis s = Linear(std::sin(b[0]), std::sin(b[1]));
+    Tri tr(s[0]);
+    double t2 = Tri(b);
+    s.push_back(Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr));
+
+    t2 *= t2;
+    for(int i = 0; i < k; i++) {
+        Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],
+                  -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);
+        bo -= s[i]*(t2/(i+1));
+
+
+        s.push_back(bo/double(i+2));
+    }
+
+    return s;
+}
+
+SBasis cos(Linear bo, int k) {
+    return sin(Linear(bo[0] + M_PI/2,
+                      bo[1] + M_PI/2),
+               k);
+}
+
+//compute fog^-1. ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)
+//TODO: compute order according to tol?
+//TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!
+SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){
+    SBasis result; //result
+    SBasis r=f; //remainder
+    SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;
+    Pk.truncate(order);
+    Qk.truncate(order);
+    Pk.resize(order,Linear(0.));
+    Qk.resize(order,Linear(0.));
+    r.resize(order,Linear(0.));
+
+    int vs= valuation(sg,zero);
+
+    for (unsigned k=0; k<order; k+=vs){
+        double p10 = Pk.at(k)[0];// we have to solve the linear system:
+        double p01 = Pk.at(k)[1];//
+        double q10 = Qk.at(k)[0];//   p10*a + q10*b = r10
+        double q01 = Qk.at(k)[1];// &
+        double r10 =  r.at(k)[0];//   p01*a + q01*b = r01
+        double r01 =  r.at(k)[1];//
+        double a,b;
+        double det = p10*q01-p01*q10;
+
+        //TODO: handle det~0!!
+        if (fabs(det)<zero){
+            det = zero;
+            a=b=0;
+        }else{
+            a=( q01*r10-q10*r01)/det;
+            b=(-p01*r10+p10*r01)/det;
+        }
+        result.push_back(Linear(a,b));
+        r=r-Pk*a-Qk*b;
+
+        Pk=Pk*sg;
+        Qk=Qk*sg;
+        Pk.truncate(order);
+        Qk.truncate(order);
+        r.truncate(order);
+    }
+    result.normalize();
+    return result;
+}
+
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 4329bc559..08674ab2f 100644
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
@@ -1,104 +1,104 @@
-#include "sweep.h"

-

-#include <algorithm>

-

-namespace Geom {

-

-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> rs) {

-    std::vector<Event> events; events.reserve(rs.size()*2);

-    std::vector<std::vector<unsigned> > pairs(rs.size());

-    

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

-        events.push_back(Event(rs[i].left(), i, false));

-        events.push_back(Event(rs[i].right(), i, true));

-    }

-    std::sort(events.begin(), events.end());

-

-    std::vector<unsigned> open;

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

-        unsigned ix = events[i].ix;

-        if(events[i].closing) {

-            std::vector<unsigned>::iterator iter = std::find(open.begin(), open.end(), ix);

-            //if(iter != open.end())

-            open.erase(iter);

-        } else {

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

-                unsigned jx = open[j];

-                if(rs[jx][Y].intersects(rs[ix][Y])) {

-                    pairs[jx].push_back(ix);

-                }

-            }

-            open.push_back(ix);

-        }

-    }

-    return pairs;

-}

-

-std::vector<std::vector<unsigned> > sweep_bounds(std::vector<Rect> a, std::vector<Rect> b) {

-    std::vector<std::vector<unsigned> > pairs(a.size());

-    if(a.empty() || b.empty()) return pairs;

-    std::vector<Event> events[2];

-    events[0].reserve(a.size()*2);

-    events[1].reserve(b.size()*2);

-    

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

-        unsigned sz = n ? b.size() : a.size();

-        events[n].reserve(sz*2);

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

-            events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));

-            events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));

-        }

-        std::sort(events[n].begin(), events[n].end());

-    }

-

-    std::vector<unsigned> open[2];

-    bool n = events[1].front() < events[0].front();

-    for(unsigned i[] = {0,0}; i[n] < events[n].size();) {

-        unsigned ix = events[n][i[n]].ix;

-        bool closing = events[n][i[n]].closing;

-        //std::cout << n << "[" << ix << "] - " << (closing ? "closer" : "opener") << "\n";

-        if(closing) {

-            open[n].erase(std::find(open[n].begin(), open[n].end(), ix));

-        } else {

-            if(n) {

-                //n = 1

-                //opening a B, add to all open a

-                for(unsigned j = 0; j < open[0].size(); j++) {

-                    unsigned jx = open[0][j];

-                    if(a[jx][Y].intersects(b[ix][Y])) {

-                        pairs[jx].push_back(ix);

-                    }

-                }

-            } else {

-                //n = 0

-                //opening an A, add all open b

-                for(unsigned j = 0; j < open[1].size(); j++) {

-                    unsigned jx = open[1][j];

-                    if(b[jx][Y].intersects(a[ix][Y])) {

-                        pairs[ix].push_back(jx);

-                    }

-                }

-            }

-            open[n].push_back(ix);

-        }

-        i[n]++;

-        n = (events[!n][i[!n]] < events[n][i[n]]) ? !n : n;

-    }

-    return pairs;

-}

-

-//Fake cull, until the switch to the real sweep is made.

-std::vector<std::vector<unsigned> > fake_cull(unsigned a, unsigned b) {

-    std::vector<std::vector<unsigned> > ret;

-    

-    std::vector<unsigned> all;

-    for(unsigned j = 0; j < b; j++)

-        all.push_back(j);

-    

-    for(unsigned i = 0; i < a; i++)

-        ret.push_back(all);

-    

-    return ret;

-}

-

-}

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