clipper.cc

Go to the documentation of this file.

// -*- mode: c++; fill-column: 78; coding: latin-1 -*-
/** \mainpage

  This is an implementation of three root finding algorithms for polynomials
  using geometric clipping in the Bezier representation. It was written for a
  lab exercise course in numerical mathematics in winter semester 2011/12 at
  the FernUniversity of Hagen, Germany.

  For more information see the website:
  http://idlebox.net/2012/0320-Finding-Roots-of-Polynomials-by-Clipping/

  The code implements the three algorithms BezierClip, QuadClip [1] and
  CubeClip [2] (both in one class KCLip) described in the referenced
  papers. The program follows the unusual method to printing out debug
  information as a LaTeX document, which itself must be compiled by pdflatex
  for comfortable viewing. This technique allows the algorithm code to be
  interleaved with ltx(...) debug lines which are can contain plots and
  equations.

  There are six main parts in this source:
  - floating-point Numeric classes
  - polynomial classes for monomial and Bezier representation
  - implementation of the degree reduction matrices
  - the BezierClip and KClip algorithms
  - a test suite to check above algorithms
  - a collection of demos and speedtests

  Each part is separately documented below (or in Related Pages).

  References:

  [1] Barton, M. and Juettler, B. "Computing roots of polynomials by quadratic
  clipping", Computer Aided Geometric Design 24, (2007) p. 125-141.

  [2] Liu, L., Zhang, L., Lin, B., Wang, G. "Fast approach for computing roots
  of polynomials using cubic clipping", Computer Aided Geometric Design 26,
  (2009), p. 524-599.

\verbatim
  ----------------------------------------------------------------------------

  Copyright (c) 2011-2012, Timo Bingmann
  All rights reserved. Permissions granted under the BSD 2-Clause License:

  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that the following conditions are met:

  * Redistributions of source code must retain the above copyright notice,
    this list of conditions and the following disclaimer.

  * Redistributions in binary form must reproduce the above copyright notice,
    this list of conditions and the following disclaimer in the documentation
    and/or other materials provided with the distribution.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  POSSIBILITY OF SUCH DAMAGE.

  ----------------------------------------------------------------------------
\endverbatim
*/

#include <stdlib.h>
#include <assert.h>
#include <getopt.h>
#include <sys/time.h>

#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <sstream>
#include <stdexcept>
#include <cmath>
#include <limits>

#include <mpfr.h>

unsigned int output_latex = false;      // output latex or text debug

#if !SPEEDTEST

#define dbg(X)  do { if (!output_latex) std::cout << X; } while(0)
#define ltx(X)  do { if (output_latex) std::cout << X; } while(0)
#define spd(X)  dbg(X)

#else

// all output functions are disabled during speedtest except spd(...)
#define dbg(X)  do { } while(0)
#define ltx(X)  do { } while(0)
#define spd(X)  do { std::cout << X; } while(0)

#endif

bool debug_more = false;                // be more verbose in some parts
unsigned int print_precision = 6;       // precision of mpfr numbers in output

// ****************************************************************************
// *** Auxiliary Functions                                                  ***
// ****************************************************************************

/// format many different types as strings
template <typename Type>
std::string S(Type i)
{
    std::ostringstream oss;
    oss << i;
    return oss.str();
}

/// format a number converting "e+10" to latex math strings
std::string ltxNum(std::string os)
{
    int i = os.size();
    while(--i > 0)
    {
        if (os[i] == 'e')
        {
            if (os[i+1] == '+')
                os.erase(os.begin()+i, os.begin()+i+1);

            while (os[i+1] == '0')
                os.erase(os.begin()+i, os.begin()+i+1);

            os.replace(i,1," \\!\\cdot\\! 10^{");
            os += "}";
            return os;
        }
    }
    return os;
}

/// format a number converting "e+10" to latex math strings
template <typename Numeric>
std::string ltxNum(Numeric n)
{
    return ltxNum(n.toString());
}

/// format interval in latex headers
template <typename Numeric>
std::string ltxInterval(Numeric l, Numeric r)
{
    std::ostringstream oss;
    oss << "\\texorpdfstring{$[" << l << "," << r << "]$}{" << l << "," << r << "}";
    return oss.str();
}

/// time measuring using gettimeofday().
inline double timestamp()
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec * 1e-6;
}

/** \page numeric Floating Point Numeric Classes

*******************************************************************************
*** Floating Point Numeric Classes                                          ***
*******************************************************************************

 This part contains three classes implementing floating point numbers with
 different precisions and properties. These are later used for all
 calculations requiring a "Numeric" type and define many common math functions
 like sqrt() and also most C++ arithmetic operators.

 The first class Fraction keeps numerator and denominator of a rational number
 separately, and supports pretty-printed using LaTeX's \frac{} command. The
 purpose of this class is to output nicely formatted fractions and it should
 not be used for serious calculations.

 The second class Double takes a standard hardware floating-point type as
 template parameter: "float", "double" and "long double" are valid choices. It
 defines all later required operations using the standard libc implements of
 e.g. sqrt(), cbrt(), sin() and so on. Most of this is done using <cmath>, but
 some explicit template specializations were necessary to correctly switch
 between other libc functions.

 The final class MpfrFloat utilizes the GNU multi-precision float library with
 rounding to create a "Numeric" implementation with a fixed arbitrary
 precision. All necessary C++ arithmetic operators are overloaded in
 implemented using appropriate mpfr_functions.

******************************************************************************/

/// Contains a pair of doubles signifying the numerator and denominator of a
/// fraction. The fraction can be printed using LaTeX syntax and many common
/// operators are defined.
class Fraction
{
private:
    double      m_numerator;
    double      m_denominator;

    /// Calculate the greatest common divisor using Euclid's algorithm.
    static double gcd(double a, double b)
    {
        while (b != 0) {
            double t = b; b = fmod(a,b); a = t;
        }
        return a;
    }

    void reduce()
    {
        double g = gcd(m_numerator, m_denominator);
        m_numerator /= g;
        m_denominator /= g;

        if (m_denominator < 0 && m_numerator > 0) {
            m_denominator *= -1;
            m_numerator *= -1;
        }
    }

public:
    Fraction(double numerator = 0, double denominator = 1)
        : m_numerator(numerator), m_denominator(denominator)
    {
        if (denominator == 0) throw(std::runtime_error("Divide by zero!"));
        reduce();
    }

    double numerator() const
    {
        return m_numerator;
    }

    double denominator() const
    {
        return m_denominator;
    }

    std::string toString() const
    {
        std::ostringstream oss;
        if (m_denominator == 1)
            oss << ltxNum(S(m_numerator));
        else
            oss << "\\frac{" << ltxNum(S(m_numerator)) << "}"
                "{" << ltxNum(S(m_denominator)) << "}";
        return oss.str();
    }

    friend std::ostream& operator<< (std::ostream& os, const Fraction& f)
    {
        return os << f.toString();
    }

    double toDouble() const
    {
        return m_numerator / m_denominator;
    }

    Fraction operator+ (const Fraction& b) const
    {
        return Fraction( m_numerator * b.m_denominator + b.m_numerator * m_denominator, m_denominator * b.m_denominator );
    }

    Fraction operator- (const Fraction& b) const
    {
        return Fraction( m_numerator * b.m_denominator - b.m_numerator * m_denominator, m_denominator * b.m_denominator );
    }

    Fraction operator* (const Fraction& b) const
    {
        return Fraction( m_numerator * b.m_numerator, m_denominator * b.m_denominator );
    }

    Fraction operator/ (const Fraction& b) const
    {
        return Fraction( m_numerator * b.m_denominator, m_denominator * b.m_numerator );
    }

    Fraction& operator+= (const Fraction& b)
    {
        m_numerator = m_numerator * b.m_denominator + b.m_numerator * m_denominator;
        m_denominator = m_denominator * b.m_denominator;
        reduce();
        return *this;
    }

    /// simple method to calculate binomial coefficients for small (n,k)
    static Fraction binomial(unsigned int n, unsigned int k)
    {
        static const double factorial[21] = {
            1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800,
            39916800, 479001600, 6227020800, 87178291200, 1307674368000,
            20922789888000, 355687428096000, 6402373705728000,
            121645100408832000, 2432902008176640000
        };

        assert( n < sizeof(factorial) / sizeof(factorial[0]) );
        assert( k < sizeof(factorial) / sizeof(factorial[0]) );

        if (k > n) return 0;

        return Fraction(factorial[n], factorial[k] * factorial[n-k]);
    }
};

/// Implements a wrapper around one of the hardware floating-point data types
/// so that these can be used as a "Numeric" template parameter.
template <typename Type = double>
class Double
{
private:
    Type        m_v;    // the value

public:
    inline Double(const Type& v = 0)
        : m_v(v)
    {
    }

    // multiplied to the number of iterations in a speed test
    static const unsigned int IterationScale    = 20;

    static std::string getName();

    static inline Double MINUS_INFINITY()
    {
        return Double( - std::numeric_limits<Type>::infinity() );
    }

    static inline Double binomial(unsigned int n, unsigned int k);

    inline const Type& get() const
    {
        return m_v;
    }

    inline std::string toString() const
    {
        return S(m_v);
    }

    inline friend std::ostream& operator<< (std::ostream& os, const Double<Type>& d)
    {
        return os << d.toString();
    }

    inline Double operator- () const
    {
        return Double( - m_v );
    }

    inline Double operator+ (const Double& o) const
    {
        return Double( m_v + o.m_v );
    }

    inline Double operator- (const Double& o) const
    {
        return Double( m_v - o.m_v );
    }

    inline Double operator* (const Double& o) const
    {
        return Double( m_v * o.m_v );
    }

    inline Double operator/ (const Double& o) const
    {
        return Double( m_v / o.m_v );
    }

    inline Double& operator+= (const Double& o)
    {
        m_v += o.m_v;
        return *this;
    }

    inline Double& operator-= (const Double& o)
    {
        m_v -= o.m_v;
        return *this;
    }

    inline Double& operator*= (const Double& o)
    {
        m_v *= o.m_v;
        return *this;
    }

    inline Double& operator/= (const Double& o)
    {
        m_v /= o.m_v;
        return *this;
    }

    inline bool operator== (const Double &o) const
    {
        return (m_v == o.m_v);
    }

    inline bool operator!= (const Double &o) const
    {
        return (m_v != o.m_v);
    }

    inline bool operator< (const Double &o) const
    {
        return (m_v < o.m_v);
    }

    inline bool operator> (const Double &o) const
    {
        return (m_v > o.m_v);
    }

    inline bool operator>= (const Double &o) const
    {
        return (m_v >= o.m_v);
    }

    inline bool operator<= (const Double &o) const
    {
        return (m_v <= o.m_v);
    }
};

template <>
std::string Double<float>::getName()
{
    return "flt";
}

template <>
std::string Double<double>::getName()
{
    return "dbl";
}

template <>
std::string Double<long double>::getName()
{
    return "ldbl";
}

template <typename Type>
inline Double<Type> operator- (const double& a, const Double<Type>& b)
{
    return Double<Type>( a - b.get() );
}

template <typename Type>
inline Double<Type> operator* (const double& a, const Double<Type>& b)
{
    return Double<Type>( a * b.get() );
}

template <typename Type>
inline Double<Type> abs(const Double<Type>& v)
{
    return Double<Type>( std::abs(v.get()) );
}

template <typename Type>
inline Double<Type> sqrt(const Double<Type>& v)
{
    return Double<Type>( std::sqrt(v.get()) );
}

inline Double<float> cbrt(const Double<float>& v)
{
    return Double<float>( cbrtf(v.get()) );
}

inline Double<double> cbrt(const Double<double>& v)
{
    return Double<double>( cbrt(v.get()) );
}

inline Double<long double> cbrt(const Double<long double>& v)
{
    return Double<long double>( cbrtl(v.get()) );
}

template <typename Type>
inline Double<Type> cos(const Double<Type>& v)
{
    return Double<Type>( std::cos(v.get()) );
}

template <typename Type>
inline Double<Type> acos(const Double<Type>& v)
{
    return Double<Type>( std::acos(v.get()) );
}

template <typename Type, typename Exponent>
inline Double<Type> pow(const Double<Type>& v, Exponent e)
{
    return Double<Type>( std::pow(v.get(), e) );
}

template <>
inline Double<float> Double<float>::binomial(unsigned int n, unsigned int k)
{
    return Double<float>( std::exp(lgammaf(n+1) - lgammaf(k+1) - lgammaf(n-k+1)) );
}

template <>
inline Double<double> Double<double>::binomial(unsigned int n, unsigned int k)
{
    return Double<double>( std::exp(lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1)) );
}

template <>
inline Double<long double> Double<long double>::binomial(unsigned int n, unsigned int k)
{
    return Double<long double>( std::exp(lgammal(n+1) - lgammal(k+1) - lgammal(n-k+1)) );
}

template <typename Type>
inline int isnormal(const Double<Type>& v)
{
    return std::isnormal(v.get());
}

/// Wraps a mpfr_t object of the multiple-precision floating-pointer number
/// library mpfr into a Numeric object and provides all necessary operators
/// and functions for further calculations.
class MpfrFloat
{
private:
    mpfr_t      m_v;    // the value

public:
    inline MpfrFloat()
    {
        mpfr_init(m_v);
    }

    inline MpfrFloat(const double& d)
    {
        mpfr_init_set_d(m_v, d, MPFR_RNDN);
    }

    /// copy constructor
    inline MpfrFloat(const MpfrFloat& d)
    {
        mpfr_init_set(m_v, d.m_v, MPFR_RNDN);
    }

    inline MpfrFloat& operator= (const MpfrFloat& d)
    {
        mpfr_set(m_v, d.m_v, MPFR_RNDN);
        return *this;
    }

    ~MpfrFloat()
    {
        mpfr_clear(m_v);
    }

    // multiplied to the number of iterations in a speed test
    static const unsigned int IterationScale    = 1;

    static std::string getName()
    {
        return "mpfr" + S(mpfr_get_default_prec());
    }

    static inline MpfrFloat MINUS_INFINITY()
    {
        MpfrFloat a; mpfr_set_inf(a.get(), -1); return a;
    }

    static inline MpfrFloat binomial(unsigned int n, unsigned int k)
    {
        MpfrFloat a, b;

        mpfr_fac_ui(a.get(), n, MPFR_RNDN);

        mpfr_fac_ui(b.get(), k, MPFR_RNDN);
        mpfr_div(a.get(), a.get(), b.get(), MPFR_RNDN);

        mpfr_fac_ui(b.get(), n - k, MPFR_RNDN);
        mpfr_div(a.get(), a.get(), b.get(), MPFR_RNDN);

        return a;
    }

    inline const mpfr_t& get() const
    {
        return m_v;
    }

    inline mpfr_t& get()
    {
        return m_v;
    }

    inline std::string toString() const
    {
        std::string str; str.resize(32);
        int r = mpfr_snprintf((char*)str.data(), str.size(), "%.*Rg", print_precision, m_v);
        str.resize(r);
        return str;
    }

    inline friend std::ostream& operator<< (std::ostream& os, const MpfrFloat& a)
    {
        return os << a.toString();
    }

    inline MpfrFloat operator- () const
    {
        MpfrFloat a; mpfr_neg(a.m_v, m_v, MPFR_RNDN); return a;
    }

    inline MpfrFloat operator+ (const MpfrFloat& o) const
    {
        MpfrFloat a; mpfr_add(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
    }

    inline MpfrFloat operator- (const MpfrFloat& o) const
    {
        MpfrFloat a; mpfr_sub(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
    }

    inline MpfrFloat operator* (const MpfrFloat& o) const
    {
        MpfrFloat a; mpfr_mul(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
    }

    inline MpfrFloat operator/ (const MpfrFloat& o) const
    {
        MpfrFloat a; mpfr_div(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
    }

    inline MpfrFloat& operator+= (const MpfrFloat& o)
    {
        mpfr_add(m_v, m_v, o.m_v, MPFR_RNDN);
        return *this;
    }

    inline MpfrFloat& operator-= (const MpfrFloat& o)
    {
        mpfr_sub(m_v, m_v, o.m_v, MPFR_RNDN);
        return *this;
    }

    inline MpfrFloat& operator*= (const MpfrFloat& o)
    {
        mpfr_mul(m_v, m_v, o.m_v, MPFR_RNDN);
        return *this;
    }

    inline MpfrFloat& operator/= (const MpfrFloat& o)
    {
        mpfr_div(m_v, m_v, o.m_v, MPFR_RNDN);
        return *this;
    }

    inline bool operator== (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) == 0;
    }

    inline bool operator!= (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) != 0;
    }

    inline bool operator< (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) < 0;
    }

    inline bool operator> (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) > 0;
    }

    inline bool operator<= (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) <= 0;
    }

    inline bool operator>= (const MpfrFloat &o) const
    {
        return mpfr_cmp(m_v, o.m_v) >= 0;
    }
};

inline MpfrFloat operator- (const double& a, const MpfrFloat& b)
{
    MpfrFloat x; mpfr_d_sub(x.get(), a, b.get(), MPFR_RNDN); return x;
}

inline MpfrFloat operator* (const double& a, const MpfrFloat& b)
{
    MpfrFloat x; mpfr_mul_d(x.get(), b.get(), a, MPFR_RNDN); return x;
}

inline MpfrFloat abs(const MpfrFloat& v)
{
    MpfrFloat x; mpfr_abs(x.get(), v.get(), MPFR_RNDN); return x;
}

inline MpfrFloat sqrt(const MpfrFloat& v)
{
    MpfrFloat x; mpfr_sqrt(x.get(), v.get(), MPFR_RNDN); return x;
}

inline MpfrFloat cbrt(const MpfrFloat& v)
{
    MpfrFloat x; mpfr_cbrt(x.get(), v.get(), MPFR_RNDN); return x;
}

inline MpfrFloat cos(const MpfrFloat& v)
{
    MpfrFloat x; mpfr_cos(x.get(), v.get(), MPFR_RNDN); return x;
}

inline MpfrFloat acos(const MpfrFloat& v)
{
    MpfrFloat x; mpfr_acos(x.get(), v.get(), MPFR_RNDN); return x;
}

inline MpfrFloat pow(const MpfrFloat& v, long int e)
{
    MpfrFloat x; mpfr_pow_si(x.get(), v.get(), e, MPFR_RNDN); return x;
}

inline int isnormal(const MpfrFloat& v)
{
    return mpfr_number_p(v.get());
}

/** \page polynomials Polynomial Classes in Monomial and Bezier Bases

*******************************************************************************
*** Polynomial Classes in Monomial and Bezier Bases                         ***
*******************************************************************************

 This section of the implementation utilizes the basic arithmetic operations
 of a "Numeric" implementation to represent polynomials. Two representations
 are available: the standard monomial base and the Bezier
 representation. These are implemented in the classes PolynomialStandard and
 PolynomialBezier, respectively.

 Both classes implement many smaller algorithms later needed for implementing
 the three root finding algorithms. Among the base algorithms are:

 - conversion between the two polynomial representations
 - constructing a polynomial from given roots
 - finding the roots of quadratic and cubic polynomial using the direct
   formulas.
 - evaluation of the polynomial at a give position (using Horner or de
   Casteljau)
 - splitting a Bezier representation using de Castljau's algorithm
 - calculating the convex hull of the Bezier polygon using Jarvis's march

******************************************************************************/

/// Base class for following polynomials in different bases. This could be
/// used for some common virtual functions.
template <typename Numeric>
class Polynomial
{
};

// *** Prototype declarations
template <typename Numeric>
class PolynomialBezier;

template <typename Numeric>
class Matrix;

/// Represents a polynomial in standard monomial base.
template <typename Numeric>
class PolynomialStandard : public Polynomial<Numeric>
{
private:

    /// Coefficients of the polynomial. Degree is m_coefficients.size()-1
    std::vector<Numeric>        m_coefficient;

public:
    PolynomialStandard()
    {
    }

    explicit PolynomialStandard(const std::vector<Numeric>& coefficient)
        : m_coefficient(coefficient)
    {
    }

    /// Construct polynomial from roots and a scaling multiplier
    PolynomialStandard(Numeric factor, const std::vector<Numeric>& roots)
        : m_coefficient(roots.size()+1, 0)
    {
        std::vector<bool> enumerate(roots.size(), 0);

        // iterate over all 2^n combinations to select a set of roots
        while(1)
        {
            // increment enumeration, calculate produkt and count factors
            Numeric a = 1.0;
            unsigned int afact = 0;

            unsigned int i = 0;
            for (; i < roots.size(); ++i)
            {
                if (enumerate[i] == false) {
                    enumerate[i] = true;
                    a *= (- roots[i]);
                    ++afact;
                    break;
                }
                else {
                    enumerate[i] = false;
                }
            }

            if (i == roots.size()) break; // all 2^n combinations done

            // multiply with remaining factors
            for (++i; i < roots.size(); ++i)
            {
                if (enumerate[i] == true) {
                    a *= (- roots[i]);
                    ++afact;
                }
            }

            // add to correct summand determinted by factor number
            m_coefficient[roots.size() - afact] += a;
        }

        m_coefficient[roots.size()] = 1.0;

        for (unsigned int j = 0; j <= roots.size(); ++j)
            m_coefficient[j] *= factor;
    }

    /// Return degree of the polynomial.
    unsigned int degree() const
    {
        return m_coefficient.size()-1;
    }

    /// Return a monomial coefficient.
    Numeric coefficient(unsigned int i) const
    {
        assert(i <= degree());
        return m_coefficient[i];
    }

    /// Return polynomial as latex math string.
    std::string toString(const std::string& X = "X") const
    {
        std::ostringstream oss;
        for (int i = m_coefficient.size()-1; i >= 0; --i)
        {
            if (m_coefficient[i] != 0) {
                if (oss.str().size())
                    oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
                else if (m_coefficient[i] < 0)
                    oss << "- ";

                if (i >= 10)
                    oss << ltxNum(abs(m_coefficient[i])) << " " << X << "^{" << i << "}";
                else if (i > 1)
                    oss << ltxNum(abs(m_coefficient[i])) << " " << X << "^" << i;
                else if (i == 1)
                    oss << ltxNum(abs(m_coefficient[i])) << " " << X;
                else
                    oss << ltxNum(abs(m_coefficient[i]));
            }
        }
        return oss.str();
    }

    /// Return polynomial as string for plotting with gnuplot.
    std::string toStringPlot() const
    {
        std::ostringstream oss;
        for (int i = m_coefficient.size()-1; i >= 0; --i)
        {
            if (m_coefficient[i] != 0) {
                if (oss.str().size())
                    oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
                else if (m_coefficient[i] < 0)
                    oss << "- ";

                if (i > 1)
                    oss << abs(m_coefficient[i]) << " * x**" << i;
                else if (i == 1)
                    oss << abs(m_coefficient[i]) << " * x";
                else
                    oss << abs(m_coefficient[i]);
            }
        }
        return oss.str();
    }

    /// Calculate polynomial's value at the position X using Horner's schema.
    Numeric evaluate(Numeric X) const
    {
        Numeric val = 0;
        for (int i = m_coefficient.size()-1; i >= 0; --i)
        {
            val *= X;
            val += m_coefficient[i];
        }
        return val;
    }

    /// Calculate all real roots for a quadratic or cubic polynomial using the
    /// p-q- or Cardano's formulas.
    std::vector<Numeric> findRoots() const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::vector<Numeric> roots;

        if (N == 0
            || (N == 1 && m_coefficient[1] == 0)
            || (N == 2 && m_coefficient[1] == 0 && m_coefficient[2] == 0)
            || (N == 3 && m_coefficient[1] == 0 && m_coefficient[2] == 0 && m_coefficient[3] == 0)
            )
        {
            return roots;
        }
        else if (N == 1
                 || (N == 2 && m_coefficient[2] == 0)
                 || (N == 3 && m_coefficient[2] == 0 && m_coefficient[3] == 0)
            )
        {
            roots.push_back( - m_coefficient[0] / m_coefficient[1] );
        }
        else if (N == 2
                 || (N == 3 && m_coefficient[3] == 0)
            )
        {
            Numeric D = (m_coefficient[1] * m_coefficient[1]) - (4.0 * m_coefficient[2] * m_coefficient[0]);

            if (D == 0)                 // one double root
            {
                roots.push_back( - m_coefficient[1] / (2 * m_coefficient[0]) );
            }
            else if (D > 0)             // two real roots
            {
                D = sqrt(D);
                roots.push_back( (- m_coefficient[1] - D) / (2.0 * m_coefficient[2]) );
                roots.push_back( (- m_coefficient[1] + D) / (2.0 * m_coefficient[2]) );

                if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            }
            // else (D < 0) only imaginary roots
            return roots;
        }
        else if (N == 3)
        {
            const std::vector<Numeric>& c = m_coefficient;

            // Cardano's formulas
            Numeric p = ( c[1] - (c[2] * c[2]) / (3 * c[3]) ) / c[3];
            Numeric q = ( ( 2.0 * c[2] * c[2] * c[2] / (27.0 * c[3])  - (c[2] * c[1] / 3.0) )
                         / c[3] + c[0] ) / c[3];

            //std::cout << "p = " << p << " - q = " << q << "\n";

            if (p == 0)                 // simple case: no quadratic term
            {
                roots.push_back( cbrt(q) - c[2] / 3.0 / c[3] );
                return roots;
            }

            Numeric D = p * p * p / 27.0 + q * q / 4.0;

            //std::cout << D << " - D\n";

            if (D > 0)                  // one real root and two complex ones
            {
                D = sqrt(D);

                Numeric u = cbrt( - q / 2.0 + D );
                Numeric v = cbrt( - q / 2.0 - D );

                roots.push_back( u + v - c[2] / 3.0 / c[3] );
            }
            else if (D == 0)            // three real roots: one single and one double
            {
                Numeric u = cbrt( - q / 2.0 );
                Numeric v = cbrt( - q / 2.0 );

                roots.push_back( u + v - c[2] / 3.0 / c[3] );
                roots.push_back( -(u + v) / 2.0 - c[2] / 3.0 / c[3] );

                if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            }
            else                        // three real roots
            {
                Numeric phi = acos( - q / (2.0 * sqrt( - p * p * p / 27.0 )) );
                Numeric w = 2 * sqrt( - p / 3.0 );

                roots.push_back( w * cos( (phi + 2.0 * M_PI) / 3.0 ) );
                roots.push_back( w * cos( (phi + 4.0 * M_PI) / 3.0 ) );
                roots.push_back( w * cos( phi / 3.0 ) );

                roots[0] = roots[0] - (c[2] / 3.0 / c[3]);
                roots[1] = roots[1] - (c[2] / 3.0 / c[3]);
                roots[2] = roots[2] - (c[2] / 3.0 / c[3]);

                if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            }
        }
        else {
            assert("Here roots are only calculable for quadratic and cubic polynomials.");
        }

        return roots;
    }

    /// Apply the linear transformation "X * (beta - alpha) + alpha" to the
    /// polynomial. This remaps the interval [alpha,beta] to [0,1].
    PolynomialStandard scale(Numeric alpha, Numeric beta) const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::vector<Numeric> coeff(N+1);

        for (unsigned int i = 0; i <= N; ++i)
        {
            coeff[i] = 0;
            for (unsigned int j = i; j <= N; ++j)
            {
                coeff[i] += m_coefficient[j] * pow(alpha, j-i) * Numeric::binomial(j,i);
            }
            coeff[i] *= pow((beta - alpha), i);
        }

        return PolynomialStandard(coeff);
    }

    /// Convert the polynomial from standard base to the Bernstein-Bezier
    /// basis in the unit interval. Use scale() to transform the coefficients
    /// before applying this base conversion.
    class PolynomialBezier<Numeric> toBezier() const;
};

/// Represents a polynomial in the Bernstein-Bezier basis using
/// Bezier coefficients.
template <typename Numeric>
class PolynomialBezier : public Polynomial<Numeric>
{
private:

    /// The n+1 Bezier coefficients to the Bernstein polynomials basis
    /// (\f$B_{k,n}\f$)
    std::vector<Numeric>         m_coefficient;

    // Thus the polynomial's degree n is m_coefficient.size()-1.

public:
    explicit PolynomialBezier(const std::vector<Numeric>& coefficient)
        : m_coefficient(coefficient)
    {
    }

    /// Return degree of the polynomial.
    unsigned int degree() const
    {
        return m_coefficient.size()-1;
    }

    /// Return a Bezier coefficient.
    Numeric coefficient(unsigned int i) const
    {
        assert(i <= degree());
        return m_coefficient[i];
    }

