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. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. 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.
 */
 
/**
 * @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 oopse {
  
  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 oopse
#endif //MATH_POLYNOMIAL_HPP
Revision:	1362
Committed:	Tue Sep 8 15:48:33 2009 UTC (16 years, 1 month ago) by cli2
File size:	20174 byte(s)
Log Message:	Fixed bug in Polynomial
#	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. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
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	#include <complex>
57	#include "config.h"
58	#include "math/Eigenvalue.hpp"
59
60	namespace oopse {
61
62	template<typename Real> Real fastpow(Real x, int N) {
63	Real result(1); //or 1.0?
64
65	for (int i = 0; i < N; ++i) {
66	result *= x;
67	}
68
69	return result;
70	}
71
72	/**
73	* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
74	* A generic Polynomial class
75	*/
76	template<typename Real>
77	class Polynomial {
78
79	public:
80	typedef Polynomial<Real> PolynomialType;
81	typedef int ExponentType;
82	typedef Real CoefficientType;
83	typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
84	typedef typename PolynomialPairMap::iterator iterator;
85	typedef typename PolynomialPairMap::const_iterator const_iterator;
86
87	Polynomial() {}
88	Polynomial(Real v) {setCoefficient(0, v);}
89	/**
90	* Calculates the value of this Polynomial evaluated at the given x value.
91	* @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	*/
95	Real evaluate(const Real& x) {
96	Real result = Real();
97	ExponentType exponent;
98	CoefficientType coefficient;
99
100	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
101	exponent = i->first;
102	coefficient = i->second;
103	result += fastpow(x, exponent) * coefficient;
104	}
105
106	return result;
107	}
108
109	/**
110	* Returns the first derivative of this polynomial.
111	* @return the first derivative of this polynomial
112	* @param x
113	*/
114	Real evaluateDerivative(const Real& x) {
115	Real result = Real();
116	ExponentType exponent;
117	CoefficientType coefficient;
118
119	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
120	exponent = i->first;
121	coefficient = i->second;
122	result += fastpow(x, exponent - 1) * coefficient * exponent;
123	}
124
125	return result;
126	}
127
128
129	/**
130	* Set the coefficent of the specified exponent, if the
131	* coefficient is already there, it will be overwritten.
132	* @param exponent exponent of a term in this Polynomial
133	* @param coefficient multiplier of a term in this Polynomial
134	*/
135	void setCoefficient(int exponent, const Real& coefficient) {
136	polyPairMap_[exponent] = coefficient;
137	}
138
139	/**
140	* 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	* @param exponent exponent of a term in this Polynomial
144	* @param coefficient multiplier of a term in this Polynomial
145	*/
146	void addCoefficient(int exponent, const Real& coefficient) {
147	iterator i = polyPairMap_.find(exponent);
148
149	if (i != end()) {
150	i->second += coefficient;
151	} else {
152	setCoefficient(exponent, coefficient);
153	}
154	}
155
156	/**
157	* 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	* @exponent exponent of any term in this Polynomial
162	*/
163	Real getCoefficient(ExponentType exponent) {
164	iterator i = polyPairMap_.find(exponent);
165
166	if (i != end()) {
167	return i->second;
168	} else {
169	return Real(0);
170	}
171	}
172
173	iterator begin() {
174	return polyPairMap_.begin();
175	}
176
177	const_iterator begin() const{
178	return polyPairMap_.begin();
179	}
180
181	iterator end() {
182	return polyPairMap_.end();
183	}
184
185	const_iterator end() const{
186	return polyPairMap_.end();
187	}
188
189	iterator find(ExponentType exponent) {
190	return polyPairMap_.find(exponent);
191	}
192
193	size_t size() {
194	return polyPairMap_.size();
195	}
196
197	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	PolynomialType& operator = (const PolynomialType& p) {
207
208	if (this != &p) // protect against invalid self-assignment
209	{
210	typename Polynomial<Real>::const_iterator i;
211
212	polyPairMap_.clear(); // clear out the old map
213
214	for (i = p.begin(); i != p.end(); ++i) {
215	this->setCoefficient(i->first, i->second);
216	}
217	}
218	// by convention, always return *this
219	return *this;
220	}
221
222	PolynomialType& operator += (const PolynomialType& p) {
223	typename Polynomial<Real>::const_iterator i;
224
225	for (i = p.begin(); i != p.end(); ++i) {
226	this->addCoefficient(i->first, i->second);
227	}
228
229	return *this;
230	}
231
232	PolynomialType& operator -= (const PolynomialType& p) {
233	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	}
239
240	PolynomialType& operator *= (const PolynomialType& p) {
241	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	}
251	return *this;
252	}
253
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	}
265
266	PolynomialType& operator += (const Real v) {
267	this->addCoefficient( 0, v);
268	return *this;
269	}
270
271	/**
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	}
290
291	// 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	}
307
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	roots.push_back(complex<Real>(reals(i), imags(i)));
320	}
321
322	return roots;
323	}
324
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
470	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	private:
540
541	PolynomialPairMap polyPairMap_;
542	};
543
544
545	/**
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	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
555	for (i = p1.begin(); i !=p1.end(); ++i) {
556	for (j = p2.begin(); j !=p2.end(); ++j) {
557	p.addCoefficient( i->first + j->first, i->second * j->second);
558	}
559	}
560
561	return p;
562	}
563
564	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
569	for (i = p.begin(); i !=p.end(); ++i) {
570	result.setCoefficient( i->first , i->second * v);
571	}
572
573	return result;
574	}
575
576	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
581	for (i = p.begin(); i !=p.end(); ++i) {
582	result.setCoefficient( i->first , i->second * v);
583	}
584
585	return result;
586	}
587
588	/**
589	* Generates and returns the sum of two given Polynomials.
590	* @param p1 the first polynomial
591	* @param p2 the second polynomial
592	*/
593	template<typename Real>
594	Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
595	Polynomial<Real> p(p1);
596
597	typename Polynomial<Real>::const_iterator i;
598
599	for (i = p2.begin(); i != p2.end(); ++i) {
600	p.addCoefficient(i->first, i->second);
601	}
602
603	return p;
604
605	}
606
607	/**
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	template<typename Real>
614	Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
615	Polynomial<Real> p(p1);
616
617	typename Polynomial<Real>::const_iterator i;
618
619	for (i = p2.begin(); i != p2.end(); ++i) {
620	p.addCoefficient(i->first, -i->second);
621	}
622
623	return p;
624
625	}
626
627	/**
628	* 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	int exponent;
637	Real coefficient;
638
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	* Tests if two polynomial have the same exponents
650	* @return true if all of the exponents in these Polynomial are identical
651	* @param p1 the first polynomial
652	* @param p2 the second polynomial
653	* @note this function does not compare the coefficient
654	*/
655	template<typename Real>
656	bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
657
658	typename Polynomial<Real>::const_iterator i;
659	typename Polynomial<Real>::const_iterator j;
660
661	if (p1.size() != p2.size() ) {
662	return false;
663	}
664
665	for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
666	if (i->first != j->first) {
667	return false;
668	}
669	}
670
671	return true;
672	}
673
674
675
676	typedef Polynomial<RealType> DoublePolynomial;
677
678	} //end namespace oopse
679	#endif //MATH_POLYNOMIAL_HPP