clipper
clipper.cc
Go to the documentation of this file.
00001 // -*- mode: c++; fill-column: 78; coding: latin-1 -*-
00002 /** \mainpage
00003 
00004   This is an implementation of three root finding algorithms for polynomials
00005   using geometric clipping in the Bezier representation. It was written for a
00006   lab exercise course in numerical mathematics in winter semester 2011/12 at
00007   the FernUniversity of Hagen, Germany.
00008 
00009   For more information see the website:
00010   http://idlebox.net/2012/0320-Finding-Roots-of-Polynomials-by-Clipping/
00011 
00012   The code implements the three algorithms BezierClip, QuadClip [1] and
00013   CubeClip [2] (both in one class KCLip) described in the referenced
00014   papers. The program follows the unusual method to printing out debug
00015   information as a LaTeX document, which itself must be compiled by pdflatex
00016   for comfortable viewing. This technique allows the algorithm code to be
00017   interleaved with ltx(...) debug lines which are can contain plots and
00018   equations.
00019 
00020   There are six main parts in this source:
00021   - floating-point Numeric classes
00022   - polynomial classes for monomial and Bezier representation
00023   - implementation of the degree reduction matrices
00024   - the BezierClip and KClip algorithms
00025   - a test suite to check above algorithms
00026   - a collection of demos and speedtests
00027 
00028   Each part is separately documented below (or in Related Pages).
00029 
00030   References:
00031 
00032   [1] Barton, M. and Juettler, B. "Computing roots of polynomials by quadratic
00033   clipping", Computer Aided Geometric Design 24, (2007) p. 125-141.
00034 
00035   [2] Liu, L., Zhang, L., Lin, B., Wang, G. "Fast approach for computing roots
00036   of polynomials using cubic clipping", Computer Aided Geometric Design 26,
00037   (2009), p. 524-599.
00038 
00039 \verbatim
00040   ----------------------------------------------------------------------------
00041 
00042   Copyright (c) 2011-2012, Timo Bingmann
00043   All rights reserved. Permissions granted under the BSD 2-Clause License:
00044 
00045   Redistribution and use in source and binary forms, with or without
00046   modification, are permitted provided that the following conditions are met:
00047 
00048   * Redistributions of source code must retain the above copyright notice,
00049     this list of conditions and the following disclaimer.
00050 
00051   * Redistributions in binary form must reproduce the above copyright notice,
00052     this list of conditions and the following disclaimer in the documentation
00053     and/or other materials provided with the distribution.
00054 
00055   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00056   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00057   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00058   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00059   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00060   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
00061   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
00062   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
00063   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
00064   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00065   POSSIBILITY OF SUCH DAMAGE.
00066 
00067   ----------------------------------------------------------------------------
00068 \endverbatim
00069 */
00070 
00071 #include <stdlib.h>
00072 #include <assert.h>
00073 #include <getopt.h>
00074 #include <sys/time.h>
00075 
00076 #include <iostream>
00077 #include <iomanip>
00078 #include <vector>
00079 #include <string>
00080 #include <sstream>
00081 #include <stdexcept>
00082 #include <cmath>
00083 #include <limits>
00084 
00085 #include <mpfr.h>
00086 
00087 unsigned int output_latex = false;      // output latex or text debug
00088 
00089 #if !SPEEDTEST
00090 
00091 #define dbg(X)  do { if (!output_latex) std::cout << X; } while(0)
00092 #define ltx(X)  do { if (output_latex) std::cout << X; } while(0)
00093 #define spd(X)  dbg(X)
00094 
00095 #else
00096 
00097 // all output functions are disabled during speedtest except spd(...)
00098 #define dbg(X)  do { } while(0)
00099 #define ltx(X)  do { } while(0)
00100 #define spd(X)  do { std::cout << X; } while(0)
00101 
00102 #endif
00103 
00104 bool debug_more = false;                // be more verbose in some parts
00105 unsigned int print_precision = 6;       // precision of mpfr numbers in output
00106 
00107 // ****************************************************************************
00108 // *** Auxiliary Functions                                                  ***
00109 // ****************************************************************************
00110 
00111 /// format many different types as strings
00112 template <typename Type>
00113 std::string S(Type i)
00114 {
00115     std::ostringstream oss;
00116     oss << i;
00117     return oss.str();
00118 }
00119 
00120 /// format a number converting "e+10" to latex math strings
00121 std::string ltxNum(std::string os)
00122 {
00123     int i = os.size();
00124     while(--i > 0)
00125     {
00126         if (os[i] == 'e')
00127         {
00128             if (os[i+1] == '+')
00129                 os.erase(os.begin()+i, os.begin()+i+1);
00130 
00131             while (os[i+1] == '0')
00132                 os.erase(os.begin()+i, os.begin()+i+1);
00133 
00134             os.replace(i,1," \\!\\cdot\\! 10^{");
00135             os += "}";
00136             return os;
00137         }
00138     }
00139     return os;
00140 }
00141 
00142 /// format a number converting "e+10" to latex math strings
00143 template <typename Numeric>
00144 std::string ltxNum(Numeric n)
00145 {
00146     return ltxNum(n.toString());
00147 }
00148 
00149 /// format interval in latex headers
00150 template <typename Numeric>
00151 std::string ltxInterval(Numeric l, Numeric r)
00152 {
00153     std::ostringstream oss;
00154     oss << "\\texorpdfstring{$[" << l << "," << r << "]$}{" << l << "," << r << "}";
00155     return oss.str();
00156 }
00157 
00158 /// time measuring using gettimeofday().
00159 inline double timestamp()
00160 {
00161     struct timeval tv;
00162     gettimeofday(&tv, NULL);
00163     return tv.tv_sec + tv.tv_usec * 1e-6;
00164 }
00165 
00166 /** \page numeric Floating Point Numeric Classes
00167 
00168 *******************************************************************************
00169 *** Floating Point Numeric Classes                                          ***
00170 *******************************************************************************
00171 
00172  This part contains three classes implementing floating point numbers with
00173  different precisions and properties. These are later used for all
00174  calculations requiring a "Numeric" type and define many common math functions
00175  like sqrt() and also most C++ arithmetic operators.
00176 
00177  The first class Fraction keeps numerator and denominator of a rational number
00178  separately, and supports pretty-printed using LaTeX's \frac{} command. The
00179  purpose of this class is to output nicely formatted fractions and it should
00180  not be used for serious calculations.
00181 
00182  The second class Double takes a standard hardware floating-point type as
00183  template parameter: "float", "double" and "long double" are valid choices. It
00184  defines all later required operations using the standard libc implements of
00185  e.g. sqrt(), cbrt(), sin() and so on. Most of this is done using <cmath>, but
00186  some explicit template specializations were necessary to correctly switch
00187  between other libc functions.
00188 
00189  The final class MpfrFloat utilizes the GNU multi-precision float library with
00190  rounding to create a "Numeric" implementation with a fixed arbitrary
00191  precision. All necessary C++ arithmetic operators are overloaded in
00192  implemented using appropriate mpfr_functions.
00193 
00194 ******************************************************************************/
00195 
00196 /// Contains a pair of doubles signifying the numerator and denominator of a
00197 /// fraction. The fraction can be printed using LaTeX syntax and many common
00198 /// operators are defined.
00199 class Fraction
00200 {
00201 private:
00202     double      m_numerator;
00203     double      m_denominator;
00204 
00205     /// Calculate the greatest common divisor using Euclid's algorithm.
00206     static double gcd(double a, double b)
00207     {
00208         while (b != 0) {
00209             double t = b; b = fmod(a,b); a = t;
00210         }
00211         return a;
00212     }
00213 
00214     void reduce()
00215     {
00216         double g = gcd(m_numerator, m_denominator);
00217         m_numerator /= g;
00218         m_denominator /= g;
00219 
00220         if (m_denominator < 0 && m_numerator > 0) {
00221             m_denominator *= -1;
00222             m_numerator *= -1;
00223         }
00224     }
00225 
00226 public:
00227     Fraction(double numerator = 0, double denominator = 1)
00228         : m_numerator(numerator), m_denominator(denominator)
00229     {
00230         if (denominator == 0) throw(std::runtime_error("Divide by zero!"));
00231         reduce();
00232     }
00233 
00234     double numerator() const
00235     {
00236         return m_numerator;
00237     }
00238 
00239     double denominator() const
00240     {
00241         return m_denominator;
00242     }
00243 
00244     std::string toString() const
00245     {
00246         std::ostringstream oss;
00247         if (m_denominator == 1)
00248             oss << ltxNum(S(m_numerator));
00249         else
00250             oss << "\\frac{" << ltxNum(S(m_numerator)) << "}"
00251                 "{" << ltxNum(S(m_denominator)) << "}";
00252         return oss.str();
00253     }
00254 
00255     friend std::ostream& operator<< (std::ostream& os, const Fraction& f)
00256     {
00257         return os << f.toString();
00258     }
00259 
00260     double toDouble() const
00261     {
00262         return m_numerator / m_denominator;
00263     }
00264 
00265     Fraction operator+ (const Fraction& b) const
00266     {
00267         return Fraction( m_numerator * b.m_denominator + b.m_numerator * m_denominator, m_denominator * b.m_denominator );
00268     }
00269 
00270     Fraction operator- (const Fraction& b) const
00271     {
00272         return Fraction( m_numerator * b.m_denominator - b.m_numerator * m_denominator, m_denominator * b.m_denominator );
00273     }
00274 
00275     Fraction operator* (const Fraction& b) const
00276     {
00277         return Fraction( m_numerator * b.m_numerator, m_denominator * b.m_denominator );
00278     }
00279 
00280     Fraction operator/ (const Fraction& b) const
00281     {
00282         return Fraction( m_numerator * b.m_denominator, m_denominator * b.m_numerator );
00283     }
00284 
00285     Fraction& operator+= (const Fraction& b)
00286     {
00287         m_numerator = m_numerator * b.m_denominator + b.m_numerator * m_denominator;
00288         m_denominator = m_denominator * b.m_denominator;
00289         reduce();
00290         return *this;
00291     }
00292 
00293     /// simple method to calculate binomial coefficients for small (n,k)
00294     static Fraction binomial(unsigned int n, unsigned int k)
00295     {
00296         static const double factorial[21] = {
00297             1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800,
00298             39916800, 479001600, 6227020800, 87178291200, 1307674368000,
00299             20922789888000, 355687428096000, 6402373705728000,
00300             121645100408832000, 2432902008176640000
00301         };
00302 
00303         assert( n < sizeof(factorial) / sizeof(factorial[0]) );
00304         assert( k < sizeof(factorial) / sizeof(factorial[0]) );
00305 
00306         if (k > n) return 0;
00307 
00308         return Fraction(factorial[n], factorial[k] * factorial[n-k]);
00309     }
00310 };
00311 
00312 /// Implements a wrapper around one of the hardware floating-point data types
00313 /// so that these can be used as a "Numeric" template parameter.
00314 template <typename Type = double>
00315 class Double
00316 {
00317 private:
00318     Type        m_v;    // the value
00319 
00320 public:
00321     inline Double(const Type& v = 0)
00322         : m_v(v)
00323     {
00324     }
00325 
00326     // multiplied to the number of iterations in a speed test
00327     static const unsigned int IterationScale    = 20;
00328 
00329     static std::string getName();
00330 
00331     static inline Double MINUS_INFINITY()
00332     {
00333         return Double( - std::numeric_limits<Type>::infinity() );
00334     }
00335 
00336     static inline Double binomial(unsigned int n, unsigned int k);
00337 
00338     inline const Type& get() const
00339     {
00340         return m_v;
00341     }
00342 
00343     inline std::string toString() const
00344     {
00345         return S(m_v);
00346     }
00347 
00348     inline friend std::ostream& operator<< (std::ostream& os, const Double<Type>& d)
00349     {
00350         return os << d.toString();
00351     }
00352 
00353     inline Double operator- () const
00354     {
00355         return Double( - m_v );
00356     }
00357 
00358     inline Double operator+ (const Double& o) const
00359     {
00360         return Double( m_v + o.m_v );
00361     }
00362 
00363     inline Double operator- (const Double& o) const
00364     {
00365         return Double( m_v - o.m_v );
00366     }
00367 
00368     inline Double operator* (const Double& o) const
00369     {
00370         return Double( m_v * o.m_v );
00371     }
00372 
00373     inline Double operator/ (const Double& o) const
00374     {
00375         return Double( m_v / o.m_v );
00376     }
00377 
00378     inline Double& operator+= (const Double& o)
00379     {
00380         m_v += o.m_v;
00381         return *this;
00382     }
00383 
00384     inline Double& operator-= (const Double& o)
00385     {
00386         m_v -= o.m_v;
00387         return *this;
00388     }
00389 
00390     inline Double& operator*= (const Double& o)
00391     {
00392         m_v *= o.m_v;
00393         return *this;
00394     }
00395 
00396     inline Double& operator/= (const Double& o)
00397     {
00398         m_v /= o.m_v;
00399         return *this;
00400     }
00401 
00402     inline bool operator== (const Double &o) const
00403     {
00404         return (m_v == o.m_v);
00405     }
00406 
00407     inline bool operator!= (const Double &o) const
00408     {
00409         return (m_v != o.m_v);
00410     }
00411 
00412     inline bool operator< (const Double &o) const
00413     {
00414         return (m_v < o.m_v);
00415     }
00416 
00417     inline bool operator> (const Double &o) const
00418     {
00419         return (m_v > o.m_v);
00420     }
00421 
00422     inline bool operator>= (const Double &o) const
00423     {
00424         return (m_v >= o.m_v);
00425     }
00426 
00427     inline bool operator<= (const Double &o) const
00428     {
00429         return (m_v <= o.m_v);
00430     }
00431 };
00432 
00433 template <>
00434 std::string Double<float>::getName()
00435 {
00436     return "flt";
00437 }
00438 
00439 template <>
00440 std::string Double<double>::getName()
00441 {
00442     return "dbl";
00443 }
00444 
00445 template <>
00446 std::string Double<long double>::getName()
00447 {
00448     return "ldbl";
00449 }
00450 
00451 template <typename Type>
00452 inline Double<Type> operator- (const double& a, const Double<Type>& b)
00453 {
00454     return Double<Type>( a - b.get() );
00455 }
00456 
00457 template <typename Type>
00458 inline Double<Type> operator* (const double& a, const Double<Type>& b)
00459 {
00460     return Double<Type>( a * b.get() );
00461 }
00462 
00463 template <typename Type>
00464 inline Double<Type> abs(const Double<Type>& v)
00465 {
00466     return Double<Type>( std::abs(v.get()) );
00467 }
00468 
00469 template <typename Type>
00470 inline Double<Type> sqrt(const Double<Type>& v)
00471 {
00472     return Double<Type>( std::sqrt(v.get()) );
00473 }
00474 
00475 inline Double<float> cbrt(const Double<float>& v)
00476 {
00477     return Double<float>( cbrtf(v.get()) );
00478 }
00479 
00480 inline Double<double> cbrt(const Double<double>& v)
00481 {
00482     return Double<double>( cbrt(v.get()) );
00483 }
00484 
00485 inline Double<long double> cbrt(const Double<long double>& v)
00486 {
00487     return Double<long double>( cbrtl(v.get()) );
00488 }
00489 
00490 template <typename Type>
00491 inline Double<Type> cos(const Double<Type>& v)
00492 {
00493     return Double<Type>( std::cos(v.get()) );
00494 }
00495 
00496 template <typename Type>
00497 inline Double<Type> acos(const Double<Type>& v)
00498 {
00499     return Double<Type>( std::acos(v.get()) );
00500 }
00501 
00502 template <typename Type, typename Exponent>
00503 inline Double<Type> pow(const Double<Type>& v, Exponent e)
00504 {
00505     return Double<Type>( std::pow(v.get(), e) );
00506 }
00507 
00508 template <>
00509 inline Double<float> Double<float>::binomial(unsigned int n, unsigned int k)
00510 {
00511     return Double<float>( std::exp(lgammaf(n+1) - lgammaf(k+1) - lgammaf(n-k+1)) );
00512 }
00513 
00514 template <>
00515 inline Double<double> Double<double>::binomial(unsigned int n, unsigned int k)
00516 {
00517     return Double<double>( std::exp(lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1)) );
00518 }
00519 
00520 template <>
00521 inline Double<long double> Double<long double>::binomial(unsigned int n, unsigned int k)
00522 {
00523     return Double<long double>( std::exp(lgammal(n+1) - lgammal(k+1) - lgammal(n-k+1)) );
00524 }
00525 
00526 template <typename Type>
00527 inline int isnormal(const Double<Type>& v)
00528 {
00529     return std::isnormal(v.get());
00530 }
00531 
00532 /// Wraps a mpfr_t object of the multiple-precision floating-pointer number
00533 /// library mpfr into a Numeric object and provides all necessary operators
00534 /// and functions for further calculations.
00535 class MpfrFloat
00536 {
00537 private:
00538     mpfr_t      m_v;    // the value
00539 
00540 public:
00541     inline MpfrFloat()
00542     {
00543         mpfr_init(m_v);
00544     }
00545 
00546     inline MpfrFloat(const double& d)
00547     {
00548         mpfr_init_set_d(m_v, d, MPFR_RNDN);
00549     }
00550 
00551     /// copy constructor
00552     inline MpfrFloat(const MpfrFloat& d)
00553     {
00554         mpfr_init_set(m_v, d.m_v, MPFR_RNDN);
00555     }
00556 
00557     inline MpfrFloat& operator= (const MpfrFloat& d)
00558     {
00559         mpfr_set(m_v, d.m_v, MPFR_RNDN);
00560         return *this;
00561     }
00562 
00563     ~MpfrFloat()
00564     {
00565         mpfr_clear(m_v);
00566     }
00567 
00568     // multiplied to the number of iterations in a speed test
00569     static const unsigned int IterationScale    = 1;
00570 
00571     static std::string getName()
00572     {
00573         return "mpfr" + S(mpfr_get_default_prec());
00574     }
00575 
00576     static inline MpfrFloat MINUS_INFINITY()
00577     {
00578         MpfrFloat a; mpfr_set_inf(a.get(), -1); return a;
00579     }
00580 
00581     static inline MpfrFloat binomial(unsigned int n, unsigned int k)
00582     {
00583         MpfrFloat a, b;
00584 
00585         mpfr_fac_ui(a.get(), n, MPFR_RNDN);
00586 
00587         mpfr_fac_ui(b.get(), k, MPFR_RNDN);
00588         mpfr_div(a.get(), a.get(), b.get(), MPFR_RNDN);
00589 
00590         mpfr_fac_ui(b.get(), n - k, MPFR_RNDN);
00591         mpfr_div(a.get(), a.get(), b.get(), MPFR_RNDN);
00592 
00593         return a;
00594     }
00595 
00596     inline const mpfr_t& get() const
00597     {
00598         return m_v;
00599     }
00600 
00601     inline mpfr_t& get()
00602     {
00603         return m_v;
00604     }
00605 
00606     inline std::string toString() const
00607     {
00608         std::string str; str.resize(32);
00609         int r = mpfr_snprintf((char*)str.data(), str.size(), "%.*Rg", print_precision, m_v);
00610         str.resize(r);
00611         return str;
00612     }
00613 
00614     inline friend std::ostream& operator<< (std::ostream& os, const MpfrFloat& a)
00615     {
00616         return os << a.toString();
00617     }
00618 
00619     inline MpfrFloat operator- () const
00620     {
00621         MpfrFloat a; mpfr_neg(a.m_v, m_v, MPFR_RNDN); return a;
00622     }
00623 
00624     inline MpfrFloat operator+ (const MpfrFloat& o) const
00625     {
00626         MpfrFloat a; mpfr_add(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
00627     }
00628 
00629     inline MpfrFloat operator- (const MpfrFloat& o) const
00630     {
00631         MpfrFloat a; mpfr_sub(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
00632     }
00633 
00634     inline MpfrFloat operator* (const MpfrFloat& o) const
00635     {
00636         MpfrFloat a; mpfr_mul(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
00637     }
00638 
00639     inline MpfrFloat operator/ (const MpfrFloat& o) const
00640     {
00641         MpfrFloat a; mpfr_div(a.m_v, m_v, o.m_v, MPFR_RNDN); return a;
00642     }
00643 
00644     inline MpfrFloat& operator+= (const MpfrFloat& o)
00645     {
00646         mpfr_add(m_v, m_v, o.m_v, MPFR_RNDN);
00647         return *this;
00648     }
00649 
00650     inline MpfrFloat& operator-= (const MpfrFloat& o)
00651     {
00652         mpfr_sub(m_v, m_v, o.m_v, MPFR_RNDN);
00653         return *this;
00654     }
00655 
00656     inline MpfrFloat& operator*= (const MpfrFloat& o)
00657     {
00658         mpfr_mul(m_v, m_v, o.m_v, MPFR_RNDN);
00659         return *this;
00660     }
00661 
00662     inline MpfrFloat& operator/= (const MpfrFloat& o)
00663     {
00664         mpfr_div(m_v, m_v, o.m_v, MPFR_RNDN);
00665         return *this;
00666     }
00667 
00668     inline bool operator== (const MpfrFloat &o) const
00669     {
00670         return mpfr_cmp(m_v, o.m_v) == 0;
00671     }
00672 
00673     inline bool operator!= (const MpfrFloat &o) const
00674     {
00675         return mpfr_cmp(m_v, o.m_v) != 0;
00676     }
00677 
00678     inline bool operator< (const MpfrFloat &o) const
00679     {
00680         return mpfr_cmp(m_v, o.m_v) < 0;
00681     }
00682 
00683     inline bool operator> (const MpfrFloat &o) const
00684     {
00685         return mpfr_cmp(m_v, o.m_v) > 0;
00686     }
00687 
00688     inline bool operator<= (const MpfrFloat &o) const
00689     {
00690         return mpfr_cmp(m_v, o.m_v) <= 0;
00691     }
00692 
00693     inline bool operator>= (const MpfrFloat &o) const
00694     {
00695         return mpfr_cmp(m_v, o.m_v) >= 0;
00696     }
00697 };
00698 
00699 inline MpfrFloat operator- (const double& a, const MpfrFloat& b)
00700 {
00701     MpfrFloat x; mpfr_d_sub(x.get(), a, b.get(), MPFR_RNDN); return x;
00702 }
00703 
00704 inline MpfrFloat operator* (const double& a, const MpfrFloat& b)
00705 {
00706     MpfrFloat x; mpfr_mul_d(x.get(), b.get(), a, MPFR_RNDN); return x;
00707 }
00708 
00709 inline MpfrFloat abs(const MpfrFloat& v)
00710 {
00711     MpfrFloat x; mpfr_abs(x.get(), v.get(), MPFR_RNDN); return x;
00712 }
00713 
00714 inline MpfrFloat sqrt(const MpfrFloat& v)
00715 {
00716     MpfrFloat x; mpfr_sqrt(x.get(), v.get(), MPFR_RNDN); return x;
00717 }
00718 
00719 inline MpfrFloat cbrt(const MpfrFloat& v)
00720 {
00721     MpfrFloat x; mpfr_cbrt(x.get(), v.get(), MPFR_RNDN); return x;
00722 }
00723 
00724 inline MpfrFloat cos(const MpfrFloat& v)
00725 {
00726     MpfrFloat x; mpfr_cos(x.get(), v.get(), MPFR_RNDN); return x;
00727 }
00728 
00729 inline MpfrFloat acos(const MpfrFloat& v)
00730 {
00731     MpfrFloat x; mpfr_acos(x.get(), v.get(), MPFR_RNDN); return x;
00732 }
00733 
00734 inline MpfrFloat pow(const MpfrFloat& v, long int e)
00735 {
00736     MpfrFloat x; mpfr_pow_si(x.get(), v.get(), e, MPFR_RNDN); return x;
00737 }
00738 
00739 inline int isnormal(const MpfrFloat& v)
00740 {
00741     return mpfr_number_p(v.get());
00742 }
00743 
00744 /** \page polynomials Polynomial Classes in Monomial and Bezier Bases
00745 
00746 *******************************************************************************
00747 *** Polynomial Classes in Monomial and Bezier Bases                         ***
00748 *******************************************************************************
00749 
00750  This section of the implementation utilizes the basic arithmetic operations
00751  of a "Numeric" implementation to represent polynomials. Two representations
00752  are available: the standard monomial base and the Bezier
00753  representation. These are implemented in the classes PolynomialStandard and
00754  PolynomialBezier, respectively.
00755 
00756  Both classes implement many smaller algorithms later needed for implementing
00757  the three root finding algorithms. Among the base algorithms are:
00758 
00759  - conversion between the two polynomial representations
00760  - constructing a polynomial from given roots
00761  - finding the roots of quadratic and cubic polynomial using the direct
00762    formulas.
00763  - evaluation of the polynomial at a give position (using Horner or de
00764    Casteljau)
00765  - splitting a Bezier representation using de Castljau's algorithm
00766  - calculating the convex hull of the Bezier polygon using Jarvis's march
00767 
00768 ******************************************************************************/
00769 
00770 /// Base class for following polynomials in different bases. This could be
00771 /// used for some common virtual functions.
00772 template <typename Numeric>
00773 class Polynomial
00774 {
00775 };
00776 
00777 // *** Prototype declarations
00778 template <typename Numeric>
00779 class PolynomialBezier;
00780 
00781 template <typename Numeric>
00782 class Matrix;
00783 
00784 /// Represents a polynomial in standard monomial base.
00785 template <typename Numeric>
00786 class PolynomialStandard : public Polynomial<Numeric>
00787 {
00788 private:
00789 
00790     /// Coefficients of the polynomial. Degree is m_coefficients.size()-1
00791     std::vector<Numeric>        m_coefficient;
00792 
00793 public:
00794     PolynomialStandard()
00795     {
00796     }
00797 
00798     explicit PolynomialStandard(const std::vector<Numeric>& coefficient)
00799         : m_coefficient(coefficient)
00800     {
00801     }
00802 
00803     /// Construct polynomial from roots and a scaling multiplier
00804     PolynomialStandard(Numeric factor, const std::vector<Numeric>& roots)
00805         : m_coefficient(roots.size()+1, 0)
00806     {
00807         std::vector<bool> enumerate(roots.size(), 0);
00808 
00809         // iterate over all 2^n combinations to select a set of roots
00810         while(1)
00811         {
00812             // increment enumeration, calculate produkt and count factors
00813             Numeric a = 1.0;
00814             unsigned int afact = 0;
00815 
00816             unsigned int i = 0;
00817             for (; i < roots.size(); ++i)
00818             {
00819                 if (enumerate[i] == false) {
00820                     enumerate[i] = true;
00821                     a *= (- roots[i]);
00822                     ++afact;
00823                     break;
00824                 }
00825                 else {
00826                     enumerate[i] = false;
00827                 }
00828             }
00829 
00830             if (i == roots.size()) break; // all 2^n combinations done
00831 
00832             // multiply with remaining factors
00833             for (++i; i < roots.size(); ++i)
00834             {
00835                 if (enumerate[i] == true) {
00836                     a *= (- roots[i]);
00837                     ++afact;
00838                 }
00839             }
00840 
00841             // add to correct summand determinted by factor number
00842             m_coefficient[roots.size() - afact] += a;
00843         }
00844 
00845         m_coefficient[roots.size()] = 1.0;
00846 
00847         for (unsigned int j = 0; j <= roots.size(); ++j)
00848             m_coefficient[j] *= factor;
00849     }
00850 
00851     /// Return degree of the polynomial.
00852     unsigned int degree() const
00853     {
00854         return m_coefficient.size()-1;
00855     }
00856 
00857     /// Return a monomial coefficient.
00858     Numeric coefficient(unsigned int i) const
00859     {
00860         assert(i <= degree());
00861         return m_coefficient[i];
00862     }
00863 
00864     /// Return polynomial as latex math string.
00865     std::string toString(const std::string& X = "X") const
00866     {
00867         std::ostringstream oss;
00868         for (int i = m_coefficient.size()-1; i >= 0; --i)
00869         {
00870             if (m_coefficient[i] != 0) {
00871                 if (oss.str().size())
00872                     oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
00873                 else if (m_coefficient[i] < 0)
00874                     oss << "- ";
00875 
00876                 if (i >= 10)
00877                     oss << ltxNum(abs(m_coefficient[i])) << " " << X << "^{" << i << "}";
00878                 else if (i > 1)
00879                     oss << ltxNum(abs(m_coefficient[i])) << " " << X << "^" << i;
00880                 else if (i == 1)
00881                     oss << ltxNum(abs(m_coefficient[i])) << " " << X;
00882                 else
00883                     oss << ltxNum(abs(m_coefficient[i]));
00884             }
00885         }
00886         return oss.str();
00887     }
00888 
00889     /// Return polynomial as string for plotting with gnuplot.
00890     std::string toStringPlot() const
00891     {
00892         std::ostringstream oss;
00893         for (int i = m_coefficient.size()-1; i >= 0; --i)
00894         {
00895             if (m_coefficient[i] != 0) {
00896                 if (oss.str().size())
00897                     oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
00898                 else if (m_coefficient[i] < 0)
00899                     oss << "- ";
00900 
00901                 if (i > 1)
00902                     oss << abs(m_coefficient[i]) << " * x**" << i;
00903                 else if (i == 1)
00904                     oss << abs(m_coefficient[i]) << " * x";
00905                 else
00906                     oss << abs(m_coefficient[i]);
00907             }
00908         }
00909         return oss.str();
00910     }
00911 
00912     /// Calculate polynomial's value at the position X using Horner's schema.
00913     Numeric evaluate(Numeric X) const
00914     {
00915         Numeric val = 0;
00916         for (int i = m_coefficient.size()-1; i >= 0; --i)
00917         {
00918             val *= X;
00919             val += m_coefficient[i];
00920         }
00921         return val;
00922     }
00923 
00924     /// Calculate all real roots for a quadratic or cubic polynomial using the
00925     /// p-q- or Cardano's formulas.
00926     std::vector<Numeric> findRoots() const
00927     {
00928         const unsigned int N = m_coefficient.size()-1;
00929 
00930         std::vector<Numeric> roots;
00931 
00932         if (N == 0
00933             || (N == 1 && m_coefficient[1] == 0)
00934             || (N == 2 && m_coefficient[1] == 0 && m_coefficient[2] == 0)
00935             || (N == 3 && m_coefficient[1] == 0 && m_coefficient[2] == 0 && m_coefficient[3] == 0)
00936             )
00937         {
00938             return roots;
00939         }
00940         else if (N == 1
00941                  || (N == 2 && m_coefficient[2] == 0)
00942                  || (N == 3 && m_coefficient[2] == 0 && m_coefficient[3] == 0)
00943             )
00944         {
00945             roots.push_back( - m_coefficient[0] / m_coefficient[1] );
00946         }
00947         else if (N == 2
00948                  || (N == 3 && m_coefficient[3] == 0)
00949             )
00950         {
00951             Numeric D = (m_coefficient[1] * m_coefficient[1]) - (4.0 * m_coefficient[2] * m_coefficient[0]);
00952 
00953             if (D == 0)                 // one double root
00954             {
00955                 roots.push_back( - m_coefficient[1] / (2 * m_coefficient[0]) );
00956             }
00957             else if (D > 0)             // two real roots
00958             {
00959                 D = sqrt(D);
00960                 roots.push_back( (- m_coefficient[1] - D) / (2.0 * m_coefficient[2]) );
00961                 roots.push_back( (- m_coefficient[1] + D) / (2.0 * m_coefficient[2]) );
00962 
00963                 if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
00964             }
00965             // else (D < 0) only imaginary roots
00966             return roots;
00967         }
00968         else if (N == 3)
00969         {
00970             const std::vector<Numeric>& c = m_coefficient;
00971 
00972             // Cardano's formulas
00973             Numeric p = ( c[1] - (c[2] * c[2]) / (3 * c[3]) ) / c[3];
00974             Numeric q = ( ( 2.0 * c[2] * c[2] * c[2] / (27.0 * c[3])  - (c[2] * c[1] / 3.0) )
00975                          / c[3] + c[0] ) / c[3];
00976 
00977             //std::cout << "p = " << p << " - q = " << q << "\n";
00978 
00979             if (p == 0)                 // simple case: no quadratic term
00980             {
00981                 roots.push_back( cbrt(q) - c[2] / 3.0 / c[3] );
00982                 return roots;
00983             }
00984 
00985             Numeric D = p * p * p / 27.0 + q * q / 4.0;
00986 
00987             //std::cout << D << " - D\n";
00988 
00989             if (D > 0)                  // one real root and two complex ones
00990             {
00991                 D = sqrt(D);
00992 
00993                 Numeric u = cbrt( - q / 2.0 + D );
00994                 Numeric v = cbrt( - q / 2.0 - D );
00995 
00996                 roots.push_back( u + v - c[2] / 3.0 / c[3] );
00997             }
00998             else if (D == 0)            // three real roots: one single and one double
00999             {
01000                 Numeric u = cbrt( - q / 2.0 );
01001                 Numeric v = cbrt( - q / 2.0 );
01002 
01003                 roots.push_back( u + v - c[2] / 3.0 / c[3] );
01004                 roots.push_back( -(u + v) / 2.0 - c[2] / 3.0 / c[3] );
01005 
01006                 if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
01007             }
01008             else                        // three real roots
01009             {
01010                 Numeric phi = acos( - q / (2.0 * sqrt( - p * p * p / 27.0 )) );
01011                 Numeric w = 2 * sqrt( - p / 3.0 );
01012 
01013                 roots.push_back( w * cos( (phi + 2.0 * M_PI) / 3.0 ) );
01014                 roots.push_back( w * cos( (phi + 4.0 * M_PI) / 3.0 ) );
01015                 roots.push_back( w * cos( phi / 3.0 ) );
01016 
01017                 roots[0] = roots[0] - (c[2] / 3.0 / c[3]);
01018                 roots[1] = roots[1] - (c[2] / 3.0 / c[3]);
01019                 roots[2] = roots[2] - (c[2] / 3.0 / c[3]);
01020 
01021                 if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
01022             }
01023         }
01024         else {
01025             assert("Here roots are only calculable for quadratic and cubic polynomials.");
01026         }
01027 
01028         return roots;
01029     }
01030 
01031     /// Apply the linear transformation "X * (beta - alpha) + alpha" to the
01032     /// polynomial. This remaps the interval [alpha,beta] to [0,1].
01033     PolynomialStandard scale(Numeric alpha, Numeric beta) const
01034     {
01035         const unsigned int N = m_coefficient.size()-1;
01036 
01037         std::vector<Numeric> coeff(N+1);
01038 
01039         for (unsigned int i = 0; i <= N; ++i)
01040         {
01041             coeff[i] = 0;
01042             for (unsigned int j = i; j <= N; ++j)
01043             {
01044                 coeff[i] += m_coefficient[j] * pow(alpha, j-i) * Numeric::binomial(j,i);
01045             }
01046             coeff[i] *= pow((beta - alpha), i);
01047         }
01048 
01049         return PolynomialStandard(coeff);
01050     }
01051 
01052     /// Convert the polynomial from standard base to the Bernstein-Bezier
01053     /// basis in the unit interval. Use scale() to transform the coefficients
01054     /// before applying this base conversion.
01055     class PolynomialBezier<Numeric> toBezier() const;
01056 };
01057 
01058 /// Represents a polynomial in the Bernstein-Bezier basis using
01059 /// Bezier coefficients.
01060 template <typename Numeric>
01061 class PolynomialBezier : public Polynomial<Numeric>
01062 {
01063 private:
01064 
01065     /// The n+1 Bezier coefficients to the Bernstein polynomials basis
01066     /// (\f$B_{k,n}\f$)
01067     std::vector<Numeric>         m_coefficient;
01068 
01069     // Thus the polynomial's degree n is m_coefficient.size()-1.
01070 
01071 public:
01072     explicit PolynomialBezier(const std::vector<Numeric>& coefficient)
01073         : m_coefficient(coefficient)
01074     {
01075     }
01076 
01077     /// Return degree of the polynomial.
01078     unsigned int degree() const
01079     {
01080         return m_coefficient.size()-1;
01081     }
01082 
01083     /// Return a Bezier coefficient.
01084     Numeric coefficient(unsigned int i) const
01085     {
01086         assert(i <= degree());
01087         return m_coefficient[i];
01088     }
01089 
01090     /// Return polynomial as latex math string.
01091     std::string toString(const std::string& T = "") const
01092     {
01093         const unsigned int N = m_coefficient.size()-1;
01094 
01095         std::ostringstream oss;
01096         for (unsigned int i = 0; i <= N; ++i)
01097         {
01098             if (m_coefficient[i] != 0 || 1) {
01099                 if (oss.str().size())
01100                     oss << ((m_coefficient[i] >= 0) ? " + " : " - ");
01101                 else if (m_coefficient[i] < 0)
01102                     oss << "- ";
01103 
01104                 oss << ltxNum(abs(m_coefficient[i])) << " B_{" << i << "," << N << "}" << T;
01105             }
01106         }
01107 
01108         return oss.str();
01109     }
01110 
01111     /// Return Bezier coefficients as a polygon.
01112     std::string toPolygon() const
01113     {
01114         const unsigned int N = m_coefficient.size()-1;
01115 
01116         std::ostringstream oss;
01117         for (unsigned int i = 0; i <= N; ++i)
01118         {
01119             oss << "(" << Numeric( (double)i / N ) << "," << m_coefficient[i] << ") ";
01120         }
01121         return oss.str();
01122     }
01123 
01124     /// Return the convex hull of the Bezier coefficients calculated using
01125     /// Jarvis's march.
01126     std::vector< std::pair<Numeric,Numeric> > getConvexHull() const
01127     {
01128         std::vector< std::pair<Numeric,Numeric> > ch;
01129 
01130         const unsigned int N = m_coefficient.size()-1;
01131 
01132         unsigned int i = 0;     // current point in convex hull;
01133         ch.push_back( std::make_pair(0, m_coefficient[0]) );
01134 
01135         while( i < N ) // search right for point with highest angle
01136         {
01137             Numeric amax = Numeric::MINUS_INFINITY();
01138             unsigned int imax = i;
01139 
01140             for (unsigned int j = i+1; j <= N; ++j)
01141             {
01142                 Numeric a = (m_coefficient[j] - m_coefficient[i]) / (j - i);
01143 
01144                 if (amax < a) {
01145                     amax = a;
01146                     imax = j;
01147                 }
01148             }
01149 
01150             ch.push_back( std::make_pair(Numeric( (double)imax / N ), m_coefficient[imax]) );
01151             i = imax;
01152         }
01153 
01154         while ( i > 0 ) // search backwards for point with lowest angle
01155         {
01156             Numeric amax = Numeric::MINUS_INFINITY();
01157             unsigned int imax = i;
01158 
01159             for (unsigned int j = 0; j < i; ++j)
01160             {
01161                 Numeric a = (m_coefficient[i] - m_coefficient[j]) / (i - j);
01162 
01163                 if (amax < a) {
01164                     amax = a;
01165                     imax = j;
01166                 }
01167             }
01168 
01169             ch.push_back( std::make_pair(Numeric( (double)imax / N ), m_coefficient[imax]) );
01170             i = imax;
01171         }
01172 
01173         return ch;
01174     }
01175 
01176     /// Return the convex hull of the Bezier coefficients calculated using
01177     /// Jarvis's march.
01178     std::vector<Numeric> getConvexHullIntersection() const
01179     {
01180         const unsigned int N = m_coefficient.size()-1;
01181 
01182         unsigned int i = 0;             // current point in convex hull
01183         Numeric v = m_coefficient[0];   // current y position
01184 
01185         std::vector<Numeric> zerocross; // zero crossings
01186 
01187         while( i < N ) // search right for point with highest angle
01188         {
01189             Numeric amax = Numeric::MINUS_INFINITY();
01190             unsigned int imax = i;
01191 
01192             for (unsigned int j = i+1; j <= N; ++j)
01193             {
01194                 Numeric a = (m_coefficient[j] - m_coefficient[i]) / (j - i);
01195 
01196                 if (amax < a) {
01197                     amax = a;
01198                     imax = j;
01199                 }
01200             }
01201 
01202             if (m_coefficient[imax] * v < 0 || m_coefficient[imax] == 0) {
01203                 //std::cout << "amax " << amax << "\n";
01204                 //std::cout << "df " << v << "\n";
01205                 //std::cout << "df " << v / amax << "\n";
01206                 zerocross.push_back( Numeric( (double)i / N ) - v / (amax * N) );
01207             }
01208 
01209             v = m_coefficient[imax];
01210             i = imax;
01211         }
01212 
01213         while ( i > 0 ) // search backwards for point with lowest angle
01214         {
01215             Numeric amax = Numeric::MINUS_INFINITY();
01216             unsigned int imax = i;
01217 
01218             for (unsigned int j = 0; j < i; ++j)
01219             {
01220                 Numeric a = (m_coefficient[i] - m_coefficient[j]) / (i - j);
01221 
01222                 if (amax < a) {
01223                     amax = a;
01224                     imax = j;
01225                 }
01226             }
01227 
01228             if (m_coefficient[imax] * v < 0 || m_coefficient[imax] == 0) {
01229                 //std::cout << "amax " << amax << "\n";
01230                 //std::cout << "df " << v << "\n";
01231                 //std::cout << "df " << v / amax << "\n";
01232                 zerocross.push_back( Numeric( (double)i / N ) - v / (amax * N) );
01233             }
01234 
01235             v = m_coefficient[imax];
01236             i = imax;
01237         }
01238 
01239 #if 0
01240         for (unsigned int i = 0; i < zerocross.size(); ++i)
01241             std::cout << "zerocross[" << i << "] = " << zerocross[i] << "\n";
01242 #endif
01243 
01244         for (unsigned int i = 0; i < zerocross.size(); ++i)
01245             if (!isnormal(zerocross[i])) zerocross.clear();
01246 
01247         assert(zerocross.size() == 0 || zerocross.size() == 2);
01248 
01249         if (zerocross.size() == 2) {
01250             if (zerocross[0] > zerocross[1])
01251                 std::swap(zerocross[0], zerocross[1]);
01252         }
01253 
01254         return zerocross;
01255     }
01256 
01257     /// Return the convex hull of the Bezier coefficients for plotting.
01258     std::string toConvexHullString() const
01259     {
01260         std::vector< std::pair<Numeric,Numeric> > ch = getConvexHull();
01261 
01262         std::ostringstream os;
01263         for (unsigned int i = 0; i < ch.size(); ++i)
01264         {
01265             if (i != 0) os << " ";
01266             os << "(" << ch[i].first << "," << ch[i].second << ")";
01267         }
01268         return os.str();
01269     }
01270 
01271     /// Transform the polynomial from Bernstein-Bezier basis back to the
01272     /// standard monomial basis.
01273     class PolynomialStandard<Numeric> toStandard() const;
01274 
01275     /// Calculate the polynomial's value at the position X using the schema of
01276     /// de Casteljau.
01277     Numeric evaluate(Numeric X) const
01278     {
01279         const unsigned int N = m_coefficient.size()-1;
01280 
01281         std::vector<Numeric> coeff = m_coefficient;
01282 
01283         for (unsigned int i = 1; i <= N; ++i)
01284         {
01285             Numeric save = coeff[0];
01286 
01287             for (unsigned int j = 0; j <= N-i; ++j)
01288             {
01289                 coeff[j] = X * coeff[j+1] + (1.0 - X) * save;
01290                 save = coeff[j+1];
01291             }
01292         }
01293 
01294         return coeff[0];
01295     }
01296 
01297     /// Similar to evaluate: run de Casteljau's algorithm and output the
01298     /// Bezier coefficients of interval [0,t0] and [t0,1] as separate
01299     /// polynomial objects. This is later used to branch when finding roots.
01300     std::pair<PolynomialBezier,PolynomialBezier> deCasteljauSplit(Numeric t0) const
01301     {
01302         const unsigned int N = m_coefficient.size()-1;
01303 
01304         std::vector<Numeric> coeff = m_coefficient;
01305 
01306         std::vector<Numeric> coeff_left(N+1), coeff_right(N+1);
01307 
01308         coeff_left[0] = coeff[0];
01309         coeff_right[N] = coeff[N];
01310 
01311         for (unsigned int i = 1; i <= N; ++i)
01312         {
01313             Numeric save = coeff[0];
01314 
01315             for (unsigned int j = 0; j <= N-i; ++j)
01316             {
01317                 coeff[j] = t0 * coeff[j+1] + (1.0 - t0) * save;
01318                 save = coeff[j+1];
01319             }
01320 
01321             coeff_left[i] = coeff[0];
01322             coeff_right[N-i] = coeff[N-i];
01323         }
01324 
01325         return std::make_pair( PolynomialBezier(coeff_left), PolynomialBezier(coeff_right) );
01326     }
01327 
01328     /// Calculate all real roots for a quadratic or cubic polynomial in Bezier
01329     /// representation (according to Remark 3 in [1]). Imaginary roots are
01330     /// ignored.
01331     std::vector<Numeric> findRoots() const
01332     {
01333         const unsigned int N = m_coefficient.size()-1;
01334 
01335         std::vector<Numeric> roots;
01336 
01337         if (N == 2)
01338         {
01339             Numeric D = (m_coefficient[1] * m_coefficient[1]) - (m_coefficient[0] * m_coefficient[2]);
01340 
01341             // numerator of all roots
01342             Numeric N = m_coefficient[2] - 2 * m_coefficient[1] + m_coefficient[0];
01343             if (N == 0) return roots;
01344 
01345             if (D == 0)         // one double root
01346             {
01347                 roots.push_back( (m_coefficient[0] - m_coefficient[1]) / N );
01348             }
01349             else if (D > 0)             // two roots
01350             {
01351                 D = sqrt(D);
01352 
01353                 roots.push_back( (m_coefficient[0] - m_coefficient[1] - D) / N );
01354 
01355                 roots.push_back( (m_coefficient[0] - m_coefficient[1] + D) / N );
01356 
01357                 if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
01358             }
01359             // else (D < 0) only imaginary roots
01360             return roots;
01361         }
01362         else if (N == 3)
01363         {
01364             const std::vector<Numeric>& c = m_coefficient;
01365 
01366             /// Instead of using the Cardano's formulas in Bezier form as in
01367             /// [2], we return to the monomial base case for calculating the
01368             /// roots
01369             std::vector<Numeric> newCoeff(4);
01370 
01371             newCoeff[0] = c[0];
01372             newCoeff[1] = 3 * c[1] - 3 * c[0];
01373             newCoeff[2] = 3 * c[0] - 6 * c[1] + 3 * c[2];
01374             newCoeff[3] = c[3] - 3 * c[2] + 3 * c[1] - c[0];
01375 
01376             return PolynomialStandard<Numeric>(newCoeff).findRoots();
01377         }
01378         else {
01379             assert("Here roots are only calculable for quadratic and cubic polynomials.");
01380         }
01381 
01382         return roots;
01383     }
01384 
01385     /// This function is the heart of the root finding algorithm: reduce or
01386     /// raise the degree of this polynomial to k.
01387     PolynomialBezier    interpolateDegree(int k) const;
01388 
01389     /// This function is the heart of the root finding algorithm: reduce or
01390     /// raise the degree of this polynomial by applying a matrix.
01391     PolynomialBezier    interpolateDegree(const class Matrix<Numeric> &m) const;
01392 
01393     /// Compare two Bezier polynomials of same degree and return the maximum
01394     /// difference between pairs of coefficients.
01395     Numeric maxCoefficientDeltaTo(const PolynomialBezier& b) const
01396     {
01397         assert(m_coefficient.size() == b.m_coefficient.size());
01398         Numeric delta = abs(m_coefficient[0] - b.m_coefficient[0]);
01399         for (unsigned int i = 1; i < m_coefficient.size(); ++i)
01400             delta = std::max<Numeric>(delta, abs(m_coefficient[i] - b.m_coefficient[i]));
01401         return delta;
01402     }
01403 
01404     /// Create a new polynomial shifted upwards by adding v to all
01405     /// coefficients.
01406     PolynomialBezier operator+ (const Numeric& v) const
01407     {
01408         std::vector<Numeric> newCoeff(m_coefficient.begin(), m_coefficient.end());
01409         for (unsigned int i = 0; i < newCoeff.size(); ++i)
01410             newCoeff[i] += v;
01411         return PolynomialBezier(newCoeff);
01412     }
01413 
01414     /// Create a new polynomial shifted downwards by subtracting v from all
01415     /// coefficients.
01416     PolynomialBezier operator- (const Numeric& v) const
01417     {
01418         std::vector<Numeric> newCoeff(m_coefficient.begin(), m_coefficient.end());
01419         for (unsigned int i = 0; i < newCoeff.size(); ++i)
01420             newCoeff[i] -= v;
01421         return PolynomialBezier(newCoeff);
01422     }
01423 };
01424 
01425 /// Convert the polynomial from monomial representation to the
01426 /// Bernstein-Bezier basis in the unit interval. Use scale() to transform the
01427 /// coefficients before applying this base conversion.
01428 template <typename Numeric>
01429 PolynomialBezier<Numeric> PolynomialStandard<Numeric>::toBezier() const
01430 {
01431     const unsigned int N = m_coefficient.size()-1;
01432 
01433     std::vector<Numeric> coeff(N+1);
01434 
01435     for (unsigned int i = 0; i <= N; ++i)
01436     {
01437         coeff[i] = 0;
01438 
01439         for (unsigned int j = 0; j <= i; ++j)
01440         {
01441             coeff[i] += Numeric::binomial(i, j) / Numeric::binomial(N, i - j) * m_coefficient[i - j];
01442         }
01443     }
01444 
01445     return PolynomialBezier<Numeric>(coeff);
01446 }
01447 
01448 /// Transform the polynomial from Bernstein-Bezier basis back to the standard
01449 /// monomial basis.
01450 template <typename Numeric>
01451 PolynomialStandard<Numeric> PolynomialBezier<Numeric>::toStandard() const
01452 {
01453     const unsigned int N = m_coefficient.size()-1;
01454 
01455     std::vector<Numeric> coeff(N+1);
01456 
01457     for (unsigned int i = 0; i <= N; ++i)
01458     {
01459         coeff[i] = 0;
01460 
01461         for (unsigned int j = 0; j <= i; ++j)
01462         {
01463             coeff[i] += (j % 2 ? -1 : +1) * Numeric::binomial(i,j) * m_coefficient[i - j];
01464         }
01465 
01466         coeff[i] *= Numeric::binomial(N,i);
01467     }
01468 
01469     return PolynomialStandard<Numeric>(coeff);
01470 }
01471 
01472 /** \page matrixdegree Polynomial Degree Reduction and Raising
01473 
01474 *******************************************************************************
01475 *** Polynomial Degree Reduction and Raising                                 ***
01476 *******************************************************************************
01477 
01478  The most complex step in the KClip algorithms is reduction of the degree of
01479  the input polynomial. For this a Matrix can be precomputed using equations
01480  found in [1] and later only matrix multiplication is needed for reducing or
01481  raising a polynomial in Bezier representation.
01482 
01483  The actual formulas are implemented in "DegreeInterpolationGenerator" which
01484  can be used to fill() a Matrix object of the corresponding size. The
01485  Generator object is implemented as a functional which calculates a matrix
01486  entry. See the code of KClip and testsuite() for details.
01487 
01488 ******************************************************************************/
01489 
01490 /// Implements a simple matrix that can be filled using a functional object.
01491 template <typename Numeric>
01492 class Matrix
01493 {
01494 private:
01495 
01496     /// number of columns in the matrix
01497     unsigned int        m_cols;
01498 
01499     /// matrix data
01500     std::vector<Numeric> m_data;
01501 
01502 public:
01503 
01504     Matrix(unsigned int rows, unsigned int cols)
01505         : m_cols(cols)
01506     {
01507         m_data.resize(rows * cols);
01508     }
01509 
01510     unsigned int cols() const
01511     {
01512         return m_cols;
01513     }
01514 
01515     unsigned int rows() const
01516     {
01517         return m_data.size() / m_cols;
01518     }
01519 
01520     /// Return value at a position of the matrix.
01521     inline Numeric& at(unsigned int row, unsigned int col)
01522     {
01523         assert( row < rows() );
01524         assert( col < cols() );
01525         return m_data[row * m_cols + col];
01526     }
01527 
01528     /// Return value at a position of the matrix.
01529     inline const Numeric& at(unsigned int row, unsigned int col) const
01530     {
01531         assert( row < rows() );
01532         assert( col < cols() );
01533         return m_data[row * m_cols + col];
01534     }
01535 
01536     /// Fill in whole matrix from a function calculating values.
01537     template <typename Functional>
01538     Matrix& fill(Functional func)
01539     {
01540         for (unsigned int i = 0; i < rows(); ++i)
01541             for (unsigned int j = 0; j < cols(); ++j)
01542                 at(i,j) = func(i,j);
01543 
01544         return *this;
01545     }
01546 
01547     /// Return the matrix as a latex pmatrix math string.
01548     std::string toString() const
01549     {
01550         std::ostringstream oss;
01551         oss << "\\begin{pmatrix}\n";
01552         for (unsigned int i = 0; i < rows(); ++i)
01553         {
01554             for (unsigned int j = 0; j < cols(); ++j)
01555             {
01556                 if (j != 0) oss << " & ";
01557                 oss << ltxNum(at(i,j));
01558             }
01559             oss << " \\\\\n";
01560         }
01561         oss << "\\end{pmatrix}";
01562         return oss.str();
01563     }
01564 };
01565 
01566 /// Functional class to calculate the entries of the degree reducing or
01567 /// raising matrices.
01568 template <typename Numeric>
01569 class DegreeInterpolationGenerator
01570 {
01571 public:
01572 
01573     /// Calculate the scalar produkt of the two Bernstein polynomials
01574     /// B_{i,m} and B_{j,n} in the unit interval according to (13) in [1]
01575     static Numeric scalarProduktBernsteinPolynomial(int i, int m, int j, int n)
01576     {
01577         assert(i <= m); assert(j <= n);
01578         return Numeric( Numeric::binomial(m,i) * Numeric::binomial(n,j) ) / Numeric( Numeric(m + n + 1) * Numeric::binomial(m + n, i + j) );
01579     }
01580 
01581     /// Calculate the coefficients c_{p,q} of the dual basis \f$D_{p}^k\f$
01582     /// according to (10) in [1].
01583     static Numeric CoefficientDualBasis(int p, int q, int k)
01584     {
01585         Numeric a = Numeric( ((p+q) % 2 ? -1 : +1) ) / Numeric( Numeric::binomial(k,p) * Numeric::binomial(k,q) );
01586         Numeric b = 0;
01587 
01588         for (int j = 0; j <= std::min(p,q); ++j)
01589             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);
01590 
01591         return a * b;
01592     }
01593 
01594     /// Calculate the coefficients \f$\beta_{i,j}^{n,k}\f$ of the
01595     /// interpolation matrix according to (12) in [1].
01596     static Numeric CoefficientBeta(int i, int j, int n, int k)
01597     {
01598         Numeric a = 0;
01599         for (int l = 0; l <= k; ++l)
01600         {
01601             a += CoefficientDualBasis(j,l,k) * scalarProduktBernsteinPolynomial(i,n,l,k);
01602         }
01603         return a;
01604     }
01605 
01606 public:
01607     /// Dimensions of degree interpolation
01608     unsigned int        m_n, m_k;
01609 
01610     /// Constructor saving the degree parameters.
01611     DegreeInterpolationGenerator(unsigned int n, unsigned int k)
01612         : m_n(n), m_k(k)
01613     {
01614     }
01615 
01616     /// Functional operator for Matrix::fill().
01617     Numeric operator()(unsigned int i, unsigned int j) const
01618     {
01619         return CoefficientBeta(i,j,m_n,m_k);
01620     }
01621 };
01622 
01623 /// Apply a degree raising of reducing matrix to the given Bezier polynomial
01624 /// via simple matrix multiplication.
01625 template <typename Numeric>
01626 PolynomialBezier<Numeric> PolynomialBezier<Numeric>::interpolateDegree(const Matrix<Numeric> &m) const
01627 {
01628     const unsigned int N = m_coefficient.size()-1;
01629 
01630     assert(m.rows() == N+1);
01631     const unsigned int k = m.cols()-1;
01632 
01633     std::vector<Numeric> newCoeff (k+1);
01634 
01635     for (unsigned int i = 0; i <= k; ++i)
01636     {
01637         newCoeff[i] = 0;
01638 
01639         for (unsigned int j = 0; j <= N; ++j)
01640         {
01641             newCoeff[i] += m_coefficient[j] * m.at(j,i);
01642         }
01643     }
01644 
01645     return PolynomialBezier(newCoeff);
01646 }
01647 
01648 /// Apply a degree rasing of reducing matrix, where the matrix has to be
01649 /// constructed before application.
01650 template <typename Numeric>
01651 PolynomialBezier<Numeric> PolynomialBezier<Numeric>::interpolateDegree(int k) const
01652 {
01653     const int N = m_coefficient.size()-1;
01654 
01655     Matrix<Numeric> m(N+1, k+1);
01656     m.fill( DegreeInterpolationGenerator<Numeric>(N,k) );
01657 
01658     return interpolateDegree(m);
01659 }
01660 
01661 /** \page bezclip The Bezier Clipping Algorithm
01662 
01663 *******************************************************************************
01664 *** The Bezier Clipping Algorithm                                           ***
01665 *******************************************************************************
01666 
01667  Using the base algorithms implemented in PolynomialBezier the class
01668  BezierClip implements root finding by intersecting the convex hull of the
01669  Bezier polygon with the x-axis.
01670 
01671 ******************************************************************************/
01672 
01673 /// Implements the BezierClip algorithm
01674 template <typename Numeric>
01675 class BezierClip
01676 {
01677 private:
01678 
01679     /// roots found
01680     std::vector< std::pair<Numeric,Numeric> >   m_roots;
01681 
01682     /// epsilon used
01683     Numeric             m_epsilon;
01684 
01685     /// prefix for latex labels
01686     std::string         m_markprefix;
01687 
01688     /// maximum recursion depth
01689     unsigned int        m_maxdepth;
01690 
01691     typedef PolynomialStandard<Numeric> NPolynomialStandard;
01692     typedef PolynomialBezier<Numeric>   NPolynomialBezier;
01693 
01694 public:
01695 
01696     BezierClip(Numeric epsilon, const std::string& markprefix = "r")
01697         : m_epsilon(epsilon),
01698           m_markprefix(markprefix)
01699     {
01700     }
01701 
01702     unsigned int maxdepth() const
01703     { return m_maxdepth; }
01704 
01705     std::vector< std::pair<Numeric,Numeric> > findRoots(const NPolynomialStandard& p1, Numeric left, Numeric right)
01706     {
01707         m_roots.clear();
01708         m_maxdepth = 0;
01709 
01710         dbg("Running BezClip on p = " << p1.toString() << " within [" << left << "," << right << "]\n");
01711         ltx("\\begin{center}\\bf\\Large $" << p1.toString() << "$\\end{center}\n");
01712 
01713         // *** Input Polynomial
01714 
01715         ltx("\\textbf{Called \\texttt{BezClip} with input polynomial on interval " << ltxInterval(left,right) << ":}\n"
01716             "\\begin{dmath*}\n"
01717             "p = " << p1.toString() << "\n"
01718             "\\end{dmath*}\n\n");
01719 
01720         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
01721             "\\begin{axis}[\n"
01722             "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
01723             "       samples=100]\n");
01724 
01725         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");
01726 
01727         ltx("\\end{axis}\n"
01728             "\\end{tikzpicture}\n\\fi\\end{center}\n");
01729 
01730         ltx("\\subsection{Recursion Branch 1 for Input Interval " << ltxInterval(left,right) << "}\n\n");
01731 
01732         NPolynomialStandard p2 = p1.scale(left, right);
01733 
01734         findRootsRecursive(p2.toBezier(), left, right, "1", 1);
01735 
01736         // *** Result
01737 
01738         ltx("\\clearpage\n");
01739 
01740         ltx("\\subsection{Result: " << m_roots.size() << " Root Intervals}\n");
01741         ltx("\\label{" << m_markprefix << "result}\n");
01742 
01743         ltx("\\textbf{Input Polynomial on Interval $[" << left << "," << right << "]$}\n"
01744             "\\begin{dmath*}\n"
01745             "p = " << p1.toString() << "\n"
01746             "\\end{dmath*}\n\n");
01747 
01748         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
01749             "\\begin{axis}[\n"
01750             "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
01751             "       samples=100]\n");
01752 
01753         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");
01754 
01755         ltx("\\end{axis}\n"
01756             "\\end{tikzpicture}\n\\fi\\end{center}\n");
01757 
01758         ltx("\\textbf{Result: Root Intervals}\n");
01759 
01760         ltx("\\begin{center}\n");
01761         for (unsigned int i = 0; i < m_roots.size(); ++i)
01762         {
01763             if (i != 0) ltx(", ");
01764             ltx("$[" << m_roots[i].first << "," << m_roots[i].second << "]$");
01765         }
01766         ltx("\\end{center}\n");
01767 
01768         ltx("with precision $\\varepsilon = " << ltxNum(m_epsilon) << "$.\n\n");
01769 
01770         return m_roots;
01771     }
01772 
01773     // The polynomial p1 is here already normed to [0:1], left and right are
01774     // the _original_ interval and only needed for scaling the found roots.
01775     void findRootsRecursive(const NPolynomialBezier& p1, Numeric left, Numeric right, const std::string& mark, unsigned int depth)
01776     {
01777         if (right - left < m_epsilon)
01778         {
01779             if (p1.coefficient(0) == p1.coefficient(p1.degree()) ||
01780                 (p1.coefficient(0) > 0.0 && p1.coefficient(p1.degree()) < 0.0) ||
01781                 (p1.coefficient(0) < 0.0 && p1.coefficient(p1.degree()) > 0.0) ||
01782                 (abs(p1.coefficient(0)) < m_epsilon && abs(p1.coefficient(p1.degree())) < m_epsilon)
01783                 )
01784             {
01785                 dbg("Found root in interval [" << left << "," << right << "] at recursion depth " << depth << "!\n");
01786                 ltx("Found root in interval $[" << left << "," << right << "]$ at recursion depth " << depth << "!\n");
01787                 m_roots.push_back( std::make_pair(left,right) );
01788             }
01789             else {
01790                 dbg("Reached interval [" << left << "," << right << "] without sign change at depth " << depth << "!\n");
01791                 ltx("Reached interval $[" << left << "," << right << "]$ \\textbf{without sign change} at depth " << depth << "!\n\n");
01792             }
01793             return;
01794         }
01795 
01796         m_maxdepth = std::max(m_maxdepth, depth);
01797 
01798         // *** Normalization on [0:1] ***
01799 
01800         ltx("\\textbf{Normalized monomial und B\\'ezier representations and the B\\'ezier polygon:}\n"
01801             "\\begin{dmath*}\n"
01802             "p = " << p1.toStandard().toString() << "\n"
01803             "  = " << p1.toString("(X)") << "\n"
01804             "\\end{dmath*}\n\n");
01805 
01806         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
01807             "\\begin{axis}\n");
01808 
01809         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
01810 
01811         ltx("\\addplot[blue,dashed,mark=*] coordinates { " << p1.toPolygon() << " };\n");
01812 
01813         ltx("\\addplot[blue,fill=blue!50!white,opacity=0.2] coordinates { " << p1.toConvexHullString() << " };\n");
01814 
01815         std::vector<Numeric> chi = p1.getConvexHullIntersection();
01816         if (chi.size() == 2) {
01817             ltx("\\addplot[red,ultra thick] coordinates { (" << chi[0] << ",0) (" << chi[1] << ",0) };\n");
01818         }
01819 
01820         ltx("\\addplot[forget plot] coordinates { (0,0) };\n");
01821 
01822         ltx("\\end{axis}\n"
01823             "\\end{tikzpicture}\n\\fi\\end{center}\n");
01824 
01825         // *** Intersection of the Bezier hull ***
01826 
01827         ltx("\\textbf{Intersection of the convex hull with the $x$ axis:}\n");
01828 
01829         ltx("\\begin{dmath*}\n \\{ ");
01830         for (unsigned int i = 0; i < chi.size(); ++i) {
01831             if (i != 0) ltx(", "); ltx(chi[i]);
01832             if (i != 0) assert(chi[i-1] < chi[i]);
01833         }
01834         ltx("\\} \n\\end{dmath*}\n");
01835 
01836         std::vector< std::pair<Numeric,Numeric> > intervals;
01837 
01838         if (chi.size() == 2)
01839             intervals.push_back( std::make_pair(chi[0], chi[1]) );
01840 
01841         ltx("\\textbf{Intersection intervals with the $x$ axis:}\n\n");
01842 
01843         if (intervals.size() == 0) {
01844             ltx("No intersection with the $x$ axis. Done.\n\n");
01845             return;
01846         }
01847 
01848         ltx("\\begin{dmath*}\n");
01849         for (unsigned int i = 0; i < intervals.size(); ++i)
01850         {
01851             if (i != 0) ltx(" && ");
01852             ltx("[" << intervals[i].first << "," << intervals[i].second << "] ");
01853         }
01854         ltx("\\end{dmath*}\n");
01855 
01856         // find maximum interval size
01857 
01858         Numeric maxInterval = intervals[0].second - intervals[0].first;
01859         for (unsigned int i = 1; i < intervals.size(); ++i)
01860         {
01861             maxInterval = std::max(maxInterval, intervals[i].second - intervals[i].first);
01862         }
01863 
01864         ltx("Longest intersection interval: $" << ltxNum(maxInterval) << "$\n\n");
01865 
01866         if (maxInterval > 0.5)
01867         {
01868             Numeric middle = left + (right - left) / 2.0;
01869 
01870             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");
01871 
01872             std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(0.5);
01873 
01874             ltx("\\subsection{Recursion Branch " << mark << " 1 on the First Half " << ltxInterval(left,middle) << "}\n");
01875             ltx("\\label{" << m_markprefix << mark << " 1}\n");
01876 
01877             findRootsRecursive(polySplit.first, left, middle, mark + " 1", depth+1);
01878 
01879             ltx("\\subsection{Recursion Branch " << mark << " 2 on the Second Half " << ltxInterval(middle,right) << "}\n");
01880             ltx("\\label{" << m_markprefix << mark << " 2}\n");
01881 
01882             findRootsRecursive(polySplit.second, middle, right, mark + " 2", depth+1);
01883         }
01884         else
01885         {
01886             ltx(" $\\Longrightarrow$ Selective recursion: ");
01887             for (unsigned int i = 0; i < intervals.size(); ++i)
01888             {
01889 #if !SPEEDTEST
01890                 Numeric ileft = left + (right - left) * intervals[i].first;
01891                 Numeric iright = left + (right - left) * intervals[i].second;
01892 
01893                 ltx("\\hyperref[" << m_markprefix << mark << " " << i+1 << "]{interval " << i+1 << ": " << ltxInterval(ileft, iright) << "}, ");
01894 #endif
01895             }
01896             ltx("\n\n");
01897 
01898             for (unsigned int i = 0; i < intervals.size(); ++i)
01899             {
01900                 Numeric ileft = left + (right - left) * intervals[i].first;
01901                 Numeric iright = left + (right - left) * intervals[i].second;
01902 
01903                 if (intervals[i].first == 0)
01904                 {
01905                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].second);
01906 
01907                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
01908                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
01909 
01910                     findRootsRecursive(polySplit.first, ileft, iright, mark + " " + S(i+1), depth+1);
01911                 }
01912                 else if (intervals[i].second == 1.0)
01913                 {
01914                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].first);
01915 
01916                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
01917                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
01918 
01919                     findRootsRecursive(polySplit.second, ileft, iright, mark + " " + S(i+1), depth+1);
01920                 }
01921                 else
01922                 {
01923                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit1 = p1.deCasteljauSplit(intervals[i].first);
01924 
01925                     Numeric split2 = (intervals[i].second - intervals[i].first) / (1.0 - intervals[i].first);
01926 
01927                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit2 = polySplit1.second.deCasteljauSplit(split2);
01928 
01929                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
01930                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
01931 
01932                     findRootsRecursive(polySplit2.first, ileft, iright, mark + " " + S(i+1), depth+1);
01933                 }
01934             }
01935         }
01936     }
01937 };
01938 
01939 /** \page kclip The Quadclip and K-Clip Algorithms
01940 
01941 *******************************************************************************
01942 *** The Quadclip and K-Clip Algorithms                                      ***
01943 *******************************************************************************
01944 
01945  Using the base algorithms implemented in PolynomialBezier and the degree
01946  reduction and raising classes KClip implements the algorithms QuadClip and
01947  CubeClip (for parameters k=2 and k=3).
01948 
01949  See [1] and [2] for a description and analysis of the algorithms.
01950 
01951 ******************************************************************************/
01952 
01953 /// Implements the K-Clip algorithm
01954 template <typename Numeric>
01955 class KClip
01956 {
01957 private:
01958 
01959     /// degree of polynomial which this class processes
01960     unsigned int        m_degree;
01961 
01962     /// degree of reduced polynomial and clipping process
01963     unsigned int        m_clip;
01964 
01965     /// constant degree reducing matrix
01966     Matrix<Numeric>     m_degreeReducing;
01967 
01968     /// constant degree raising matrix
01969     Matrix<Numeric>     m_degreeRaising;
01970 
01971     /// roots found
01972     std::vector< std::pair<Numeric,Numeric> >   m_roots;
01973 
01974     /// epsilon used
01975     Numeric             m_epsilon;
01976 
01977     /// prefix for latex labels
01978     std::string         m_markprefix;
01979 
01980     /// maximum recursion depth
01981     unsigned int        m_maxdepth;
01982 
01983     typedef PolynomialStandard<Numeric> NPolynomialStandard;
01984     typedef PolynomialBezier<Numeric>   NPolynomialBezier;
01985 
01986 public:
01987 
01988     KClip(unsigned int degree, unsigned int clip = 2, Numeric epsilon = 0.001, const std::string& markprefix = "r")
01989         : m_degree(degree),
01990           m_clip(clip),
01991           m_degreeReducing(degree+1, clip+1),
01992           m_degreeRaising(clip+1, degree+1),
01993           m_epsilon(epsilon),
01994           m_markprefix(markprefix)
01995     {
01996         // construct constant matrices
01997         m_degreeReducing.fill( DegreeInterpolationGenerator<Numeric>(degree,clip) );
01998 
01999         m_degreeRaising.fill( DegreeInterpolationGenerator<Numeric>(clip,degree) );
02000     }
02001 
02002     unsigned int maxdepth() const
02003     { return m_maxdepth; }
02004 
02005     std::vector< std::pair<Numeric,Numeric> > findRoots(NPolynomialStandard p1, Numeric left, Numeric right)
02006     {
02007         m_roots.clear();
02008         m_maxdepth = 0;
02009 
02010         dbg("Running " << (m_clip == 2 ? "QuadClip" : "CubeClip") << " on p = " << p1.toString() << " in interval [" << left << "," << right << "]\n");
02011         ltx("\\begin{center}\\bf\\Large $" << p1.toString() << "$\\end{center}\n");
02012 
02013         // *** Input Polynom
02014 
02015         ltx("\\textbf{Called \\texttt{" << (m_clip == 2 ? "QuadClip" : "CubeClip") << "} with input polynomial on interval " << ltxInterval(left,right) << ":}\n"
02016             "\\begin{dmath*}\n"
02017             "p = " << p1.toString() << "\n"
02018             "\\end{dmath*}\n\n");
02019 
02020         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02021             "\\begin{axis}[\n"
02022             "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
02023             "       samples=100]\n");
02024 
02025         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");
02026 
02027         ltx("\\end{axis}\n"
02028             "\\end{tikzpicture}\n\\fi\\end{center}\n");
02029 
02030         ltx("\\subsection{Recursion Branch 1 for Input Interval " << ltxInterval(left,right) << "}\n\n");
02031 
02032         if (left != 0.0 || right != 0.0) {
02033             p1 = p1.scale(left, right);
02034         }
02035 
02036         findRootsRecursive(p1.toBezier(), left, right, "1", 1);
02037 
02038         // *** Result
02039 
02040         ltx("\\clearpage\n");
02041 
02042         ltx("\\subsection{Result: " << m_roots.size() << " Root Intervals}\n");
02043         ltx("\\label{" << m_markprefix << "result}\n");
02044 
02045         ltx("\\textbf{Input Polynomial on Interval $[" << left << "," << right << "]$}\n"
02046             "\\begin{dmath*}\n"
02047             "p = " << p1.toString() << "\n"
02048             "\\end{dmath*}\n\n");
02049 
02050         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02051             "\\begin{axis}[\n"
02052             "       domain=" << left << ":" << right << ",xmin=" << left << ",xmax=" << right << ",\n"
02053             "       samples=100]\n");
02054 
02055         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStringPlot() << " };\n");
02056 
02057         ltx("\\end{axis}\n"
02058             "\\end{tikzpicture}\n\\fi\\end{center}\n");
02059 
02060         ltx("\\textbf{Result: Root Intervals}\n");
02061 
02062         ltx("\\begin{center}\n");
02063         for (unsigned int i = 0; i < m_roots.size(); ++i)
02064         {
02065             if (i != 0) ltx(", ");
02066             ltx("$[" << m_roots[i].first << "," << m_roots[i].second << "]$");
02067         }
02068         ltx("\\end{center}\n");
02069 
02070         ltx("with precision $\\varepsilon = " << ltxNum(m_epsilon) << "$.\n\n");
02071 
02072         return m_roots;
02073     }
02074 
02075     // the polynomial p1 is here already normed to [0:1], left and right are the
02076     // _original_ interval and only needed for scaling the found roots.
02077     void findRootsRecursive(const NPolynomialBezier& p1, Numeric left, Numeric right, const std::string& mark, unsigned int depth)
02078     {
02079         if (right - left < m_epsilon)
02080         {
02081             if (p1.coefficient(0) == p1.coefficient(p1.degree()) ||
02082                 (p1.coefficient(0) > 0.0 && p1.coefficient(p1.degree()) < 0.0) ||
02083                 (p1.coefficient(0) < 0.0 && p1.coefficient(p1.degree()) > 0.0) ||
02084                 (abs(p1.coefficient(0)) < m_epsilon && abs(p1.coefficient(p1.degree())) < m_epsilon)
02085                 )
02086             {
02087                 dbg("Found root in interval [" << left << "," << right << "] at recursion depth " << depth << "!\n");
02088                 ltx("Found root in interval $[" << left << "," << right << "]$ at recursion depth " << depth << "!\n");
02089                 m_roots.push_back( std::make_pair(left,right) );
02090             }
02091             else {
02092                 dbg("Reached interval [" << left << "," << right << "] without sign change at depth " << depth << "!\n");
02093                 ltx("Reached interval $[" << left << "," << right << "]$ \\textbf{without sign change} at depth " << depth << "!\n\n");
02094 
02095                 ltx("p(0) = " << p1.coefficient(0) << " - p(1) " << p1.coefficient(p1.degree()) << "\n\n");
02096 
02097             }
02098             return;
02099         }
02100 
02101         m_maxdepth = std::max(m_maxdepth, depth);
02102 
02103         // *** Normalization on [0:1] ***
02104 
02105         ltx("\\textbf{Normalized monomial und B\\'ezier representations and the B\\'ezier polygon:}\n"
02106             "\\begin{dmath*}\n"
02107             "p = " << p1.toStandard().toString() << "\n"
02108             "  = " << p1.toString("(X)") << "\n"
02109             "\\end{dmath*}\n\n");
02110 
02111         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02112             "\\begin{axis}\n");
02113 
02114         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
02115 
02116         ltx("\\addplot[blue,dashed,mark=*] coordinates { " << p1.toPolygon() << " };\n");
02117 
02118         ltx("\\addplot[forget plot] coordinates { (0,0) };\n");
02119 
02120         ltx("\\end{axis}\n"
02121             "\\end{tikzpicture}\n\\fi\\end{center}\n");
02122 
02123 #if !SPEEDTEST
02124         if (debug_more)
02125         {
02126             // *** Reduced Polynomials ***
02127 
02128             ltx("\\textbf{Best approximation polynomials of degree 0, 1, 2, 3 and 4:}\n");
02129 
02130             ltx("\\begin{dgroup*}\n"
02131                 "\\begin{dmath*}\n q_0 = " << p1.interpolateDegree(0).toStandard().toString() << " = " << p1.interpolateDegree(0).toString() << " \\end{dmath*}\n"
02132                 "\\begin{dmath*}\n q_1 = " << p1.interpolateDegree(1).toStandard().toString() << " = " << p1.interpolateDegree(1).toString() << " \\end{dmath*}\n"
02133                 "\\begin{dmath*}\n q_2 = " << p1.interpolateDegree(2).toStandard().toString() << " = " << p1.interpolateDegree(2).toString() << " \\end{dmath*}\n"
02134                 "\\begin{dmath*}\n q_3 = " << p1.interpolateDegree(3).toStandard().toString() << " = " << p1.interpolateDegree(3).toString() << " \\end{dmath*}\n"
02135                 "\\begin{dmath*}\n q_4 = " << p1.interpolateDegree(4).toStandard().toString() << " = " << p1.interpolateDegree(4).toString() << " \\end{dmath*}\n"
02136                 "\\end{dgroup*}\n\n");
02137 
02138             ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02139                 "\\begin{axis}\n");
02140 
02141             ltx("\\addplot[purple] gnuplot { " << p1.interpolateDegree(0).toStandard().toStringPlot() << " };\n");
02142             ltx("\\addplot[cyan] gnuplot { " << p1.interpolateDegree(1).toStandard().toStringPlot() << " };\n");
02143             ltx("\\addplot[red] gnuplot { " << p1.interpolateDegree(2).toStandard().toStringPlot() << " };\n");
02144             ltx("\\addplot[green!60!black] gnuplot { " << p1.interpolateDegree(3).toStandard().toStringPlot() << " };\n");
02145             ltx("\\addplot[orange] gnuplot { " << p1.interpolateDegree(4).toStandard().toStringPlot() << " };\n");
02146 
02147             ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
02148 
02149             ltx("\\addplot[forget plot] coordinates { (0,0) };\n");
02150 
02151             ltx("\\legend{$q_0$,$q_1$,$q_2$,$q_3$,$q_4$,$p$};\n");
02152 
02153             ltx("\\end{axis}\n"
02154                 "\\end{tikzpicture}\n\\fi\\end{center}\n");
02155         }
02156 
02157         if (debug_more)
02158         {
02159             // *** Matrices ***
02160 
02161             ltx("\\textbf{Degree reduction and raising matrices:}\n");
02162 
02163             ltx("\\begin{align*}\n"
02164                 "M_{" << m_degree << "," << m_clip << "} = " << m_degreeReducing.toString() << " && "
02165                 "M_{" << m_clip << "," << m_degree << "} = " << m_degreeRaising.toString() << "\n"
02166                 "\\end{align*}\n");
02167         }
02168 #endif
02169 
02170         // *** Interpolation ***
02171 
02172         ltx("\\textbf{Degree reduction and raising:}\n");
02173 
02174         NPolynomialBezier p2 = p1.interpolateDegree(m_degreeReducing);
02175 
02176         NPolynomialBezier p3 = p2.interpolateDegree(m_degreeRaising);
02177 
02178         ltx("\\begin{dgroup*}\n"
02179             "\\begin{dmath*}\n q_" << m_clip << " = " << p2.toStandard().toString() << "\\\\\n"
02180             "       = " << p2.toString() << " \\end{dmath*}\n"
02181             "\\begin{dmath*}\n \\widetilde{q_" << m_clip << "} = " << p3.toStandard().toString() << "\\\\\n"
02182             "       = " << p3.toString() << " \\end{dmath*}\n"
02183             "\\end{dgroup*}\n\n");
02184 
02185         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02186             "\\begin{axis}\n");
02187 
02188         ltx("\\addplot[blue,no markers,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
02189         ltx("\\addplot[blue,dashed,mark=*,forget plot] coordinates { " << p1.toPolygon() << " };\n");
02190 
02191         ltx("\\addplot[orange,thick] gnuplot { " << p2.toStandard().toStringPlot() << " };\n");
02192         ltx("\\addplot[red,dashed,mark=*,forget plot] coordinates { " << p3.toPolygon() << " };\n");
02193         ltx("\\addplot[orange,dashed,mark=*] coordinates { " << p2.toPolygon() << " };\n");
02194 
02195         for (unsigned int i = 0; i <= m_degree; ++i)
02196         {
02197             ltx("\\addplot[green!60!black,forget plot,ultra thick] coordinates { "
02198                 "(" << i*(1.0/m_degree) << "," << p1.coefficient(i) << ") "
02199                 "(" << i*(1.0/m_degree) << "," << p3.coefficient(i) << ") };\n");
02200         }
02201 
02202         ltx("\\addplot[forget plot] coordinates { (0,0) };\n");
02203 
02204         ltx("\\legend{$p$,$q_" << m_clip << "$,$\\widetilde{q_" << m_clip << "}$};\n");
02205 
02206         ltx("\\end{axis}\n"
02207             "\\end{tikzpicture}\n\\fi\\end{center}\n");
02208 
02209         Numeric delta = p3.maxCoefficientDeltaTo(p1);
02210 
02211         ltx("The maximum difference of the B\\'ezier coefficients is $\\delta = " << ltxNum(delta) << "$.\n\n");
02212 
02213         // *** Bounding Polynomials M and m ***
02214 
02215         ltx("\\textbf{Bounding polynomials $M$ and $m$:}\n");
02216 
02217         NPolynomialBezier p2M = p2 + delta;
02218         NPolynomialBezier p2m = p2 - delta;
02219 
02220         ltx("\\begin{dgroup*}\n"
02221             "\\begin{dmath*}\n M = " << p2M.toStandard().toString() << " \\end{dmath*}\n"
02222             "\\begin{dmath*}\n m = " << p2m.toStandard().toString() << " \\end{dmath*}\n"
02223             "\\end{dgroup*}\n");
02224 
02225         // *** Root ***
02226 
02227         ltx("\\textbf{Root of $M$ and $m$:}\n");
02228 
02229         std::vector<Numeric> rootsM = p2M.findRoots();
02230         std::vector<Numeric> rootsm = p2m.findRoots();
02231 
02232         ltx("\\begin{align*}\n N(M) = \\{ ");
02233         for (unsigned int i = 0; i < rootsM.size(); ++i) {
02234             if (i != 0) ltx(", ");
02235             ltx(ltxNum(rootsM[i]));
02236             if (i != 0) assert(rootsM[i-1] <= rootsM[i]);
02237         }
02238         ltx("\\} &&");
02239 
02240         ltx("N(m) = \\{ ");
02241         for (unsigned int i = 0; i < rootsm.size(); ++i) {
02242             if (i != 0) ltx(", ");
02243             ltx(ltxNum(rootsm[i]));
02244             if (i != 0) assert(rootsm[i-1] <= rootsm[i]);
02245         }
02246         ltx("\\}\n\\end{align*}\n");
02247 
02248         std::vector< std::pair<Numeric,Numeric> > intervals;
02249 
02250         if (m_clip == 2)
02251         {
02252             // *** Quadratic Clipping
02253 
02254             if (rootsM.size() == 0 && rootsm.size() == 0)
02255             {
02256                 // no intervals
02257             }
02258             else if (rootsM.size() == 1 && rootsm.size() == 0)
02259             {
02260                 // no intervals -- but a double root?
02261                 assert(0);
02262             }
02263             else if (rootsM.size() == 2 && rootsm.size() == 0)
02264             {
02265                 intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
02266                                                     std::min<Numeric>(1,rootsM[1])) );
02267             }
02268             else if (rootsM.size() == 0 && rootsm.size() == 1)
02269             {
02270                 // no intervals -- but a double root?
02271                 assert(0);
02272             }
02273             else if (rootsM.size() == 0 && rootsm.size() == 2)
02274             {
02275                 intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
02276                                                     std::min<Numeric>(1,rootsm[1])) );
02277             }
02278             else if (rootsM.size() == 1 && rootsm.size() == 1)
02279             {
02280                 // impossible.
02281                 assert(0);
02282             }
02283             else if (rootsM.size() == 2 && rootsm.size() == 1)
02284             {
02285                 intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
02286                                                     std::min<Numeric>(1,rootsM[1])) );
02287             }
02288             else if (rootsM.size() == 1 && rootsm.size() == 2)
02289             {
02290                 intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
02291                                                     std::min<Numeric>(1,rootsm[1])) );
02292             }
02293             else if (rootsM.size() == 2 && rootsm.size() == 2)
02294             {
02295                 if (rootsM[0] > rootsm[0])
02296                 {
02297                     if (rootsM[0] > 0 && rootsm[0] < 1) {
02298                         intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[0]),
02299                                                             std::min<Numeric>(1,rootsM[0])) );
02300                     }
02301                     if (rootsM[1] < 1 && rootsm[1] > 0) {
02302                         intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[1]),
02303                                                             std::min<Numeric>(1,rootsm[1])) );
02304                     }
02305                 }
02306                 else
02307                 {
02308                     if (rootsm[0] > 0 && rootsM[0] < 1) {
02309                         intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsM[0]),
02310                                                             std::min<Numeric>(1,rootsm[0])) );
02311                     }
02312                     if (rootsm[1] < 1 && rootsM[1] > 0) {
02313                         intervals.push_back( std::make_pair(std::max<Numeric>(0,rootsm[1]),
02314                                                             std::min<Numeric>(1,rootsM[1])) );
02315                     }
02316                 }
02317             }
02318             else {
02319                 assert("Invalid combination of roots");
02320             }
02321         }
02322         else if (m_clip == 3)
02323         {
02324             // *** Cubic Clipping
02325 
02326             if (rootsM.size() == 1 && rootsm.size() == 1)
02327             {
02328                 if (rootsM[0] < rootsm[0])
02329                 {
02330                     intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
02331                 }
02332                 else // rootsM[0] > rootsm[0]
02333                 {
02334                     intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
02335                 }
02336             }
02337             else if (rootsM.size() == 1 && rootsm.size() == 2)
02338             {
02339                 assert(0); // unknown
02340             }
02341             else if (rootsM.size() == 2 && rootsm.size() == 1)
02342             {
02343                 assert(0); // unknown
02344             }
02345             else if (rootsM.size() == 2 && rootsm.size() == 2)
02346             {
02347                 if (rootsM[0] < rootsm[0])
02348                 {
02349                     intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
02350                     intervals.push_back( std::make_pair(rootsM[0], rootsm[1]) );
02351                     intervals.push_back( std::make_pair(rootsm[1], rootsM[1]) );
02352                 }
02353                 else // rootsM[0] > rootsm[0]
02354                 {
02355                     intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
02356                     intervals.push_back( std::make_pair(rootsm[0], rootsM[1]) );
02357                     intervals.push_back( std::make_pair(rootsM[1], rootsm[1]) );
02358                 }
02359             }
02360             else if (rootsM.size() == 1 && rootsm.size() == 3)
02361             {
02362                 if (rootsM[0] < rootsm[0])
02363                 {
02364                     intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
02365                     intervals.push_back( std::make_pair(rootsm[1], rootsm[2]) );
02366                 }
02367                 else // rootsM[0] > rootsm[0]
02368                 {
02369                     intervals.push_back( std::make_pair(rootsm[0], rootsm[1]) );
02370                     intervals.push_back( std::make_pair(rootsm[2], rootsM[0]) );
02371                 }
02372             }
02373             else if (rootsM.size() == 3 && rootsm.size() == 1)
02374             {
02375                 if (rootsM[0] < rootsm[0])
02376                 {
02377                     intervals.push_back( std::make_pair(rootsM[0], rootsM[1]) );
02378                     intervals.push_back( std::make_pair(rootsM[2], rootsm[0]) );
02379                 }
02380                 else // rootsM[0] > rootsm[0]
02381                 {
02382                     intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
02383                     intervals.push_back( std::make_pair(rootsM[1], rootsM[2]) );
02384                 }
02385             }
02386             else if (rootsM.size() == 2 && rootsm.size() == 3)
02387             {
02388                 assert(0); // unknown
02389             }
02390             else if (rootsM.size() == 3 && rootsm.size() == 2)
02391             {
02392                 assert(0); // unknown
02393             }
02394             else if (rootsM.size() == 3 && rootsm.size() == 3)
02395             {
02396                 if (rootsM[0] < rootsm[0])
02397                 {
02398                     intervals.push_back( std::make_pair(rootsM[0], rootsm[0]) );
02399                     intervals.push_back( std::make_pair(rootsm[1], rootsM[1]) );
02400                     intervals.push_back( std::make_pair(rootsM[2], rootsm[2]) );
02401                 }
02402                 else // rootsM[0] > rootsm[0]
02403                 {
02404                     intervals.push_back( std::make_pair(rootsm[0], rootsM[0]) );
02405                     intervals.push_back( std::make_pair(rootsM[1], rootsm[1]) );
02406                     intervals.push_back( std::make_pair(rootsm[2], rootsM[2]) );
02407                 }
02408             }
02409             else {
02410                 assert(0); // invalid combination
02411             }
02412 
02413             std::vector< std::pair<Numeric,Numeric> > newintervals;
02414 
02415             for (unsigned int i = 0; i < intervals.size(); ++i)
02416             {
02417                 if (intervals[i].first > intervals[i].second)
02418                 {
02419                     std::swap(intervals[i].first, intervals[i].second);
02420                 }
02421                 assert( intervals[i].first <= intervals[i].second);
02422 
02423                 if (intervals[i].second < 0.0) continue;
02424                 if (intervals[i].first > 1.0) continue;
02425 
02426                 newintervals.push_back( std::make_pair(std::max<Numeric>(0.0, intervals[i].first),
02427                                                        std::min<Numeric>(1.0, intervals[i].second)) );
02428             }
02429 
02430             std::swap(intervals, newintervals);
02431         }
02432         else {
02433             assert("Invalid combination of roots");
02434         }
02435 
02436         ltx("\\textbf{Intersection intervals:}\n");
02437 
02438         ltx("\\begin{center}\\iffinal\\begin{tikzpicture}\n"
02439             "\\begin{axis}\n");
02440 
02441         ltx("\\addplot[blue,very thick] gnuplot { " << p1.toStandard().toStringPlot() << " };\n");
02442 
02443         ltx("\\addplot[green!60!black] gnuplot { " << p2M.toStandard().toStringPlot() << " };\n");
02444         ltx("\\addplot[green!60!black] gnuplot { " << p2m.toStandard().toStringPlot() << " };\n");
02445 
02446         for (unsigned int i = 0; i < intervals.size(); ++i)
02447             ltx("\\addplot[red,ultra thick,forget plot] coordinates { (" << intervals[i].first << ",0) (" << intervals[i].second << ",0) };\n");
02448 
02449         ltx("\\addplot[forget plot] coordinates { (0,0) };\n");
02450 
02451         //ltx("\\legend{$p$,$M$,$m$};\n");
02452 
02453         ltx("\\end{axis}\n"
02454             "\\end{tikzpicture}\n\\fi\\end{center}\n");
02455 
02456         if (intervals.size() == 0) {
02457             ltx("No intersection intervals with the $x$ axis.\n\n");
02458             return;
02459         }
02460 
02461         ltx("\\begin{center}\n");
02462         for (unsigned int i = 0; i < intervals.size(); ++i)
02463         {
02464             if (i != 0) ltx(", ");
02465             ltx("$[" << intervals[i].first << "," << intervals[i].second << "]$");
02466         }
02467         ltx("\\end{center}\n");
02468 
02469         // find maximum interval size
02470 
02471         Numeric maxInterval = intervals[0].second - intervals[0].first;
02472         for (unsigned int i = 1; i < intervals.size(); ++i)
02473         {
02474             maxInterval = std::max(maxInterval, intervals[i].second - intervals[i].first);
02475         }
02476 
02477         ltx("Longest intersection interval: $" << ltxNum(maxInterval) << "$\n\n");
02478 
02479         if (maxInterval > 0.5)
02480         {
02481             Numeric middle = left + (right - left) / 2.0;
02482 
02483             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");
02484 
02485             std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(0.5);
02486 
02487             if (polySplit.second.coefficient(0) < m_epsilon)
02488             {
02489                 ltx("Bisection point is very near to a root?!?\n\n");
02490             }
02491 
02492             ltx("\\subsection{Recursion Branch " << mark << " 1 on the First Half " << ltxInterval(left,middle) << "}\n");
02493             ltx("\\label{" << m_markprefix << mark << " 1}\n");
02494 
02495             findRootsRecursive(polySplit.first, left, middle, mark + " 1", depth+1);
02496 
02497             ltx("\\subsection{Recursion Branch " << mark << " 2 on the Second Half " << ltxInterval(middle,right) << "}\n");
02498             ltx("\\label{" << m_markprefix << mark << " 2}\n");
02499 
02500             findRootsRecursive(polySplit.second, middle, right, mark + " 2", depth+1);
02501         }
02502         else
02503         {
02504             ltx(" $\\Longrightarrow$ Selective recursion: ");
02505             for (unsigned int i = 0; i < intervals.size(); ++i)
02506             {
02507 #if !SPEEDTEST
02508                 Numeric ileft = left + (right - left) * intervals[i].first;
02509                 Numeric iright = left + (right - left) * intervals[i].second;
02510 
02511                 ltx("\\hyperref[" << m_markprefix << mark << " " << i+1 << "]{interval " << i+1 << ": " << ltxInterval(ileft, iright) << "}, ");
02512 #endif
02513             }
02514             ltx("\n\n");
02515 
02516             for (unsigned int i = 0; i < intervals.size(); ++i)
02517             {
02518                 Numeric ileft = left + (right - left) * intervals[i].first;
02519                 Numeric iright = left + (right - left) * intervals[i].second;
02520 
02521                 if (intervals[i].first == 0) // only one de Casteljau split required
02522                 {
02523                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].second);
02524 
02525                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
02526                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
02527 
02528                     findRootsRecursive(polySplit.first, ileft, iright, mark + " " + S(i+1), depth+1);
02529                 }
02530                 else if (intervals[i].second == 1.0) // only one de Casteljau split required
02531                 {
02532                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit = p1.deCasteljauSplit(intervals[i].first);
02533 
02534                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
02535                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
02536 
02537                     findRootsRecursive(polySplit.second, ileft, iright, mark + " " + S(i+1), depth+1);
02538                 }
02539                 else // requires two de Casteljau splits
02540                 {
02541                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit1 = p1.deCasteljauSplit(intervals[i].first);
02542 
02543                     Numeric split2 = (intervals[i].second - intervals[i].first) / (1.0 - intervals[i].first);
02544 
02545                     std::pair<NPolynomialBezier,NPolynomialBezier> polySplit2 = polySplit1.second.deCasteljauSplit(split2);
02546 
02547                     ltx("\\subsection{Recursion Branch " << mark << " " << i+1 << " in Interval " << i+1 << ": " << ltxInterval(ileft,iright) << "}\n");
02548                     ltx("\\label{" << m_markprefix << mark << " " << i+1 << "}\n");
02549 
02550                     findRootsRecursive(polySplit2.first, ileft, iright, mark + " " + S(i+1), depth+1);
02551                 }
02552             }
02553         }
02554     }
02555 };
02556 
02557 /** \page testsuite A Suite of Test Cases
02558 
02559 *******************************************************************************
02560 *** A Suite of Test Cases                                                   ***
02561 *******************************************************************************
02562 
02563  In the testsuite() function many of the basic algorithm implemented above are
02564  checked against precomputed values.
02565 
02566 ******************************************************************************/
02567 
02568 #define ASSERT(x)                                                \
02569     do { if (!(x)) {                                             \
02570         std::cout << "Testsuite assertion \"" #x "\" on line "   \
02571                   << __LINE__ << " failed!\n";                   \
02572         abort();                                                 \
02573         }                                                        \
02574     } while(0)
02575 
02576 #define CHKSTR(p,str)                                           \
02577     do { if (p.toString() != str) {                             \
02578         std::cout << "Testsuite error: "                        \
02579                   << p.toString() << " != " << str << "\n";     \
02580         ASSERT(p.toString() == str);                            \
02581         }                                                       \
02582     } while(0)
02583 
02584 template <typename Numeric>
02585 void testsuite()
02586 {
02587     typedef PolynomialStandard<Numeric> NPolynomialStandard;
02588     typedef PolynomialBezier<Numeric> NPolynomialBezier;
02589 
02590     { // test Bezier transformation, evaluation and inverse transformation
02591 
02592         Numeric coeff[6] = { 1, 0, 0, 0, -16, 16 };
02593         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[6]));
02594 
02595         CHKSTR( p1, "16 X^5 - 16 X^4 + 1" );
02596 
02597         ASSERT( abs(p1.evaluate(0.0) - 1.0) < 0.00001 );
02598         ASSERT( abs(p1.evaluate(0.25) - 0.953125) < 0.00001 );
02599         ASSERT( abs(p1.evaluate(0.5) - 0.5) < 0.00001 );
02600         ASSERT( abs(p1.evaluate(0.75) + 0.265625) < 0.00001 );
02601         ASSERT( abs(p1.evaluate(1.0) - 1.0) < 0.00001 );
02602 
02603         NPolynomialBezier p2 = p1.toBezier();
02604 
02605         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}" );
02606 
02607         ASSERT( abs(p2.evaluate(0.0) - 1.0) < 0.00001 );
02608         ASSERT( abs(p2.evaluate(0.25) - 0.953125) < 0.00001 );
02609         ASSERT( abs(p2.evaluate(0.5) - 0.5) < 0.00001 );
02610         ASSERT( abs(p2.evaluate(0.75) + 0.265625) < 0.00001 );
02611         ASSERT( abs(p2.evaluate(1.0) - 1.0) < 0.00001 );
02612 
02613         NPolynomialStandard p3 = p2.toStandard();
02614 
02615         ASSERT( p3.toString() == "16 X^5 - 16 X^4 + 1" );
02616     }
02617 
02618     if (0) { // test roots calculation using Cardano's formulas: template
02619 
02620         Numeric coeff[4] = { -3.2, 20, -40, 25 };
02621         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02622 
02623         std::cout << "p = " << p1.toString() << "\n";
02624 
02625         std::vector<Numeric> roots = p1.findRoots();
02626 
02627         for (unsigned int i = 0; i < roots.size(); ++i)
02628             std::cout << "root[" << i << "] = " << roots[i] << "\n";
02629     }
02630 
02631     { // test roots calculation using Cardano's formulas: three real roots
02632 
02633         Numeric coeff[4] = { -3, 20, -40, 25 };
02634         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02635 
02636         CHKSTR( p1, "25 X^3 - 40 X^2 + 20 X - 3" );
02637 
02638         std::vector<Numeric> roots = p1.findRoots();
02639 
02640         ASSERT(roots.size() == 3);
02641         ASSERT( abs(roots[0] - 0.276393) < 0.000001 );
02642         ASSERT( abs(roots[1] - 0.6) < 0.000001 );
02643         ASSERT( abs(roots[2] - 0.723607) < 0.000001 );
02644     }
02645 
02646     { // test roots calculation using Cardano's formulas: one single real root
02647 
02648         Numeric coeff[4] = { -0.5, 4, -8, 5 };
02649         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02650 
02651         CHKSTR( p1, "5 X^3 - 8 X^2 + 4 X - 0.5" );
02652 
02653         std::vector<Numeric> roots = p1.findRoots();
02654 
02655         ASSERT(roots.size() == 1);
02656         ASSERT( abs(roots[0] - 0.186385) < 0.000001 );
02657     }
02658 
02659     { // test roots calculation using Cardano's formulas: case of a double root?
02660 
02661         Numeric coeff[4] = { -12, 80, -175, 125 };
02662         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02663 
02664         CHKSTR( p1, "125 X^3 - 175 X^2 + 80 X - 12" );
02665 
02666         std::vector<Numeric> roots1 = p1.findRoots();
02667 
02668 #if 0
02669         ASSERT( roots1.size() == 3 || roots1.size() == 2 );
02670         if (roots1.size() == 3) {
02671             ASSERT( abs(roots1[0] - 0.4) < 0.000001 );
02672             ASSERT( abs(roots1[1] - 0.4) < 0.000001 );
02673             ASSERT( abs(roots1[2] - 0.6) < 0.000001 );
02674         }
02675         else {
02676             ASSERT( abs(roots1[0] - 0.4) < 0.000001 );
02677             ASSERT( abs(roots1[1] - 0.6) < 0.000001 );
02678         }
02679 #endif
02680     }
02681 
02682     { // test roots calculation using Cardano's formulas: case of a triple root
02683 
02684         Numeric coeff[4] = { -8, 60, -150, 125 };
02685         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02686 
02687         CHKSTR( p1, "125 X^3 - 150 X^2 + 60 X - 8" );
02688 
02689         std::vector<Numeric> roots = p1.findRoots();
02690 
02691         ASSERT(roots.size() == 1);
02692         ASSERT( abs(roots[0] - 0.4) < 0.000001 );
02693     }
02694 
02695     { // test roots of quadratic Bezier and monomial polynomials
02696 
02697         Numeric coeff[3] = { 0.1, -1, 2 };
02698         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[3]));
02699 
02700         CHKSTR( p1, "2 X^2 - 1 X + 0.1" );
02701 
02702         std::vector<Numeric> roots1 = p1.findRoots();
02703 
02704         ASSERT( roots1.size() == 2 );
02705         ASSERT( abs(roots1[0] - 0.138197) < 0.000001 );
02706         ASSERT( abs(roots1[1] - 0.361803) < 0.000001 );
02707 
02708         NPolynomialBezier p2 = p1.toBezier();
02709 
02710         CHKSTR( p2, "0.1 B_{0,2} - 0.4 B_{1,2} + 1.1 B_{2,2}" );
02711 
02712         std::vector<Numeric> roots2 = p2.findRoots();
02713 
02714         ASSERT( roots2.size() == 2 );
02715         ASSERT( abs(roots2[0] - 0.138197) < 0.000001 );
02716         ASSERT( abs(roots2[1] - 0.361803) < 0.000001 );
02717     }
02718 
02719     { // test roots of cubic monomial and Bezier polynomials
02720 
02721         Numeric coeff[4] = { -2, 36, -96, 64 };
02722         NPolynomialStandard p1 (std::vector<Numeric>(&coeff[0], &coeff[4]));
02723 
02724         CHKSTR( p1, "64 X^3 - 96 X^2 + 36 X - 2" );
02725 
02726         std::vector<Numeric> roots1 = p1.findRoots();
02727 
02728         ASSERT( roots1.size() == 3 );
02729         ASSERT( abs(roots1[0] - 0.066987) < 0.000001 );
02730         ASSERT( abs(roots1[1] - 0.5) < 0.000001 );
02731         ASSERT( abs(roots1[2] - 0.933012) < 0.000001 );
02732 
02733         NPolynomialBezier p2 = p1.toBezier();
02734 
02735         CHKSTR( p2, "- 2 B_{0,3} + 10 B_{1,3} - 10 B_{2,3} + 2 B_{3,3}" );
02736 
02737         std::vector<Numeric> roots2 = p2.findRoots();
02738 
02739         ASSERT( roots2.size() == 3 );
02740         ASSERT( abs(roots2[0] - 0.066987) < 0.000001 );
02741         ASSERT( abs(roots2[1] - 0.5) < 0.000001 );
02742         ASSERT( abs(roots2[2] - 0.933012) < 0.000001 );
02743     }
02744 
02745     { // interpolation coefficients: scalar product of Bernstein polynomials.
02746 
02747         Numeric testValues[6][3] = {
02748             { 0.125,      0.0357143, 0.00595238 },
02749             { 0.0892857,  0.0595238, 0.0178571 },
02750             { 0.0595238,  0.0714286, 0.0357143 },
02751             { 0.0357143,  0.0714286, 0.0595238 },
02752             { 0.0178571,  0.0595238, 0.0892857 },
02753             { 0.00595238, 0.0357143, 0.125 }
02754         };
02755 
02756         for (unsigned int i = 0; i < 6; ++i) {
02757             for (unsigned int j = 0; j < 3; ++j) {
02758                 ASSERT( abs(DegreeInterpolationGenerator<Numeric>::scalarProduktBernsteinPolynomial(i, 5, j, 2)
02759                             - testValues[i][j]) < 0.000001 );
02760             }
02761         }
02762 
02763         // debug code to print out a matrix of scalar products:
02764         //std::cout << Matrix<Numeric>(6,3).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::scalarProduktBernsteinPolynomial(i, 5, j, 2); }).toString();
02765     }
02766 
02767     { // interpolation coefficients: check coefficients of dual basis.
02768 
02769         Numeric testValues[3][3] = {
02770             {  16, -24,         16 },
02771             { -24,  69.333333, -57.333333 },
02772             {  16, -57.333333,  69.333333 }
02773         };
02774 
02775         for (unsigned int i = 0; i < 3; ++i) {
02776             for (unsigned int j = 0; j < 3; ++j) {
02777                 ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientDualBasis(i, j, 3)
02778                             - testValues[i][j]) < 0.000001 );
02779             }
02780         }
02781 
02782         // debug code to print out a matrix of dual basis coefficients:
02783         //std::cout << Matrix<Numeric>(4,4).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientDualBasis(i, j, 3); }).toString();
02784     }
02785 
02786     { // interpolation coefficients: check beta coefficients of degree
02787       // reducing matrix. The values used here are the same as in Example 1
02788       // and Example 2 of [1].
02789 
02790         Numeric testValues1[6][3] = {
02791             {  0.821429, -0.428571,  0.107143 },
02792             {  0.321429,  0.285714, -0.107143 },
02793             {  0,         0.642857, -0.142857 },
02794             { -0.142857,  0.642857,  0 },
02795             { -0.107143,  0.285714,  0.321429 },
02796             {  0.107143, -0.428571,  0.821429 }
02797         };
02798 
02799         for (unsigned int i = 0; i < 3; ++i) {
02800             for (unsigned int j = 0; j < 3; ++j) {
02801                 ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 5, 2)
02802                             - testValues1[i][j]) < 0.000001 );
02803             }
02804         }
02805 
02806         Numeric testValues2[3][6] = {
02807             { 1, 0.6, 0.3, 0.1, 0,   0 },
02808             { 0, 0.4, 0.6, 0.6, 0.4, 0 },
02809             { 0, 0,   0.1, 0.3, 0.6, 1 }
02810         };
02811 
02812         for (unsigned int i = 0; i < 3; ++i) {
02813             for (unsigned int j = 0; j < 3; ++j) {
02814                 ASSERT( abs(DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 2, 5)
02815                             - testValues2[i][j]) < 0.000001 );
02816             }
02817         }
02818 
02819         // debug code to print out a matrix of beta coefficients:
02820         //std::cout << Matrix<Numeric>(6,3).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 5, 2); }).toString();
02821         //std::cout << Matrix<Numeric>(3,6).fill([](int i, int j) { return DegreeInterpolationGenerator<Numeric>::CoefficientBeta(i, j, 2, 5); }).toString();
02822     }
02823 
02824     { // test polynomial generator from roots
02825 
02826         Numeric coeff2[2] = { 1.0/3.0, 3 };
02827         NPolynomialStandard p2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));
02828 
02829         Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
02830         NPolynomialStandard p4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));
02831 
02832         Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
02833         NPolynomialStandard p8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));
02834 
02835         Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
02836         NPolynomialStandard p16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));
02837 
02838         ASSERT( p2.toString() == "- 1 X^2 + 3.33333 X - 1" );
02839 
02840         ASSERT( p4.toString() == "- 1 X^4 - 7.66667 X^3 - 2.33333 X^2 + 51.6667 X - 16.6667" );
02841 
02842         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" );
02843 
02844         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}" );
02845     }
02846 }
02847 
02848 /** \page demos Demos and Speedtests
02849 
02850 *******************************************************************************
02851 *** Demos and Speedtests                                                    ***
02852 *******************************************************************************
02853 
02854  To demonstrate and evaluate the root finding implementations the program
02855  contains 8 predefined demos, which contain different tests and runs of groups
02856  of polynomials.
02857 
02858  Demos 1-5 are simple examples used to illustrate the workings of the
02859  algorithms.
02860 
02861  Demos 6-8 are speedtests which will repeat the algorithms a large number of
02862  times and calculate their average speed for three different sets of
02863  polynomials.
02864 
02865 ******************************************************************************/
02866 
02867 template <typename Clipper, typename Numeric>
02868 double runSpeedtest(unsigned int iterations, Clipper& clipper, const PolynomialStandard<Numeric>& p, double left, double right)
02869 {
02870 #if !SPEEDTEST
02871     iterations = 1;
02872 #endif
02873 
02874     double ts1 = timestamp();
02875 
02876     for (unsigned int i = 0; i < iterations; ++i)
02877     {
02878         clipper.findRoots(p, left, right);
02879     }
02880 
02881     double ts2 = timestamp();
02882     return ts2 - ts1;
02883 }
02884 
02885 template <typename Numeric>
02886 void runAllClippers(const std::string& name, unsigned int iterations, const PolynomialStandard<Numeric>& p, int epsexp, double left = 0.0, double right = 1.0)
02887 {
02888     Numeric epsilon(10);
02889     epsilon = pow(epsilon, -epsexp);
02890 
02891     std::string ltxlabel = name + "-" + Numeric::getName() + "-" + S(epsexp);
02892 
02893     spd(std::setprecision(12));
02894 
02895     BezierClip<Numeric> c1(epsilon, ltxlabel + "b:");
02896     ltx("\\section{Running \\texttt{BezClip} on " << name << " with epsilon " << epsexp << "}\n\n");
02897     double spd1 = runSpeedtest(iterations, c1, p, left, right);
02898     spd("algo=BezClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
02899         " iterations=" << iterations << " maxdepth=" << c1.maxdepth() <<
02900         " totaltime=" << spd1 << " speed=" << (spd1 / iterations) << "\n");
02901     ltx("\\clearpage\n");
02902 
02903     KClip<Numeric> c2(p.degree(), 2, epsilon, ltxlabel + "q:");
02904     ltx("\\section{Running \\texttt{QuadClip} on " << name << " with epsilon " << epsexp << "}\n\n");
02905     double spd2 = runSpeedtest(iterations, c2, p, left, right);
02906     spd("algo=QuadClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
02907         " iterations=" << iterations << " maxdepth=" << c2.maxdepth() <<
02908         " totaltime=" << spd2 << " speed=" << (spd2 / iterations) << "\n");
02909     ltx("\\clearpage\n");
02910 
02911     KClip<Numeric> c3(p.degree(), 3, epsilon, ltxlabel + "c:");
02912     ltx("\\section{Running \\texttt{CubeClip} on " << name << " with epsilon " << epsexp << "}\n\n");
02913     double spd3 = runSpeedtest(iterations, c3, p, left, right);
02914     spd("algo=CubeClip input=" << name << " numeric=" << Numeric::getName() << " epsilon=" << epsexp <<
02915         " iterations=" << iterations << " maxdepth=" << c3.maxdepth() <<
02916         " totaltime=" << spd3 << " speed=" << (spd3 / iterations) << "\n");
02917     ltx("\\clearpage\n");
02918 }
02919 
02920 template <typename Numeric>
02921 void runDemo6()
02922 {
02923     typedef PolynomialStandard<Numeric> NPolynomialStandard;
02924 
02925     spd("Single root functions from Barton and Juettler (Example 9):\n");
02926     spd("Numeric = " << Numeric::getName() << "\n");
02927 
02928     Numeric coeff2[2] = { 1.0/3.0, 3 };
02929     NPolynomialStandard f2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));
02930 
02931     Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
02932     NPolynomialStandard f4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));
02933 
02934     Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
02935     NPolynomialStandard f8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));
02936 
02937     Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
02938     NPolynomialStandard f16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));
02939 
02940     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02941         runAllClippers("$f_2$", 50000 * Numeric::IterationScale, f2, epsexp);
02942 
02943     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02944         runAllClippers("$f_4$", 20000 * Numeric::IterationScale, f4, epsexp);
02945 
02946     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02947         runAllClippers("$f_8$", 10000 * Numeric::IterationScale, f8, epsexp);
02948 
02949     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02950         runAllClippers("$f_{16}$", 1000 * Numeric::IterationScale, f16, epsexp);
02951 }
02952 
02953 template <typename Numeric>
02954 void runDemo7()
02955 {
02956     typedef PolynomialStandard<Numeric> NPolynomialStandard;
02957 
02958     spd("Triple root functions from Liu et al. (Example 3):\n");
02959     spd("Numeric = " << Numeric::getName() << "\n");
02960 
02961     Numeric coeff4[4] = { 1.0/3.0, 1.0/3.0, 1.0/3.0, 5 };
02962     NPolynomialStandard g4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));
02963 
02964     Numeric coeff8[8] = { 1.0/3.0, 1.0/3.0, 1.0/3.0, -2, -2, -2, 5, 5 };
02965     NPolynomialStandard g8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));
02966 
02967     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 };
02968     NPolynomialStandard g16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));
02969 
02970     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02971         runAllClippers("$g_4$", 20000 * Numeric::IterationScale, g4, epsexp);
02972 
02973     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02974         runAllClippers("$g_8$", 10000 * Numeric::IterationScale, g8, epsexp);
02975 
02976     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
02977         runAllClippers("$g_{16}$", 1000 * Numeric::IterationScale, g16, epsexp);
02978 }
02979 
02980 template <typename Numeric>
02981 void runDemo8()
02982 {
02983     typedef PolynomialStandard<Numeric> NPolynomialStandard;
02984 
02985     spd("Near double root functions from Barton and Juettler (Example 11):\n");
02986     spd("Numeric = " << Numeric::getName() << "\n");
02987 
02988     print_precision = 16;
02989 
02990     Numeric coeff2[2] = { 0.56, 0.57 };
02991     NPolynomialStandard h2 (1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));
02992 
02993     Numeric coeff4[4] = { 0.4, 0.40000001, -1, 2 };
02994     NPolynomialStandard h4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));
02995 
02996     Numeric coeff8[8] = { 0.50000002, 0.50000003, -5, -5, -5, -7, -7, -7 };
02997     NPolynomialStandard h8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));
02998 
02999     Numeric coeff16[16] = { 0.30000008, 0.30000009, 6, 6, 6, 6, 6, 6, 6, -5, -5, -5, -5, -5, -5, -7 };
03000     NPolynomialStandard h16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));
03001 
03002     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
03003         runAllClippers("$h_2$", 50000 * Numeric::IterationScale, h2, epsexp);
03004 
03005     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
03006         runAllClippers("$h_4$", 20000 * Numeric::IterationScale, h4, epsexp);
03007 
03008     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
03009         runAllClippers("$h_8$", 10000 * Numeric::IterationScale, h8, epsexp);
03010 
03011     for (unsigned int epsexp = 2; epsexp <= 128; epsexp *= 2)
03012         runAllClippers("$h_{16}$", 1000 * Numeric::IterationScale, h16, epsexp);
03013 }
03014 
03015 void runDemo(unsigned int demo)
03016 {
03017     ltx("{\\Large\\bf\\centerline{Demo " << demo << "}\\par}\n\\bigskip\n");
03018 
03019     if (demo == 1)
03020     {
03021         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");
03022 
03023         ltx("This demo works with the standard datatype \\texttt{double} and precision $\\varepsilon = 0.001$.\n");
03024 
03025         ltx("\\tableofcontents\n\n"
03026             "\\clearpage\n");
03027 
03028         typedef Double<double>  Numeric;
03029         Numeric epsilon = 0.001;
03030 
03031         //Numeric coeff1[6] = { 1, -2, -1, 2, 0.5, 2 };
03032         Numeric coeff1[6] = { 1, -2, -1, 2.5, 0, 1 };
03033         PolynomialBezier<Numeric> p1 (std::vector<Numeric>(&coeff1[0], &coeff1[6]));
03034 
03035         debug_more = true;
03036 
03037         ltx("\\section{\\texttt{BezClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
03038         BezierClip<Numeric>(epsilon, "p1b:").findRoots(p1.toStandard(), 0.0, 1.0);
03039         ltx("\\clearpage\n");
03040 
03041         ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
03042         KClip<Numeric>(p1.degree(), 2, epsilon, "p1q:").findRoots(p1.toStandard(), 0.0, 1.0);
03043         ltx("\\clearpage\n");
03044 
03045         ltx("\\section{\\texttt{CubeClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
03046         KClip<Numeric>(p1.degree(), 3, epsilon, "p1c:").findRoots(p1.toStandard(), 0.0, 1.0);
03047         ltx("\\clearpage\n");
03048     }
03049     else if (demo == 2)
03050     {
03051         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");
03052 
03053         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");
03054 
03055         ltx("\\tableofcontents\n\n"
03056             "\\clearpage\n");
03057 
03058         typedef Double<long double>     Numeric;
03059         Numeric epsilon = 0.001;
03060 
03061         Numeric coeff2[9] = { 1, -5, 8, -8, 8, -10, 9, -4, 1 };
03062         PolynomialBezier<Numeric> p8 (std::vector<Numeric>(&coeff2[0], &coeff2[9]));
03063 
03064         ltx("\\section{\\texttt{BezClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
03065         BezierClip<Numeric>(epsilon, "p8b:").findRoots(p8.toStandard(), 0.0, 1.0);
03066         ltx("\\clearpage\n");
03067 
03068         ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
03069         KClip<Numeric>(p8.degree(), 2, epsilon, "p8q:").findRoots(p8.toStandard(), 0.0, 1.0);
03070         ltx("\\clearpage\n");
03071 
03072         ltx("\\section{\\texttt{CubeClip} Applied to a Polynomial of 8th Degree with Six Roots}\n\n");
03073         KClip<Numeric>(p8.degree(), 3, epsilon, "p8c:").findRoots(p8.toStandard(), 0.0, 1.0);
03074         ltx("\\clearpage\n");
03075     }
03076     else if (demo == 3)
03077     {
03078         typedef Double<double>  Numeric;
03079         Numeric epsilon = 0.001;
03080 
03081         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");
03082 
03083         ltx("\\tableofcontents\n\n"
03084             "\\clearpage\n");
03085 
03086         Numeric coeff20[20] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20 };
03087         PolynomialStandard<Numeric> p20 (1.0, std::vector<Numeric>(&coeff20[0], &coeff20[20]));
03088 
03089         ltx("\\section{\\texttt{BezClip} Applied to the Wilkinson Polynomial}\n\n");
03090         BezierClip<Numeric>(epsilon, "p20b:").findRoots(p20, 0.0, 25.0);
03091         ltx("\\clearpage\n");
03092 
03093         ltx("\\section{\\texttt{QuadClip} Applied to the Wilkinson Polynomial}\n\n");
03094         KClip<Numeric>(p20.degree(), 2, epsilon, "p20q:").findRoots(p20, 0.0, 25.0);
03095         ltx("\\clearpage\n");
03096 
03097         ltx("\\section{\\texttt{CubeClip} Applied to the Wilkinson Polynomial}\n\n");
03098         KClip<Numeric>(p20.degree(), 3, epsilon, "p20c:").findRoots(p20, 0.0, 25.0);
03099         ltx("\\clearpage\n");
03100     }
03101     else if (demo == 4)
03102     {
03103         typedef Double<long double>             Numeric;
03104         typedef PolynomialStandard<Numeric> NPolynomialStandard;
03105 
03106         ltx("This demo entry is used to test out further polynomials.\n\n");
03107 
03108         ltx("\\tableofcontents\n\n"
03109             "\\clearpage\n");
03110 
03111         Numeric coeff2[2] = { 1.0/3.0, 3 };
03112         NPolynomialStandard p2 (-1.0, std::vector<Numeric>(&coeff2[0], &coeff2[2]));
03113 
03114         Numeric coeff4[4] = { 1.0/3.0, 2, -5, -5 };
03115         NPolynomialStandard p4 (-1.0, std::vector<Numeric>(&coeff4[0], &coeff4[4]));
03116 
03117         Numeric coeff8[8] = { 1.0/3.0, 2, 2, 2, -5, -5, -5, -5 };
03118         NPolynomialStandard p8 (-1.0, std::vector<Numeric>(&coeff8[0], &coeff8[8]));
03119 
03120         Numeric coeff16[16] = { 1.0/3.0, 2, 2, 2, 2, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5 };
03121         NPolynomialStandard p16 (-1.0, std::vector<Numeric>(&coeff16[0], &coeff16[16]));
03122 
03123         NPolynomialStandard p = p16;
03124 
03125         Numeric epsilon(10);
03126         epsilon = pow(epsilon, -32);
03127 
03128         ltx("\\section{\\texttt{BezClip} Applied to the Example Polynomial}\n\n");
03129         BezierClip<Numeric>(epsilon, "b:").findRoots(p, 0.0, 1.0);
03130         ltx("\\clearpage\n");
03131 
03132         ltx("\\section{\\texttt{QuadClip} Applied to the Example Polynomial}\n\n");
03133         KClip<Numeric>(p.degree(), 2, epsilon, "q:").findRoots(p, 0.0, 1.0);
03134         ltx("\\clearpage\n");
03135 
03136         ltx("\\section{\\texttt{CubeClip} Applied to the Example Polynomial}\n\n");
03137         KClip<Numeric>(p.degree(), 3, epsilon, "c:").findRoots(p, 0.0, 1.0);
03138         ltx("\\clearpage\n");
03139     }
03140     else if (demo == 5)
03141     {
03142         typedef Double<double>  Numeric;
03143         Numeric epsilon = 0.001;
03144 
03145         ltx("This demo contains the calculation procedures to generate the results shown in Appendix A of the lab report.\n\n");
03146 
03147         ltx("\\section{Degree Reduction and Raising Matrices for Degree $5$}\n\n");
03148         for (int i = 4; i >= 0; --i)
03149         {
03150             Matrix<Fraction> M1(6,i+1);
03151             M1.fill( DegreeInterpolationGenerator<Fraction>(5,i) );
03152 
03153             Matrix<Fraction> M2(i+1,6);
03154             M2.fill( DegreeInterpolationGenerator<Fraction>(i,5) );
03155 
03156             ltx("\\begin{align*}\n"
03157                 "M_{5," << i << "} = " << M1.toString() << " &&\n"
03158                 "M_{" << i << ",5} = " << M2.toString() << "\n"
03159                 "\\end{align*}\n");
03160         }
03161         ltx("\\clearpage\n");
03162 
03163         Numeric coeff1[6] = { 1, -2, -1, 2.5, 0, 1 };
03164         PolynomialBezier<Numeric> p1 (std::vector<Numeric>(&coeff1[0], &coeff1[6]));
03165 
03166         debug_more = true;
03167 
03168         ltx("\\section{\\texttt{QuadClip} Applied to a Polynomial of 5th Degree with Two Roots}\n\n");
03169         KClip<Numeric>(p1.degree(), 2, epsilon, "p1q:").findRoots(p1.toStandard(), 0.0, 1.0);
03170         ltx("\\clearpage\n");
03171     }
03172     else if (demo == 6)
03173     {
03174         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");
03175 
03176         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"
03177             "\\begin{align*}\n"
03178             "f_2 &:= (t - \\tfrac{1}{3}) (3 - t) \\\\\n"
03179             "f_4 &:= (t - \\tfrac{1}{3}) (2 - t) (t + 5)^2 \\\\\n"
03180             "f_8 &:= (t - \\tfrac{1}{3}) (2 - t)^3 (t + 5)^4 \\\\\n"
03181             "f_{16} &:= (t - \\tfrac{1}{3}) (2 - t)^5 (t + 5)^{10}\n"
03182             "\\end{align*}\n\n");
03183 
03184         ltx("\\tableofcontents\n\n"
03185             "\\clearpage\n");
03186 
03187         ltx("\\part{Numeric = double}\n");
03188         runDemo6< Double<double> >();
03189 
03190         ltx("\\part{Numeric = long double}\n");
03191         runDemo6< Double<long double> >();
03192 
03193         ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
03194         mpfr_set_default_prec(1024);
03195         runDemo6< MpfrFloat >();
03196     }
03197     else if (demo == 7)
03198     {
03199         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");
03200 
03201         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"
03202             "\\begin{align*}\n"
03203             "g_4 &:= (t - \\tfrac{1}{3})^3 (t - 5) \\\\\n"
03204             "g_8 &:= (t - \\tfrac{1}{3})^3 (2 + t)^3 (t - 5)^2 \\\\\n"
03205             "g_{16} &:= (t - \\tfrac{1}{3})^3 (2 + t)^2 (t - 5)^7 (t + 7)^4\n"
03206             "\\end{align*}\n\n");
03207 
03208         ltx("\\tableofcontents\n\n"
03209             "\\clearpage\n");
03210 
03211         //runDemo7< Double<double> >();
03212         //runDemo7< Double<long double> >();
03213 
03214         ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
03215         mpfr_set_default_prec(1024);
03216         runDemo7< MpfrFloat >();
03217     }
03218     else if (demo == 8)
03219     {
03220         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");
03221 
03222         ltx("In this demo or speed test the following polynomials with two near roots are used (Example 3 from Liu et al.):\n"
03223             "\\begin{align*}\n"
03224             "h_2 &:= (t - 0.56)(t - 0.57) \\\\\n"
03225             "h_4 &:= (t - 0.4) (t - 0.40000001) (t+1) (2-t) \\\\\n"
03226             "h_8 &:= (t - 0.50000002) (t - 0.50000003) (t+5)^3 (t+7)^3 \\\\\n"
03227             "h_{16} &:= (t - 0.30000008) (t - 0.30000009) (6-t)^7 (t+5)^6 (t+7)\n"
03228             "\\end{align*}\n\n");
03229 
03230         ltx("\\tableofcontents\n\n"
03231             "\\clearpage\n");
03232 
03233         //runDemo7< Double<double> >();
03234         //runDemo7< Double<long double> >();
03235 
03236         ltx("\\part{Numeric = MpfrFloat with precision 1024}\n");
03237         mpfr_set_default_prec(1024);
03238         runDemo8< MpfrFloat >();
03239     }
03240     else if (demo == 9)
03241     {
03242         typedef Double<long double>             Numeric;
03243         typedef PolynomialStandard<Numeric> NPolynomialStandard;
03244 
03245         ltx("This demo entry is used to test out further polynomials.\n\n");
03246 
03247         ltx("\\tableofcontents\n\n"
03248             "\\clearpage\n");
03249 
03250         std::vector<NPolynomialStandard> polys;
03251 
03252         Numeric coeff3[4] = { -1.0/16, +5.0/8.0, -3.0/2.0, 1.0 };
03253         polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff3[0], &coeff3[4])));
03254 
03255         Numeric coeff4[5] = { +5.0/256.0, -5.0/16.0, +21.0/16.0, -2.0, 1.0 };
03256         polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff4[0], &coeff4[5])));
03257 
03258         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 };
03259         polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff5[0], &coeff5[6])));
03260 
03261         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 };
03262         polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff6[0], &coeff6[7])));
03263 
03264         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 };
03265         polys.push_back(NPolynomialStandard(std::vector<Numeric>(&coeff7[0], &coeff7[8])));
03266 
03267         for (unsigned int i = 0; i < polys.size(); ++i)
03268         {
03269             runAllClippers("p"+S(i+3), 1000 * Numeric::IterationScale, polys[i], 6);
03270         }
03271     }
03272 }
03273 
03274 // ****************************************************************************
03275 // *** main() with Command Line Parser                                      ***
03276 // ****************************************************************************
03277 
03278 void printLatexPreamble()
03279 {
03280     ltx("\\documentclass[a4paper,12pt]{article}\n"
03281 
03282         "\\usepackage[latin1]{inputenc}\n"
03283         "\\usepackage[T1]{fontenc}\n"
03284 
03285         "\\usepackage{lmodern}\n"
03286         "\\usepackage[english]{babel}\n"
03287 
03288         "\\usepackage[tmargin=20mm,bmargin=20mm,lmargin=15mm,rmargin=15mm]{geometry}\n"
03289 
03290         "\\parskip=0pt\n"
03291         "\\parindent=0pt\n"
03292         "\\hbadness=10000\n"
03293 
03294         "\\usepackage{amsmath,amssymb,breqn}\n"
03295 
03296         "\\usepackage[final,colorlinks=true,linkcolor=blue!60!black]{hyperref}\n"
03297 
03298         "\\usepackage{pgfplots}\n"
03299 
03300         "\\pgfplotsset{"
03301         "  width=16cm,height=6cm,\n"
03302         "  axis x line=middle,axis y line=left,\n"
03303         "  domain=0:1,xmin=0,xmax=1,\n"
03304         "  samples=50,\n"
03305         "  legend style={at={(0.5,-0.1)}, anchor=north}, legend columns=-1}\n"
03306 
03307         "\\def\\arraystretch{1.2}\n"
03308 
03309         "% special \\iffinal construct to exclude large drawings\n"
03310         "\\newif\\iffinal\n"
03311         "\\finaltrue\n"
03312 
03313         "% enlarge space for numbering in toc\n"
03314         "\\makeatletter"
03315         "\\renewcommand\\@pnumwidth{1.85em}\n"
03316         "\\def\\sectionnumwidth#1{"
03317         "\\gdef\\numberline##1{\\hb@xt@ #1{##1\\hfil}\\hskip 1ex}"
03318         "}"
03319         "\\sectionnumwidth{6ex}\n"
03320         "\\makeatother\n"
03321 
03322         "\\begin{document}\n");
03323 }
03324 
03325 void printUsage(char* argv[])
03326 {
03327     std::cerr <<
03328         "Usage: " << argv[0] << " [options] <coefficients>\n"
03329         "Options:\n"
03330         "  -B, --bezier         Coefficients on command line are for Bezier representation.\n"
03331         "  -D, --demo <n>       Run demo of different polynomials.\n"
03332         "  -M, --monomial       Coefficients on command line are for monomial representation.\n"
03333         "  -R, --roots <Scale>  Construct a polynomial with roots given as parameters and scale.\n"
03334         "  -a, --algo <Algo>    Select algorithm to run: Bezclip, Quadclip or Cubeclip.\n"
03335         "  -e, --epsilon <e>    Calculate until precision epsilon is reached.\n"
03336         "  -f, --fragment       Suppress printing the LaTeX preamble.\n"
03337         "  -h, --help           The help on command line parameters you are reading.\n"
03338         "  -i, --interval <a,b> Calculate roots in interval [a,b].\n"
03339         "  -l, --latex          Output LaTeX code of calculations (if compiled in).\n"
03340         "  -t, --testsuite      Run testsuite of checks.\n"
03341         ;
03342 }
03343 
03344 bool isCoefficient(const char* s)
03345 {
03346     if (!*s) return false;
03347 
03348     if (*s == '-' || *s == '+') ++s;    // scan sign
03349 
03350     if (*s == '-') return false;        // -e and --epsilon parameters
03351     if (*s == 'e') return false;
03352 
03353     while (*s)
03354     {
03355         if (!isdigit(*s) &&
03356             *s != '-' && *s != '+' && *s != 'e' && *s != '.')
03357         {
03358             return false;
03359         }
03360         ++s;
03361     }
03362     return true;
03363 }
03364 
03365 int main(int argc, char* argv[])
03366 {
03367     mpfr_set_default_prec (1024);
03368 
03369     char* endptr;
03370     bool badinput = false;
03371 
03372     // *** process command line parameters
03373 
03374     enum { ALGO_BEZCLIP, ALGO_QUADCLIP, ALGO_CUBECLIP } opt_algo = ALGO_QUADCLIP;
03375 
03376     enum { INPUT_BEZIER, INPUT_MONOMIAL, INPUT_ROOTS } opt_input = INPUT_BEZIER;
03377 
03378     bool opt_output_latexpreamble = true;
03379 
03380     unsigned int opt_rundemo = 0;
03381 
03382     typedef Double<long double> Numeric;
03383 
03384     Numeric opt_input_roots_scale = 1.0;
03385 
03386     Numeric opt_interval_left = 0.0;
03387     Numeric opt_interval_right = 1.0;
03388 
03389     Numeric opt_epsilon = 0.001;
03390 
03391     std::vector<Numeric> coefficients;
03392 
03393     while (1) {
03394         static struct option long_options[] = {
03395             { "algo",     required_argument,    0,  'a' },
03396             { "bezier",   no_argument,          0,  'B' },
03397             { "demo",     required_argument,    0,  'D' },
03398             { "epsilon",  required_argument,    0,  'e' },
03399             { "fragment", no_argument,          0,  'f' },
03400             { "help",     no_argument,          0,  'h' },
03401             { "interval", required_argument,    0,  'i' },
03402             { "latex",    no_argument,          0,  'l' },
03403             { "monomial", no_argument,          0,  'M' },
03404             { "roots",    required_argument,    0,  'R' },
03405             { "testsuite",no_argument,          0,  't' },
03406             { 0,          0,                    0,  0 }
03407         };
03408 
03409         // check parameter for coefficent syntax before getopt_long can misparse
03410         // negative numbers
03411 
03412         while ( optind < argc && isCoefficient(argv[optind]) )
03413         {
03414             Numeric coeff = strtod(argv[optind], &endptr);
03415 
03416             if (endptr && *endptr == 0)
03417                 coefficients.push_back(coeff);
03418             else
03419             {
03420                 std::cerr << "Invalid coefficient \"" << argv[optind] << "\"\n";
03421                 badinput = true;
03422             }
03423 
03424             optind++;
03425         }
03426 
03427         int option_index = 0;
03428         int c = getopt_long(argc, argv, "h?lMBR:a:fe:i:Dt",
03429                             long_options, &option_index);
03430 
03431         if (c == -1) break;
03432 
03433         switch (c) {
03434         case 'l':
03435             output_latex = true;
03436             break;
03437 
03438         case 'f':
03439             opt_output_latexpreamble = false;
03440             break;
03441 
03442         case 'D':
03443             opt_rundemo = strtoul(optarg, &endptr, 10);
03444             if (!endptr || *endptr != 0) {
03445                 std::cout << "Invalid demo number \"" << optarg << "\"\n";
03446                 badinput = true;
03447             }
03448             break;
03449 
03450         case 'a':
03451             if (optarg[0] == 'b') {
03452                 opt_algo = ALGO_BEZCLIP;
03453             }
03454             else if (optarg[0] == 'q') {
03455                 opt_algo = ALGO_QUADCLIP;
03456             }
03457             else if (optarg[0] == 'c') {
03458                 opt_algo = ALGO_CUBECLIP;
03459             }
03460             else {
03461                 std::cerr << "Unknown algorithm \"" << optarg << "\" requested.\n";
03462             }
03463             break;
03464 
03465         case 'M':
03466             opt_input = INPUT_MONOMIAL;
03467             break;
03468 
03469         case 'B':
03470             opt_input = INPUT_BEZIER;
03471             break;
03472 
03473         case 'R':
03474             opt_input = INPUT_ROOTS;
03475 
03476             opt_input_roots_scale = strtod(optarg, &endptr);
03477             if (!endptr || *endptr != 0) {
03478                 std::cout << "Invalid scale parameter \"" << optarg << "\"\n";
03479                 badinput = true;
03480             }
03481             break;
03482 
03483         case 'i':
03484             opt_interval_left = strtod(optarg, &endptr);
03485             if (!endptr || *endptr != ',') {
03486                 std::cout << "Invalid interval parameter \"" << optarg << "\"\n";
03487                 badinput = true;
03488                 break;
03489             }
03490 
03491             opt_interval_right = strtod(endptr+1, &endptr);
03492             if (!endptr || *endptr != 0) {
03493                 std::cout << "Invalid interval parameter \"" << optarg << "\"\n";
03494                 badinput = true;
03495             }
03496             break;
03497 
03498         case 'e':
03499             opt_epsilon = strtod(optarg, &endptr);
03500             if (!endptr || *endptr != 0) {
03501                 std::cout << "Invalid epsilon parameter \"" << optarg << "\"\n";
03502                 badinput = true;
03503             }
03504             break;
03505 
03506         case 't':
03507             (std::cout << "Running test suite: ").flush();
03508 
03509             testsuite< Double<double> >();
03510             testsuite< Double<long double> >();
03511             testsuite< MpfrFloat >();
03512             std::cout << "done.\n";
03513 
03514             return 0;
03515 
03516         case 'h':
03517         case '?':
03518             printUsage(argv);
03519             return 0;
03520 
03521         default:
03522             std::cerr << "?? getopt returned character code " << c << " ???\n";;
03523         }
03524     }
03525 
03526     if (badinput) return 0;
03527 
03528     if (coefficients.size() == 0 && !opt_rundemo)
03529     {
03530         printUsage(argv);
03531         return 0;
03532     }
03533 
03534     // *** main program directed by command line input
03535 
03536     if (opt_output_latexpreamble)
03537         printLatexPreamble();
03538 
03539     if (opt_rundemo) {
03540         runDemo(opt_rundemo);
03541 
03542         if (opt_output_latexpreamble)
03543             ltx("\\end{document}\n");
03544         return 0;
03545     }
03546 
03547     PolynomialStandard<Numeric> p;
03548 
03549     // create polynomial depending on input parameters
03550     if (opt_input == INPUT_BEZIER)
03551     {
03552         PolynomialBezier<Numeric> pB (coefficients);
03553         p = pB.toStandard();
03554     }
03555     else if (opt_input == INPUT_MONOMIAL)
03556     {
03557         p = PolynomialStandard<Numeric>(coefficients);
03558     }
03559     else if (opt_input == INPUT_ROOTS)
03560     {
03561         p = PolynomialStandard<Numeric>(opt_input_roots_scale, coefficients);
03562     }
03563 
03564     // run selected algorithm on polynomial
03565     if (opt_algo == ALGO_BEZCLIP)
03566     {
03567         BezierClip<Numeric>(opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
03568     }
03569     else if (opt_algo == ALGO_QUADCLIP)
03570     {
03571         KClip<Numeric>(p.degree(), 2, opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
03572     }
03573     else if (opt_algo == ALGO_CUBECLIP)
03574     {
03575         KClip<Numeric>(p.degree(), 3, opt_epsilon).findRoots(p, opt_interval_left, opt_interval_right);
03576     }
03577 
03578     if (opt_output_latexpreamble)
03579         ltx("\\end{document}\n");
03580 
03581     return 0;
03582 }
03583 
03584 /*****************************************************************************/
 All Classes Files Functions Variables Typedefs Friends Defines