ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/branches/development/src/math/Polynomial.hpp

Comparing:
trunk/src/math/Polynomial.hpp (file contents), Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
branches/development/src/math/Polynomial.hpp (file contents), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC

#	Line 6 | Line 6
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
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
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.
#	Line 37 | Line 28
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		/**
#	Line 53 | Line 54
54		#include <list>
55		#include <map>
56		#include <utility>
57	+	#include <complex>
58		#include "config.h"
59	<	namespace oopse {
59	>	#include "math/Eigenvalue.hpp"
60
61	<	template<typename ElemType> ElemType pow(ElemType x, int N) {
62	<	ElemType result(1);
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;
#	Line 70 | Line 74 | namespace oopse {
74		* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
75		* A generic Polynomial class
76		*/
77	<	template<typename ElemType>
77	>	template<typename Real>
78		class Polynomial {
79
80		public:
81	<	typedef Polynomial<ElemType> PolynomialType;
81	>	typedef Polynomial<Real> PolynomialType;
82		typedef int ExponentType;
83	<	typedef ElemType CoefficientType;
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(ElemType v) {setCoefficient(0, v);}
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 Polynomial function
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	<	ElemType evaluate(const ElemType& x) {
97	<	ElemType result = ElemType();
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  += pow(x, exponent) * coefficient;
104	>	result  += fastpow(x, exponent) * coefficient;
105		}
106
107		return result;
#	Line 107 | Line 112 | namespace oopse {
112		* @return the first derivative of this polynomial
113		* @param x
114		*/
115	<	ElemType evaluateDerivative(const ElemType& x) {
116	<	ElemType result = ElemType();
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  += pow(x, exponent - 1) * coefficient * exponent;
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 coefficient is already there, it
132	<	* will be overwritten.
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	<
131	<	void setCoefficient(int exponent, const ElemType& coefficient) {
135	>	*/
136	>	void setCoefficient(int exponent, const Real& coefficient) {
137		polyPairMap_[exponent] = coefficient;
138		}
139	<
139	>
140		/**
141	<	* Set the coefficent of the specified exponent. If the coefficient is already there,  just add the
142	<	* new coefficient to the old one, otherwise,  just call setCoefficent
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	<
142	<	void addCoefficient(int exponent, const ElemType& coefficient) {
146	>	*/
147	>	void addCoefficient(int exponent, const Real& coefficient) {
148		iterator i = polyPairMap_.find(exponent);
149
150		if (i != end()) {
#	Line 150 | Line 155 | namespace oopse {
155		}
156
157		/**
158	<	* Returns the coefficient associated with the given power for this Polynomial.
159	<	* @return the coefficient associated with the given power for this Polynomial
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		* @exponent exponent of any term in this Polynomial
163		*/
164	<	ElemType getCoefficient(ExponentType exponent) {
164	>	Real getCoefficient(ExponentType exponent) {
165		iterator i = polyPairMap_.find(exponent);
166
167		if (i != end()) {
168		return i->second;
169		} else {
170	<	return ElemType(0);
170	>	return Real(0);
171		}
172		}
173
#	Line 188 | Line 195 | namespace oopse {
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<ElemType>::const_iterator i;
210	>	{
211	>	typename Polynomial<Real>::const_iterator i;
212
213	<	polyPairMap_.clear();  // clear out the old map
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);
215	>	for (i =  p.begin(); i != p.end(); ++i) {
216	>	this->setCoefficient(i->first, i->second);
217	>	}
218		}
202	–	}
219		// by convention, always return *this
220		return *this;
221		}
222
223		PolynomialType& operator += (const PolynomialType& p) {
224	<	typename Polynomial<ElemType>::const_iterator i;
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	<	}
226	>	for (i =  p.begin(); i  != p.end(); ++i) {
227	>	this->addCoefficient(i->first, i->second);
228	>	}
229
230	<	return *this;
230	>	return *this;
231		}
232
233		PolynomialType& operator -= (const PolynomialType& p) {
234	<	typename Polynomial<ElemType>::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<ElemType>::const_iterator i;
243	<	typename Polynomial<ElemType>::const_iterator j;
244	<	Polynomial<ElemType> 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);
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	<	return *this;
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 ElemType v)
268	<	PolynomialType& operator *= (const ElemType v) {
269	<	typename Polynomial<ElemType>::const_iterator i;
270	<	//Polynomial<ElemType> result;
267	>	PolynomialType& operator += (const Real v) {
268	>	this->addCoefficient( 0, v);
269	>	return *this;
270	>	}
271
272	<	for (i = this->begin(); i != this->end(); ++i) {
273	<	this->setCoefficient( i->first, i->second*v);
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	<	return *this;
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	<	PolynomialType& operator += (const ElemType v) {
252	<	this->addCoefficient( 0, v);
253	<	return *this;
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	+	break;
342	+	case 2: {
343	+	Real fC2 = getCoefficient(2);
344	+	Real fC1 = getCoefficient(1);
345	+	Real fC0 = getCoefficient(0);
346	+	Real fDiscr = fC1*fC1 - 4.0*fC0*fC2;
347	+	if (abs(fDiscr) <= fEpsilon) {
348	+	fDiscr = (Real)0.0;
349	+	}
350	+
351	+	if (fDiscr < (Real)0.0) {  // complex roots only
352	+	return roots;
353	+	}
354	+
355	+	Real fTmp = ((Real)0.5)/fC2;
356	+
357	+	if (fDiscr > (Real)0.0) { // 2 real roots
358	+	fDiscr = sqrt(fDiscr);
359	+	roots.push_back(fTmp*(-fC1 - fDiscr));
360	+	roots.push_back(fTmp*(-fC1 + fDiscr));
361	+	} else {
362	+	roots.push_back(-fTmp * fC1);  // 1 real root
363	+	}
364	+	}
365	+	return roots;
366	+	break;
367	+
368	+	case 3: {
369	+	Real fC3 = getCoefficient(3);
370	+	Real fC2 = getCoefficient(2);
371	+	Real fC1 = getCoefficient(1);
372	+	Real fC0 = getCoefficient(0);
373	+
374	+	// make polynomial monic, x^3+c2*x^2+c1*x+c0
375	+	Real fInvC3 = ((Real)1.0)/fC3;
376	+	fC0 *= fInvC3;
377	+	fC1 *= fInvC3;
378	+	fC2 *= fInvC3;
379	+
380	+	// convert to y^3+a*y+b = 0 by x = y-c2/3
381	+	const Real fThird = (Real)1.0/(Real)3.0;
382	+	const Real fTwentySeventh = (Real)1.0/(Real)27.0;
383	+	Real fOffset = fThird*fC2;
384	+	Real fA = fC1 - fC2*fOffset;
385	+	Real fB = fC0+fC2*(((Real)2.0)*fC2*fC2-((Real)9.0)*fC1)*fTwentySeventh;
386	+	Real fHalfB = ((Real)0.5)*fB;
387	+
388	+	Real fDiscr = fHalfB*fHalfB + fA*fA*fA*fTwentySeventh;
389	+	if (fabs(fDiscr) <= fEpsilon) {
390	+	fDiscr = (Real)0.0;
391	+	}
392	+
393	+	if (fDiscr > (Real)0.0) {  // 1 real, 2 complex roots
394	+
395	+	fDiscr = sqrt(fDiscr);
396	+	Real fTemp = -fHalfB + fDiscr;
397	+	Real root;
398	+	if (fTemp >= (Real)0.0) {
399	+	root = pow(fTemp,fThird);
400	+	} else {
401	+	root = -pow(-fTemp,fThird);
402	+	}
403	+	fTemp = -fHalfB - fDiscr;
404	+	if ( fTemp >= (Real)0.0 ) {
405	+	root += pow(fTemp,fThird);
406	+	} else {
407	+	root -= pow(-fTemp,fThird);
408	+	}
409	+	root -= fOffset;
410	+
411	+	roots.push_back(root);
412	+	} else if (fDiscr < (Real)0.0) {
413	+	const Real fSqrt3 = sqrt((Real)3.0);
414	+	Real fDist = sqrt(-fThird*fA);
415	+	Real fAngle = fThird*atan2(sqrt(-fDiscr), -fHalfB);
416	+	Real fCos = cos(fAngle);
417	+	Real fSin = sin(fAngle);
418	+	roots.push_back(((Real)2.0)*fDist*fCos-fOffset);
419	+	roots.push_back(-fDist*(fCos+fSqrt3*fSin)-fOffset);
420	+	roots.push_back(-fDist*(fCos-fSqrt3*fSin)-fOffset);
421	+	} else {
422	+	Real fTemp;
423	+	if (fHalfB >= (Real)0.0) {
424	+	fTemp = -pow(fHalfB,fThird);
425	+	} else {
426	+	fTemp = pow(-fHalfB,fThird);
427	+	}
428	+	roots.push_back(((Real)2.0)*fTemp-fOffset);
429	+	roots.push_back(-fTemp-fOffset);
430	+	roots.push_back(-fTemp-fOffset);
431	+	}
432	+	}
433	+	return roots;
434	+	break;
435	+	case 4: {
436	+	Real fC4 = getCoefficient(4);
437	+	Real fC3 = getCoefficient(3);
438	+	Real fC2 = getCoefficient(2);
439	+	Real fC1 = getCoefficient(1);
440	+	Real fC0 = getCoefficient(0);
441	+
442	+	// make polynomial monic, x^4+c3*x^3+c2*x^2+c1*x+c0
443	+	Real fInvC4 = ((Real)1.0)/fC4;
444	+	fC0 *= fInvC4;
445	+	fC1 *= fInvC4;
446	+	fC2 *= fInvC4;
447	+	fC3 *= fInvC4;
448	+
449	+	// reduction to resolvent cubic polynomial y^3+r2*y^2+r1*y+r0 = 0
450	+	Real fR0 = -fC3*fC3*fC0 + ((Real)4.0)*fC2*fC0 - fC1*fC1;
451	+	Real fR1 = fC3*fC1 - ((Real)4.0)*fC0;
452	+	Real fR2 = -fC2;
453	+	Polynomial<Real> tempCubic;
454	+	tempCubic.setCoefficient(0, fR0);
455	+	tempCubic.setCoefficient(1, fR1);
456	+	tempCubic.setCoefficient(2, fR2);
457	+	tempCubic.setCoefficient(3, 1.0);
458	+	std::vector<Real> cubeRoots = tempCubic.FindRealRoots(); // always
459	+	// produces
460	+	// at
461	+	// least
462	+	// one
463	+	// root
464	+	Real fY = cubeRoots[0];
465	+
466	+	Real fDiscr = ((Real)0.25)*fC3*fC3 - fC2 + fY;
467	+	if (fabs(fDiscr) <= fEpsilon) {
468	+	fDiscr = (Real)0.0;
469	+	}
470
471	+	if (fDiscr > (Real)0.0) {
472	+	Real fR = sqrt(fDiscr);
473	+	Real fT1 = ((Real)0.75)*fC3*fC3 - fR*fR - ((Real)2.0)*fC2;
474	+	Real fT2 = (((Real)4.0)*fC3*fC2 - ((Real)8.0)*fC1 - fC3*fC3*fC3) /
475	+	(((Real)4.0)*fR);
476	+
477	+	Real fTplus = fT1+fT2;
478	+	Real fTminus = fT1-fT2;
479	+	if (fabs(fTplus) <= fEpsilon) {
480	+	fTplus = (Real)0.0;
481	+	}
482	+	if (fabs(fTminus) <= fEpsilon) {
483	+	fTminus = (Real)0.0;
484	+	}
485	+
486	+	if (fTplus >= (Real)0.0) {
487	+	Real fD = sqrt(fTplus);
488	+	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR+fD));
489	+	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR-fD));
490	+	}
491	+	if (fTminus >= (Real)0.0) {
492	+	Real fE = sqrt(fTminus);
493	+	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fE-fR));
494	+	roots.push_back(-((Real)0.25)*fC3-((Real)0.5)*(fE+fR));
495	+	}
496	+	} else if (fDiscr < (Real)0.0) {
497	+	//roots.clear();
498	+	} else {
499	+	Real fT2 = fY*fY-((Real)4.0)*fC0;
500	+	if (fT2 >= -fEpsilon) {
501	+	if (fT2 < (Real)0.0) { // round to zero
502	+	fT2 = (Real)0.0;
503	+	}
504	+	fT2 = ((Real)2.0)*sqrt(fT2);
505	+	Real fT1 = ((Real)0.75)*fC3*fC3 - ((Real)2.0)*fC2;
506	+	if (fT1+fT2 >= fEpsilon) {
507	+	Real fD = sqrt(fT1+fT2);
508	+	roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fD);
509	+	roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fD);
510	+	}
511	+	if (fT1-fT2 >= fEpsilon) {
512	+	Real fE = sqrt(fT1-fT2);
513	+	roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fE);
514	+	roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fE);
515	+	}
516	+	}
517	+	}
518	+	}
519	+	return roots;
520	+	break;
521	+	default: {
522	+	DynamicRectMatrix<Real> companion = CreateCompanion();
523	+	JAMA::Eigenvalue<Real> eig(companion);
524	+	DynamicVector<Real> reals, imags;
525	+	eig.getRealEigenvalues(reals);
526	+	eig.getImagEigenvalues(imags);
527	+
528	+	for (int i = 0; i < deg; i++) {
529	+	if (fabs(imags(i)) < fEpsilon)
530	+	roots.push_back(reals(i));
531	+	}
532	+	}
533	+	return roots;
534	+	break;
535	+	}
536	+
537	+	return roots; // should be empty if you got here
538	+	}
539	+
540		private:
541
542		PolynomialPairMap polyPairMap_;
543		};
544
545	<
545	>
546		/**
547		* Generates and returns the product of two given Polynomials.
548		* @return A Polynomial containing the product of the two given Polynomial parameters
549		*/
550	<	template<typename ElemType>
551	<	Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
552	<	typename Polynomial<ElemType>::const_iterator i;
553	<	typename Polynomial<ElemType>::const_iterator j;
554	<	Polynomial<ElemType> p;
550	>	template<typename Real>
551	>	Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
552	>	typename Polynomial<Real>::const_iterator i;
553	>	typename Polynomial<Real>::const_iterator j;
554	>	Polynomial<Real> p;
555
556		for (i = p1.begin(); i !=p1.end(); ++i) {
557		for (j = p2.begin(); j !=p2.end(); ++j) {
#	Line 278 | Line 562 | namespace oopse {
562		return p;
563		}
564
565	<	template<typename ElemType>
566	<	Polynomial<ElemType> operator *(const Polynomial<ElemType>& p, const ElemType v) {
567	<	typename Polynomial<ElemType>::const_iterator i;
568	<	Polynomial<ElemType> result;
565	>	template<typename Real>
566	>	Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
567	>	typename Polynomial<Real>::const_iterator i;
568	>	Polynomial<Real> result;
569
570		for (i = p.begin(); i !=p.end(); ++i) {
571		result.setCoefficient( i->first , i->second * v);
#	Line 290 | Line 574 | namespace oopse {
574		return result;
575		}
576
577	<	template<typename ElemType>
578	<	Polynomial<ElemType> operator *( const ElemType v, const Polynomial<ElemType>& p) {
579	<	typename Polynomial<ElemType>::const_iterator i;
580	<	Polynomial<ElemType> result;
577	>	template<typename Real>
578	>	Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
579	>	typename Polynomial<Real>::const_iterator i;
580	>	Polynomial<Real> result;
581
582		for (i = p.begin(); i !=p.end(); ++i) {
583		result.setCoefficient( i->first , i->second * v);
#	Line 307 | Line 591 | namespace oopse {
591		* @param p1 the first polynomial
592		* @param p2 the second polynomial
593		*/
594	<	template<typename ElemType>
595	<	Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
596	<	Polynomial<ElemType> p(p1);
594	>	template<typename Real>
595	>	Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
596	>	Polynomial<Real> p(p1);
597
598	<	typename Polynomial<ElemType>::const_iterator i;
598	>	typename Polynomial<Real>::const_iterator i;
599
600		for (i =  p2.begin(); i  != p2.end(); ++i) {
601		p.addCoefficient(i->first, i->second);
#	Line 327 | Line 611 | namespace oopse {
611		* @param p1 the first polynomial
612		* @param p2 the second polynomial
613		*/
614	<	template<typename ElemType>
615	<	Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
616	<	Polynomial<ElemType> p(p1);
614	>	template<typename Real>
615	>	Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
616	>	Polynomial<Real> p(p1);
617
618	<	typename Polynomial<ElemType>::const_iterator i;
618	>	typename Polynomial<Real>::const_iterator i;
619
620		for (i =  p2.begin(); i  != p2.end(); ++i) {
621		p.addCoefficient(i->first, -i->second);
#	Line 342 | Line 626 | namespace oopse {
626		}
627
628		/**
629	+	* Returns the first derivative of this polynomial.
630	+	* @return the first derivative of this polynomial
631	+	*/
632	+	template<typename Real>
633	+	Polynomial<Real> getDerivative(const Polynomial<Real>& p1) {
634	+	Polynomial<Real> p;
635	+
636	+	typename Polynomial<Real>::const_iterator i;
637	+	int exponent;
638	+	Real coefficient;
639	+
640	+	for (i =  p1.begin(); i  != p1.end(); ++i) {
641	+	exponent = i->first;
642	+	coefficient = i->second;
643	+	p.setCoefficient(exponent-1, coefficient * exponent);
644	+	}
645	+
646	+	return p;
647	+	}
648	+
649	+	/**
650		* Tests if two polynomial have the same exponents
651		* @return true if all of the exponents in these Polynomial are identical
652		* @param p1 the first polynomial
653		* @param p2 the second polynomial
654		* @note this function does not compare the coefficient
655		*/
656	<	template<typename ElemType>
657	<	bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
656	>	template<typename Real>
657	>	bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
658
659	<	typename Polynomial<ElemType>::const_iterator i;
660	<	typename Polynomial<ElemType>::const_iterator j;
659	>	typename Polynomial<Real>::const_iterator i;
660	>	typename Polynomial<Real>::const_iterator j;
661
662		if (p1.size() != p2.size() ) {
663		return false;
#	Line 367 | Line 672 | namespace oopse {
672		return true;
673		}
674
675	+
676	+
677		typedef Polynomial<RealType> DoublePolynomial;
678
679	<	} //end namespace oopse
679	>	} //end namespace OpenMD
680		#endif //MATH_POLYNOMIAL_HPP

Comparing:
trunk/src/math/Polynomial.hpp (property svn:keywords), Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
branches/development/src/math/Polynomial.hpp (property svn:keywords), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines