src/math/Polynomial.hpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @file Polynomial.hpp
 * @author    teng lin
 * @date  11/16/2004
 * @version 1.0
 */ 

#ifndef MATH_POLYNOMIAL_HPP
#define MATH_POLYNOMIAL_HPP

#include <list>
#include <utility>

namespace oopse {

template<typename ElemType> pow(ElemType x, int N) {
    ElemType result(1);

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

    return result;
}

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

    public:
        
        typedef int ExponentType;
        typedef ElemType CoefficientType;
        typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
        typedef PolynomialPairMap::iterator PolynomialIterator;

        /** 
         * 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
         */
        ElemType evaluate(const ElemType& x) {
            ElemType result;
            double exponent;
            double coefficient;
            
            for (PolynomialIterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
                exponent = i->first;
                coefficient = i->second;
                result  += pow(x, exponent) * coefficient;
            }

            return result;
        }

        /**
         * Returns the first derivative of this polynomial.
         * @return the first derivative of this polynomial
         * @param x
         */
        ElemType evaluateDerivative(const ElemType& x) {
            ElemType result;
            double exponent;
            double coefficient;
            
            for (PolynomialIterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
                exponent = i->first;
                coefficient = i->second;
                result  += pow(x, exponent - 1) * coefficient * exponent;
            }

            return result;
        }

        /**
         * @param exponent exponent of a term in this Polynomial 
         * @param coefficient multiplier of a term in this Polynomial 
         */
        
