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 <iostream>
#include <list>
#include <map>
#include <utility>

namespace oopse {

template<typename ElemType> 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 iterator;
        typedef PolynomialPairMap::const_iterator const_iterator;
        /** 
         * 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 = ElemType();
            double exponent;
            double coefficient;
            
            for (iterator 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 = ElemType();
            double exponent;
            double coefficient;
            
            for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
                exponent = i->first;
                coefficient = i->second;
                result  += pow(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 ElemType& coefficient) {
            polyPairMap_.insert(PolynomialPairMap::value_type(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 ElemType& 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
         */
        ElemType getCoefficient(ExponentType exponent) {
            iterator i = polyPairMap_.find(exponent);

            if (i != end()) {
                return i->second;
            } else {
                return ElemType(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();
        }
        
    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>::const_iterator i;
    typename Polynomial<ElemType>::const_iterator j;
    Polynomial<ElemType> 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;
}

/**
 * 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>::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 ElemType>
Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
    Polynomial<ElemType> p(p1);

    typename Polynomial<ElemType>::const_iterator i;

    for (i =  p2.begin(); i  != p2.end(); ++i) {
        p.addCoefficient(i->first, -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>::const_iterator i;
    typename Polynomial<ElemType>::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<double> DoublePolynomial;

} //end namespace oopse
#endif //MATH_POLYNOMIAL_HPP
Revision:	1796
Committed:	Mon Nov 29 21:42:08 2004 UTC (19 years, 9 months ago) by tim
File size:	7990 byte(s)
Log Message:	Polynomial passes unit test
#	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	tim	1796	#include <iostream>
37	tim	1744	#include <list>
38	tim	1796	#include <map>
39	tim	1744	#include <utility>
40
41			namespace oopse {
42
43	tim	1796	template<typename ElemType> ElemType pow(ElemType x, int N) {
44	tim	1744	ElemType result(1);
45
46			for (int i = 0; i < N; ++i) {
47			result *= x;
48			}
49
50			return result;
51			}
52
53	tim	1746	/**
54			* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
55			* A generic Polynomial class
56			*/
57	tim	1744	template<typename ElemType>
58			class Polynomial {
59
60			public:
61
62			typedef int ExponentType;
63	tim	1746	typedef ElemType CoefficientType;
64			typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
65	tim	1796	typedef PolynomialPairMap::iterator iterator;
66			typedef PolynomialPairMap::const_iterator const_iterator;
67	tim	1744	/**
68	tim	1746	* Calculates the value of this Polynomial evaluated at the given x value.
69			* @return The value of this Polynomial evaluates at the given x value
70			* @param x the value of the independent variable for this Polynomial function
71	tim	1744	*/
72			ElemType evaluate(const ElemType& x) {
73	tim	1796	ElemType result = ElemType();
74	tim	1744	double exponent;
75	tim	1746	double coefficient;
76	tim	1744
77	tim	1796	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
78	tim	1744	exponent = i->first;
79	tim	1746	coefficient = i->second;
80			result += pow(x, exponent) * coefficient;
81	tim	1744	}
82
83			return result;
84			}
85
86	tim	1746	/**
87			* Returns the first derivative of this polynomial.
88			* @return the first derivative of this polynomial
89			* @param x
90			*/
91			ElemType evaluateDerivative(const ElemType& x) {
92	tim	1796	ElemType result = ElemType();
93	tim	1744	double exponent;
94	tim	1746	double coefficient;
95	tim	1744
96	tim	1796	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
97	tim	1744	exponent = i->first;
98	tim	1746	coefficient = i->second;
99			result += pow(x, exponent - 1) * coefficient * exponent;
100	tim	1744	}
101
102			return result;
103			}
104
105	tim	1746	/**
106	tim	1747	* Set the coefficent of the specified exponent, if the coefficient is already there, it
107			* will be overwritten.
108	tim	1746	* @param exponent exponent of a term in this Polynomial
109			* @param coefficient multiplier of a term in this Polynomial
110			*/
111
112			void setCoefficient(int exponent, const ElemType& coefficient) {
113			polyPairMap_.insert(PolynomialPairMap::value_type(exponent, coefficient));
114	tim	1744	}
115
116	tim	1746	/**
117	tim	1747	* Set the coefficent of the specified exponent. If the coefficient is already there, just add the
118			* new coefficient to the old one, otherwise, just call setCoefficent
119			* @param exponent exponent of a term in this Polynomial
120			* @param coefficient multiplier of a term in this Polynomial
121			*/
122
123			void addCoefficient(int exponent, const ElemType& coefficient) {
124	tim	1796	iterator i = polyPairMap_.find(exponent);
125	tim	1747
126			if (i != end()) {
127			i->second += coefficient;
128			} else {
129			setCoefficient(exponent, coefficient);
130			}
131			}
132
133
134			/**
135	tim	1746	* Returns the coefficient associated with the given power for this Polynomial.
136			* @return the coefficient associated with the given power for this Polynomial
137			* @exponent exponent of any term in this Polynomial
138			*/
139			ElemType getCoefficient(ExponentType exponent) {
140	tim	1796	iterator i = polyPairMap_.find(exponent);
141	tim	1744
142	tim	1746	if (i != end()) {
143			return i->second;
144			} else {
145			return ElemType(0);
146			}
147	tim	1744	}
148
149	tim	1796	iterator begin() {
150	tim	1744	return polyPairMap_.begin();
151			}
152	tim	1796
153			const_iterator begin() const{
154			return polyPairMap_.begin();
155			}
156	tim	1744
157	tim	1796	iterator end() {
158	tim	1744	return polyPairMap_.end();
159			}
160
161	tim	1796	const_iterator end() const{
162			return polyPairMap_.end();
163	tim	1744	}
164	tim	1746
165	tim	1796	iterator find(ExponentType exponent) {
166			return polyPairMap_.find(exponent);
167			}
168
169	tim	1746	size_t size() {
170			return polyPairMap_.size();
171			}
172
173	tim	1744	private:
174
175			PolynomialPairMap polyPairMap_;
176			};
177
178	tim	1746
179			/**
180			* Generates and returns the product of two given Polynomials.
181			* @return A Polynomial containing the product of the two given Polynomial parameters
182			*/
183			template<typename ElemType>
184			Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
185	tim	1796	typename Polynomial<ElemType>::const_iterator i;
186			typename Polynomial<ElemType>::const_iterator j;
187	tim	1746	Polynomial<ElemType> p;
188
189			for (i = p1.begin(); i !=p1.end(); ++i) {
190	tim	1796	for (j = p2.begin(); j !=p2.end(); ++j) {
191			p.addCoefficient( i->first + j->first, i->second * j->second);
192	tim	1746	}
193			}
194	tim	1796
195			return p;
196	tim	1746	}
197
198			/**
199			* Generates and returns the sum of two given Polynomials.
200			* @param p1 the first polynomial
201			* @param p2 the second polynomial
202			*/
203			template<typename ElemType>
204			Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
205			Polynomial<ElemType> p(p1);
206
207	tim	1796	typename Polynomial<ElemType>::const_iterator i;
208	tim	1746
209			for (i = p2.begin(); i != p2.end(); ++i) {
210	tim	1796	p.addCoefficient(i->first, i->second);
211	tim	1746	}
212
213			return p;
214
215			}
216
217			/**
218			* Generates and returns the difference of two given Polynomials.
219			* @return
220			* @param p1 the first polynomial
221			* @param p2 the second polynomial
222			*/
223			template<typename ElemType>
224			Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
225			Polynomial<ElemType> p(p1);
226
227	tim	1796	typename Polynomial<ElemType>::const_iterator i;
228	tim	1746
229			for (i = p2.begin(); i != p2.end(); ++i) {
230	tim	1796	p.addCoefficient(i->first, -i->second);
231	tim	1746	}
232
233			return p;
234
235			}
236
237			/**
238			* Tests if two polynomial have the same exponents
239			* @return true if these all of the exponents in these Polynomial are identical
240			* @param p1 the first polynomial
241			* @param p2 the second polynomial
242			* @note this function does not compare the coefficient
243			*/
244			template<typename ElemType>
245			bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
246
247	tim	1796	typename Polynomial<ElemType>::const_iterator i;
248			typename Polynomial<ElemType>::const_iterator j;
249	tim	1746
250			if (p1.size() !== p2.size() ) {
251			return false;
252			}
253
254			for (i = p1.begin(), j = p2.begin(); i != p1.end() && j != p2.end(); ++i, ++j) {
255			if (i->first != j->first) {
256			return false;
257			}
258			}
259
260			return true;
261			}
262
263			typedef Polynomial<double> DoublePolynomial;
264
265	tim	1744	} //end namespace oopse
266			#endif //MATH_POLYNOMIAL_HPP