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