1 /****************************************************************************** 2 * Author: Laurent Kneip * 3 * Contact: kneip.laurent@gmail.com * 4 * License: Copyright (c) 2013 Laurent Kneip, ANU. All rights reserved. * 5 * * 6 * Redistribution and use in source and binary forms, with or without * 7 * modification, are permitted provided that the following conditions * 8 * are met: * 9 * * Redistributions of source code must retain the above copyright * 10 * notice, this list of conditions and the following disclaimer. * 11 * * Redistributions in binary form must reproduce the above copyright * 12 * notice, this list of conditions and the following disclaimer in the * 13 * documentation and/or other materials provided with the distribution. * 14 * * Neither the name of ANU nor the names of its contributors may be * 15 * used to endorse or promote products derived from this software without * 16 * specific prior written permission. * 17 * * 18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"* 19 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * 20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * 21 * ARE DISCLAIMED. IN NO EVENT SHALL ANU OR THE CONTRIBUTORS BE LIABLE * 22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * 23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR * 24 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER * 25 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * 26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * 27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * 28 * SUCH DAMAGE. * 29 ******************************************************************************/ 30 31 /** 32 * \file Sturm.hpp 33 * \brief Class for evaluating Sturm chains on polynomials, and bracketing the 34 * real roots. 35 */ 36 37 #ifndef OPENGV_STURM_HPP_ 38 #define OPENGV_STURM_HPP_ 39 40 #include <memory> 41 #include <stdlib.h> 42 #include <vector> 43 #include <list> 44 #include <Eigen/Eigen> 45 #include <Eigen/src/Core/util/DisableStupidWarnings.h> 46 47 /** 48 * \brief The namespace of this library. 49 */ 50 namespace opengv 51 { 52 /** 53 * \brief The namespace of the math tools. 54 */ 55 namespace math 56 { 57 58 class Bracket 59 { 60 public: 61 typedef std::shared_ptr<Bracket> Ptr; 62 typedef std::shared_ptr<const Bracket> ConstPtr; 63 64 Bracket( double lowerBound, double upperBound ); 65 Bracket( double lowerBound, double upperBound, size_t changes, bool setUpperBoundChanges ); 66 virtual ~Bracket(); 67 68 bool dividable( double eps ) const; 69 void divide( std::list<Ptr> & brackets ) const; 70 double lowerBound() const; 71 double upperBound() const; 72 bool lowerBoundChangesComputed() const; 73 bool upperBoundChangesComputed() const; 74 void setLowerBoundChanges( size_t changes ); 75 void setUpperBoundChanges( size_t changes ); 76 size_t numberRoots() const; 77 78 private: 79 double _lowerBound; 80 double _upperBound; 81 bool _lowerBoundChangesComputed; 82 bool _upperBoundChangesComputed; 83 size_t _lowerBoundChanges; 84 size_t _upperBoundChanges; 85 }; 86 87 /** 88 * Sturm is initialized over polynomials of arbitrary order, and used to compute 89 * the real roots of the polynomial. 90 */ 91 class Sturm 92 { 93 94 public: 95 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 96 /** A pair of values bracketing a real root */ 97 typedef std::pair<double,double> bracket_t; 98 99 /** 100 * \brief Contructor. 101 * \param[in] p The polynomial coefficients (poly = p(0,0)*x^n + p(0,1)*x^(n-1) ...). 102 */ 103 Sturm( const Eigen::MatrixXd & p ); 104 /** 105 * \brief Contructor. 106 * \param[in] p The polynomial coefficients (poly = p[0]*x^n + p[1]*x^(n-1) ...). 107 */ 108 Sturm( const std::vector<double> & p ); 109 /** 110 * \brief Destructor. 111 */ 112 virtual ~Sturm(); 113 114 void findRoots2( std::vector<double> & roots, double eps_x = 0.001, double eps_val = 0.001 ); 115 /** 116 * \brief Finds the roots of the polynomial. 117 * \return An array with the real roots of the polynomial. 118 */ 119 std::vector<double> findRoots(); 120 /** 121 * \brief Finds brackets for the real roots of the polynomial. 122 * \return A list of brackets for the real roots of the polynomial. 123 */ 124 void bracketRoots( std::vector<double> & roots, double eps = -1.0 ); 125 /** 126 * \brief Evaluates the Sturm chain at a single bound. 127 * \param[in] bound The bound. 128 * \return The number of sign changes on the bound. 129 */ 130 size_t evaluateChain( double bound ); 131 /** 132 * \brief Evaluates the Sturm chain at a single bound. 133 * \param[in] bound The bound. 134 * \return The number of sign changes on the bound. 135 */ 136 size_t evaluateChain2( double bound ); 137 /** 138 * \brief Composes an initial bracket for all the roots of the polynomial. 139 * \return The maximum of the absolute values of the bracket-values (That's 140 * what the Lagrangian bound is able to find). 141 */ 142 double computeLagrangianBound(); 143 144 private: 145 /** 146 * \brief Internal function used for composing the Sturm chain 147 * \param[in] p1 First polynomial. 148 * \param[in] p2 Second polynomial. 149 * \param[out] r The negated remainder of the polynomial division p1/p2. 150 */ 151 void computeNegatedRemainder( 152 const Eigen::MatrixXd & p1, 153 const Eigen::MatrixXd & p2, 154 Eigen::MatrixXd & r ); 155 156 /** A matrix containing the coefficients of the Sturm-chain of the polynomial */ 157 Eigen::MatrixXd _C; 158 /** The dimension _C, which corresponds to (polynomial order+1) */ 159 size_t _dimension; 160 }; 161 162 } 163 } 164 165 #endif /* OPENGV_STURM_HPP_ */ 166