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, 24107 (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
     * @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;
      
      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;
      }
        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<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;
        break;
      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;
        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<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;
    
    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:	1782
Committed:	Wed Aug 22 02:28:28 2012 UTC (13 years, 2 months ago) by gezelter
File size:	20279 byte(s)
Log Message:	MERGE OpenMD development branch 1465:1781 into trunk
#	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			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1782	* [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	507	* @exponent exponent of any term in this Polynomial
163			*/
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			PolynomialType & getDerivative() {
277	gezelter	1782	Polynomial<Real> p;
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			p.setCoefficient(exponent-1, coefficient * exponent);
287			}
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			break;
342			case 2: {
343			Real fC2 = getCoefficient(2);
344			Real fC1 = getCoefficient(1);
345			Real fC0 = getCoefficient(0);
346			Real fDiscr = fC1fC1 - 4.0fC0*fC2;
347			if (abs(fDiscr) <= fEpsilon) {
348			fDiscr = (Real)0.0;
349			}
350
351			if (fDiscr < (Real)0.0) { // complex roots only
352			return roots;
353			}
354
355			Real fTmp = ((Real)0.5)/fC2;
356
357			if (fDiscr > (Real)0.0) { // 2 real roots
358			fDiscr = sqrt(fDiscr);
359			roots.push_back(fTmp*(-fC1 - fDiscr));
360			roots.push_back(fTmp*(-fC1 + fDiscr));
361			} else {
362			roots.push_back(-fTmp * fC1); // 1 real root
363			}
364			}
365			return roots;
366			break;
367
368			case 3: {
369			Real fC3 = getCoefficient(3);
370			Real fC2 = getCoefficient(2);
371			Real fC1 = getCoefficient(1);
372			Real fC0 = getCoefficient(0);
373
374			// make polynomial monic, x^3+c2x^2+c1x+c0
375			Real fInvC3 = ((Real)1.0)/fC3;
376			fC0 *= fInvC3;
377			fC1 *= fInvC3;
378			fC2 *= fInvC3;
379
380			// convert to y^3+a*y+b = 0 by x = y-c2/3
381			const Real fThird = (Real)1.0/(Real)3.0;
382			const Real fTwentySeventh = (Real)1.0/(Real)27.0;
383			Real fOffset = fThird*fC2;
384			Real fA = fC1 - fC2*fOffset;
385			Real fB = fC0+fC2(((Real)2.0)fC2fC2-((Real)9.0)fC1)*fTwentySeventh;
386			Real fHalfB = ((Real)0.5)*fB;
387
388			Real fDiscr = fHalfBfHalfB + fAfAfAfTwentySeventh;
389			if (fabs(fDiscr) <= fEpsilon) {
390			fDiscr = (Real)0.0;
391			}
392
393			if (fDiscr > (Real)0.0) { // 1 real, 2 complex roots
394
395			fDiscr = sqrt(fDiscr);
396			Real fTemp = -fHalfB + fDiscr;
397			Real root;
398			if (fTemp >= (Real)0.0) {
399			root = pow(fTemp,fThird);
400			} else {
401			root = -pow(-fTemp,fThird);
402			}
403			fTemp = -fHalfB - fDiscr;
404			if ( fTemp >= (Real)0.0 ) {
405			root += pow(fTemp,fThird);
406			} else {
407			root -= pow(-fTemp,fThird);
408			}
409			root -= fOffset;
410
411			roots.push_back(root);
412			} else if (fDiscr < (Real)0.0) {
413			const Real fSqrt3 = sqrt((Real)3.0);
414			Real fDist = sqrt(-fThird*fA);
415			Real fAngle = fThird*atan2(sqrt(-fDiscr), -fHalfB);
416			Real fCos = cos(fAngle);
417			Real fSin = sin(fAngle);
418			roots.push_back(((Real)2.0)fDistfCos-fOffset);
419			roots.push_back(-fDist(fCos+fSqrt3fSin)-fOffset);
420			roots.push_back(-fDist(fCos-fSqrt3fSin)-fOffset);
421			} else {
422			Real fTemp;
423			if (fHalfB >= (Real)0.0) {
424			fTemp = -pow(fHalfB,fThird);
425			} else {
426			fTemp = pow(-fHalfB,fThird);
427			}
428			roots.push_back(((Real)2.0)*fTemp-fOffset);
429			roots.push_back(-fTemp-fOffset);
430			roots.push_back(-fTemp-fOffset);
431			}
432			}
433			return roots;
434			break;
435			case 4: {
436			Real fC4 = getCoefficient(4);
437			Real fC3 = getCoefficient(3);
438			Real fC2 = getCoefficient(2);
439			Real fC1 = getCoefficient(1);
440			Real fC0 = getCoefficient(0);
441
442			// make polynomial monic, x^4+c3x^3+c2x^2+c1*x+c0
443			Real fInvC4 = ((Real)1.0)/fC4;
444			fC0 *= fInvC4;
445			fC1 *= fInvC4;
446			fC2 *= fInvC4;
447			fC3 *= fInvC4;
448
449			// reduction to resolvent cubic polynomial y^3+r2y^2+r1y+r0 = 0
450			Real fR0 = -fC3fC3fC0 + ((Real)4.0)fC2fC0 - fC1*fC1;
451			Real fR1 = fC3fC1 - ((Real)4.0)fC0;
452			Real fR2 = -fC2;
453			Polynomial<Real> tempCubic;
454			tempCubic.setCoefficient(0, fR0);
455			tempCubic.setCoefficient(1, fR1);
456			tempCubic.setCoefficient(2, fR2);
457			tempCubic.setCoefficient(3, 1.0);
458			std::vector<Real> cubeRoots = tempCubic.FindRealRoots(); // always
459			// produces
460			// at
461			// least
462			// one
463			// root
464			Real fY = cubeRoots[0];
465
466			Real fDiscr = ((Real)0.25)fC3fC3 - fC2 + fY;
467			if (fabs(fDiscr) <= fEpsilon) {
468			fDiscr = (Real)0.0;
469			}
470	tim	749
471	gezelter	1358	if (fDiscr > (Real)0.0) {
472			Real fR = sqrt(fDiscr);
473			Real fT1 = ((Real)0.75)fC3fC3 - fRfR - ((Real)2.0)fC2;
474			Real fT2 = (((Real)4.0)fC3fC2 - ((Real)8.0)fC1 - fC3fC3*fC3) /
475			(((Real)4.0)*fR);
476
477			Real fTplus = fT1+fT2;
478			Real fTminus = fT1-fT2;
479			if (fabs(fTplus) <= fEpsilon) {
480			fTplus = (Real)0.0;
481			}
482			if (fabs(fTminus) <= fEpsilon) {
483			fTminus = (Real)0.0;
484			}
485
486			if (fTplus >= (Real)0.0) {
487			Real fD = sqrt(fTplus);
488			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fR+fD));
489			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fR-fD));
490			}
491			if (fTminus >= (Real)0.0) {
492			Real fE = sqrt(fTminus);
493			roots.push_back(-((Real)0.25)fC3+((Real)0.5)(fE-fR));
494			roots.push_back(-((Real)0.25)fC3-((Real)0.5)(fE+fR));
495			}
496			} else if (fDiscr < (Real)0.0) {
497			//roots.clear();
498			} else {
499			Real fT2 = fYfY-((Real)4.0)fC0;
500			if (fT2 >= -fEpsilon) {
501			if (fT2 < (Real)0.0) { // round to zero
502			fT2 = (Real)0.0;
503			}
504			fT2 = ((Real)2.0)*sqrt(fT2);
505			Real fT1 = ((Real)0.75)fC3fC3 - ((Real)2.0)*fC2;
506			if (fT1+fT2 >= fEpsilon) {
507			Real fD = sqrt(fT1+fT2);
508			roots.push_back( -((Real)0.25)fC3+((Real)0.5)fD);
509			roots.push_back( -((Real)0.25)fC3-((Real)0.5)fD);
510			}
511			if (fT1-fT2 >= fEpsilon) {
512			Real fE = sqrt(fT1-fT2);
513			roots.push_back( -((Real)0.25)fC3+((Real)0.5)fE);
514			roots.push_back( -((Real)0.25)fC3-((Real)0.5)fE);
515			}
516			}
517			}
518			}
519			return roots;
520			break;
521			default: {
522			DynamicRectMatrix<Real> companion = CreateCompanion();
523			JAMA::Eigenvalue<Real> eig(companion);
524			DynamicVector<Real> reals, imags;
525			eig.getRealEigenvalues(reals);
526			eig.getImagEigenvalues(imags);
527
528			for (int i = 0; i < deg; i++) {
529			if (fabs(imags(i)) < fEpsilon)
530			roots.push_back(reals(i));
531			}
532			}
533			return roots;
534			break;
535			}
536
537			return roots; // should be empty if you got here
538			}
539
540	gezelter	507	private:
541	gezelter	246
542	gezelter	507	PolynomialPairMap polyPairMap_;
543			};
544	gezelter	246
545	gezelter	1358
546	gezelter	507	/**
547			* Generates and returns the product of two given Polynomials.
548			* @return A Polynomial containing the product of the two given Polynomial parameters
549			*/
550	gezelter	1358	template<typename Real>
551			Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
552			typename Polynomial<Real>::const_iterator i;
553			typename Polynomial<Real>::const_iterator j;
554			Polynomial<Real> p;
555	gezelter	246
556			for (i = p1.begin(); i !=p1.end(); ++i) {
557	gezelter	507	for (j = p2.begin(); j !=p2.end(); ++j) {
558			p.addCoefficient( i->first + j->first, i->second * j->second);
559			}
560	gezelter	246	}
561
562			return p;
563	gezelter	507	}
564	gezelter	246
565	gezelter	1358	template<typename Real>
566			Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
567			typename Polynomial<Real>::const_iterator i;
568			Polynomial<Real> result;
569	tim	876
570			for (i = p.begin(); i !=p.end(); ++i) {
571	cli2	1290	result.setCoefficient( i->first , i->second * v);
572	tim	876	}
573
574			return result;
575			}
576
577	gezelter	1358	template<typename Real>
578			Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
579			typename Polynomial<Real>::const_iterator i;
580			Polynomial<Real> result;
581	tim	876
582			for (i = p.begin(); i !=p.end(); ++i) {
583	cli2	1290	result.setCoefficient( i->first , i->second * v);
584	tim	876	}
585
586			return result;
587			}
588
589	gezelter	507	/**
590			* Generates and returns the sum of two given Polynomials.
591			* @param p1 the first polynomial
592			* @param p2 the second polynomial
593			*/
594	gezelter	1358	template<typename Real>
595			Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
596			Polynomial<Real> p(p1);
597	gezelter	246
598	gezelter	1358	typename Polynomial<Real>::const_iterator i;
599	gezelter	246
600			for (i = p2.begin(); i != p2.end(); ++i) {
601	gezelter	507	p.addCoefficient(i->first, i->second);
602	gezelter	246	}
603
604			return p;
605
606	gezelter	507	}
607	gezelter	246
608	gezelter	507	/**
609			* Generates and returns the difference of two given Polynomials.
610			* @return
611			* @param p1 the first polynomial
612			* @param p2 the second polynomial
613			*/
614	gezelter	1358	template<typename Real>
615			Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
616			Polynomial<Real> p(p1);
617	gezelter	246
618	gezelter	1358	typename Polynomial<Real>::const_iterator i;
619	gezelter	246
620			for (i = p2.begin(); i != p2.end(); ++i) {
621	gezelter	507	p.addCoefficient(i->first, -i->second);
622	gezelter	246	}
623
624			return p;
625
626	gezelter	507	}
627	gezelter	246
628	gezelter	507	/**
629	gezelter	1358	* Returns the first derivative of this polynomial.
630			* @return the first derivative of this polynomial
631			*/
632			template<typename Real>
633			Polynomial<Real> getDerivative(const Polynomial<Real>& p1) {
634	gezelter	1782	Polynomial<Real> p;
635	gezelter	1358
636			typename Polynomial<Real>::const_iterator i;
637	cli2	1362	int exponent;
638			Real coefficient;
639	gezelter	1358
640			for (i = p1.begin(); i != p1.end(); ++i) {
641			exponent = i->first;
642			coefficient = i->second;
643			p.setCoefficient(exponent-1, coefficient * exponent);
644			}
645
646			return p;
647			}
648
649			/**
650	gezelter	507	* Tests if two polynomial have the same exponents
651	tim	883	* @return true if all of the exponents in these Polynomial are identical
652	gezelter	507	* @param p1 the first polynomial
653			* @param p2 the second polynomial
654			* @note this function does not compare the coefficient
655			*/
656	gezelter	1358	template<typename Real>
657			bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
658	gezelter	246
659	gezelter	1358	typename Polynomial<Real>::const_iterator i;
660			typename Polynomial<Real>::const_iterator j;
661	gezelter	246
662			if (p1.size() != p2.size() ) {
663	gezelter	507	return false;
664	gezelter	246	}
665
666			for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
667	gezelter	507	if (i->first != j->first) {
668			return false;
669			}
670	gezelter	246	}
671
672			return true;
673	gezelter	507	}
674	gezelter	246
675	gezelter	1358
676
677	tim	963	typedef Polynomial<RealType> DoublePolynomial;
678	gezelter	246
679	gezelter	1390	} //end namespace OpenMD
680	gezelter	246	#endif //MATH_POLYNOMIAL_HPP