ViewVC Help

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

root/OpenMD/trunk/src/math/Polynomial.hpp

Comparing trunk/src/math/Polynomial.hpp (file contents):
Revision 877 by gezelter, Wed Feb 1 17:26:43 2006 UTC vs.
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 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]  Vardeman & Gezelter, in progress (2009).
40		*/
41
42		/**
#	Line 53 | Line 53
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 {
60	>	namespace OpenMD {
61	>
62	>	template<typename Real> Real fastpow(Real x, int N) {
63	>	Real result(1); //or 1.0?
64
59	–	template<typename ElemType> ElemType pow(ElemType x, int N) {
60	–	ElemType result(1);
61	–
65		for (int i = 0; i < N; ++i) {
66		result *= x;
67		}
#	Line 70 | Line 73 | namespace oopse {
73		* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
74		* A generic Polynomial class
75		*/
76	<	template<typename ElemType>
76	>	template<typename Real>
77		class Polynomial {
78
79		public:
80	<	typedef Polynomial<ElemType> PolynomialType;
80	>	typedef Polynomial<Real> PolynomialType;
81		typedef int ExponentType;
82	<	typedef ElemType CoefficientType;
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(ElemType v) {setCoefficient(0, v);}
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 Polynomial function
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	<	ElemType evaluate(const ElemType& x) {
96	<	ElemType result = ElemType();
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  += pow(x, exponent) * coefficient;
103	>	result  += fastpow(x, exponent) * coefficient;
104		}
105
106		return result;
#	Line 107 | Line 111 | namespace oopse {
111		* @return the first derivative of this polynomial
112		* @param x
113		*/
114	<	ElemType evaluateDerivative(const ElemType& x) {
115	<	ElemType result = ElemType();
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  += pow(x, exponent - 1) * coefficient * exponent;
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 coefficient is already there, it
131	<	* will be overwritten.
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	<
136	<	void setCoefficient(int exponent, const ElemType& coefficient) {
132	<	polyPairMap_.insert(typename PolynomialPairMap::value_type(exponent, coefficient));
134	>	*/
135	>	void setCoefficient(int exponent, const Real& coefficient) {
136	>	polyPairMap_[exponent] = coefficient;
137		}
138	<
138	>
139		/**
140	<	* Set the coefficent of the specified exponent. If the coefficient is already there,  just add the
141	<	* new coefficient to the old one, otherwise,  just call setCoefficent
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	<
142	<	void addCoefficient(int exponent, const ElemType& coefficient) {
145	>	*/
146	>	void addCoefficient(int exponent, const Real& coefficient) {
147		iterator i = polyPairMap_.find(exponent);
148
149		if (i != end()) {
#	Line 149 | Line 153 | namespace oopse {
153		}
154		}
155
152	–
156		/**
157	<	* Returns the coefficient associated with the given power for this Polynomial.
158	<	* @return the coefficient associated with the given power for this Polynomial
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	<	ElemType getCoefficient(ExponentType exponent) {
163	>	Real getCoefficient(ExponentType exponent) {
164		iterator i = polyPairMap_.find(exponent);
165
166		if (i != end()) {
167		return i->second;
168		} else {
169	<	return ElemType(0);
169	>	return Real(0);
170		}
171		}
172
#	Line 189 | Line 194 | namespace oopse {
194		return polyPairMap_.size();
195		}
196
197	<	PolynomialType& operator += (const PolynomialType& p) {
198	<	typename Polynomial<ElemType>::const_iterator i;
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	<	for (i =  p.begin(); i  != p.end(); ++i) {
196	<	this->addCoefficient(i->first, i->second);
197	<	}
206	>	PolynomialType& operator = (const PolynomialType& p) {
207
208	<	return *this;
209	<	}
208	>	if (this != &p)  // protect against invalid self-assignment
209	>	{
210	>	typename Polynomial<Real>::const_iterator i;
211
212	<	PolynomialType& operator -= (const PolynomialType& p) {
213	<	typename Polynomial<ElemType>::const_iterator i;
214	<	for (i =  p.begin(); i  != p.end(); ++i) {
215	<	this->addCoefficient(i->first, -i->second);
216	<	}
217	<	return *this;
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<ElemType>::const_iterator i;
224	<	typename Polynomial<ElemType>::const_iterator j;
225	<
226	<	for (i = this->begin(); i !=this->end(); ++i) {
215	<	for (j = p.begin(); j !=p.end(); ++j) {
216	<	this->addCoefficient( i->first + j->first, i->second * j->second);
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	<	return *this;
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 = fC1*fC1 - 4.0*fC0*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+c2*x^2+c1*x+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)*fC2*fC2-((Real)9.0)*fC1)*fTwentySeventh;
385	+	Real fHalfB = ((Real)0.5)*fB;
386	+
387	+	Real fDiscr = fHalfB*fHalfB + fA*fA*fA*fTwentySeventh;
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)*fDist*fCos-fOffset);
418	+	roots.push_back(-fDist*(fCos+fSqrt3*fSin)-fOffset);
419	+	roots.push_back(-fDist*(fCos-fSqrt3*fSin)-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+c3*x^3+c2*x^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+r2*y^2+r1*y+r0 = 0
449	+	Real fR0 = -fC3*fC3*fC0 + ((Real)4.0)*fC2*fC0 - fC1*fC1;
450	+	Real fR1 = fC3*fC1 - ((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)*fC3*fC3 - 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)*fC3*fC3 - fR*fR - ((Real)2.0)*fC2;
473	+	Real fT2 = (((Real)4.0)*fC3*fC2 - ((Real)8.0)*fC1 - fC3*fC3*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 = fY*fY-((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)*fC3*fC3 - ((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	<
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 ElemType>
550	<	Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
551	<	typename Polynomial<ElemType>::const_iterator i;
552	<	typename Polynomial<ElemType>::const_iterator j;
553	<	Polynomial<ElemType> p;
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) {
#	Line 246 | Line 561 | namespace oopse {
561		return p;
562		}
563
564	<	template<typename ElemType>
565	<	Polynomial<ElemType> operator *(const Polynomial<ElemType>& p, const ElemType v) {
566	<	typename Polynomial<ElemType>::const_iterator i;
567	<	Polynomial<ElemType> result;
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.addCoefficient( i->first , i->second * v);
570	>	result.setCoefficient( i->first , i->second * v);
571		}
572
573		return result;
574		}
575
576	<	template<typename ElemType>
577	<	Polynomial<ElemType> operator *( const ElemType v, const Polynomial<ElemType>& p) {
578	<	typename Polynomial<ElemType>::const_iterator i;
579	<	Polynomial<ElemType> result;
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.addCoefficient( i->first , i->second * v);
582	>	result.setCoefficient( i->first , i->second * v);
583		}
584
585		return result;
#	Line 275 | Line 590 | namespace oopse {
590		* @param p1 the first polynomial
591		* @param p2 the second polynomial
592		*/
593	<	template<typename ElemType>
594	<	Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
595	<	Polynomial<ElemType> p(p1);
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<ElemType>::const_iterator i;
597	>	typename Polynomial<Real>::const_iterator i;
598
599		for (i =  p2.begin(); i  != p2.end(); ++i) {
600		p.addCoefficient(i->first, i->second);
#	Line 295 | Line 610 | namespace oopse {
610		* @param p1 the first polynomial
611		* @param p2 the second polynomial
612		*/
613	<	template<typename ElemType>
614	<	Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
615	<	Polynomial<ElemType> p(p1);
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<ElemType>::const_iterator i;
617	>	typename Polynomial<Real>::const_iterator i;
618
619		for (i =  p2.begin(); i  != p2.end(); ++i) {
620		p.addCoefficient(i->first, -i->second);
#	Line 310 | Line 625 | namespace oopse {
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 these all of the exponents in these Polynomial are identical
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 ElemType>
656	<	bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
655	>	template<typename Real>
656	>	bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
657
658	<	typename Polynomial<ElemType>::const_iterator i;
659	<	typename Polynomial<ElemType>::const_iterator j;
658	>	typename Polynomial<Real>::const_iterator i;
659	>	typename Polynomial<Real>::const_iterator j;
660
661		if (p1.size() != p2.size() ) {
662		return false;
#	Line 335 | Line 671 | namespace oopse {
671		return true;
672		}
673
338	–	typedef Polynomial<double> DoublePolynomial;
674
675	<	} //end namespace oopse
675	>
676	>	typedef Polynomial<RealType> DoublePolynomial;
677	>
678	>	} //end namespace OpenMD
679		#endif //MATH_POLYNOMIAL_HPP

Comparing trunk/src/math/Polynomial.hpp (property svn:keywords):
Revision 877 by gezelter, Wed Feb 1 17:26:43 2006 UTC vs.
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC

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

Diff Legend

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