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
     * @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;
      
      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;
      }

      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:	1827
Committed:	Wed Jan 9 19:49:07 2013 UTC (12 years, 3 months ago) by gezelter
File size:	20210 byte(s)
Log Message:	More cleaning, removing non-C++ MPI calls
#	Content
1	/*
2	* 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	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* 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	*
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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
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	#include <complex>
58	#include "config.h"
59	#include "math/Eigenvalue.hpp"
60
61	namespace OpenMD {
62
63	template<typename Real> Real fastpow(Real x, int N) {
64	Real result(1); //or 1.0?
65
66	for (int i = 0; i < N; ++i) {
67	result *= x;
68	}
69
70	return result;
71	}
72
73	/**
74	* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
75	* A generic Polynomial class
76	*/
77	template<typename Real>
78	class Polynomial {
79
80	public:
81	typedef Polynomial<Real> PolynomialType;
82	typedef int ExponentType;
83	typedef Real CoefficientType;
84	typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
85	typedef typename PolynomialPairMap::iterator iterator;
86	typedef typename PolynomialPairMap::const_iterator const_iterator;
87
88	Polynomial() {}
89	Polynomial(Real v) {setCoefficient(0, v);}
90	/**
91	* Calculates the value of this Polynomial evaluated at the given x value.
92	* @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	*/
96	Real evaluate(const Real& x) {
97	Real result = Real();
98	ExponentType exponent;
99	CoefficientType coefficient;
100
101	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
102	exponent = i->first;
103	coefficient = i->second;
104	result += fastpow(x, exponent) * coefficient;
105	}
106
107	return result;
108	}
109
110	/**
111	* Returns the first derivative of this polynomial.
112	* @return the first derivative of this polynomial
113	* @param x
114	*/
115	Real evaluateDerivative(const Real& x) {
116	Real result = Real();
117	ExponentType exponent;
118	CoefficientType coefficient;
119
120	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
121	exponent = i->first;
122	coefficient = i->second;
123	result += fastpow(x, exponent - 1) * coefficient * exponent;
124	}
125
126	return result;
127	}
128
129
130	/**
131	* Set the coefficent of the specified exponent, if the
132	* coefficient is already there, it will be overwritten.
133	* @param exponent exponent of a term in this Polynomial
134	* @param coefficient multiplier of a term in this Polynomial
135	*/
136	void setCoefficient(int exponent, const Real& coefficient) {
137	polyPairMap_[exponent] = coefficient;
138	}
139
140	/**
141	* 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	* @param exponent exponent of a term in this Polynomial
145	* @param coefficient multiplier of a term in this Polynomial
146	*/
147	void addCoefficient(int exponent, const Real& coefficient) {
148	iterator i = polyPairMap_.find(exponent);
149
150	if (i != end()) {
151	i->second += coefficient;
152	} else {
153	setCoefficient(exponent, coefficient);
154	}
155	}
156
157	/**
158	* 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	* @param exponent exponent of any term in this Polynomial
163	*/
164	Real getCoefficient(ExponentType exponent) {
165	iterator i = polyPairMap_.find(exponent);
166
167	if (i != end()) {
168	return i->second;
169	} else {
170	return Real(0);
171	}
172	}
173
174	iterator begin() {
175	return polyPairMap_.begin();
176	}
177
178	const_iterator begin() const{
179	return polyPairMap_.begin();
180	}
181
182	iterator end() {
183	return polyPairMap_.end();
184	}
185
186	const_iterator end() const{
187	return polyPairMap_.end();
188	}
189
190	iterator find(ExponentType exponent) {
191	return polyPairMap_.find(exponent);
192	}
193
194	size_t size() {
195	return polyPairMap_.size();
196	}
197
198	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	PolynomialType& operator = (const PolynomialType& p) {
208
209	if (this != &p) // protect against invalid self-assignment
210	{
211	typename Polynomial<Real>::const_iterator i;
212
213	polyPairMap_.clear(); // clear out the old map
214
215	for (i = p.begin(); i != p.end(); ++i) {
216	this->setCoefficient(i->first, i->second);
217	}
218	}
219	// by convention, always return *this
220	return *this;
221	}
222
223	PolynomialType& operator += (const PolynomialType& p) {
224	typename Polynomial<Real>::const_iterator i;
225
226	for (i = p.begin(); i != p.end(); ++i) {
227	this->addCoefficient(i->first, i->second);
228	}
229
230	return *this;
231	}
232
233	PolynomialType& operator -= (const PolynomialType& p) {
234	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	}
240
241	PolynomialType& operator *= (const PolynomialType& p) {
242	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	}
252	return *this;
253	}
254
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	}
266
267	PolynomialType& operator += (const Real v) {
268	this->addCoefficient( 0, v);
269	return *this;
270	}
271
272	/**
273	* Returns the first derivative of this polynomial.
274	* @return the first derivative of this polynomial
275	*/
276	PolynomialType & getDerivative() {
277	Polynomial<Real> p;
278
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	}
291
292	// 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	}
308
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	roots.push_back(complex<Real>(reals(i), imags(i)));
321	}
322
323	return roots;
324	}
325
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
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	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
467	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	return roots; // should be empty if you got here
532	}
533
534	private:
535
536	PolynomialPairMap polyPairMap_;
537	};
538
539
540	/**
541	* Generates and returns the product of two given Polynomials.
542	* @return A Polynomial containing the product of the two given Polynomial parameters
543	*/
544	template<typename Real>
545	Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
546	typename Polynomial<Real>::const_iterator i;
547	typename Polynomial<Real>::const_iterator j;
548	Polynomial<Real> p;
549
550	for (i = p1.begin(); i !=p1.end(); ++i) {
551	for (j = p2.begin(); j !=p2.end(); ++j) {
552	p.addCoefficient( i->first + j->first, i->second * j->second);
553	}
554	}
555
556	return p;
557	}
558
559	template<typename Real>
560	Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
561	typename Polynomial<Real>::const_iterator i;
562	Polynomial<Real> result;
563
564	for (i = p.begin(); i !=p.end(); ++i) {
565	result.setCoefficient( i->first , i->second * v);
566	}
567
568	return result;
569	}
570
571	template<typename Real>
572	Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
573	typename Polynomial<Real>::const_iterator i;
574	Polynomial<Real> result;
575
576	for (i = p.begin(); i !=p.end(); ++i) {
577	result.setCoefficient( i->first , i->second * v);
578	}
579
580	return result;
581	}
582
583	/**
584	* Generates and returns the sum of two given Polynomials.
585	* @param p1 the first polynomial
586	* @param p2 the second polynomial
587	*/
588	template<typename Real>
589	Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
590	Polynomial<Real> p(p1);
591
592	typename Polynomial<Real>::const_iterator i;
593
594	for (i = p2.begin(); i != p2.end(); ++i) {
595	p.addCoefficient(i->first, i->second);
596	}
597
598	return p;
599
600	}
601
602	/**
603	* Generates and returns the difference of two given Polynomials.
604	* @return
605	* @param p1 the first polynomial
606	* @param p2 the second polynomial
607	*/
608	template<typename Real>
609	Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
610	Polynomial<Real> p(p1);
611
612	typename Polynomial<Real>::const_iterator i;
613
614	for (i = p2.begin(); i != p2.end(); ++i) {
615	p.addCoefficient(i->first, -i->second);
616	}
617
618	return p;
619
620	}
621
622	/**
623	* Returns the first derivative of this polynomial.
624	* @return the first derivative of this polynomial
625	*/
626	template<typename Real>
627	Polynomial<Real> getDerivative(const Polynomial<Real>& p1) {
628	Polynomial<Real> p;
629
630	typename Polynomial<Real>::const_iterator i;
631	int exponent;
632	Real coefficient;
633
634	for (i = p1.begin(); i != p1.end(); ++i) {
635	exponent = i->first;
636	coefficient = i->second;
637	p.setCoefficient(exponent-1, coefficient * exponent);
638	}
639
640	return p;
641	}
642
643	/**
644	* Tests if two polynomial have the same exponents
645	* @return true if all of the exponents in these Polynomial are identical
646	* @param p1 the first polynomial
647	* @param p2 the second polynomial
648	* @note this function does not compare the coefficient
649	*/
650	template<typename Real>
651	bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
652
653	typename Polynomial<Real>::const_iterator i;
654	typename Polynomial<Real>::const_iterator j;
655
656	if (p1.size() != p2.size() ) {
657	return false;
658	}
659
660	for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
661	if (i->first != j->first) {
662	return false;
663	}
664	}
665
666	return true;
667	}
668
669
670
671	typedef Polynomial<RealType> DoublePolynomial;
672
673	} //end namespace OpenMD
674	#endif //MATH_POLYNOMIAL_HPP