    /// Return polynomial as latex math string.
    std::string toString(const std::string& T = "") const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::ostringstream oss;
        for (unsigned int i = 0; i <= N; ++i)
        {
            if (m_coefficient[i] != 0 || 1) {
                if (oss.str().size())
                    oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
                else if (m_coefficient[i] < 0)
                    oss << "- ";

                oss << ltxNum(abs(m_coefficient[i])) << " B_{" << i << "," << N << "}" << T;
            }
        }

        return oss.str();
    }

    /// Return Bezier coefficients as a polygon.
    std::string toPolygon() const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::ostringstream oss;
        for (unsigned int i = 0; i <= N; ++i)
        {
            oss << "(" << Numeric( (double)i / N ) << "," << m_coefficient[i] << ") ";
        }
        return oss.str();
    }

    /// Return the convex hull of the Bezier coefficients calculated using
    /// Jarvis's march.
    std::vector< std::pair<Numeric,Numeric> > getConvexHull() const
    {
        std::vector< std::pair<Numeric,Numeric> > ch;

        const unsigned int N = m_coefficient.size()-1;

        unsigned int i = 0;     // current point in convex hull;
        ch.push_back( std::make_pair(0, m_coefficient[0]) );

        while( i < N ) // search right for point with highest angle
        {
            Numeric amax = Numeric::MINUS_INFINITY();
            unsigned int imax = i;

            for (unsigned int j = i+1; j <= N; ++j)
            {
                Numeric a = (m_coefficient[j] - m_coefficient[i]) / (j - i);

                if (amax < a) {
                    amax = a;
                    imax = j;
                }
            }

            ch.push_back( std::make_pair(Numeric( (double)imax / N ), m_coefficient[imax]) );
            i = imax;
        }

        while ( i > 0 ) // search backwards for point with lowest angle
        {
            Numeric amax = Numeric::MINUS_INFINITY();
            unsigned int imax = i;

            for (unsigned int j = 0; j < i; ++j)
            {
                Numeric a = (m_coefficient[i] - m_coefficient[j]) / (i - j);

                if (amax < a) {
                    amax = a;
                    imax = j;
                }
            }

            ch.push_back( std::make_pair(Numeric( (double)imax / N ), m_coefficient[imax]) );
            i = imax;
        }

        return ch;
    }

    /// Return the convex hull of the Bezier coefficients calculated using
    /// Jarvis's march.
    std::vector<Numeric> getConvexHullIntersection() const
    {
        const unsigned int N = m_coefficient.size()-1;

        unsigned int i = 0;             // current point in convex hull
        Numeric v = m_coefficient[0];   // current y position

        std::vector<Numeric> zerocross; // zero crossings

        while( i < N ) // search right for point with highest angle
        {
            Numeric amax = Numeric::MINUS_INFINITY();
            unsigned int imax = i;

            for (unsigned int j = i+1; j <= N; ++j)
            {
                Numeric a = (m_coefficient[j] - m_coefficient[i]) / (j - i);

                if (amax < a) {
                    amax = a;
                    imax = j;
                }
            }

            if (m_coefficient[imax] * v < 0 || m_coefficient[imax] == 0) {
                //std::cout << "amax " << amax << "\n";
                //std::cout << "df " << v << "\n";
                //std::cout << "df " << v / amax << "\n";
                zerocross.push_back( Numeric( (double)i / N ) - v / (amax * N) );
            }

            v = m_coefficient[imax];
            i = imax;
        }

        while ( i > 0 ) // search backwards for point with lowest angle
        {
            Numeric amax = Numeric::MINUS_INFINITY();
            unsigned int imax = i;

            for (unsigned int j = 0; j < i; ++j)
            {
                Numeric a = (m_coefficient[i] - m_coefficient[j]) / (i - j);

                if (amax < a) {
                    amax = a;
                    imax = j;
                }
            }

            if (m_coefficient[imax] * v < 0 || m_coefficient[imax] == 0) {
                //std::cout << "amax " << amax << "\n";
                //std::cout << "df " << v << "\n";
                //std::cout << "df " << v / amax << "\n";
                zerocross.push_back( Numeric( (double)i / N ) - v / (amax * N) );
            }

            v = m_coefficient[imax];
            i = imax;
        }

#if 0
        for (unsigned int i = 0; i < zerocross.size(); ++i)
            std::cout << "zerocross[" << i << "] = " << zerocross[i] << "\n";
#endif

        for (unsigned int i = 0; i < zerocross.size(); ++i)
            if (!isnormal(zerocross[i])) zerocross.clear();

        assert(zerocross.size() == 0 || zerocross.size() == 2);

        if (zerocross.size() == 2) {
            if (zerocross[0] > zerocross[1])
                std::swap(zerocross[0], zerocross[1]);
        }

        return zerocross;
    }

    /// Return the convex hull of the Bezier coefficients for plotting.
    std::string toConvexHullString() const
    {
        std::vector< std::pair<Numeric,Numeric> > ch = getConvexHull();

        std::ostringstream os;
        for (unsigned int i = 0; i < ch.size(); ++i)
        {
            if (i != 0) os << " ";
            os << "(" << ch[i].first << "," << ch[i].second << ")";
        }
        return os.str();
    }

    /// Transform the polynomial from Bernstein-Bezier basis back to the
    /// standard monomial basis.
    class PolynomialStandard<Numeric> toStandard() const;

    /// Calculate the polynomial's value at the position X using the schema of
    /// de Casteljau.
    Numeric evaluate(Numeric X) const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::vector<Numeric> coeff = m_coefficient;

        for (unsigned int i = 1; i <= N; ++i)
        {
            Numeric save = coeff[0];

            for (unsigned int j = 0; j <= N-i; ++j)
            {
                coeff[j] = X * coeff[j+1] + (1.0 - X) * save;
                save = coeff[j+1];
            }
        }

        return coeff[0];
    }

    /// Similar to evaluate: run de Casteljau's algorithm and output the
    /// Bezier coefficients of interval [0,t0] and [t0,1] as separate
    /// polynomial objects. This is later used to branch when finding roots.
    std::pair<PolynomialBezier,PolynomialBezier> deCasteljauSplit(Numeric t0) const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::vector<Numeric> coeff = m_coefficient;

        std::vector<Numeric> coeff_left(N+1), coeff_right(N+1);

        coeff_left[0] = coeff[0];
        coeff_right[N] = coeff[N];

        for (unsigned int i = 1; i <= N; ++i)
        {
            Numeric save = coeff[0];

            for (unsigned int j = 0; j <= N-i; ++j)
            {
                coeff[j] = t0 * coeff[j+1] + (1.0 - t0) * save;
                save = coeff[j+1];
            }

            coeff_left[i] = coeff[0];
            coeff_right[N-i] = coeff[N-i];
        }

        return std::make_pair( PolynomialBezier(coeff_left), PolynomialBezier(coeff_right) );
    }

    /// Calculate all real roots for a quadratic or cubic polynomial in Bezier
    /// representation (according to Remark 3 in [1]). Imaginary roots are
    /// ignored.
    std::vector<Numeric> findRoots() const
    {
        const unsigned int N = m_coefficient.size()-1;

        std::vector<Numeric> roots;

        if (N == 2)
        {
            Numeric D = (m_coefficient[1] * m_coefficient[1]) - (m_coefficient[0] * m_coefficient[2]);

            // numerator of all roots
            Numeric N = m_coefficient[2] - 2 * m_coefficient[1] + m_coefficient[0];
            if (N == 0) return roots;

            if (D == 0)         // one double root
            {
                roots.push_back( (m_coefficient[0] - m_coefficient[1]) / N );
            }
            else if (D > 0)             // two roots
            {
                D = sqrt(D);

                roots.push_back( (m_coefficient[0] - m_coefficient[1] - D) / N );

                roots.push_back( (m_coefficient[0] - m_coefficient[1] + D) / N );

                if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
            }
            // else (D < 0) only imaginary roots
            return roots;
        }
        else if (N == 3)
        {
            const std::vector<Numeric>& c = m_coefficient;

            /// Instead of using the Cardano's formulas in Bezier form as in
            /// [2], we return to the monomial base case for calculating the
            /// roots
            std::vector<Numeric> newCoeff(4);

            newCoeff[0] = c[0];
            newCoeff[1] = 3 * c[1] - 3 * c[0];
            newCoeff[2] = 3 * c[0] - 6 * c[1] + 3 * c[2];
            newCoeff[3] = c[3] - 3 * c[2] + 3 * c[1] - c[0];

            return PolynomialStandard<Numeric>(newCoeff).findRoots();
        }
        else {
            assert("Here roots are only calculable for quadratic and cubic polynomials.");
        }

        return roots;
    }

    /// This function is the heart of the root finding algorithm: reduce or
    /// raise the degree of this polynomial to k.
    PolynomialBezier    interpolateDegree(int k) const;

    /// This function is the heart of the root finding algorithm: reduce or
    /// raise the degree of this polynomial by applying a matrix.
    PolynomialBezier    interpolateDegree(const class Matrix<Numeric> &m) const;

    /// Compare two Bezier polynomials of same degree and return the maximum
    /// difference between pairs of coefficients.
    Numeric maxCoefficientDeltaTo(const PolynomialBezier& b) const
    {
        assert(m_coefficient.size() == b.m_coefficient.size());
        Numeric delta = abs(m_coefficient[0] - b.m_coefficient[0]);
        for (unsigned int i = 1; i < m_coefficient.size(); ++i)
            delta = std::max<Numeric>(delta, abs(m_coefficient[i] - b.m_coefficient[i]));
        return delta;
    }

    /// Create a new polynomial shifted upwards by adding v to all
    /// coefficients.
    PolynomialBezier operator+ (const Numeric& v) const
    {
        std::vector<Numeric> newCoeff(m_coefficient.begin(), m_coefficient.end());
        for (unsigned int i = 0; i < newCoeff.size(); ++i)
            newCoeff[i] += v;
        return PolynomialBezier(newCoeff);
    }

    /// Create a new polynomial shifted downwards by subtracting v from all
    /// coefficients.
    PolynomialBezier operator- (const Numeric& v) const
    {
        std::vector<Numeric> newCoeff(m_coefficient.begin(), m_coefficient.end());
        for (unsigned int i = 0; i < newCoeff.size(); ++i)
            newCoeff[i] -= v;
        return PolynomialBezier(newCoeff);
    }
};

/// Convert the polynomial from monomial representation to the
/// Bernstein-Bezier basis in the unit interval. Use scale() to transform the
/// coefficients before applying this base conversion.
template <typename Numeric>
PolynomialBezier<Numeric> PolynomialStandard<Numeric>::toBezier() const
{
    const unsigned int N = m_coefficient.size()-1;

    std::vector<Numeric> coeff(N+1);

    for (unsigned int i = 0; i <= N; ++i)
    {
        coeff[i] = 0;

        for (unsigned int j = 0; j <= i; ++j)
        {
            coeff[i] += Numeric::binomial(i, j) / Numeric::binomial(N, i - j) * m_coefficient[i - j];
        }
    }

    return PolynomialBezier<Numeric>(coeff);
}

/// Transform the polynomial from Bernstein-Bezier basis back to the standard
/// monomial basis.
template <typename Numeric>
PolynomialStandard<Numeric> PolynomialBezier<Numeric>::toStandard() const
{
    const unsigned int N = m_coefficient.size()-1;

    std::vector<Numeric> coeff(N+1);

    for (unsigned int i = 0; i <= N; ++i)
    {
        coeff[i] = 0;

        for (unsigned int j = 0; j <= i; ++j)
        {
            coeff[i] += (j % 2 ? -1 : +1) * Numeric::binomial(i,j) * m_coefficient[i - j];
        }

        coeff[i] *= Numeric::binomial(N,i);
    }

    return PolynomialStandard<Numeric>(coeff);
}

/** \page matrixdegree Polynomial Degree Reduction and Raising

*******************************************************************************
*** Polynomial Degree Reduction and Raising                                 ***
*******************************************************************************

 The most complex step in the KClip algorithms is reduction of the degree of
 the input polynomial. For this a Matrix can be precomputed using equations
 found in [1] and later only matrix multiplication is needed for reducing or
 raising a polynomial in Bezier representation.

 The actual formulas are implemented in "DegreeInterpolationGenerator" which
 can be used to fill() a Matrix object of the corresponding size. The
 Generator object is implemented as a functional which calculates a matrix
 entry. See the code of KClip and testsuite() for details.

******************************************************************************/

/// Implements a simple matrix that can be filled using a functional object.
template <typename Numeric>
class Matrix
{
private:

    /// number of columns in the matrix
    unsigned int        m_cols;

    /// matrix data
    std::vector<Numeric> m_data;

public:

    Matrix(unsigned int rows, unsigned int cols)
        : m_cols(cols)
    {
        m_data.resize(rows * cols);
    }

    unsigned int cols() const
    {
        return m_cols;
    }

    unsigned int rows() const
    {
        return m_data.size() / m_cols;
    }

    /// Return value at a position of the matrix.
    inline Numeric& at(unsigned int row, unsigned int col)
    {
        assert( row < rows() );
        assert( col < cols() );
        return m_data[row * m_cols + col];
    }

    /// Return value at a position of the matrix.
    inline const Numeric& at(unsigned int row, unsigned int col) const
    {
        assert( row < rows() );
        assert( col < cols() );
        return m_data[row * m_cols + col];
    }

    /// Fill in whole matrix from a function calculating values.
    template <typename Functional>
    Matrix& fill(Functional func)
    {
        for (unsigned int i = 0; i < rows(); ++i)
            for (unsigned int j = 0; j < cols(); ++j)
                at(i,j) = func(i,j);

        return *this;
    }

    /// Return the matrix as a latex pmatrix math string.
    std::string toString() const
    {
        std::ostringstream oss;
        oss << "\\begin{pmatrix}\n";
        for (unsigned int i = 0; i < rows(); ++i)
        {
            for (unsigned int j = 0; j < cols(); ++j)
            {
                if (j != 0) oss << " & ";
                oss << ltxNum(at(i,j));
            }
            oss << " \\\\\n";
        }
        oss << "\\end{pmatrix}";
        return oss.str();
    }
};

/// Functional class to calculate the entries of the degree reducing or
/// raising matrices.
template <typename Numeric>
class DegreeInterpolationGenerator
{
public:

