src/math/Polynomial.hpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 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.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * 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, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
/**
 * @file Polynomial.hpp
 * @author    teng lin
 * @date  11/16/2004
 * @version 1.0
 */ 

#ifndef MATH_POLYNOMIAL_HPP
#define MATH_POLYNOMIAL_HPP

#include <iostream>
#include <list>
#include <map>
#include <utility>
#include <complex>
#include "config.h"
#include "math/Eigenvalue.hpp"

namespace OpenMD {
  
  template<typename Real> Real fastpow(Real x, int N) {
    Real result(1); //or 1.0?

    for (int i = 0; i < N; ++i) {
      result *= x;
    }

    return result;
  }

  /**
   * @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
   * A generic Polynomial class
   */
  template<typename Real>
  class Polynomial {

  public:
    typedef Polynomial<Real> PolynomialType;    
    typedef int ExponentType;
    typedef Real CoefficientType;
    typedef std::map<ExponentType, CoefficientType> 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  += fastpow(x, exponent) * coefficient;
      }

      return result;
    }

    /**
     * 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  += fastpow(x, exponent - 1) * coefficient * exponent;
      }

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

    /**
     * Returns the coefficient associated with the given power for
     * this Polynomial.
     * @return the coefficient associated with the given power for
     * this Polynomial
     * @param exponent exponent of any term in this Polynomial
     */
    Real getCoefficient(ExponentType exponent) {
      iterator i = polyPairMap_.find(exponent);

      if (i != end()) {
        return i->second;
      } else {
        return Real(0);
      }
    }

    iterator begin() {
      return polyPairMap_.begin();
    }

    const_iterator begin() const{
      return polyPairMap_.begin();
    }
        
    iterator end() {
      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<Real>::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; 
    }

    PolynomialType& operator += (const PolynomialType& p) {
      typename Polynomial<Real>::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<Real>::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<Real>::const_iterator i;
      typename Polynomial<Real>::const_iterator j;
      Polynomial<Real> 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;
    }

    //PolynomialType& operator *= (const Real v)
    PolynomialType& operator *= (const Real v) {
      typename Polynomial<Real>::const_iterator i;
      //Polynomial<Real> 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<Real>* p = new Polynomial<Real>();
      
      typename Polynomial<Real>::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<Real> CreateCompanion() {
      int rank = degree();
      DynamicRectMatrix<Real> 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;
          }
        }
      }
      return mat;
    }
    
    // Find the Roots of a given polynomial
    std::vector<complex<Real> > FindRoots() {
      int rank = degree();
      DynamicRectMatrix<Real> companion = CreateCompanion();
      JAMA::Eigenvalue<Real> eig(companion);
      DynamicVector<Real> reals, imags;
      eig.getRealEigenvalues(reals);
      eig.getImagEigenvalues(imags);
      
      std::vector<complex<Real> > roots;
      for (int i = 0; i < rank; i++) {
        roots.push_back(complex<Real>(reals(i), imags(i)));
      }

      return roots;
    }

    std::vector<Real> FindRealRoots() {
      
      const Real fEpsilon = 1.0e-8;
      std::vector<Real> 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;
      }
      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;      
      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;
      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<Real> tempCubic;
        tempCubic.setCoefficient(0, fR0);
        tempCubic.setCoefficient(1, fR1);
        tempCubic.setCoefficient(2, fR2);
        tempCubic.setCoefficient(3, 1.0);
        std::vector<Real> 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;
      default: {
        DynamicRectMatrix<Real> companion = CreateCompanion();
        JAMA::Eigenvalue<Real> eig(companion);
        DynamicVector<Real> 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;
      }
    }
    
  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<typename Real>
  Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
    typename Polynomial<Real>::const_iterator i;
    typename Polynomial<Real>::const_iterator j;
    Polynomial<Real> 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);
      }
    }

    return p;
  }

  template<typename Real>
  Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
    typename Polynomial<Real>::const_iterator i;
    Polynomial<Real> result;
    
    for (i = p.begin(); i !=p.end(); ++i) {
        result.setCoefficient( i->first , i->second * v);
    }

    return result;
  }

  template<typename Real>
  Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
    typename Polynomial<Real>::const_iterator i;
    Polynomial<Real> 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<typename Real>
  Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
    Polynomial<Real> p(p1);

    typename Polynomial<Real>::const_iterator i;

    for (i =  p2.begin(); i  != p2.end(); ++i) {
      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<typename Real>
  Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
    Polynomial<Real> p(p1);

    typename Polynomial<Real>::const_iterator i;

    for (i =  p2.begin(); i  != p2.end(); ++i) {
      p.addCoefficient(i->first, -i->second);
    }

    return p;

  }

  /**
   * Returns the first derivative of this polynomial.
   * @return the first derivative of this polynomial
   */
  template<typename Real>
  Polynomial<Real> * getDerivative(const Polynomial<Real>& p1) {
    Polynomial<Real> * p = new Polynomial<Real>();
    
    typename Polynomial<Real>::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;
  }

  /**
   * 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<typename Real>
  bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {

    typename Polynomial<Real>::const_iterator i;
    typename Polynomial<Real>::const_iterator j;

    if (p1.size() != p2.size() ) {
      return false;
    }
    
    for (i =  p1.begin(), j = p2.begin(); i  != p1.end() && j != p2.end(); ++i, ++j) {
      if (i->first != j->first) {
        return false;
      }
    }

    return true;
  }


  typedef Polynomial<RealType> DoublePolynomial;

} //end namespace OpenMD
#endif //MATH_POLYNOMIAL_HPP
Revision:	1874
Committed:	Wed May 15 15:09:35 2013 UTC (11 years, 11 months ago) by gezelter
File size:	20211 byte(s)
Log Message:	Fixed a bunch of cppcheck warnings.
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	gezelter	1850	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43			/**
44			* @file Polynomial.hpp
45			* @author teng lin
46			* @date 11/16/2004
47			* @version 1.0
48			*/
49
50			#ifndef MATH_POLYNOMIAL_HPP
51			#define MATH_POLYNOMIAL_HPP
52
53			#include <iostream>
54			#include <list>
55			#include <map>
56			#include <utility>
57	gezelter	1358	#include <complex>
58	tim	963	#include "config.h"
59	gezelter	1358	#include "math/Eigenvalue.hpp"
60
61	gezelter	1390	namespace OpenMD {
62	gezelter	1358
63			template<typename Real> Real fastpow(Real x, int N) {
64			Real result(1); //or 1.0?
65	gezelter	246
66			for (int i = 0; i < N; ++i) {
67	gezelter	507	result *= x;
68	gezelter	246	}
69
70			return result;
71	gezelter	507	}
72	gezelter	246
73	gezelter	507	/**
74			* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
75			* A generic Polynomial class
76			*/
77	gezelter	1358	template<typename Real>
78	gezelter	507	class Polynomial {
79	gezelter	246
80	gezelter	507	public:
81	gezelter	1358	typedef Polynomial<Real> PolynomialType;
82	gezelter	507	typedef int ExponentType;
83	gezelter	1358	typedef Real CoefficientType;
84	gezelter	507	typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
85			typedef typename PolynomialPairMap::iterator iterator;
86			typedef typename PolynomialPairMap::const_iterator const_iterator;
87	tim	749
88			Polynomial() {}
89	gezelter	1358	Polynomial(Real v) {setCoefficient(0, v);}
90	gezelter	507	/**
91			* Calculates the value of this Polynomial evaluated at the given x value.
92	gezelter	1358	* @return The value of this Polynomial evaluates at the given x value
93			* @param x the value of the independent variable for this
94			* Polynomial function
95	gezelter	507	*/
96	gezelter	1358	Real evaluate(const Real& x) {
97			Real result = Real();
98	gezelter	507	ExponentType exponent;
99			CoefficientType coefficient;
100	gezelter	246
101	gezelter	507	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
102			exponent = i->first;
103			coefficient = i->second;
104	gezelter	1358	result += fastpow(x, exponent) * coefficient;
105	gezelter	507	}
106	gezelter	246
107	gezelter	507	return result;
108			}
109	gezelter	246
110	gezelter	507	/**
111			* Returns the first derivative of this polynomial.
112			* @return the first derivative of this polynomial
113			* @param x
114			*/
115	gezelter	1358	Real evaluateDerivative(const Real& x) {
116			Real result = Real();
117	gezelter	507	ExponentType exponent;
118			CoefficientType coefficient;
119	gezelter	246
120	gezelter	507	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
121			exponent = i->first;
122			coefficient = i->second;
123	gezelter	1358	result += fastpow(x, exponent - 1) * coefficient * exponent;
124	gezelter	507	}
125	gezelter	246
126	gezelter	507	return result;
127			}
128	gezelter	246
129	gezelter	1358
130	gezelter	507	/**
131	gezelter	1358	* Set the coefficent of the specified exponent, if the
132			* coefficient is already there, it will be overwritten.
133	gezelter	507	* @param exponent exponent of a term in this Polynomial
134			* @param coefficient multiplier of a term in this Polynomial
135	gezelter	1358	*/
136			void setCoefficient(int exponent, const Real& coefficient) {
137	cli2	1290	polyPairMap_[exponent] = coefficient;
138	gezelter	507	}
139	gezelter	1358
140	gezelter	507	/**
141	gezelter	1358	* Set the coefficent of the specified exponent. If the
142			* coefficient is already there, just add the new coefficient to
143			* the old one, otherwise, just call setCoefficent
144	gezelter	507	* @param exponent exponent of a term in this Polynomial
145			* @param coefficient multiplier of a term in this Polynomial
146	gezelter	1358	*/
147			void addCoefficient(int exponent, const Real& coefficient) {
148	gezelter	507	iterator i = polyPairMap_.find(exponent);
149	gezelter	246
150	gezelter	507	if (i != end()) {
151			i->second += coefficient;
152			} else {
153			setCoefficient(exponent, coefficient);
154			}
155			}
156	gezelter	246
157	gezelter	507	/**
158	gezelter	1358	* Returns the coefficient associated with the given power for
159			* this Polynomial.
160			* @return the coefficient associated with the given power for
161			* this Polynomial
162	gezelter	1808	* @param exponent exponent of any term in this Polynomial
163	gezelter	507	*/
164	gezelter	1358	Real getCoefficient(ExponentType exponent) {
165	gezelter	507	iterator i = polyPairMap_.find(exponent);
166	gezelter	246
167	gezelter	507	if (i != end()) {
168			return i->second;
169			} else {
170	gezelter	1358	return Real(0);
171	gezelter	507	}
172			}
173	gezelter	246
174	gezelter	507	iterator begin() {
175			return polyPairMap_.begin();
176			}
177	gezelter	246
178	gezelter	507	const_iterator begin() const{
179			return polyPairMap_.begin();
180			}
181	gezelter	246
182	gezelter	507	iterator end() {
183			return polyPairMap_.end();
184			}
185	gezelter	246
186	gezelter	507	const_iterator end() const{
187			return polyPairMap_.end();
188			}
189	gezelter	246
190	gezelter	507	iterator find(ExponentType exponent) {
191			return polyPairMap_.find(exponent);
192			}
193	gezelter	246
194	gezelter	507	size_t size() {
195			return polyPairMap_.size();
196			}
197	tim	749
198	gezelter	1358	int degree() {
199			int deg = 0;
200			for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
201			if (i->first > deg)
202			deg = i->first;
203			}
204			return deg;
205			}
206
207	cpuglis	1230	PolynomialType& operator = (const PolynomialType& p) {
208
209			if (this != &p) // protect against invalid self-assignment
210	gezelter	1358	{
211			typename Polynomial<Real>::const_iterator i;
212	cpuglis	1230
213	gezelter	1358	polyPairMap_.clear(); // clear out the old map
214	cpuglis	1230
215	gezelter	1358	for (i = p.begin(); i != p.end(); ++i) {
216			this->setCoefficient(i->first, i->second);
217			}
218	cpuglis	1230	}
219			// by convention, always return *this
220			return *this;
221			}
222
223	tim	749	PolynomialType& operator += (const PolynomialType& p) {
224	gezelter	1358	typename Polynomial<Real>::const_iterator i;
225	tim	749
226	gezelter	1358	for (i = p.begin(); i != p.end(); ++i) {
227			this->addCoefficient(i->first, i->second);
228			}
229	tim	749
230	gezelter	1358	return *this;
231	tim	749	}
232
233			PolynomialType& operator -= (const PolynomialType& p) {
234	gezelter	1358	typename Polynomial<Real>::const_iterator i;
235			for (i = p.begin(); i != p.end(); ++i) {
236			this->addCoefficient(i->first, -i->second);
237			}
238			return *this;
239	tim	749	}
240	gezelter	1358
241	tim	749	PolynomialType& operator *= (const PolynomialType& p) {
242	gezelter	1358	typename Polynomial<Real>::const_iterator i;
243			typename Polynomial<Real>::const_iterator j;
244			Polynomial<Real> p2(*this);
245
246			polyPairMap_.clear(); // clear out old map
247			for (i = p2.begin(); i !=p2.end(); ++i) {
248			for (j = p.begin(); j !=p.end(); ++j) {
249			this->addCoefficient( i->first + j->first, i->second * j->second);
250			}
251	tim	749	}
252	gezelter	1358	return *this;
253	tim	749	}
254	gezelter	1358
255			//PolynomialType& operator *= (const Real v)
256			PolynomialType& operator *= (const Real v) {
257			typename Polynomial<Real>::const_iterator i;
258			//Polynomial<Real> result;
259
260			for (i = this->begin(); i != this->end(); ++i) {
261			this->setCoefficient( i->first, i->second*v);
262			}
263
264			return *this;
265	cli2	1290	}
266	tim	749
267	gezelter	1358	PolynomialType& operator += (const Real v) {
268			this->addCoefficient( 0, v);
269			return *this;
270			}
271	cli2	1290
272	gezelter	1358	/**
273			* Returns the first derivative of this polynomial.
274			* @return the first derivative of this polynomial
275			*/
276	gezelter	1874	PolynomialType* getDerivative() {
277			Polynomial<Real>* p = new Polynomial<Real>();
278	gezelter	1358
279			typename Polynomial<Real>::const_iterator i;
280			ExponentType exponent;
281			CoefficientType coefficient;
282
283			for (i = this->begin(); i != this->end(); ++i) {
284			exponent = i->first;
285			coefficient = i->second;
286	gezelter	1874	p->setCoefficient(exponent-1, coefficient * exponent);
287	gezelter	1358	}
288
289			return p;
290	cli2	1290	}
291
292	gezelter	1358	// Creates the Companion matrix for a given polynomial
293			DynamicRectMatrix<Real> CreateCompanion() {
294			int rank = degree();
295			DynamicRectMatrix<Real> mat(rank, rank);
296			Real majorCoeff = getCoefficient(rank);
297			for(int i = 0; i < rank; ++i) {
298			for(int j = 0; j < rank; ++j) {
299			if(i - j == 1) {
300			mat(i, j) = 1;
301			} else if(j == rank-1) {
302			mat(i, j) = -1 * getCoefficient(i) / majorCoeff;
303			}
304			}
305			}
306			return mat;
307	tim	749	}
308	gezelter	1358
309			// Find the Roots of a given polynomial
310			std::vector<complex<Real> > FindRoots() {
311			int rank = degree();
312			DynamicRectMatrix<Real> companion = CreateCompanion();
313			JAMA::Eigenvalue<Real> eig(companion);
314			DynamicVector<Real> reals, imags;
315			eig.getRealEigenvalues(reals);
316			eig.getImagEigenvalues(imags);
317
318			std::vector<complex<Real> > roots;
319			for (int i = 0; i < rank; i++) {
320	cli2	1362	roots.push_back(complex<Real>(reals(i), imags(i)));
321	gezelter	1358	}
322	tim	749
323	gezelter	1358	return roots;
324	cli2	1290	}
325	gezelter	1358
326			std::vector<Real> FindRealRoots() {
327
328			const Real fEpsilon = 1.0e-8;
329			std::vector<Real> roots;
330			roots.clear();
331
332			const int deg = degree();
333
334			switch (deg) {
335			case 1: {
336			Real fC1 = getCoefficient(1);
337			Real fC0 = getCoefficient(0);
338			roots.push_back( -fC0 / fC1);
339			return roots;
340			}
341			case 2: {
342			Real fC2 = getCoefficient(2);
343			Real fC1 = getCoefficient(1);
344			Real fC0 = getCoefficient(0);
345			Real fDiscr = fC1fC1 - 4.0fC0*fC2;
346			if (abs(fDiscr) <= fEpsilon) {
347			fDiscr = (Real)0.0;
348			}
349	gezelter	1826
350	gezelter	1358	if (fDiscr < (Real)0.0) { // complex roots only
351			return roots;
352			}
353	gezelter	1826
354	gezelter	1358	Real fTmp = ((Real)0.5)/fC2;
355	gezelter	1826
356	gezelter	1358	if (fDiscr > (Real)0.0) { // 2 real roots
357			fDiscr = sqrt(fDiscr);
358			roots.push_back(fTmp*(-fC1 - fDiscr));
359			roots.push_back(fTmp*(-fC1 + fDiscr));
360			} else {
361			roots.push_back(-fTmp * fC1); // 1 real root
362			}
363			}
364	gezelter	1826	return roots;
365	gezelter	1358	case 3: {
366			Real fC3 = getCoefficient(3);
367			Real fC2 = getCoefficient(2);
368			Real fC1 = getCoefficient(1);
369			Real fC0 = getCoefficient(0);
370
371			// make polynomial monic, x^3+c2x^2+c1x+c0
372			Real fInvC3 = ((Real)1.0)/fC3;
373			fC0 *= fInvC3;
374			fC1 *= fInvC3;
375			fC2 *= fInvC3;
376
377			// convert to y^3+a*y+b = 0 by x = y-c2/3
378			const Real fThird = (Real)1.0/(Real)3.0;
379			const Real fTwentySeventh = (Real)1.0/(Real)27.0;
380			Real fOffset = fThird*fC2;
381			Real fA = fC1 - fC2*fOffset;
382			Real fB = fC0+fC2(((Real)2.0)fC2fC2-((Real)9.0)fC1)*fTwentySeventh;
383			Real fHalfB = ((Real)0.5)*fB;
384
385			Real fDiscr = fHalfBfHalfB + fAfAfAfTwentySeventh;
386			if (fabs(fDiscr) <= fEpsilon) {
387			fDiscr = (Real)0.0;
388			}
389
390			if (fDiscr > (Real)0.0) { // 1 real, 2 complex roots
391
392			fDiscr = sqrt(fDiscr);
393			Real fTemp = -fHalfB + fDiscr;
394			Real root;
395			if (fTemp >= (Real)0.0) {
396			root = pow(fTemp,fThird);
397			} else {
398			root = -pow(-fTemp,fThird);
399			}
400			fTemp = -fHalfB - fDiscr;
401			if ( fTemp >= (Real)0.0 ) {
402			root += pow(fTemp,fThird);
403			} else {
404			root -= pow(-fTemp,fThird);
405			}
406			root -= fOffset;
407
408			roots.push_back(root);
409			} else if (fDiscr < (Real)0.0) {
410			const Real fSqrt3 = sqrt((Real)3.0);
411			Real fDist = sqrt(-fThird*fA);
412			Real fAngle = fThird*atan2(sqrt(-fDiscr), -fHalfB);
413			Real fCos = cos(fAngle);
414			Real fSin = sin(fAngle);
415			roots.push_back(((Real)2.0)fDistfCos-fOffset);
416			roots.push_back(-fDist(fCos+fSqrt3fSin)-fOffset);
417			roots.push_back(-fDist(fCos-fSqrt3fSin)-fOffset);
418			} else {
419			Real fTemp;
420			if (fHalfB >= (Real)0.0) {
421			fTemp = -pow(fHalfB,fThird);
422			} else {
423			fTemp = pow(-fHalfB,fThird);
424			}
425			roots.push_back(((Real)2.0)*fTemp-fOffset);
426			roots.push_back(-fTemp-fOffset);
427			roots.push_back(-fTemp-fOffset);
428			}
429			}
430			return roots;
431			case 4: {
432			Real fC4 = getCoefficient(4);
433			Real fC3 = getCoefficient(3);
434			Real fC2 = getCoefficient(2);
435			Real fC1 = getCoefficient(1);
436			Real fC0 = getCoefficient(0);
437
438			// make polynomial monic, x^4+c3x^3+c2x^2+c1*x+c0
439			Real fInvC4 = ((Real)1.0)/fC4;
440			fC0 *= fInvC4;
441			fC1 *= fInvC4;
442			fC2 *= fInvC4;
443			fC3 *= fInvC4;
444
445			// reduction to resolvent cubic polynomial y^3+r2y^2+r1y+r0 = 0
446			Real fR0 = -fC3fC3fC0 + ((Real)4.0)fC2fC0 - fC1*fC1;
447			Real fR1 = fC3fC1 - ((Real)4.0)fC0;
448			Real fR2 = -fC2;
449			Polynomial<Real> tempCubic;
450			tempCubic.setCoefficient(0, fR0);
451			tempCubic.setCoefficient(1, fR1);
452			tempCubic.setCoefficient(2, fR2);
453			tempCubic.setCoefficient(3, 1.0);
454			std::vector<Real> cubeRoots = tempCubic.FindRealRoots(); // always
455			// produces
456			// at
457			// least
458			// one
459			// root
460			Real fY = cubeRoots[0];
461
462			Real fDiscr = ((Real)0.25)fC3fC3 - fC2 + fY;
463			if (fabs(fDiscr) <= fEpsilon) {
464			fDiscr = (Real)0.0;
465			}
466	tim	749
467	gezelter	1358	if (fDiscr > (Real)0.0) {
468			Real fR = sqrt(fDiscr);
469			Real fT1 = ((Real)0.75)fC3fC3 - fRfR - ((Real)2.0)fC2;
470			Real fT2 = (((Real)4.0)fC3fC2 - ((Real)8.0)fC1 - fC3fC3*fC3) /
471			(((Real)4.0)*fR);
472
473			Real fTplus = fT1+fT2;
474			Real fTminus = fT1-fT2;
475			if (fabs(fTplus) <= fEpsilon) {
476			fTplus = (Real)0.0;
477			}
478			if (fabs(fTminus) <= fEpsilon) {
479			fTminus = (Real)0.0;
480			}
481
482			if (fTplus >= (Real)0.0) {
483			Real fD = sqrt(fTplus);
484			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fR+fD));
485			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fR-fD));
486			}
487			if (fTminus >= (Real)0.0) {
488			Real fE = sqrt(fTminus);
489			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fE-fR));
490			roots.push_back(-((Real)0.25)fC3-((Real)0.5)(fE+fR));
491			}
492			} else if (fDiscr < (Real)0.0) {
493			//roots.clear();
494			} else {
495			Real fT2 = fYfY-((Real)4.0)fC0;
496			if (fT2 >= -fEpsilon) {
497			if (fT2 < (Real)0.0) { // round to zero
498			fT2 = (Real)0.0;
499			}
500			fT2 = ((Real)2.0)*sqrt(fT2);
501			Real fT1 = ((Real)0.75)fC3fC3 - ((Real)2.0)*fC2;
502			if (fT1+fT2 >= fEpsilon) {
503			Real fD = sqrt(fT1+fT2);
504			roots.push_back( -((Real)0.25)fC3+((Real)0.5)fD);
505			roots.push_back( -((Real)0.25)fC3-((Real)0.5)fD);
506			}
507			if (fT1-fT2 >= fEpsilon) {
508			Real fE = sqrt(fT1-fT2);
509			roots.push_back( -((Real)0.25)fC3+((Real)0.5)fE);
510			roots.push_back( -((Real)0.25)fC3-((Real)0.5)fE);
511			}
512			}
513			}
514			}
515			return roots;
516			default: {
517			DynamicRectMatrix<Real> companion = CreateCompanion();
518			JAMA::Eigenvalue<Real> eig(companion);
519			DynamicVector<Real> reals, imags;
520			eig.getRealEigenvalues(reals);
521			eig.getImagEigenvalues(imags);
522
523			for (int i = 0; i < deg; i++) {
524			if (fabs(imags(i)) < fEpsilon)
525			roots.push_back(reals(i));
526			}
527			}
528			return roots;
529			}
530			}
531	gezelter	1826
532	gezelter	507	private:
533	gezelter	246
534	gezelter	507	PolynomialPairMap polyPairMap_;
535			};
536	gezelter	246
537	gezelter	1358
538	gezelter	507	/**
539			* Generates and returns the product of two given Polynomials.
540			* @return A Polynomial containing the product of the two given Polynomial parameters
541			*/
542	gezelter	1358	template<typename Real>
543			Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
544			typename Polynomial<Real>::const_iterator i;
545			typename Polynomial<Real>::const_iterator j;
546			Polynomial<Real> p;
547	gezelter	246
548			for (i = p1.begin(); i !=p1.end(); ++i) {
549	gezelter	507	for (j = p2.begin(); j !=p2.end(); ++j) {
550			p.addCoefficient( i->first + j->first, i->second * j->second);
551			}
552	gezelter	246	}
553
554			return p;
555	gezelter	507	}
556	gezelter	246
557	gezelter	1358	template<typename Real>
558			Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
559			typename Polynomial<Real>::const_iterator i;
560			Polynomial<Real> result;
561	tim	876
562			for (i = p.begin(); i !=p.end(); ++i) {
563	cli2	1290	result.setCoefficient( i->first , i->second * v);
564	tim	876	}
565
566			return result;
567			}
568
569	gezelter	1358	template<typename Real>
570			Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
571			typename Polynomial<Real>::const_iterator i;
572			Polynomial<Real> result;
573	tim	876
574			for (i = p.begin(); i !=p.end(); ++i) {
575	cli2	1290	result.setCoefficient( i->first , i->second * v);
576	tim	876	}
577
578			return result;
579			}
580
581	gezelter	507	/**
582			* Generates and returns the sum of two given Polynomials.
583			* @param p1 the first polynomial
584			* @param p2 the second polynomial
585			*/
586	gezelter	1358	template<typename Real>
587			Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
588			Polynomial<Real> p(p1);
589	gezelter	246
590	gezelter	1358	typename Polynomial<Real>::const_iterator i;
591	gezelter	246
592			for (i = p2.begin(); i != p2.end(); ++i) {
593	gezelter	507	p.addCoefficient(i->first, i->second);
594	gezelter	246	}
595
596			return p;
597
598	gezelter	507	}
599	gezelter	246
600	gezelter	507	/**
601			* Generates and returns the difference of two given Polynomials.
602			* @return
603			* @param p1 the first polynomial
604			* @param p2 the second polynomial
605			*/
606	gezelter	1358	template<typename Real>
607			Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
608			Polynomial<Real> p(p1);
609	gezelter	246
610	gezelter	1358	typename Polynomial<Real>::const_iterator i;
611	gezelter	246
612			for (i = p2.begin(); i != p2.end(); ++i) {
613	gezelter	507	p.addCoefficient(i->first, -i->second);
614	gezelter	246	}
615
616			return p;
617
618	gezelter	507	}
619	gezelter	246
620	gezelter	507	/**
621	gezelter	1358	* Returns the first derivative of this polynomial.
622			* @return the first derivative of this polynomial
623			*/
624			template<typename Real>
625	gezelter	1874	Polynomial<Real> * getDerivative(const Polynomial<Real>& p1) {
626			Polynomial<Real> * p = new Polynomial<Real>();
627	gezelter	1358
628			typename Polynomial<Real>::const_iterator i;
629	cli2	1362	int exponent;
630			Real coefficient;
631	gezelter	1358
632			for (i = p1.begin(); i != p1.end(); ++i) {
633			exponent = i->first;
634			coefficient = i->second;
635	gezelter	1874	p->setCoefficient(exponent-1, coefficient * exponent);
636	gezelter	1358	}
637
638			return p;
639			}
640
641			/**
642	gezelter	507	* Tests if two polynomial have the same exponents
643	tim	883	* @return true if all of the exponents in these Polynomial are identical
644	gezelter	507	* @param p1 the first polynomial
645			* @param p2 the second polynomial
646			* @note this function does not compare the coefficient
647			*/
648	gezelter	1358	template<typename Real>
649			bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
650	gezelter	246
651	gezelter	1358	typename Polynomial<Real>::const_iterator i;
652			typename Polynomial<Real>::const_iterator j;
653	gezelter	246
654			if (p1.size() != p2.size() ) {
655	gezelter	507	return false;
656	gezelter	246	}
657
658			for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
659	gezelter	507	if (i->first != j->first) {
660			return false;
661			}
662	gezelter	246	}
663
664			return true;
665	gezelter	507	}
666	gezelter	246
667	gezelter	1358
668
669	tim	963	typedef Polynomial<RealType> DoublePolynomial;
670	gezelter	246
671	gezelter	1390	} //end namespace OpenMD
672	gezelter	246	#endif //MATH_POLYNOMIAL_HPP