        void setCoefficient(int exponent, const ElemType& coefficient) {
            polyPairMap_.insert(PolynomialPairMap::value_type(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
         */
        ElemType getCoefficient(ExponentType exponent) {
            PolynomialIterator i = find(exponent);

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

        PolynomialIterator begin() {
            return polyPairMap_.begin();
        }
        
        PolynomialIterator end() {
            return polyPairMap_.end();
        }

        PolynomialIterator find(ExponentType exponent) {
            return polyPairMap_.find();
        }

        size_t size() {
            return polyPairMap_.size();
        }
        
    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 ElemType>
Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
    typename Polynomial<ElemType>::PolynomialIterator i;
    typename Polynomial<ElemType>::PolynomialIterator j;
    typename Polynomial<ElemType>::PolynomialIterator k;
    Polynomial<ElemType> p;
    int exponent;
    int coefficient;
    
    for (i = p1.begin(); i !=p1.end(); ++i) {
        for (j = p1.begjn(); j !=p1.end(); ++j) {
            exponent = i->first + j->first;            
            coefficient = i->second * j->second;
            k = p->find(exponent);

            if (k != p.end()) {
                p[exponent] = coefficient;
            } else {
                k->second += coefficient;
            }
        }
    }
}

/**
 * Generates and returns the sum of two given Polynomials.
 * @param p1 the first polynomial
 * @param p2 the second polynomial
 */
template<typename ElemType>
Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
    Polynomial<ElemType> p(p1);

    typename Polynomial<ElemType>::PolynomialIterator i;
    typename Polynomial<ElemType>::PolynomialIterator j;

    for (i =  p2.begin(); i  != p2.end(); ++i) {
        j = p.find(i->first);
        if (j == p.end()) {
            p[j] = i->second;
        } else {
            j->second += 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 ElemType>
Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
    Polynomial<ElemType> p(p1);

    typename Polynomial<ElemType>::PolynomialIterator i;
    typename Polynomial<ElemType>::PolynomialIterator j;

    for (i =  p2.begin(); i  != p2.end(); ++i) {
        j = p.find(i->first);
        if (j == p.end()) {
            p[j] = -i->second;
        } else {
            j->second -= i->second;
        }
    }

    return p;

}

/**
 * Tests if two polynomial have the same exponents
 * @return true if these 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 ElemType>
bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {

    typename Polynomial<ElemType>::PolynomialIterator i;
    typename Polynomial<ElemType>::PolynomialIterator 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<double> DoublePolynomial;

} //end namespace oopse
#endif //MATH_POLYNOMIAL_HPP
Revision:	1746
Committed:	Wed Nov 17 06:37:56 2004 UTC (19 years, 9 months ago) by tim
File size:	7585 byte(s)
Log Message:	add PolynomialBondType, PolynomialBendType, PolynomialTorsionType, HarmonicBendType and CharmmTorsionType. Need to refine the design and add document for them
#	User	Rev	Content
1	tim	1744	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file Polynomial.hpp
28	tim	1746	* @author teng lin
29			* @date 11/16/2004
30	tim	1744	* @version 1.0
31			*/
32
33			#ifndef MATH_POLYNOMIAL_HPP
34			#define MATH_POLYNOMIAL_HPP
35
36			#include <list>
37			#include <utility>
38
39			namespace oopse {
40
41	tim	1746	template<typename ElemType> pow(ElemType x, int N) {
42	tim	1744	ElemType result(1);
43
44			for (int i = 0; i < N; ++i) {
45			result *= x;
46			}
47
48			return result;
49			}
50
51	tim	1746	/**
52			* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
53			* A generic Polynomial class
54			*/
55	tim	1744	template<typename ElemType>
56			class Polynomial {
57
58			public:
59
60			typedef int ExponentType;
61	tim	1746	typedef ElemType CoefficientType;
62			typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
63	tim	1744	typedef PolynomialPairMap::iterator PolynomialIterator;
64
65			/**
66	tim	1746	* Calculates the value of this Polynomial evaluated at the given x value.
67			* @return The value of this Polynomial evaluates at the given x value
68			* @param x the value of the independent variable for this Polynomial function
69	tim	1744	*/
70			ElemType evaluate(const ElemType& x) {
71			ElemType result;
72			double exponent;
73	tim	1746	double coefficient;
74	tim	1744
75			for (PolynomialIterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
76			exponent = i->first;
77	tim	1746	coefficient = i->second;
78			result += pow(x, exponent) * coefficient;
79	tim	1744	}
80
81			return result;
82			}
83
84	tim	1746	/**
85			* Returns the first derivative of this polynomial.
86			* @return the first derivative of this polynomial
87			* @param x
88			*/
89			ElemType evaluateDerivative(const ElemType& x) {
90	tim	1744	ElemType result;
91			double exponent;
92	tim	1746	double coefficient;
93	tim	1744
94			for (PolynomialIterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
95			exponent = i->first;
96	tim	1746	coefficient = i->second;
97			result += pow(x, exponent - 1) * coefficient * exponent;
98	tim	1744	}
99
100			return result;
101			}
102
103	tim	1746	/**
104			* @param exponent exponent of a term in this Polynomial
105			* @param coefficient multiplier of a term in this Polynomial
106			*/
107
108			void setCoefficient(int exponent, const ElemType& coefficient) {
109			polyPairMap_.insert(PolynomialPairMap::value_type(exponent, coefficient));
110	tim	1744	}
111
112	tim	1746	/**
113			* Returns the coefficient associated with the given power for this Polynomial.
114			* @return the coefficient associated with the given power for this Polynomial
115			* @exponent exponent of any term in this Polynomial
116			*/
117			ElemType getCoefficient(ExponentType exponent) {
118			PolynomialIterator i = find(exponent);
119	tim	1744
120	tim	1746	if (i != end()) {
121			return i->second;
122			} else {
123			return ElemType(0);
124			}
125	tim	1744	}
126
127			PolynomialIterator begin() {
128			return polyPairMap_.begin();
129			}
130
131			PolynomialIterator end() {
132			return polyPairMap_.end();
133			}
134
135			PolynomialIterator find(ExponentType exponent) {
136			return polyPairMap_.find();
137			}
138	tim	1746
139			size_t size() {
140			return polyPairMap_.size();
141			}
142
143	tim	1744	private:
144
145			PolynomialPairMap polyPairMap_;
146			};
147
148	tim	1746
149			/**
150			* Generates and returns the product of two given Polynomials.
151			* @return A Polynomial containing the product of the two given Polynomial parameters
152			*/
153			template<typename ElemType>
154			Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
155			typename Polynomial<ElemType>::PolynomialIterator i;
156			typename Polynomial<ElemType>::PolynomialIterator j;
157			typename Polynomial<ElemType>::PolynomialIterator k;
158			Polynomial<ElemType> p;
159			int exponent;
160			int coefficient;
161
162			for (i = p1.begin(); i !=p1.end(); ++i) {
163			for (j = p1.begjn(); j !=p1.end(); ++j) {
164			exponent = i->first + j->first;
165			coefficient = i->second * j->second;
166			k = p->find(exponent);
167
168			if (k != p.end()) {
169			p[exponent] = coefficient;
170			} else {
171			k->second += coefficient;
172			}
173			}
174			}
175			}
176
177			/**
178			* Generates and returns the sum of two given Polynomials.
179			* @param p1 the first polynomial
180			* @param p2 the second polynomial
181			*/
182			template<typename ElemType>
183			Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
184			Polynomial<ElemType> p(p1);
185
186			typename Polynomial<ElemType>::PolynomialIterator i;
187			typename Polynomial<ElemType>::PolynomialIterator j;
188
189			for (i = p2.begin(); i != p2.end(); ++i) {
190			j = p.find(i->first);
191			if (j == p.end()) {
192			p[j] = i->second;
193			} else {
194			j->second += i->second;
195			}
196			}
197
198			return p;
199
200			}
201
202			/**
203			* Generates and returns the difference of two given Polynomials.
204			* @return
205			* @param p1 the first polynomial
206			* @param p2 the second polynomial
207			*/
208			template<typename ElemType>
209			Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
210			Polynomial<ElemType> p(p1);
211
212			typename Polynomial<ElemType>::PolynomialIterator i;
213			typename Polynomial<ElemType>::PolynomialIterator j;
214
215			for (i = p2.begin(); i != p2.end(); ++i) {
216			j = p.find(i->first);
217			if (j == p.end()) {
218			p[j] = -i->second;
219			} else {
220			j->second -= i->second;
221			}
222			}
223
224			return p;
225
226			}
227
228			/**
229			* Tests if two polynomial have the same exponents
230			* @return true if these all of the exponents in these Polynomial are identical
231			* @param p1 the first polynomial
232			* @param p2 the second polynomial
233			* @note this function does not compare the coefficient
234			*/
235			template<typename ElemType>
236			bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
237
238			typename Polynomial<ElemType>::PolynomialIterator i;
239			typename Polynomial<ElemType>::PolynomialIterator j;
240
241			if (p1.size() !== p2.size() ) {
242			return false;
243			}
244
245			for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
246			if (i->first != j->first) {
247			return false;
248			}
249			}
250
251			return true;
252			}
253
254			typedef Polynomial<double> DoublePolynomial;
255
256	tim	1744	} //end namespace oopse
257			#endif //MATH_POLYNOMIAL_HPP