    /// Calculate the scalar produkt of the two Bernstein polynomials
    /// B_{i,m} and B_{j,n} in the unit interval according to (13) in [1]
    static Numeric scalarProduktBernsteinPolynomial(int i, int m, int j, int n)
    {
        assert(i <= m); assert(j <= n);
        return Numeric( Numeric::binomial(m,i) * Numeric::binomial(n,j) ) / Numeric( Numeric(m + n + 1) * Numeric::binomial(m + n, i + j) );
    }

    /// Calculate the coefficients c_{p,q} of the dual basis \f$D_{p}^k\f$
    /// according to (10) in [1].
    static Numeric CoefficientDualBasis(int p, int q, int k)
    {
        Numeric a = Numeric( ((p+q) % 2 ? -1 : +1) ) / Numeric( Numeric::binomial(k,p) * Numeric::binomial(k,q) );
        Numeric b = 0;

        for (int j = 0; j <= std::min(p,q); ++j)
            b += Numeric(2*j+1) * Numeric::binomial(k+j+1,k-p) * Numeric::binomial(k-j,k-p) * Numeric::binomial(k+j+1,k-q) * Numeric::binomial(k-j,k-q);

        return a * b;
    }

    /// Calculate the coefficients \f$\beta_{i,j}^{n,k}\f$ of the
    /// interpolation matrix according to (12) in [1].
    static Numeric CoefficientBeta(int i, int j, int n, int k)
    {
        Numeric a = 0;
        for (int l = 0; l <= k; ++l)
        {
            a += CoefficientDualBasis(j,l,k) * scalarProduktBernsteinPolynomial(i,n,l,k);
        }
        return a;
    }

public:
    /// Dimensions of degree interpolation
    unsigned int        m_n, m_k;

    /// Constructor saving the degree parameters.
    DegreeInterpolationGenerator(unsigned int n, unsigned int k)
        : m_n(n), m_k(k)
    {
    }

    /// Functional operator for Matrix::fill().
    Numeric operator()(unsigned int i, unsigned int j) const
    {
        return CoefficientBeta(i,j,m_n,m_k);
    }
};

/// Apply a degree raising of reducing matrix to the given Bezier polynomial
/// via simple matrix multiplication.
template <typename Numeric>
PolynomialBezier<Numeric> PolynomialBezier<Numeric>::interpolateDegree(const Matrix<Numeric> &m) const
{
    const unsigned int N = m_coefficient.size()-1;

    assert(m.rows() == N+1);
    const unsigned int k = m.cols()-1;

    std::vector<Numeric> newCoeff (k+1);

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

        for (unsigned int j = 0; j <= N; ++j)
        {
            newCoeff[i] += m_coefficient[j] * m.at(j,i);
        }
    }

    return PolynomialBezier(newCoeff);
}

/// Apply a degree rasing of reducing matrix, where the matrix has to be
/// constructed before application.
template <typename Numeric>
PolynomialBezier<Numeric> PolynomialBezier<Numeric>::interpolateDegree(int k) const
{
    const int N = m_coefficient.size()-1;

    Matrix<Numeric> m(N+1, k+1);
    m.fill( DegreeInterpolationGenerator<Numeric>(N,k) );

    return interpolateDegree(m);
}

/** \page bezclip The Bezier Clipping Algorithm

*******************************************************************************
*** The Bezier Clipping Algorithm                                           ***
*******************************************************************************

 Using the base algorithms implemented in PolynomialBezier the class
 BezierClip implements root finding by intersecting the convex hull of the
 Bezier polygon with the x-axis.

******************************************************************************/

/// Implements the BezierClip algorithm
template <typename Numeric>
class BezierClip
{
private:

    /// roots found
    std::vector< std::pair<Numeric,Numeric> >   m_roots;

    /// epsilon used
    Numeric             m_epsilon;

    /// prefix for latex labels
    std::string         m_markprefix;

    /// maximum recursion depth
    unsigned int        m_maxdepth;

    typedef PolynomialStandard<Numeric> NPolynomialStandard;
    typedef PolynomialBezier<Numeric>   NPolynomialBezier;

public:

    BezierClip(Numeric epsilon, const std::string& markprefix = "r")
        : m_epsilon(epsilon),
          m_markprefix(markprefix)
    {
    }

    unsigned int maxdepth() const
    { return m_maxdepth; }

    std::vector< std::pair<Numeric,Numeric> > findRoots(const NPolynomialStandard& p1, Numeric left, Numeric right)
    {
        m_roots.clear();
        m_maxdepth = 0;

        dbg("Running BezClip on p = " << p1.toString() << " within [" << left << "," << right << "]\n");
        ltx("\\begin{center}\\bf\\Large $" << p1.toString() << "$\\end{center}\n");

        // *** Input Polynomial

        ltx("\\textbf{Called \\texttt{BezClip} with input polynomial on interval " << ltxInterval(left,right) << ":}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toString() << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}[\n"
            "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
            "       samples=100]\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        ltx("\\subsection{Recursion Branch 1 for Input Interval " << ltxInterval(left,right) << "}\n\n");

        NPolynomialStandard p2 = p1.scale(left, right);

        findRootsRecursive(p2.toBezier(), left, right, "1", 1);

        // *** Result

        ltx("\\clearpage\n");

        ltx("\\subsection{Result: " << m_roots.size() << " Root Intervals}\n");
        ltx("\\label{" << m_markprefix << "result}\n");

        ltx("\\textbf{Input Polynomial on Interval $[" << left << "," << right << "]$}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toString() << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}[\n"
            "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
            "       samples=100]\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        ltx("\\textbf{Result: Root Intervals}\n");

        ltx("\\begin{center}\n");
        for (unsigned int i = 0; i < m_roots.size(); ++i)
        {
            if (i != 0) ltx(", ");
            ltx("$[" << m_roots[i].first << "," << m_roots[i].second << "]$");
        }
        ltx("\\end{center}\n");

        ltx("with precision $\\varepsilon = " << ltxNum(m_epsilon) << "$.\n\n");

        return m_roots;
    }

    // The polynomial p1 is here already normed to [0:1], left and right are
    // the _original_ interval and only needed for scaling the found roots.
    void findRootsRecursive(const NPolynomialBezier& p1, Numeric left, Numeric right, const std::string& mark, unsigned int depth)
    {
        if (right - left < m_epsilon)
        {
            if (p1.coefficient(0) == p1.coefficient(p1.degree()) ||
                (p1.coefficient(0) > 0.0 && p1.coefficient(p1.degree()) < 0.0) ||
                (p1.coefficient(0) < 0.0 && p1.coefficient(p1.degree()) > 0.0) ||
                (abs(p1.coefficient(0)) < m_epsilon && abs(p1.coefficient(p1.degree())) < m_epsilon)
                )
            {
                dbg("Found root in interval [" << left << "," << right << "] at recursion depth " << depth << "!\n");
                ltx("Found root in interval $[" << left << "," << right << "]$ at recursion depth " << depth << "!\n");
                m_roots.push_back( std::make_pair(left,right) );
            }
            else {
                dbg("Reached interval [" << left << "," << right << "] without sign change at depth " << depth << "!\n");
                ltx("Reached interval $[" << left << "," << right << "]$ \\textbf{without sign change} at depth " << depth << "!\n\n");
            }
            return;
        }

        m_maxdepth = std::max(m_maxdepth, depth);

        // *** Normalization on [0:1] ***

        ltx("\\textbf{Normalized monomial und B\\'ezier representations and the B\\'ezier polygon:}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toStandard().toString() << "\n"
            "  = " << p1.toString("(X)") << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");

        ltx("\\addplot[blue,dashed,mark=*] coordinates { " << p1.toPolygon() << " };\n");

        ltx("\\addplot[blue,fill=blue!50!white,opacity=0.2] coordinates { " << p1.toConvexHullString() << " };\n");

        std::vector<Numeric> chi = p1.getConvexHullIntersection();
        if (chi.size() == 2) {
            ltx("\\addplot[red,ultra thick] coordinates { (" << chi[0] << ",0) (" << chi[1] << ",0) };\n");
        }

        ltx("\\addplot[forget plot] coordinates { (0,0) };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        // *** Intersection of the Bezier hull ***

        ltx("\\textbf{Intersection of the convex hull with the $x$ axis:}\n");

        ltx("\\begin{dmath*}\n \\{ ");
        for (unsigned int i = 0; i < chi.size(); ++i) {
            if (i != 0) ltx(", "); ltx(chi[i]);
            if (i != 0) assert(chi[i-1] < chi[i]);
        }
        ltx("\\} \n\\end{dmath*}\n");

        std::vector< std::pair<Numeric,Numeric> > intervals;

        if (chi.size() == 2)
            intervals.push_back( std::make_pair(chi[0], chi[1]) );

        ltx("\\textbf{Intersection intervals with the $x$ axis:}\n\n");

        if (intervals.size() == 0) {
            ltx("No intersection with the $x$ axis. Done.\n\n");
            return;
        }

        ltx("\\begin{dmath*}\n");
        for (unsigned int i = 0; i < intervals.size(); ++i)
        {
            if (i != 0) ltx(" && ");
            ltx("[" << intervals[i].first << "," << intervals[i].second << "] ");
        }
        ltx("\\end{dmath*}\n");

        // find maximum interval size

        Numeric maxInterval = intervals[0].second - intervals[0].first;
        for (unsigned int i = 1; i < intervals.size(); ++i)
        {
            maxInterval = std::max(maxInterval, intervals[i].second - intervals[i].first);
        }

        ltx("Longest intersection interval: $" << ltxNum(maxInterval) << "$\n\n");

        if (maxInterval > 0.5)
        {
            Numeric middle = left + (right - left) / 2.0;

            ltx(" $\\Longrightarrow$ Bisection: \\hyperref[" << m_markprefix << mark << " 1" << "]{first half " << ltxInterval(left,middle) << "} und \\hyperref[" << m_markprefix << mark << " 2" << "]{second half " << ltxInterval(middle,right) << "}\n\n");

            std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(0.5);

            ltx("\\subsection{Recursion Branch " << mark << " 1 on the First Half " << ltxInterval(left,middle) << "}\n");
            ltx("\\label{" << m_markprefix << mark << " 1}\n");

            findRootsRecursive(polySplit.first, left, middle, mark + " 1", depth+1);

            ltx("\\subsection{Recursion Branch " << mark << " 2 on the Second Half " << ltxInterval(middle,right) << "}\n");
            ltx("\\label{" << m_markprefix << mark << " 2}\n");

            findRootsRecursive(polySplit.second, middle, right, mark + " 2", depth+1);
        }
        else
        {
            ltx(" $\\Longrightarrow$ Selective recursion: ");
            for (unsigned int i = 0; i < intervals.size(); ++i)
            {
#if !SPEEDTEST
                Numeric ileft = left + (right - left) * intervals[i].first;
                Numeric iright = left + (right - left) * intervals[i].second;

                ltx("\\hyperref[" << m_markprefix << mark << " " << i+1 << "]{interval " << i+1 << ": " << ltxInterval(ileft, iright) << "}, ");
#endif
            }
            ltx("\n\n");

            for (unsigned int i = 0; i < intervals.size(); ++i)
            {
                Numeric ileft = left + (right - left) * intervals[i].first;
                Numeric iright = left + (right - left) * intervals[i].second;

                if (intervals[i].first == 0)
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].second);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit.first, ileft, iright, mark + " " + S(i+1), depth+1);
                }
                else if (intervals[i].second == 1.0)
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].first);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit.second, ileft, iright, mark + " " + S(i+1), depth+1);
                }
                else
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit1 = p1.deCasteljauSplit(intervals[i].first);

                    Numeric split2 = (intervals[i].second - intervals[i].first) / (1.0 - intervals[i].first);

                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit2 = polySplit1.second.deCasteljauSplit(split2);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit2.first, ileft, iright, mark + " " + S(i+1), depth+1);
                }
            }
        }
    }
};

/** \page kclip The Quadclip and K-Clip Algorithms

*******************************************************************************
*** The Quadclip and K-Clip Algorithms                                      ***
*******************************************************************************

 Using the base algorithms implemented in PolynomialBezier and the degree
 reduction and raising classes KClip implements the algorithms QuadClip and
 CubeClip (for parameters k=2 and k=3).

 See [1] and [2] for a description and analysis of the algorithms.

******************************************************************************/

/// Implements the K-Clip algorithm
template <typename Numeric>
class KClip
{
private:

    /// degree of polynomial which this class processes
    unsigned int        m_degree;

    /// degree of reduced polynomial and clipping process
    unsigned int        m_clip;

    /// constant degree reducing matrix
    Matrix<Numeric>     m_degreeReducing;

    /// constant degree raising matrix
    Matrix<Numeric>     m_degreeRaising;

    /// roots found
    std::vector< std::pair<Numeric,Numeric> >   m_roots;

    /// epsilon used
    Numeric             m_epsilon;

    /// prefix for latex labels
    std::string         m_markprefix;

    /// maximum recursion depth
    unsigned int        m_maxdepth;

    typedef PolynomialStandard<Numeric> NPolynomialStandard;
    typedef PolynomialBezier<Numeric>   NPolynomialBezier;

public:

    KClip(unsigned int degree, unsigned int clip = 2, Numeric epsilon = 0.001, const std::string& markprefix = "r")
        : m_degree(degree),
          m_clip(clip),
          m_degreeReducing(degree+1, clip+1),
          m_degreeRaising(clip+1, degree+1),
          m_epsilon(epsilon),
          m_markprefix(markprefix)
    {
        // construct constant matrices
        m_degreeReducing.fill( DegreeInterpolationGenerator<Numeric>(degree,clip) );

        m_degreeRaising.fill( DegreeInterpolationGenerator<Numeric>(clip,degree) );
    }

    unsigned int maxdepth() const
    { return m_maxdepth; }

    std::vector< std::pair<Numeric,Numeric> > findRoots(NPolynomialStandard p1, Numeric left, Numeric right)
    {
        m_roots.clear();
        m_maxdepth = 0;

        dbg("Running " << (m_clip == 2 ? "QuadClip" : "CubeClip") << " on p = " << p1.toString() << " in interval [" << left << "," << right << "]\n");
        ltx("\\begin{center}\\bf\\Large $" << p1.toString() << "$\\end{center}\n");

        // *** Input Polynom

        ltx("\\textbf{Called \\texttt{" << (m_clip == 2 ? "QuadClip" : "CubeClip") << "} with input polynomial on interval " << ltxInterval(left,right) << ":}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toString() << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}[\n"
            "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
            "       samples=100]\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        ltx("\\subsection{Recursion Branch 1 for Input Interval " << ltxInterval(left,right) << "}\n\n");

        if (left != 0.0 || right != 0.0) {
            p1 = p1.scale(left, right);
        }

        findRootsRecursive(p1.toBezier(), left, right, "1", 1);

        // *** Result

        ltx("\\clearpage\n");

        ltx("\\subsection{Result: " << m_roots.size() << " Root Intervals}\n");
        ltx("\\label{" << m_markprefix << "result}\n");

        ltx("\\textbf{Input Polynomial on Interval $[" << left << "," << right << "]$}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toString() << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}[\n"
            "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
            "       samples=100]\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        ltx("\\textbf{Result: Root Intervals}\n");

        ltx("\\begin{center}\n");
        for (unsigned int i = 0; i < m_roots.size(); ++i)
        {
            if (i != 0) ltx(", ");
            ltx("$[" << m_roots[i].first << "," << m_roots[i].second << "]$");
        }
        ltx("\\end{center}\n");

        ltx("with precision $\\varepsilon = " << ltxNum(m_epsilon) << "$.\n\n");

        return m_roots;
    }

    // the polynomial p1 is here already normed to [0:1], left and right are the
    // _original_ interval and only needed for scaling the found roots.
    void findRootsRecursive(const NPolynomialBezier& p1, Numeric left, Numeric right, const std::string& mark, unsigned int depth)
    {
        if (right - left < m_epsilon)
        {
            if (p1.coefficient(0) == p1.coefficient(p1.degree()) ||
                (p1.coefficient(0) > 0.0 && p1.coefficient(p1.degree()) < 0.0) ||
                (p1.coefficient(0) < 0.0 && p1.coefficient(p1.degree()) > 0.0) ||
                (abs(p1.coefficient(0)) < m_epsilon && abs(p1.coefficient(p1.degree())) < m_epsilon)
                )
            {
                dbg("Found root in interval [" << left << "," << right << "] at recursion depth " << depth << "!\n");
                ltx("Found root in interval $[" << left << "," << right << "]$ at recursion depth " << depth << "!\n");
                m_roots.push_back( std::make_pair(left,right) );
            }
            else {
                dbg("Reached interval [" << left << "," << right << "] without sign change at depth " << depth << "!\n");
                ltx("Reached interval $[" << left << "," << right << "]$ \\textbf{without sign change} at depth " << depth << "!\n\n");

                ltx("p(0) = " << p1.coefficient(0) << " - p(1) " << p1.coefficient(p1.degree()) << "\n\n");

            }
            return;
        }

        m_maxdepth = std::max(m_maxdepth, depth);

        // *** Normalization on [0:1] ***

        ltx("\\textbf{Normalized monomial und B\\'ezier representations and the B\\'ezier polygon:}\n"
            "\\begin{dmath*}\n"
            "p = " << p1.toStandard().toString() << "\n"
            "  = " << p1.toString("(X)") << "\n"
            "\\end{dmath*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");

        ltx("\\addplot[blue,dashed,mark=*] coordinates { " << p1.toPolygon() << " };\n");

        ltx("\\addplot[forget plot] coordinates { (0,0) };\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

#if !SPEEDTEST
        if (debug_more)
        {
            // *** Reduced Polynomials ***

            ltx("\\textbf{Best approximation polynomials of degree 0, 1, 2, 3 and 4:}\n");

            ltx("\\begin{dgroup*}\n"
                "\\begin{dmath*}\n q_0 = " << p1.interpolateDegree(0).toStandard().toString() << " = " << p1.interpolateDegree(0).toString() << " \\end{dmath*}\n"
                "\\begin{dmath*}\n q_1 = " << p1.interpolateDegree(1).toStandard().toString() << " = " << p1.interpolateDegree(1).toString() << " \\end{dmath*}\n"
                "\\begin{dmath*}\n q_2 = " << p1.interpolateDegree(2).toStandard().toString() << " = " << p1.interpolateDegree(2).toString() << " \\end{dmath*}\n"
                "\\begin{dmath*}\n q_3 = " << p1.interpolateDegree(3).toStandard().toString() << " = " << p1.interpolateDegree(3).toString() << " \\end{dmath*}\n"
                "\\begin{dmath*}\n q_4 = " << p1.interpolateDegree(4).toStandard().toString() << " = " << p1.interpolateDegree(4).toString() << " \\end{dmath*}\n"
                "\\end{dgroup*}\n\n");

            ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
                "\\begin{axis}\n");

            ltx("\\addplot[purple] gnuplot { " << p1.interpolateDegree(0).toStandard().toStringPlot() << " };\n");
            ltx("\\addplot[cyan] gnuplot { " << p1.interpolateDegree(1).toStandard().toStringPlot() << " };\n");
            ltx("\\addplot[red] gnuplot { " << p1.interpolateDegree(2).toStandard().toStringPlot() << " };\n");
            ltx("\\addplot[green!60!black] gnuplot { " << p1.interpolateDegree(3).toStandard().toStringPlot() << " };\n");
            ltx("\\addplot[orange] gnuplot { " << p1.interpolateDegree(4).toStandard().toStringPlot() << " };\n");

            ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");

            ltx("\\addplot[forget plot] coordinates { (0,0) };\n");

            ltx("\\legend{$q_0$,$q_1$,$q_2$,$q_3$,$q_4$,$p$};\n");

            ltx("\\end{axis}\n"
                "\\end{tikzpicture}\n\\fi\\end{center}\n");
        }

        if (debug_more)
        {
            // *** Matrices ***

            ltx("\\textbf{Degree reduction and raising matrices:}\n");

            ltx("\\begin{align*}\n"
                "M_{" << m_degree << "," << m_clip << "} = " << m_degreeReducing.toString() << " && "
                "M_{" << m_clip << "," << m_degree << "} = " << m_degreeRaising.toString() << "\n"
                "\\end{align*}\n");
        }
#endif

        // *** Interpolation ***

        ltx("\\textbf{Degree reduction and raising:}\n");

        NPolynomialBezier p2 = p1.interpolateDegree(m_degreeReducing);

        NPolynomialBezier p3 = p2.interpolateDegree(m_degreeRaising);

        ltx("\\begin{dgroup*}\n"
            "\\begin{dmath*}\n q_" << m_clip << " = " << p2.toStandard().toString() << "\\\\\n"
            "       = " << p2.toString() << " \\end{dmath*}\n"
            "\\begin{dmath*}\n \\widetilde{q_" << m_clip << "} = " << p3.toStandard().toString() << "\\\\\n"
            "       = " << p3.toString() << " \\end{dmath*}\n"
            "\\end{dgroup*}\n\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}\n");

        ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
        ltx("\\addplot[blue,dashed,mark=*,forget plot] coordinates { " << p1.toPolygon() << " };\n");

        ltx("\\addplot[orange,thick] gnuplot { " << p2.toStandard().toStringPlot() << " };\n");
        ltx("\\addplot[red,dashed,mark=*,forget plot] coordinates { " << p3.toPolygon() << " };\n");
        ltx("\\addplot[orange,dashed,mark=*] coordinates { " << p2.toPolygon() << " };\n");

        for (unsigned int i = 0; i <= m_degree; ++i)
        {
            ltx("\\addplot[green!60!black,forget plot,ultra thick] coordinates { "
                "(" << i*(1.0/m_degree) << "," << p1.coefficient(i) << ") "
                "(" << i*(1.0/m_degree) << "," << p3.coefficient(i) << ") };\n");
        }

        ltx("\\addplot[forget plot] coordinates { (0,0) };\n");

        ltx("\\legend{$p$,$q_" << m_clip << "$,$\\widetilde{q_" << m_clip << "}$};\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        Numeric delta = p3.maxCoefficientDeltaTo(p1);

        ltx("The maximum difference of the B\\'ezier coefficients is $\\delta = " << ltxNum(delta) << "$.\n\n");

        // *** Bounding Polynomials M and m ***

        ltx("\\textbf{Bounding polynomials $M$ and $m$:}\n");

        NPolynomialBezier p2M = p2 + delta;
        NPolynomialBezier p2m = p2 - delta;

        ltx("\\begin{dgroup*}\n"
            "\\begin{dmath*}\n M = " << p2M.toStandard().toString() << " \\end{dmath*}\n"
            "\\begin{dmath*}\n m = " << p2m.toStandard().toString() << " \\end{dmath*}\n"
            "\\end{dgroup*}\n");

        // *** Root ***

        ltx("\\textbf{Root of $M$ and $m$:}\n");

        std::vector<Numeric> rootsM = p2M.findRoots();
        std::vector<Numeric> rootsm = p2m.findRoots();

        ltx("\\begin{align*}\n N(M) = \\{ ");
        for (unsigned int i = 0; i < rootsM.size(); ++i) {
            if (i != 0) ltx(", ");
            ltx(ltxNum(rootsM[i]));
            if (i != 0) assert(rootsM[i-1] <= rootsM[i]);
        }
        ltx("\\} &&");

        ltx("N(m) = \\{ ");
        for (unsigned int i = 0; i < rootsm.size(); ++i) {
            if (i != 0) ltx(", ");
            ltx(ltxNum(rootsm[i]));
            if (i != 0) assert(rootsm[i-1] <= rootsm[i]);
        }
        ltx("\\}\n\\end{align*}\n");

        std::vector< std::pair<Numeric,Numeric> > intervals;

        if (m_clip == 2)
        {
            // *** Quadratic Clipping

            if (rootsM.size() == 0 && rootsm.size() == 0)
            {
                // no intervals
            }
            else if (rootsM.size() == 1 && rootsm.size() == 0)
            {
                // no intervals -- but a double root?
                assert(0);
            }
            else if (rootsM.size() == 2 && rootsm.size() == 0)
            {
                intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
                                                    std::min<Numeric>(1,rootsM[1])) );
            }
            else if (rootsM.size() == 0 && rootsm.size() == 1)
            {
                // no intervals -- but a double root?
                assert(0);
            }
            else if (rootsM.size() == 0 && rootsm.size() == 2)
            {
                intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
                                                    std::min<Numeric>(1,rootsm[1])) );
            }
            else if (rootsM.size() == 1 && rootsm.size() == 1)
            {
                // impossible.
                assert(0);
            }
            else if (rootsM.size() == 2 && rootsm.size() == 1)
            {
                intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
                                                    std::min<Numeric>(1,rootsM[1])) );
            }
            else if (rootsM.size() == 1 && rootsm.size() == 2)
            {
                intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
                                                    std::min<Numeric>(1,rootsm[1])) );
            }
            else if (rootsM.size() == 2 && rootsm.size() == 2)
            {
                if (rootsM[0] > rootsm[0])
                {
                    if (rootsM[0] > 0 && rootsm[0] < 1) {
                        intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
                                                            std::min<Numeric>(1,rootsM[0])) );
                    }
                    if (rootsM[1] < 1 && rootsm[1] > 0) {
                        intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[1]),
                                                            std::min<Numeric>(1,rootsm[1])) );
                    }
                }
                else
                {
                    if (rootsm[0] > 0 && rootsM[0] < 1) {
                        intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
                                                            std::min<Numeric>(1,rootsm[0])) );
                    }
                    if (rootsm[1] < 1 && rootsM[1] > 0) {
                        intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[1]),
                                                            std::min<Numeric>(1,rootsM[1])) );
                    }
                }
            }
            else {
                assert("Invalid combination of roots");
            }
        }
        else if (m_clip == 3)
        {
            // *** Cubic Clipping

            if (rootsM.size() == 1 && rootsm.size() == 1)
            {
                if (rootsM[0] < rootsm[0])
                {
                    intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
                }
                else // rootsM[0] > rootsm[0]
                {
                    intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
                }
            }
            else if (rootsM.size() == 1 && rootsm.size() == 2)
            {
                assert(0); // unknown
            }
            else if (rootsM.size() == 2 && rootsm.size() == 1)
            {
                assert(0); // unknown
            }
            else if (rootsM.size() == 2 && rootsm.size() == 2)
            {
                if (rootsM[0] < rootsm[0])
                {
                    intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
                    intervals.push_back( std::make_pair(rootsM[0], rootsm[1]) );
                    intervals.push_back( std::make_pair(rootsm[1], rootsM[1]) );
                }
                else // rootsM[0] > rootsm[0]
                {
                    intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
                    intervals.push_back( std::make_pair(rootsm[0], rootsM[1]) );
                    intervals.push_back( std::make_pair(rootsM[1], rootsm[1]) );
                }
            }
            else if (rootsM.size() == 1 && rootsm.size() == 3)
            {
                if (rootsM[0] < rootsm[0])
                {
                    intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
                    intervals.push_back( std::make_pair(rootsm[1], rootsm[2]) );
                }
                else // rootsM[0] > rootsm[0]
                {
                    intervals.push_back( std::make_pair(rootsm[0], rootsm[1]) );
                    intervals.push_back( std::make_pair(rootsm[2], rootsM[0]) );
                }
            }
            else if (rootsM.size() == 3 && rootsm.size() == 1)
            {
                if (rootsM[0] < rootsm[0])
                {
                    intervals.push_back( std::make_pair(rootsM[0], rootsM[1]) );
                    intervals.push_back( std::make_pair(rootsM[2], rootsm[0]) );
                }
                else // rootsM[0] > rootsm[0]
                {
                    intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
                    intervals.push_back( std::make_pair(rootsM[1], rootsM[2]) );
                }
            }
            else if (rootsM.size() == 2 && rootsm.size() == 3)
            {
                assert(0); // unknown
            }
            else if (rootsM.size() == 3 && rootsm.size() == 2)
            {
                assert(0); // unknown
            }
            else if (rootsM.size() == 3 && rootsm.size() == 3)
            {
                if (rootsM[0] < rootsm[0])
                {
                    intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
                    intervals.push_back( std::make_pair(rootsm[1], rootsM[1]) );
                    intervals.push_back( std::make_pair(rootsM[2], rootsm[2]) );
                }
                else // rootsM[0] > rootsm[0]
                {
                    intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
                    intervals.push_back( std::make_pair(rootsM[1], rootsm[1]) );
                    intervals.push_back( std::make_pair(rootsm[2], rootsM[2]) );
                }
            }
            else {
                assert(0); // invalid combination
            }

            std::vector< std::pair<Numeric,Numeric> > newintervals;

            for (unsigned int i = 0; i < intervals.size(); ++i)
            {
                if (intervals[i].first > intervals[i].second)
                {
                    std::swap(intervals[i].first, intervals[i].second);
                }
                assert( intervals[i].first <= intervals[i].second);

                if (intervals[i].second < 0.0) continue;
                if (intervals[i].first > 1.0) continue;

                newintervals.push_back( std::make_pair(std::max<Numeric>(0.0, intervals[i].first),
                                                       std::min<Numeric>(1.0, intervals[i].second)) );
            }

            std::swap(intervals, newintervals);
        }
        else {
            assert("Invalid combination of roots");
        }

        ltx("\\textbf{Intersection intervals:}\n");

        ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
            "\\begin{axis}\n");

        ltx("\\addplot[blue,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");

        ltx("\\addplot[green!60!black] gnuplot { " << p2M.toStandard().toStringPlot() << " };\n");
        ltx("\\addplot[green!60!black] gnuplot { " << p2m.toStandard().toStringPlot() << " };\n");

        for (unsigned int i = 0; i < intervals.size(); ++i)
            ltx("\\addplot[red,ultra thick,forget plot] coordinates { (" << intervals[i].first << ",0) (" << intervals[i].second << ",0) };\n");

        ltx("\\addplot[forget plot] coordinates { (0,0) };\n");

        //ltx("\\legend{$p$,$M$,$m$};\n");

        ltx("\\end{axis}\n"
            "\\end{tikzpicture}\n\\fi\\end{center}\n");

        if (intervals.size() == 0) {
            ltx("No intersection intervals with the $x$ axis.\n\n");
            return;
        }

        ltx("\\begin{center}\n");
        for (unsigned int i = 0; i < intervals.size(); ++i)
        {
            if (i != 0) ltx(", ");
            ltx("$[" << intervals[i].first << "," << intervals[i].second << "]$");
        }
        ltx("\\end{center}\n");

        // find maximum interval size

        Numeric maxInterval = intervals[0].second - intervals[0].first;
        for (unsigned int i = 1; i < intervals.size(); ++i)
        {
            maxInterval = std::max(maxInterval, intervals[i].second - intervals[i].first);
        }

        ltx("Longest intersection interval: $" << ltxNum(maxInterval) << "$\n\n");

        if (maxInterval > 0.5)
        {
            Numeric middle = left + (right - left) / 2.0;

            ltx(" $\\Longrightarrow$ Bisection: \\hyperref[" << m_markprefix << mark << " 1" << "]{first half " << ltxInterval(left,middle) << "} und \\hyperref[" << m_markprefix << mark << " 2" << "]{second half " << ltxInterval(middle,right) << "}\n\n");

            std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(0.5);

            if (polySplit.second.coefficient(0) < m_epsilon)
            {
                ltx("Bisection point is very near to a root?!?\n\n");
            }

            ltx("\\subsection{Recursion Branch " << mark << " 1 on the First Half " << ltxInterval(left,middle) << "}\n");
            ltx("\\label{" << m_markprefix << mark << " 1}\n");

            findRootsRecursive(polySplit.first, left, middle, mark + " 1", depth+1);

            ltx("\\subsection{Recursion Branch " << mark << " 2 on the Second Half " << ltxInterval(middle,right) << "}\n");
            ltx("\\label{" << m_markprefix << mark << " 2}\n");

            findRootsRecursive(polySplit.second, middle, right, mark + " 2", depth+1);
        }
        else
        {
            ltx(" $\\Longrightarrow$ Selective recursion: ");
            for (unsigned int i = 0; i < intervals.size(); ++i)
            {
#if !SPEEDTEST
                Numeric ileft = left + (right - left) * intervals[i].first;
                Numeric iright = left + (right - left) * intervals[i].second;

                ltx("\\hyperref[" << m_markprefix << mark << " " << i+1 << "]{interval " << i+1 << ": " << ltxInterval(ileft, iright) << "}, ");
#endif
            }
            ltx("\n\n");

            for (unsigned int i = 0; i < intervals.size(); ++i)
            {
                Numeric ileft = left + (right - left) * intervals[i].first;
                Numeric iright = left + (right - left) * intervals[i].second;

                if (intervals[i].first == 0) // only one de Casteljau split required
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].second);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit.first, ileft, iright, mark + " " + S(i+1), depth+1);
                }
                else if (intervals[i].second == 1.0) // only one de Casteljau split required
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].first);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit.second, ileft, iright, mark + " " + S(i+1), depth+1);
                }
                else // requires two de Casteljau splits
                {
                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit1 = p1.deCasteljauSplit(intervals[i].first);

                    Numeric split2 = (intervals[i].second - intervals[i].first) / (1.0 - intervals[i].first);

                    std::pair<NPolynomialBezier,NPolynomialBezier> polySplit2 = polySplit1.second.deCasteljauSplit(split2);

                    ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
                    ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");

                    findRootsRecursive(polySplit2.first, ileft, iright, mark + " " + S(i+1), depth+1);
                }
            }
        }
    }
};

/** \page testsuite A Suite of Test Cases

*******************************************************************************
*** A Suite of Test Cases                                                   ***
*******************************************************************************

 In the testsuite() function many of the basic algorithm implemented above are
 checked against precomputed values.

******************************************************************************/

#define ASSERT(x)                                                \
    do { if (!(x)) {                                             \
        std::cout << "Testsuite assertion \"" #x "\" on line "   \
                  << __LINE__ << " failed!\n";                   \
        abort();                                                 \
        }                                                        \
    } while(0)

#define CHKSTR(p,str)                                           \
    do { if (p.toString() != str) {                             \
        std::cout << "Testsuite error: "                        \
                  << p.toString() << " != " << str << "\n";     \
        ASSERT(p.toString() == str);                            \
        }                                                       \
    } while(0)

template <typename Numeric>
void testsuite()
{
    typedef PolynomialStandard<Numeric> NPolynomialStandard;
    typedef PolynomialBezier<Numeric> NPolynomialBezier;

    { // test Bezier transformation, evaluation and inverse transformation

        Numeric coeff[6] = { 1, 0, 0, 0, -16, 16 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[6]));

        CHKSTR( p1, "16 X^5 - 16 X^4 + 1" );

        ASSERT( abs(p1.evaluate(0.0) - 1.0) < 0.00001 );
        ASSERT( abs(p1.evaluate(0.25) - 0.953125) < 0.00001 );
        ASSERT( abs(p1.evaluate(0.5) - 0.5) < 0.00001 );
        ASSERT( abs(p1.evaluate(0.75) + 0.265625) < 0.00001 );
        ASSERT( abs(p1.evaluate(1.0) - 1.0) < 0.00001 );

        NPolynomialBezier p2 = p1.toBezier();

        CHKSTR( p2, "1 B_{0,5} + 1 B_{1,5} + 1 B_{2,5} + 1 B_{3,5} - 2.2 B_{4,5} + 1 B_{5,5}" );

        ASSERT( abs(p2.evaluate(0.0) - 1.0) < 0.00001 );
        ASSERT( abs(p2.evaluate(0.25) - 0.953125) < 0.00001 );
        ASSERT( abs(p2.evaluate(0.5) - 0.5) < 0.00001 );
        ASSERT( abs(p2.evaluate(0.75) + 0.265625) < 0.00001 );
        ASSERT( abs(p2.evaluate(1.0) - 1.0) < 0.00001 );

        NPolynomialStandard p3 = p2.toStandard();

        ASSERT( p3.toString() == "16 X^5 - 16 X^4 + 1" );
    }

    if (0) { // test roots calculation using Cardano's formulas: template

        Numeric coeff[4] = { -3.2, 20, -40, 25 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        std::cout << "p = " << p1.toString() << "\n";

        std::vector<Numeric> roots = p1.findRoots();

        for (unsigned int i = 0; i < roots.size(); ++i)
            std::cout << "root[" << i << "] = " << roots[i] << "\n";
    }

    { // test roots calculation using Cardano's formulas: three real roots

        Numeric coeff[4] = { -3, 20, -40, 25 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        CHKSTR( p1, "25 X^3 - 40 X^2 + 20 X - 3" );

        std::vector<Numeric> roots = p1.findRoots();

        ASSERT(roots.size() == 3);
        ASSERT( abs(roots[0] - 0.276393) < 0.000001 );
        ASSERT( abs(roots[1] - 0.6) < 0.000001 );
        ASSERT( abs(roots[2] - 0.723607) < 0.000001 );
    }

    { // test roots calculation using Cardano's formulas: one single real root

        Numeric coeff[4] = { -0.5, 4, -8, 5 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        CHKSTR( p1, "5 X^3 - 8 X^2 + 4 X - 0.5" );

        std::vector<Numeric> roots = p1.findRoots();

        ASSERT(roots.size() == 1);
        ASSERT( abs(roots[0] - 0.186385) < 0.000001 );
    }

    { // test roots calculation using Cardano's formulas: case of a double root?

        Numeric coeff[4] = { -12, 80, -175, 125 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        CHKSTR( p1, "125 X^3 - 175 X^2 + 80 X - 12" );

        std::vector<Numeric> roots1 = p1.findRoots();

#if 0
        ASSERT( roots1.size() == 3 || roots1.size() == 2 );
        if (roots1.size() == 3) {
            ASSERT( abs(roots1[0] - 0.4) < 0.000001 );
            ASSERT( abs(roots1[1] - 0.4) < 0.000001 );
            ASSERT( abs(roots1[2] - 0.6) < 0.000001 );
        }
        else {
            ASSERT( abs(roots1[0] - 0.4) < 0.000001 );
            ASSERT( abs(roots1[1] - 0.6) < 0.000001 );
        }
#endif
    }

    { // test roots calculation using Cardano's formulas: case of a triple root

        Numeric coeff[4] = { -8, 60, -150, 125 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        CHKSTR( p1, "125 X^3 - 150 X^2 + 60 X - 8" );

        std::vector<Numeric> roots = p1.findRoots();

        ASSERT(roots.size() == 1);
        ASSERT( abs(roots[0] - 0.4) < 0.000001 );
    }

    { // test roots of quadratic Bezier and monomial polynomials

        Numeric coeff[3] = { 0.1, -1, 2 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[3]));

        CHKSTR( p1, "2 X^2 - 1 X + 0.1" );

        std::vector<Numeric> roots1 = p1.findRoots();

        ASSERT( roots1.size() == 2 );
        ASSERT( abs(roots1[0] - 0.138197) < 0.000001 );
        ASSERT( abs(roots1[1] - 0.361803) < 0.000001 );

        NPolynomialBezier p2 = p1.toBezier();

        CHKSTR( p2, "0.1 B_{0,2} - 0.4 B_{1,2} + 1.1 B_{2,2}" );

        std::vector<Numeric> roots2 = p2.findRoots();

        ASSERT( roots2.size() == 2 );
        ASSERT( abs(roots2[0] - 0.138197) < 0.000001 );
        ASSERT( abs(roots2[1] - 0.361803) < 0.000001 );
    }

    { // test roots of cubic monomial and Bezier polynomials

        Numeric coeff[4] = { -2, 36, -96, 64 };
        NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));

        CHKSTR( p1, "64 X^3 - 96 X^2 + 36 X - 2" );

        std::vector<Numeric> roots1 = p1.findRoots();

        ASSERT( roots1.size() == 3 );
        ASSERT( abs(roots1[0] - 0.066987) < 0.000001 );
        ASSERT( abs(roots1[1] - 0.5) < 0.000001 );
        ASSERT( abs(roots1[2] - 0.933012) < 0.000001 );

        NPolynomialBezier p2 = p1.toBezier();

        CHKSTR( p2, "- 2 B_{0,3} + 10 B_{1,3} - 10 B_{2,3} + 2 B_{3,3}" );

        std::vector<Numeric> roots2 = p2.findRoots();

        ASSERT( roots2.size() == 3 );
        ASSERT( abs(roots2[0] - 0.066987) < 0.000001 );
        ASSERT( abs(roots2[1] - 0.5) < 0.000001 );
        ASSERT( abs(roots2[2] - 0.933012) < 0.000001 );
    }

    { // interpolation coefficients: scalar product of Bernstein polynomials.

        Numeric testValues[6][3] = {
            { 0.125,      0.0357143, 0.00595238 },
            { 0.0892857,  0.0595238, 0.0178571 },
            { 0.0595238,  0.0714286, 0.0357143 },
            { 0.0357143,  0.0714286, 0.0595238 },
            { 0.0178571,  0.0595238, 0.0892857 },
            { 0.00595238, 0.0357143, 0.125 }
        };

        for (unsigned int i = 0; i < 6; ++i) {
            for (unsigned int j = 0; j < 3; ++j) {
                ASSERT( abs(DegreeInterpolationGenerator<Numeric>::scalarProduktBernsteinPolynomial(i, 5, j, 2)
                            - testValues[i][j]) < 0.000001 );
            }
        }

        // debug code to print out a matrix of scalar products:
        //std::cout << Matrix<Numeric>(6,3).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::scalarProduktBernsteinPolynomial(i, 5, j, 2); }).toString();
    }

    { // interpolation coefficients: check coefficients of dual basis.

        Numeric testValues[3][3] = {
            {  16, -24,         16 },
            { -24,  69.333333, -57.333333 },
            {  16, -57.333333,  69.333333 }
        };

        for (unsigned int i = 0; i < 3; ++i) {
            for (unsigned int j = 0; j < 3; ++j) {
                ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientDualBasis(i, j, 3)
                            - testValues[i][j]) < 0.000001 );
            }
        }

        // debug code to print out a matrix of dual basis coefficients:
        //std::cout << Matrix<Numeric>(4,4).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientDualBasis(i, j, 3); }).toString();
    }

    { // interpolation coefficients: check beta coefficients of degree
      // reducing matrix. The values used here are the same as in Example 1
      // and Example 2 of [1].

        Numeric testValues1[6][3] = {
            {  0.821429, -0.428571,  0.107143 },
            {  0.321429,  0.285714, -0.107143 },
            {  0,         0.642857, -0.142857 },
            { -0.142857,  0.642857,  0 },
            { -0.107143,  0.285714,  0.321429 },
            {  0.107143, -0.428571,  0.821429 }
        };

        for (unsigned int i = 0; i < 3; ++i) {
            for (unsigned int j = 0; j < 3; ++j) {
                ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 5, 2)
                            - testValues1[i][j]) < 0.000001 );
            }
        }

        Numeric testValues2[3][6] = {
            { 1, 0.6, 0.3, 0.1, 0,   0 },
            { 0, 0.4, 0.6, 0.6, 0.4, 0 },
            { 0, 0,   0.1, 0.3, 0.6, 1 }
        };

        for (unsigned int i = 0; i < 3; ++i) {
            for (unsigned int j = 0; j < 3; ++j) {
                ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 2, 5)
                            - testValues2[i][j]) < 0.000001 );
            }
        }

        // debug code to print out a matrix of beta coefficients:
        //std::cout << Matrix<Numeric>(6,3).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 5, 2); }).toString();
        //std::cout << Matrix<Numeric>(3,6).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 2, 5); }).toString();
    }

    { // test polynomial generator from roots

        Numeric coeff2[2] = { 1.0/3.0, 3 };
        NPolynomialStandard p2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));

        Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
        NPolynomialStandard p4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));

        Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
        NPolynomialStandard p8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));

        Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
        NPolynomialStandard p16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));

        ASSERT( p2.toString() == "- 1 X^2 + 3.33333 X - 1" );

        ASSERT( p4.toString() == "- 1 X^4 - 7.66667 X^3 - 2.33333 X^2 + 51.6667 X - 16.6667" );

        ASSERT( p8.toString() == "- 1 X^8 - 13.6667 X^7 - 37.3333 X^6 + 182 X^5 + 679 X^4 - 1295 X^3 - 3150 X^2 + 6166.67 X - 1666.67" );

        ASSERT( p16.toString() == "- 1 X^{16} - 39.6667 X^{15} - 651.667 X^{14} - 5448.33 X^{13} - 20440 X^{12} + 18475.3 X^{11} + 451673 X^{10} + 1.12172 \\!\\cdot\\! 10^{6} X^9 - 2.52262 \\!\\cdot\\! 10^{6} X^8 - 1.43506 \\!\\cdot\\! 10^{7} X^7 + 138542 X^6 + 7.92823 \\!\\cdot\\! 10^{7} X^5 + 3.97396 \\!\\cdot\\! 10^{7} X^4 - 2.40625 \\!\\cdot\\! 10^{8} X^3 - 8.33333 \\!\\cdot\\! 10^{7} X^2 + 3.64583 \\!\\cdot\\! 10^{8} X - 1.04167 \\!\\cdot\\! 10^{8}" );
    }
}

