--- trunk/src/math/Polynomial.hpp 2005/03/01 20:10:14 385 +++ trunk/src/math/Polynomial.hpp 2009/11/25 20:02:06 1390 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,15 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ /** @@ -53,230 +53,627 @@ #include #include #include +#include +#include "config.h" +#include "math/Eigenvalue.hpp" -namespace oopse { +namespace OpenMD { + + template Real fastpow(Real x, int N) { + Real result(1); //or 1.0? -template ElemType pow(ElemType x, int N) { - ElemType result(1); - for (int i = 0; i < N; ++i) { - result *= x; + result *= x; } return result; -} + } -/** - * @class Polynomial Polynomial.hpp "math/Polynomial.hpp" - * A generic Polynomial class - */ -template -class Polynomial { + /** + * @class Polynomial Polynomial.hpp "math/Polynomial.hpp" + * A generic Polynomial class + */ + template + class Polynomial { - public: - - typedef int ExponentType; - typedef ElemType CoefficientType; - typedef std::map PolynomialPairMap; - typedef typename PolynomialPairMap::iterator iterator; - typedef typename PolynomialPairMap::const_iterator const_iterator; - /** - * Calculates the value of this Polynomial evaluated at the given x value. - * @return The value of this Polynomial evaluates at the given x value - * @param x the value of the independent variable for this Polynomial function - */ - ElemType evaluate(const ElemType& x) { - ElemType result = ElemType(); - ExponentType exponent; - CoefficientType coefficient; + public: + typedef Polynomial PolynomialType; + typedef int ExponentType; + typedef Real CoefficientType; + typedef std::map PolynomialPairMap; + typedef typename PolynomialPairMap::iterator iterator; + typedef typename PolynomialPairMap::const_iterator const_iterator; + + Polynomial() {} + Polynomial(Real v) {setCoefficient(0, v);} + /** + * Calculates the value of this Polynomial evaluated at the given x value. + * @return The value of this Polynomial evaluates at the given x value + * @param x the value of the independent variable for this + * Polynomial function + */ + Real evaluate(const Real& x) { + Real result = Real(); + ExponentType exponent; + CoefficientType coefficient; - for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) { - exponent = i->first; - coefficient = i->second; - result += pow(x, exponent) * coefficient; - } + for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) { + exponent = i->first; + coefficient = i->second; + result += fastpow(x, exponent) * coefficient; + } - return result; - } + return result; + } - /** - * Returns the first derivative of this polynomial. - * @return the first derivative of this polynomial - * @param x - */ - ElemType evaluateDerivative(const ElemType& x) { - ElemType result = ElemType(); - ExponentType exponent; - CoefficientType coefficient; + /** + * Returns the first derivative of this polynomial. + * @return the first derivative of this polynomial + * @param x + */ + Real evaluateDerivative(const Real& x) { + Real result = Real(); + ExponentType exponent; + CoefficientType coefficient; - for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) { - exponent = i->first; - coefficient = i->second; - result += pow(x, exponent - 1) * coefficient * exponent; - } + for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) { + exponent = i->first; + coefficient = i->second; + result += fastpow(x, exponent - 1) * coefficient * exponent; + } - return result; - } + return result; + } - /** - * Set the coefficent of the specified exponent, if the coefficient is already there, it - * will be overwritten. - * @param exponent exponent of a term in this Polynomial - * @param coefficient multiplier of a term in this Polynomial - */ - - void setCoefficient(int exponent, const ElemType& coefficient) { - polyPairMap_.insert(typename PolynomialPairMap::value_type(exponent, coefficient)); - } - /** - * Set the coefficent of the specified exponent. If the coefficient is already there, just add the - * new coefficient to the old one, otherwise, just call setCoefficent - * @param exponent exponent of a term in this Polynomial - * @param coefficient multiplier of a term in this Polynomial - */ - - void addCoefficient(int exponent, const ElemType& coefficient) { - iterator i = polyPairMap_.find(exponent); + /** + * Set the coefficent of the specified exponent, if the + * coefficient is already there, it will be overwritten. + * @param exponent exponent of a term in this Polynomial + * @param coefficient multiplier of a term in this Polynomial + */ + void setCoefficient(int exponent, const Real& coefficient) { + polyPairMap_[exponent] = coefficient; + } + + /** + * Set the coefficent of the specified exponent. If the + * coefficient is already there, just add the new coefficient to + * the old one, otherwise, just call setCoefficent + * @param exponent exponent of a term in this Polynomial + * @param coefficient multiplier of a term in this Polynomial + */ + void addCoefficient(int exponent, const Real& coefficient) { + iterator i = polyPairMap_.find(exponent); - if (i != end()) { - i->second += coefficient; - } else { - setCoefficient(exponent, coefficient); - } - } + if (i != end()) { + i->second += coefficient; + } else { + setCoefficient(exponent, coefficient); + } + } - - /** - * Returns the coefficient associated with the given power for this Polynomial. - * @return the coefficient associated with the given power for this Polynomial - * @exponent exponent of any term in this Polynomial - */ - ElemType getCoefficient(ExponentType exponent) { - iterator i = polyPairMap_.find(exponent); - - if (i != end()) { - return i->second; - } else { - return ElemType(0); - } - } + /** + * Returns the coefficient associated with the given power for + * this Polynomial. + * @return the coefficient associated with the given power for + * this Polynomial + * @exponent exponent of any term in this Polynomial + */ + Real getCoefficient(ExponentType exponent) { + iterator i = polyPairMap_.find(exponent); - iterator begin() { - return polyPairMap_.begin(); - } + if (i != end()) { + return i->second; + } else { + return Real(0); + } + } - const_iterator begin() const{ - return polyPairMap_.begin(); - } + iterator begin() { + return polyPairMap_.begin(); + } + + const_iterator begin() const{ + return polyPairMap_.begin(); + } - iterator end() { - return polyPairMap_.end(); - } + iterator end() { + return polyPairMap_.end(); + } - const_iterator end() const{ - return polyPairMap_.end(); + const_iterator end() const{ + return polyPairMap_.end(); + } + + iterator find(ExponentType exponent) { + return polyPairMap_.find(exponent); + } + + size_t size() { + return polyPairMap_.size(); + } + + int degree() { + int deg = 0; + for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) { + if (i->first > deg) + deg = i->first; + } + return deg; + } + + PolynomialType& operator = (const PolynomialType& p) { + + if (this != &p) // protect against invalid self-assignment + { + typename Polynomial::const_iterator i; + + polyPairMap_.clear(); // clear out the old map + + for (i = p.begin(); i != p.end(); ++i) { + this->setCoefficient(i->first, i->second); + } } + // by convention, always return *this + return *this; + } - iterator find(ExponentType exponent) { - return polyPairMap_.find(exponent); + PolynomialType& operator += (const PolynomialType& p) { + typename Polynomial::const_iterator i; + + for (i = p.begin(); i != p.end(); ++i) { + this->addCoefficient(i->first, i->second); + } + + return *this; + } + + PolynomialType& operator -= (const PolynomialType& p) { + typename Polynomial::const_iterator i; + for (i = p.begin(); i != p.end(); ++i) { + this->addCoefficient(i->first, -i->second); + } + return *this; + } + + PolynomialType& operator *= (const PolynomialType& p) { + typename Polynomial::const_iterator i; + typename Polynomial::const_iterator j; + Polynomial p2(*this); + + polyPairMap_.clear(); // clear out old map + for (i = p2.begin(); i !=p2.end(); ++i) { + for (j = p.begin(); j !=p.end(); ++j) { + this->addCoefficient( i->first + j->first, i->second * j->second); } + } + return *this; + } - size_t size() { - return polyPairMap_.size(); + //PolynomialType& operator *= (const Real v) + PolynomialType& operator *= (const Real v) { + typename Polynomial::const_iterator i; + //Polynomial result; + + for (i = this->begin(); i != this->end(); ++i) { + this->setCoefficient( i->first, i->second*v); + } + + return *this; + } + + PolynomialType& operator += (const Real v) { + this->addCoefficient( 0, v); + return *this; + } + + /** + * Returns the first derivative of this polynomial. + * @return the first derivative of this polynomial + */ + PolynomialType & getDerivative() { + Polynomial p(); + + typename Polynomial::const_iterator i; + ExponentType exponent; + CoefficientType coefficient; + + for (i = this->begin(); i != this->end(); ++i) { + exponent = i->first; + coefficient = i->second; + p.setCoefficient(exponent-1, coefficient * exponent); + } + + return p; + } + + // Creates the Companion matrix for a given polynomial + DynamicRectMatrix CreateCompanion() { + int rank = degree(); + DynamicRectMatrix mat(rank, rank); + Real majorCoeff = getCoefficient(rank); + for(int i = 0; i < rank; ++i) { + for(int j = 0; j < rank; ++j) { + if(i - j == 1) { + mat(i, j) = 1; + } else if(j == rank-1) { + mat(i, j) = -1 * getCoefficient(i) / majorCoeff; + } } - - private: - - PolynomialPairMap polyPairMap_; -}; + } + return mat; + } + + // Find the Roots of a given polynomial + std::vector > FindRoots() { + int rank = degree(); + DynamicRectMatrix companion = CreateCompanion(); + JAMA::Eigenvalue eig(companion); + DynamicVector reals, imags; + eig.getRealEigenvalues(reals); + eig.getImagEigenvalues(imags); + + std::vector > roots; + for (int i = 0; i < rank; i++) { + roots.push_back(complex(reals(i), imags(i))); + } + return roots; + } -/** - * Generates and returns the product of two given Polynomials. - * @return A Polynomial containing the product of the two given Polynomial parameters - */ -template -Polynomial operator *(const Polynomial& p1, const Polynomial& p2) { - typename Polynomial::const_iterator i; - typename Polynomial::const_iterator j; - Polynomial p; + std::vector FindRealRoots() { + + const Real fEpsilon = 1.0e-8; + std::vector roots; + roots.clear(); + + const int deg = degree(); + + switch (deg) { + case 1: { + Real fC1 = getCoefficient(1); + Real fC0 = getCoefficient(0); + roots.push_back( -fC0 / fC1); + return roots; + } + break; + case 2: { + Real fC2 = getCoefficient(2); + Real fC1 = getCoefficient(1); + Real fC0 = getCoefficient(0); + Real fDiscr = fC1*fC1 - 4.0*fC0*fC2; + if (abs(fDiscr) <= fEpsilon) { + fDiscr = (Real)0.0; + } + + if (fDiscr < (Real)0.0) { // complex roots only + return roots; + } + + Real fTmp = ((Real)0.5)/fC2; + + if (fDiscr > (Real)0.0) { // 2 real roots + fDiscr = sqrt(fDiscr); + roots.push_back(fTmp*(-fC1 - fDiscr)); + roots.push_back(fTmp*(-fC1 + fDiscr)); + } else { + roots.push_back(-fTmp * fC1); // 1 real root + } + } + return roots; + break; + + case 3: { + Real fC3 = getCoefficient(3); + Real fC2 = getCoefficient(2); + Real fC1 = getCoefficient(1); + Real fC0 = getCoefficient(0); + + // make polynomial monic, x^3+c2*x^2+c1*x+c0 + Real fInvC3 = ((Real)1.0)/fC3; + fC0 *= fInvC3; + fC1 *= fInvC3; + fC2 *= fInvC3; + + // convert to y^3+a*y+b = 0 by x = y-c2/3 + const Real fThird = (Real)1.0/(Real)3.0; + const Real fTwentySeventh = (Real)1.0/(Real)27.0; + Real fOffset = fThird*fC2; + Real fA = fC1 - fC2*fOffset; + Real fB = fC0+fC2*(((Real)2.0)*fC2*fC2-((Real)9.0)*fC1)*fTwentySeventh; + Real fHalfB = ((Real)0.5)*fB; + + Real fDiscr = fHalfB*fHalfB + fA*fA*fA*fTwentySeventh; + if (fabs(fDiscr) <= fEpsilon) { + fDiscr = (Real)0.0; + } + + if (fDiscr > (Real)0.0) { // 1 real, 2 complex roots + + fDiscr = sqrt(fDiscr); + Real fTemp = -fHalfB + fDiscr; + Real root; + if (fTemp >= (Real)0.0) { + root = pow(fTemp,fThird); + } else { + root = -pow(-fTemp,fThird); + } + fTemp = -fHalfB - fDiscr; + if ( fTemp >= (Real)0.0 ) { + root += pow(fTemp,fThird); + } else { + root -= pow(-fTemp,fThird); + } + root -= fOffset; + + roots.push_back(root); + } else if (fDiscr < (Real)0.0) { + const Real fSqrt3 = sqrt((Real)3.0); + Real fDist = sqrt(-fThird*fA); + Real fAngle = fThird*atan2(sqrt(-fDiscr), -fHalfB); + Real fCos = cos(fAngle); + Real fSin = sin(fAngle); + roots.push_back(((Real)2.0)*fDist*fCos-fOffset); + roots.push_back(-fDist*(fCos+fSqrt3*fSin)-fOffset); + roots.push_back(-fDist*(fCos-fSqrt3*fSin)-fOffset); + } else { + Real fTemp; + if (fHalfB >= (Real)0.0) { + fTemp = -pow(fHalfB,fThird); + } else { + fTemp = pow(-fHalfB,fThird); + } + roots.push_back(((Real)2.0)*fTemp-fOffset); + roots.push_back(-fTemp-fOffset); + roots.push_back(-fTemp-fOffset); + } + } + return roots; + break; + case 4: { + Real fC4 = getCoefficient(4); + Real fC3 = getCoefficient(3); + Real fC2 = getCoefficient(2); + Real fC1 = getCoefficient(1); + Real fC0 = getCoefficient(0); + + // make polynomial monic, x^4+c3*x^3+c2*x^2+c1*x+c0 + Real fInvC4 = ((Real)1.0)/fC4; + fC0 *= fInvC4; + fC1 *= fInvC4; + fC2 *= fInvC4; + fC3 *= fInvC4; + + // reduction to resolvent cubic polynomial y^3+r2*y^2+r1*y+r0 = 0 + Real fR0 = -fC3*fC3*fC0 + ((Real)4.0)*fC2*fC0 - fC1*fC1; + Real fR1 = fC3*fC1 - ((Real)4.0)*fC0; + Real fR2 = -fC2; + Polynomial tempCubic; + tempCubic.setCoefficient(0, fR0); + tempCubic.setCoefficient(1, fR1); + tempCubic.setCoefficient(2, fR2); + tempCubic.setCoefficient(3, 1.0); + std::vector cubeRoots = tempCubic.FindRealRoots(); // always + // produces + // at + // least + // one + // root + Real fY = cubeRoots[0]; + + Real fDiscr = ((Real)0.25)*fC3*fC3 - fC2 + fY; + if (fabs(fDiscr) <= fEpsilon) { + fDiscr = (Real)0.0; + } + + if (fDiscr > (Real)0.0) { + Real fR = sqrt(fDiscr); + Real fT1 = ((Real)0.75)*fC3*fC3 - fR*fR - ((Real)2.0)*fC2; + Real fT2 = (((Real)4.0)*fC3*fC2 - ((Real)8.0)*fC1 - fC3*fC3*fC3) / + (((Real)4.0)*fR); + + Real fTplus = fT1+fT2; + Real fTminus = fT1-fT2; + if (fabs(fTplus) <= fEpsilon) { + fTplus = (Real)0.0; + } + if (fabs(fTminus) <= fEpsilon) { + fTminus = (Real)0.0; + } + + if (fTplus >= (Real)0.0) { + Real fD = sqrt(fTplus); + roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR+fD)); + roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR-fD)); + } + if (fTminus >= (Real)0.0) { + Real fE = sqrt(fTminus); + roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fE-fR)); + roots.push_back(-((Real)0.25)*fC3-((Real)0.5)*(fE+fR)); + } + } else if (fDiscr < (Real)0.0) { + //roots.clear(); + } else { + Real fT2 = fY*fY-((Real)4.0)*fC0; + if (fT2 >= -fEpsilon) { + if (fT2 < (Real)0.0) { // round to zero + fT2 = (Real)0.0; + } + fT2 = ((Real)2.0)*sqrt(fT2); + Real fT1 = ((Real)0.75)*fC3*fC3 - ((Real)2.0)*fC2; + if (fT1+fT2 >= fEpsilon) { + Real fD = sqrt(fT1+fT2); + roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fD); + roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fD); + } + if (fT1-fT2 >= fEpsilon) { + Real fE = sqrt(fT1-fT2); + roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fE); + roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fE); + } + } + } + } + return roots; + break; + default: { + DynamicRectMatrix companion = CreateCompanion(); + JAMA::Eigenvalue eig(companion); + DynamicVector reals, imags; + eig.getRealEigenvalues(reals); + eig.getImagEigenvalues(imags); + + for (int i = 0; i < deg; i++) { + if (fabs(imags(i)) < fEpsilon) + roots.push_back(reals(i)); + } + } + return roots; + break; + } + + return roots; // should be empty if you got here + } + + private: + + PolynomialPairMap polyPairMap_; + }; + + + /** + * Generates and returns the product of two given Polynomials. + * @return A Polynomial containing the product of the two given Polynomial parameters + */ + template + Polynomial operator *(const Polynomial& p1, const Polynomial& p2) { + typename Polynomial::const_iterator i; + typename Polynomial::const_iterator j; + Polynomial p; for (i = p1.begin(); i !=p1.end(); ++i) { - for (j = p2.begin(); j !=p2.end(); ++j) { - p.addCoefficient( i->first + j->first, i->second * j->second); - } + for (j = p2.begin(); j !=p2.end(); ++j) { + p.addCoefficient( i->first + j->first, i->second * j->second); + } } return p; -} + } -/** - * Generates and returns the sum of two given Polynomials. - * @param p1 the first polynomial - * @param p2 the second polynomial - */ -template -Polynomial operator +(const Polynomial& p1, const Polynomial& p2) { - Polynomial p(p1); + template + Polynomial operator *(const Polynomial& p, const Real v) { + typename Polynomial::const_iterator i; + Polynomial result; + + for (i = p.begin(); i !=p.end(); ++i) { + result.setCoefficient( i->first , i->second * v); + } - typename Polynomial::const_iterator i; + return result; + } + template + Polynomial operator *( const Real v, const Polynomial& p) { + typename Polynomial::const_iterator i; + Polynomial result; + + for (i = p.begin(); i !=p.end(); ++i) { + result.setCoefficient( i->first , i->second * v); + } + + return result; + } + + /** + * Generates and returns the sum of two given Polynomials. + * @param p1 the first polynomial + * @param p2 the second polynomial + */ + template + Polynomial operator +(const Polynomial& p1, const Polynomial& p2) { + Polynomial p(p1); + + typename Polynomial::const_iterator i; + for (i = p2.begin(); i != p2.end(); ++i) { - p.addCoefficient(i->first, i->second); + p.addCoefficient(i->first, i->second); } return p; -} + } -/** - * Generates and returns the difference of two given Polynomials. - * @return - * @param p1 the first polynomial - * @param p2 the second polynomial - */ -template -Polynomial operator -(const Polynomial& p1, const Polynomial& p2) { - Polynomial p(p1); + /** + * Generates and returns the difference of two given Polynomials. + * @return + * @param p1 the first polynomial + * @param p2 the second polynomial + */ + template + Polynomial operator -(const Polynomial& p1, const Polynomial& p2) { + Polynomial p(p1); - typename Polynomial::const_iterator i; + typename Polynomial::const_iterator i; for (i = p2.begin(); i != p2.end(); ++i) { - p.addCoefficient(i->first, -i->second); + p.addCoefficient(i->first, -i->second); } return p; -} + } -/** - * Tests if two polynomial have the same exponents - * @return true if these all of the exponents in these Polynomial are identical - * @param p1 the first polynomial - * @param p2 the second polynomial - * @note this function does not compare the coefficient - */ -template -bool equal(const Polynomial& p1, const Polynomial& p2) { + /** + * Returns the first derivative of this polynomial. + * @return the first derivative of this polynomial + */ + template + Polynomial getDerivative(const Polynomial& p1) { + Polynomial p(); + + typename Polynomial::const_iterator i; + int exponent; + Real coefficient; + + for (i = p1.begin(); i != p1.end(); ++i) { + exponent = i->first; + coefficient = i->second; + p.setCoefficient(exponent-1, coefficient * exponent); + } + + return p; + } - typename Polynomial::const_iterator i; - typename Polynomial::const_iterator j; + /** + * Tests if two polynomial have the same exponents + * @return true if all of the exponents in these Polynomial are identical + * @param p1 the first polynomial + * @param p2 the second polynomial + * @note this function does not compare the coefficient + */ + template + bool equal(const Polynomial& p1, const Polynomial& p2) { + typename Polynomial::const_iterator i; + typename Polynomial::const_iterator j; + if (p1.size() != p2.size() ) { - return false; + return false; } for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) { - if (i->first != j->first) { - return false; - } + if (i->first != j->first) { + return false; + } } return true; -} + } -typedef Polynomial DoublePolynomial; -} //end namespace oopse + + typedef Polynomial DoublePolynomial; + +} //end namespace OpenMD #endif //MATH_POLYNOMIAL_HPP