/** \page demos Demos and Speedtests

*******************************************************************************
*** Demos and Speedtests                                                    ***
*******************************************************************************

 To demonstrate and evaluate the root finding implementations the program
 contains 8 predefined demos, which contain different tests and runs of groups
 of polynomials.

 Demos 1-5 are simple examples used to illustrate the workings of the
 algorithms.

 Demos 6-8 are speedtests which will repeat the algorithms a large number of
 times and calculate their average speed for three different sets of
 polynomials.

******************************************************************************/

template <typename Clipper, typename Numeric>
double runSpeedtest(unsigned int iterations, Clipper& clipper, const PolynomialStandard<Numeric>& p, double left, double right)
{
#if !SPEEDTEST
    iterations = 1;
#endif

    double ts1 = timestamp();

    for (unsigned int i = 0; i < iterations; ++i)
    {
        clipper.findRoots(p, left, right);
    }

    double ts2 = timestamp();
    return ts2 - ts1;
}

template <typename Numeric>
void runAllClippers(const std::string& name, unsigned int iterations, const PolynomialStandard<Numeric>& p, int epsexp, double left = 0.0, double right = 1.0)
{
    Numeric epsilon(10);
    epsilon = pow(epsilon, -epsexp);

    std::string ltxlabel = name + "-" + Numeric::getName() + "-" + S(epsexp);

    spd(std::setprecision(12));

    BezierClip<Numeric> c1(epsilon, ltxlabel + "b:");
    ltx("\\section{Running \\texttt{BezClip} on " << name << " with epsilon " << epsexp << "}\n\n");
    double spd1 = runSpeedtest(iterations, c1, p, left, right);
    spd("algo=BezClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
        " iterations=" << iterations << " maxdepth=" << c1.maxdepth() <<
        " totaltime=" << spd1 << " speed=" << (spd1 / iterations) << "\n");
    ltx("\\clearpage\n");

    KClip<Numeric> c2(p.degree(), 2, epsilon, ltxlabel + "q:");
    ltx("\\section{Running \\texttt{QuadClip} on " << name << " with epsilon " << epsexp << "}\n\n");
    double spd2 = runSpeedtest(iterations, c2, p, left, right);
    spd("algo=QuadClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
        " iterations=" << iterations << " maxdepth=" << c2.maxdepth() <<
        " totaltime=" << spd2 << " speed=" << (spd2 / iterations) << "\n");
    ltx("\\clearpage\n");

    KClip<Numeric> c3(p.degree(), 3, epsilon, ltxlabel + "c:");
    ltx("\\section{Running \\texttt{CubeClip} on " << name << " with epsilon " << epsexp << "}\n\n");
    double spd3 = runSpeedtest(iterations, c3, p, left, right);
    spd("algo=CubeClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
        " iterations=" << iterations << " maxdepth=" << c3.maxdepth() <<
        " totaltime=" << spd3 << " speed=" << (spd3 / iterations) << "\n");
    ltx("\\clearpage\n");
}

template <typename Numeric>
void runDemo6()
{
    typedef PolynomialStandard<Numeric> NPolynomialStandard;

    spd("Single root functions from Barton and Juettler (Example 9):\n");
    spd("Numeric = " << Numeric::getName() << "\n");

    Numeric coeff2[2] = { 1.0/3.0, 3 };
    NPolynomialStandard f2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));

    Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
    NPolynomialStandard f4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));

    Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
    NPolynomialStandard f8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));

    Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
    NPolynomialStandard f16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$f_2$", 50000 * Numeric::IterationScale, f2, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$f_4$", 20000 * Numeric::IterationScale, f4, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$f_8$", 10000 * Numeric::IterationScale, f8, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$f_{16}$", 1000 * Numeric::IterationScale, f16, epsexp);
}

template <typename Numeric>
void runDemo7()
{
    typedef PolynomialStandard<Numeric> NPolynomialStandard;

    spd("Triple root functions from Liu et al. (Example 3):\n");
    spd("Numeric = " << Numeric::getName() << "\n");

    Numeric coeff4[4] = { 1.0/3.0, 1.0/3.0, 1.0/3.0, 5 };
    NPolynomialStandard g4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));

    Numeric coeff8[8] = { 1.0/3.0, 1.0/3.0, 1.0/3.0, -2, -2, -2, 5, 5 };
    NPolynomialStandard g8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));

    Numeric coeff16[16] = { 1.0/3.0, 1.0/3.0, 1.0/3.0, -2, -2, 5, 5, 5, 5, 5, 5, 5, -7, -7, -7, -7 };
    NPolynomialStandard g16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$g_4$", 20000 * Numeric::IterationScale, g4, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$g_8$", 10000 * Numeric::IterationScale, g8, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$g_{16}$", 1000 * Numeric::IterationScale, g16, epsexp);
}

template <typename Numeric>
void runDemo8()
{
    typedef PolynomialStandard<Numeric> NPolynomialStandard;

    spd("Near double root functions from Barton and Juettler (Example 11):\n");
    spd("Numeric = " << Numeric::getName() << "\n");

    print_precision = 16;

    Numeric coeff2[2] = { 0.56, 0.57 };
    NPolynomialStandard h2 (1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));

    Numeric coeff4[4] = { 0.4, 0.40000001, -1, 2 };
    NPolynomialStandard h4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));

    Numeric coeff8[8] = { 0.50000002, 0.50000003, -5, -5, -5, -7, -7, -7 };
    NPolynomialStandard h8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));

    Numeric coeff16[16] = { 0.30000008, 0.30000009, 6, 6, 6, 6, 6, 6, 6, -5, -5, -5, -5, -5, -5, -7 };
    NPolynomialStandard h16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$h_2$", 50000 * Numeric::IterationScale, h2, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$h_4$", 20000 * Numeric::IterationScale, h4, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$h_8$", 10000 * Numeric::IterationScale, h8, epsexp);

    for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
        runAllClippers("$h_{16}$", 1000 * Numeric::IterationScale, h16, epsexp);
}

void runDemo(unsigned int demo)
{
    ltx("{\\Large\\bf\\centerline{Demo " << demo << "}\\par}\n\\bigskip\n");

    if (demo == 1)
    {
        ltx("This demo exemplifies the techniques of the three clipping algorithms on a polynomial of 5th degree. The example was used in preparation of the associated talk, because most of the different cases occurring during the algorithms' execution can easily be found herein.\n\n");

        ltx("This demo works with the standard datatype \\texttt{double} and precision $\\varepsilon = 0.001$.\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        typedef Double<double>  Numeric;
        Numeric epsilon = 0.001;

        //Numeric coeff1[6] = { 1, -2, -1, 2, 0.5, 2 };
        Numeric coeff1[6] = { 1, -2, -1, 2.5, 0, 1 };
        PolynomialBezier<Numeric> p1 (std::vector<Numeric>(&coeff1[0], &coeff1[6]));

        debug_more = true;

        ltx("\\section{\\texttt{BezClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
        BezierClip<Numeric>(epsilon, "p1b:").findRoots(p1.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
        KClip<Numeric>(p1.degree(), 2, epsilon, "p1q:").findRoots(p1.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{CubeClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
        KClip<Numeric>(p1.degree(), 3, epsilon, "p1c:").findRoots(p1.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");
    }
    else if (demo == 2)
    {
        ltx("This demo visualises the techniques of the algorithms by applying them to a polynomial with many different roots. For this purpose a polynomial of 8th degree is used, that has six different roots on the unit interval. The demo works with the standard datatype \\texttt{long double} and precision $\\varepsilon = 0.001$.\n\n");

        ltx("Note that for this polynomial \\texttt{BezClip} needs 29 recursive calls with maximum depth six, \\texttt{QuadClip} needs 23 recursive calls with maximum depth five and \\texttt{CubeClip} 21 recursive calls with maximum depth five, each time to isolate all six roots.\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        typedef Double<long double>     Numeric;
        Numeric epsilon = 0.001;

        Numeric coeff2[9] = { 1, -5, 8, -8, 8, -10, 9, -4, 1 };
        PolynomialBezier<Numeric> p8 (std::vector<Numeric>(&coeff2[0], &coeff2[9]));

        ltx("\\section{\\texttt{BezClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
        BezierClip<Numeric>(epsilon, "p8b:").findRoots(p8.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
        KClip<Numeric>(p8.degree(), 2, epsilon, "p8q:").findRoots(p8.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{CubeClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
        KClip<Numeric>(p8.degree(), 3, epsilon, "p8c:").findRoots(p8.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");
    }
    else if (demo == 3)
    {
        typedef Double<double>  Numeric;
        Numeric epsilon = 0.001;

        ltx("This demonstration shows that the implementation can also successfully isolate the roots of the large Wilkinson polynomial with $n=20$. This polynomial is defined as $$p = \\prod_{i=1}^{20} (x - i)$$ and is analysed on the interval $[0,25]$. All three algorithms successfully find all roots with the standard datatype \\texttt{double} with precision $\\varepsilon = 0.001$. Regarding the results see \\autopageref{p20b:result}, \\autopageref{p20q:result} and \\autopageref{p20c:result}.\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        Numeric coeff20[20] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20 };
        PolynomialStandard<Numeric> p20 (1.0, std::vector<Numeric>(&coeff20[0], &coeff20[20]));

        ltx("\\section{\\texttt{BezClip} Applied to the Wilkinson Polynomial}\n\n");
        BezierClip<Numeric>(epsilon, "p20b:").findRoots(p20, 0.0, 25.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{QuadClip} Applied to the Wilkinson Polynomial}\n\n");
        KClip<Numeric>(p20.degree(), 2, epsilon, "p20q:").findRoots(p20, 0.0, 25.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{CubeClip} Applied to the Wilkinson Polynomial}\n\n");
        KClip<Numeric>(p20.degree(), 3, epsilon, "p20c:").findRoots(p20, 0.0, 25.0);
        ltx("\\clearpage\n");
    }
    else if (demo == 4)
    {
        typedef Double<long double>             Numeric;
        typedef PolynomialStandard<Numeric> NPolynomialStandard;

        ltx("This demo entry is used to test out further polynomials.\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        Numeric coeff2[2] = { 1.0/3.0, 3 };
        NPolynomialStandard p2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));

        Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
        NPolynomialStandard p4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));

        Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
        NPolynomialStandard p8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));

        Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
        NPolynomialStandard p16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));

        NPolynomialStandard p = p16;

        Numeric epsilon(10);
        epsilon = pow(epsilon, -32);

        ltx("\\section{\\texttt{BezClip} Applied to the Example Polynomial}\n\n");
        BezierClip<Numeric>(epsilon, "b:").findRoots(p, 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{QuadClip} Applied to the Example Polynomial}\n\n");
        KClip<Numeric>(p.degree(), 2, epsilon, "q:").findRoots(p, 0.0, 1.0);
        ltx("\\clearpage\n");

        ltx("\\section{\\texttt{CubeClip} Applied to the Example Polynomial}\n\n");
        KClip<Numeric>(p.degree(), 3, epsilon, "c:").findRoots(p, 0.0, 1.0);
        ltx("\\clearpage\n");
    }
    else if (demo == 5)
    {
        typedef Double<double>  Numeric;
        Numeric epsilon = 0.001;

        ltx("This demo contains the calculation procedures to generate the results shown in Appendix A of the lab report.\n\n");

        ltx("\\section{Degree Reduction and Raising Matrices for Degree $5$}\n\n");
        for (int i = 4; i >= 0; --i)
        {
            Matrix<Fraction> M1(6,i+1);
            M1.fill( DegreeInterpolationGenerator<Fraction>(5,i) );

            Matrix<Fraction> M2(i+1,6);
            M2.fill( DegreeInterpolationGenerator<Fraction>(i,5) );

            ltx("\\begin{align*}\n"
                "M_{5," << i << "} = " << M1.toString() << " &&\n"
                "M_{" << i << ",5} = " << M2.toString() << "\n"
                "\\end{align*}\n");
        }
        ltx("\\clearpage\n");

        Numeric coeff1[6] = { 1, -2, -1, 2.5, 0, 1 };
        PolynomialBezier<Numeric> p1 (std::vector<Numeric>(&coeff1[0], &coeff1[6]));

        debug_more = true;

        ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
        KClip<Numeric>(p1.degree(), 2, epsilon, "p1q:").findRoots(p1.toStandard(), 0.0, 1.0);
        ltx("\\clearpage\n");
    }
    else if (demo == 6)
    {
        ltx("Demo 6 is a speed test, which can be explicitly analysed using this execution protocol. While running the speed test the LaTeX output is deactivated and the root finding process iterated a large number of times to better measure the average execution time. The measurement results are saved as raw data in the file \\textt{demo6speed.txt}.\n\n");

        ltx("In this demo or speed test the following polynomials with a single root at $\\frac{1}{3}$ are used (Example 9 from Barto\\v{n} and J\\\"{u}ttler):\n"
            "\\begin{align*}\n"
            "f_2 &:= (t - \\tfrac{1}{3}) (3 - t) \\\\\n"
            "f_4 &:= (t - \\tfrac{1}{3}) (2 - t) (t + 5)^2 \\\\\n"
            "f_8 &:= (t - \\tfrac{1}{3}) (2 - t)^3 (t + 5)^4 \\\\\n"
            "f_{16} &:= (t - \\tfrac{1}{3}) (2 - t)^5 (t + 5)^{10}\n"
            "\\end{align*}\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        ltx("\\part{Numeric = double}\n");
        runDemo6< Double<double> >();

        ltx("\\part{Numeric = long double}\n");
        runDemo6< Double<long double> >();

        ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
        mpfr_set_default_prec(1024);
        runDemo6< MpfrFloat >();
    }
    else if (demo == 7)
    {
        ltx("Demo 7 is another speed test, which can be explicitly analysed using this execution protocol. While running the speed test the LaTeX output is deactivated and the root finding process iterated a large number of times to better measure the average execution time. The measurement results are saved as raw data in the file \\textt{demo7speed.txt}.\n\n");

        ltx("In this demo or speed test the following polynomials with a triple root at $\\frac{1}{3}$ as used (Example 3 from Liu et al.):\n"
            "\\begin{align*}\n"
            "g_4 &:= (t - \\tfrac{1}{3})^3 (t - 5) \\\\\n"
            "g_8 &:= (t - \\tfrac{1}{3})^3 (2 + t)^3 (t - 5)^2 \\\\\n"
            "g_{16} &:= (t - \\tfrac{1}{3})^3 (2 + t)^2 (t - 5)^7 (t + 7)^4\n"
            "\\end{align*}\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        //runDemo7< Double<double> >();
        //runDemo7< Double<long double> >();

        ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
        mpfr_set_default_prec(1024);
        runDemo7< MpfrFloat >();
    }
    else if (demo == 8)
    {
        ltx("Demo 8 is another speed test, which can be explicitly analysed using this execution protocol. While running the speed test the LaTeX output is deactivated and the root finding process iterated a large number of times to better measure the average execution time. The measurement results are saved as raw data in the file \\textt{demo8speed.txt}.\n\n");

        ltx("In this demo or speed test the following polynomials with two near roots are used (Example 3 from Liu et al.):\n"
            "\\begin{align*}\n"
            "h_2 &:= (t - 0.56)(t - 0.57) \\\\\n"
            "h_4 &:= (t - 0.4) (t - 0.40000001) (t+1) (2-t) \\\\\n"
            "h_8 &:= (t - 0.50000002) (t - 0.50000003) (t+5)^3 (t+7)^3 \\\\\n"
            "h_{16} &:= (t - 0.30000008) (t - 0.30000009) (6-t)^7 (t+5)^6 (t+7)\n"
            "\\end{align*}\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        //runDemo7< Double<double> >();
        //runDemo7< Double<long double> >();

        ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
        mpfr_set_default_prec(1024);
        runDemo8< MpfrFloat >();
    }
    else if (demo == 9)
    {
        typedef Double<long double>             Numeric;
        typedef PolynomialStandard<Numeric> NPolynomialStandard;

        ltx("This demo entry is used to test out further polynomials.\n\n");

        ltx("\\tableofcontents\n\n"
            "\\clearpage\n");

        std::vector<NPolynomialStandard> polys;

        Numeric coeff3[4] = { -1.0/16, +5.0/8.0, -3.0/2.0, 1.0 };
        polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff3[0], &coeff3[4])));

        Numeric coeff4[5] = { +5.0/256.0, -5.0/16.0, +21.0/16.0, -2.0, 1.0 };
        polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff4[0], &coeff4[5])));

        Numeric coeff5[6] = { -3.0/512.0, +35.0/256.0, -7.0/8.0, +9.0/4.0, -5.0/2.0, 1.0 };
        polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff5[0], &coeff5[6])));

        Numeric coeff6[7] = { +7.0/4096.0, -7.0/128.0, +63.0/128.0, -15.0/8.0, +55.0/16.0, -3.0, 1.0 };
        polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff6[0], &coeff6[7])));

        Numeric coeff7[8] = { -1.0/2048.0, +21.0/1024.0, -63.0/256.0, +165.0/128.0, -55.0/16.0, +39.0/8.0, -7.0/2.0, 1.0 };
        polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff7[0], &coeff7[8])));

        for (unsigned int i = 0; i < polys.size(); ++i)
        {
            runAllClippers("p"+S(i+3), 1000 * Numeric::IterationScale, polys[i], 6);
        }
    }
}

// ****************************************************************************
// *** main() with Command Line Parser                                      ***
// ****************************************************************************

void printLatexPreamble()
{
    ltx("\\documentclass[a4paper,12pt]{article}\n"

        "\\usepackage[latin1]{inputenc}\n"
        "\\usepackage[T1]{fontenc}\n"

        "\\usepackage{lmodern}\n"
        "\\usepackage[english]{babel}\n"

        "\\usepackage[tmargin=20mm,bmargin=20mm,lmargin=15mm,rmargin=15mm]{geometry}\n"

        "\\parskip=0pt\n"
        "\\parindent=0pt\n"
        "\\hbadness=10000\n"

        "\\usepackage{amsmath,amssymb,breqn}\n"

        "\\usepackage[final,colorlinks=true,linkcolor=blue!60!black]{hyperref}\n"

        "\\usepackage{pgfplots}\n"

        "\\pgfplotsset{"
        "  width=16cm,height=6cm,\n"
        "  axis x line=middle,axis y line=left,\n"
        "  domain=0:1,xmin=0,xmax=1,\n"
        "  samples=50,\n"
        "  legend style={at={(0.5,-0.1)}, anchor=north}, legend columns=-1}\n"

        "\\def\\arraystretch{1.2}\n"

        "% special \\iffinal construct to exclude large drawings\n"
        "\\newif\\iffinal\n"
        "\\finaltrue\n"

        "% enlarge space for numbering in toc\n"
        "\\makeatletter"
        "\\renewcommand\\@pnumwidth{1.85em}\n"
        "\\def\\sectionnumwidth#1{"
        "\\gdef\\numberline##1{\\hb@xt@ #1{##1\\hfil}\\hskip 1ex}"
        "}"
        "\\sectionnumwidth{6ex}\n"
        "\\makeatother\n"

        "\\begin{document}\n");
}

void printUsage(char* argv[])
{
    std::cerr <<
        "Usage: " << argv[0] << " [options] <coefficients>\n"
        "Options:\n"
        "  -B, --bezier         Coefficients on command line are for Bezier representation.\n"
        "  -D, --demo <n>       Run demo of different polynomials.\n"
        "  -M, --monomial       Coefficients on command line are for monomial representation.\n"
        "  -R, --roots <Scale>  Construct a polynomial with roots given as parameters and scale.\n"
        "  -a, --algo <Algo>    Select algorithm to run: Bezclip, Quadclip or Cubeclip.\n"
        "  -e, --epsilon <e>    Calculate until precision epsilon is reached.\n"
        "  -f, --fragment       Suppress printing the LaTeX preamble.\n"
        "  -h, --help           The help on command line parameters you are reading.\n"
        "  -i, --interval <a,b> Calculate roots in interval [a,b].\n"
        "  -l, --latex          Output LaTeX code of calculations (if compiled in).\n"
        "  -t, --testsuite      Run testsuite of checks.\n"
        ;
}

bool isCoefficient(const char* s)
{
    if (!*s) return false;

    if (*s == '-' || *s == '+') ++s;    // scan sign

    if (*s == '-') return false;        // -e and --epsilon parameters
    if (*s == 'e') return false;

    while (*s)
    {
        if (!isdigit(*s) &&
            *s != '-' && *s != '+' && *s != 'e' && *s != '.')
        {
            return false;
        }
        ++s;
    }
    return true;
}

int main(int argc, char* argv[])
{
    mpfr_set_default_prec (1024);

    char* endptr;
    bool badinput = false;

    // *** process command line parameters

    enum { ALGO_BEZCLIP, ALGO_QUADCLIP, ALGO_CUBECLIP } opt_algo = ALGO_QUADCLIP;

    enum { INPUT_BEZIER, INPUT_MONOMIAL, INPUT_ROOTS } opt_input = INPUT_BEZIER;

    bool opt_output_latexpreamble = true;

    unsigned int opt_rundemo = 0;

    typedef Double<long double> Numeric;

    Numeric opt_input_roots_scale = 1.0;

    Numeric opt_interval_left = 0.0;
    Numeric opt_interval_right = 1.0;

    Numeric opt_epsilon = 0.001;

    std::vector<Numeric> coefficients;

    while (1) {
        static struct option long_options[] = {
            { "algo",     required_argument,    0,  'a' },
            { "bezier",   no_argument,          0,  'B' },
            { "demo",     required_argument,    0,  'D' },
            { "epsilon",  required_argument,    0,  'e' },
            { "fragment", no_argument,          0,  'f' },
            { "help",     no_argument,          0,  'h' },
            { "interval", required_argument,    0,  'i' },
            { "latex",    no_argument,          0,  'l' },
            { "monomial", no_argument,          0,  'M' },
            { "roots",    required_argument,    0,  'R' },
            { "testsuite",no_argument,          0,  't' },
            { 0,          0,                    0,  0 }
        };

        // check parameter for coefficent syntax before getopt_long can misparse
        // negative numbers

        while ( optind < argc && isCoefficient(argv[optind]) )
        {
            Numeric coeff = strtod(argv[optind], &endptr);

            if (endptr && *endptr == 0)
                coefficients.push_back(coeff);
            else
            {
                std::cerr << "Invalid coefficient \"" << argv[optind] << "\"\n";
                badinput = true;
            }

            optind++;
        }

        int option_index = 0;
        int c = getopt_long(argc, argv, "h?lMBR:a:fe:i:Dt",
                            long_options, &option_index);

        if (c == -1) break;

        switch (c) {
        case 'l':
            output_latex = true;
            break;

        case 'f':
            opt_output_latexpreamble = false;
            break;

        case 'D':
            opt_rundemo = strtoul(optarg, &endptr, 10);
            if (!endptr || *endptr != 0) {
                std::cout << "Invalid demo number \"" << optarg << "\"\n";
                badinput = true;
            }
            break;

        case 'a':
            if (optarg[0] == 'b') {
                opt_algo = ALGO_BEZCLIP;
            }
            else if (optarg[0] == 'q') {
                opt_algo = ALGO_QUADCLIP;
            }
            else if (optarg[0] == 'c') {
                opt_algo = ALGO_CUBECLIP;
            }
            else {
                std::cerr << "Unknown algorithm \"" << optarg << "\" requested.\n";
            }
            break;

        case 'M':
            opt_input = INPUT_MONOMIAL;
            break;

        case 'B':
            opt_input = INPUT_BEZIER;
            break;

        case 'R':
            opt_input = INPUT_ROOTS;

            opt_input_roots_scale = strtod(optarg, &endptr);
            if (!endptr || *endptr != 0) {
                std::cout << "Invalid scale parameter \"" << optarg << "\"\n";
                badinput = true;
            }
            break;

        case 'i':
            opt_interval_left = strtod(optarg, &endptr);
            if (!endptr || *endptr != ',') {
                std::cout << "Invalid interval parameter \"" << optarg << "\"\n";
                badinput = true;
                break;
            }

            opt_interval_right = strtod(endptr+1, &endptr);
            if (!endptr || *endptr != 0) {
                std::cout << "Invalid interval parameter \"" << optarg << "\"\n";
                badinput = true;
            }
            break;

        case 'e':
            opt_epsilon = strtod(optarg, &endptr);
            if (!endptr || *endptr != 0) {
                std::cout << "Invalid epsilon parameter \"" << optarg << "\"\n";
                badinput = true;
            }
            break;

        case 't':
            (std::cout << "Running test suite: ").flush();

            testsuite< Double<double> >();
            testsuite< Double<long double> >();
            testsuite< MpfrFloat >();
            std::cout << "done.\n";

            return 0;

        case 'h':
        case '?':
            printUsage(argv);
            return 0;

        default:
            std::cerr << "?? getopt returned character code " << c << " ???\n";;
        }
    }

    if (badinput) return 0;

    if (coefficients.size() == 0 && !opt_rundemo)
    {
        printUsage(argv);
        return 0;
    }

    // *** main program directed by command line input

    if (opt_output_latexpreamble)
        printLatexPreamble();

    if (opt_rundemo) {
        runDemo(opt_rundemo);

        if (opt_output_latexpreamble)
            ltx("\\end{document}\n");
        return 0;
    }

    PolynomialStandard<Numeric> p;

    // create polynomial depending on input parameters
    if (opt_input == INPUT_BEZIER)
    {
        PolynomialBezier<Numeric> pB (coefficients);
        p = pB.toStandard();
    }
    else if (opt_input == INPUT_MONOMIAL)
    {
        p = PolynomialStandard<Numeric>(coefficients);
    }
    else if (opt_input == INPUT_ROOTS)
    {
        p = PolynomialStandard<Numeric>(opt_input_roots_scale, coefficients);
    }

    // run selected algorithm on polynomial
    if (opt_algo == ALGO_BEZCLIP)
    {
        BezierClip<Numeric>(opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
    }
    else if (opt_algo == ALGO_QUADCLIP)
    {
        KClip<Numeric>(p.degree(), 2, opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
    }
    else if (opt_algo == ALGO_CUBECLIP)
    {
        KClip<Numeric>(p.degree(), 3, opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
    }

    if (opt_output_latexpreamble)
        ltx("\\end{document}\n");

    return 0;
}

/*****************************************************************************/