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 385 by tim, Tue Mar 1 20:10:14 2005 UTC vs.
branches/development/src/math/Polynomial.hpp (file contents), Revision 1827 by gezelter, Wed Jan 9 19:49:07 2013 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	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	+	#include "math/Eigenvalue.hpp"
60
61	<	namespace oopse {
61	>	namespace OpenMD {
62	>
63	>	template<typename Real> Real fastpow(Real x, int N) {
64	>	Real result(1); //or 1.0?
65
59	–	template<typename ElemType> ElemType pow(ElemType x, int N) {
60	–	ElemType result(1);
61	–
66		for (int i = 0; i < N; ++i) {
67	<	result *= x;
67	>	result *= x;
68		}
69
70		return result;
71	<	}
71	>	}
72
73	<	/**
74	<	* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
75	<	* A generic Polynomial class
76	<	*/
77	<	template<typename ElemType>
78	<	class Polynomial {
73	>	/**
74	>	* @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
75	>	* A generic Polynomial class
76	>	*/
77	>	template<typename Real>
78	>	class Polynomial {
79
80	<	public:
81	<
82	<	typedef int ExponentType;
83	<	typedef ElemType CoefficientType;
84	<	typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
85	<	typedef typename PolynomialPairMap::iterator iterator;
86	<	typedef typename PolynomialPairMap::const_iterator const_iterator;
87	<	/**
88	<	* Calculates the value of this Polynomial evaluated at the given x value.
89	<	* @return The value of this Polynomial evaluates at the given x value
90	<	* @param x the value of the independent variable for this Polynomial function
91	<	*/
92	<	ElemType evaluate(const ElemType& x) {
93	<	ElemType result = ElemType();
94	<	ExponentType exponent;
95	<	CoefficientType coefficient;
80	>	public:
81	>	typedef Polynomial<Real> PolynomialType;
82	>	typedef int ExponentType;
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(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
94	>	* Polynomial function
95	>	*/
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;
105	<	}
101	>	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
102	>	exponent = i->first;
103	>	coefficient = i->second;
104	>	result  += fastpow(x, exponent) * coefficient;
105	>	}
106
107	<	return result;
108	<	}
107	>	return result;
108	>	}
109
110	<	/**
111	<	* Returns the first derivative of this polynomial.
112	<	* @return the first derivative of this polynomial
113	<	* @param x
114	<	*/
115	<	ElemType evaluateDerivative(const ElemType& x) {
116	<	ElemType result = ElemType();
117	<	ExponentType exponent;
118	<	CoefficientType coefficient;
110	>	/**
111	>	* Returns the first derivative of this polynomial.
112	>	* @return the first derivative of this polynomial
113	>	* @param x
114	>	*/
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;
124	<	}
120	>	for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
121	>	exponent = i->first;
122	>	coefficient = i->second;
123	>	result  += fastpow(x, exponent - 1) * coefficient * exponent;
124	>	}
125
126	<	return result;
127	<	}
126	>	return result;
127	>	}
128
121	–	/**
122	–	* Set the coefficent of the specified exponent, if the coefficient is already there, it
123	–	* will be overwritten.
124	–	* @param exponent exponent of a term in this Polynomial
125	–	* @param coefficient multiplier of a term in this Polynomial
126	–	*/
127	–
128	–	void setCoefficient(int exponent, const ElemType& coefficient) {
129	–	polyPairMap_.insert(typename PolynomialPairMap::value_type(exponent, coefficient));
130	–	}
129
130	<	/**
131	<	* Set the coefficent of the specified exponent. If the coefficient is already there,  just add the
132	<	* new coefficient to the old one, otherwise,  just call setCoefficent
133	<	* @param exponent exponent of a term in this Polynomial
134	<	* @param coefficient multiplier of a term in this Polynomial
135	<	*/
136	<
137	<	void addCoefficient(int exponent, const ElemType& coefficient) {
138	<	iterator i = polyPairMap_.find(exponent);
130	>	/**
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	>	void setCoefficient(int exponent, const Real& coefficient) {
137	>	polyPairMap_[exponent] = coefficient;
138	>	}
139	>
140	>	/**
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	>	void addCoefficient(int exponent, const Real& coefficient) {
148	>	iterator i = polyPairMap_.find(exponent);
149
150	<	if (i != end()) {
151	<	i->second += coefficient;
152	<	} else {
153	<	setCoefficient(exponent, coefficient);
154	<	}
155	<	}
150	>	if (i != end()) {
151	>	i->second += coefficient;
152	>	} else {
153	>	setCoefficient(exponent, coefficient);
154	>	}
155	>	}
156
157	<
158	<	/**
159	<	* Returns the coefficient associated with the given power for this Polynomial.
160	<	* @return the coefficient associated with the given power for this Polynomial
161	<	* @exponent exponent of any term in this Polynomial
162	<	*/
163	<	ElemType getCoefficient(ExponentType exponent) {
164	<	iterator i = polyPairMap_.find(exponent);
157	>	/**
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	>	* @param exponent exponent of any term in this Polynomial
163	>	*/
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);
171	<	}
172	<	}
167	>	if (i != end()) {
168	>	return i->second;
169	>	} else {
170	>	return Real(0);
171	>	}
172	>	}
173
174	<	iterator begin() {
175	<	return polyPairMap_.begin();
176	<	}
174	>	iterator begin() {
175	>	return polyPairMap_.begin();
176	>	}
177
178	<	const_iterator begin() const{
179	<	return polyPairMap_.begin();
180	<	}
178	>	const_iterator begin() const{
179	>	return polyPairMap_.begin();
180	>	}
181
182	<	iterator end() {
183	<	return polyPairMap_.end();
184	<	}
182	>	iterator end() {
183	>	return polyPairMap_.end();
184	>	}
185
186	<	const_iterator end() const{
187	<	return polyPairMap_.end();
186	>	const_iterator end() const{
187	>	return polyPairMap_.end();
188	>	}
189	>
190	>	iterator find(ExponentType exponent) {
191	>	return polyPairMap_.find(exponent);
192	>	}
193	>
194	>	size_t size() {
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<Real>::const_iterator i;
212	>
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);
217	>	}
218		}
219	+	// by convention, always return *this
220	+	return *this;
221	+	}
222
223	<	iterator find(ExponentType exponent) {
224	<	return polyPairMap_.find(exponent);
223	>	PolynomialType& operator += (const PolynomialType& p) {
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	>	}
229	>
230	>	return *this;
231	>	}
232	>
233	>	PolynomialType& operator -= (const PolynomialType& p) {
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
255	<	size_t size() {
256	<	return polyPairMap_.size();
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 Real v) {
268	>	this->addCoefficient( 0, v);
269	>	return *this;
270	>	}
271	>
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	>	// 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	<	private:
307	<
308	<	PolynomialPairMap polyPairMap_;
309	<	};
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	+	return roots;
324	+	}
325
326	<	/**
327	<	* Generates and returns the product of two given Polynomials.
328	<	* @return A Polynomial containing the product of the two given Polynomial parameters
329	<	*/
330	<	template<typename ElemType>
331	<	Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
332	<	typename Polynomial<ElemType>::const_iterator i;
333	<	typename Polynomial<ElemType>::const_iterator j;
334	<	Polynomial<ElemType> p;
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	>	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	>	case 3: {
366	>	Real fC3 = getCoefficient(3);
367	>	Real fC2 = getCoefficient(2);
368	>	Real fC1 = getCoefficient(1);
369	>	Real fC0 = getCoefficient(0);
370	>
371	>	// make polynomial monic, x^3+c2*x^2+c1*x+c0
372	>	Real fInvC3 = ((Real)1.0)/fC3;
373	>	fC0 *= fInvC3;
374	>	fC1 *= fInvC3;
375	>	fC2 *= fInvC3;
376	>
377	>	// convert to y^3+a*y+b = 0 by x = y-c2/3
378	>	const Real fThird = (Real)1.0/(Real)3.0;
379	>	const Real fTwentySeventh = (Real)1.0/(Real)27.0;
380	>	Real fOffset = fThird*fC2;
381	>	Real fA = fC1 - fC2*fOffset;
382	>	Real fB = fC0+fC2*(((Real)2.0)*fC2*fC2-((Real)9.0)*fC1)*fTwentySeventh;
383	>	Real fHalfB = ((Real)0.5)*fB;
384	>
385	>	Real fDiscr = fHalfB*fHalfB + fA*fA*fA*fTwentySeventh;
386	>	if (fabs(fDiscr) <= fEpsilon) {
387	>	fDiscr = (Real)0.0;
388	>	}
389	>
390	>	if (fDiscr > (Real)0.0) {  // 1 real, 2 complex roots
391	>
392	>	fDiscr = sqrt(fDiscr);
393	>	Real fTemp = -fHalfB + fDiscr;
394	>	Real root;
395	>	if (fTemp >= (Real)0.0) {
396	>	root = pow(fTemp,fThird);
397	>	} else {
398	>	root = -pow(-fTemp,fThird);
399	>	}
400	>	fTemp = -fHalfB - fDiscr;
401	>	if ( fTemp >= (Real)0.0 ) {
402	>	root += pow(fTemp,fThird);
403	>	} else {
404	>	root -= pow(-fTemp,fThird);
405	>	}
406	>	root -= fOffset;
407	>
408	>	roots.push_back(root);
409	>	} else if (fDiscr < (Real)0.0) {
410	>	const Real fSqrt3 = sqrt((Real)3.0);
411	>	Real fDist = sqrt(-fThird*fA);
412	>	Real fAngle = fThird*atan2(sqrt(-fDiscr), -fHalfB);
413	>	Real fCos = cos(fAngle);
414	>	Real fSin = sin(fAngle);
415	>	roots.push_back(((Real)2.0)*fDist*fCos-fOffset);
416	>	roots.push_back(-fDist*(fCos+fSqrt3*fSin)-fOffset);
417	>	roots.push_back(-fDist*(fCos-fSqrt3*fSin)-fOffset);
418	>	} else {
419	>	Real fTemp;
420	>	if (fHalfB >= (Real)0.0) {
421	>	fTemp = -pow(fHalfB,fThird);
422	>	} else {
423	>	fTemp = pow(-fHalfB,fThird);
424	>	}
425	>	roots.push_back(((Real)2.0)*fTemp-fOffset);
426	>	roots.push_back(-fTemp-fOffset);
427	>	roots.push_back(-fTemp-fOffset);
428	>	}
429	>	}
430	>	return roots;
431	>	case 4: {
432	>	Real fC4 = getCoefficient(4);
433	>	Real fC3 = getCoefficient(3);
434	>	Real fC2 = getCoefficient(2);
435	>	Real fC1 = getCoefficient(1);
436	>	Real fC0 = getCoefficient(0);
437	>
438	>	// make polynomial monic, x^4+c3*x^3+c2*x^2+c1*x+c0
439	>	Real fInvC4 = ((Real)1.0)/fC4;
440	>	fC0 *= fInvC4;
441	>	fC1 *= fInvC4;
442	>	fC2 *= fInvC4;
443	>	fC3 *= fInvC4;
444	>
445	>	// reduction to resolvent cubic polynomial y^3+r2*y^2+r1*y+r0 = 0
446	>	Real fR0 = -fC3*fC3*fC0 + ((Real)4.0)*fC2*fC0 - fC1*fC1;
447	>	Real fR1 = fC3*fC1 - ((Real)4.0)*fC0;
448	>	Real fR2 = -fC2;
449	>	Polynomial<Real> tempCubic;
450	>	tempCubic.setCoefficient(0, fR0);
451	>	tempCubic.setCoefficient(1, fR1);
452	>	tempCubic.setCoefficient(2, fR2);
453	>	tempCubic.setCoefficient(3, 1.0);
454	>	std::vector<Real> cubeRoots = tempCubic.FindRealRoots(); // always
455	>	// produces
456	>	// at
457	>	// least
458	>	// one
459	>	// root
460	>	Real fY = cubeRoots[0];
461	>
462	>	Real fDiscr = ((Real)0.25)*fC3*fC3 - fC2 + fY;
463	>	if (fabs(fDiscr) <= fEpsilon) {
464	>	fDiscr = (Real)0.0;
465	>	}
466	>
467	>	if (fDiscr > (Real)0.0) {
468	>	Real fR = sqrt(fDiscr);
469	>	Real fT1 = ((Real)0.75)*fC3*fC3 - fR*fR - ((Real)2.0)*fC2;
470	>	Real fT2 = (((Real)4.0)*fC3*fC2 - ((Real)8.0)*fC1 - fC3*fC3*fC3) /
471	>	(((Real)4.0)*fR);
472	>
473	>	Real fTplus = fT1+fT2;
474	>	Real fTminus = fT1-fT2;
475	>	if (fabs(fTplus) <= fEpsilon) {
476	>	fTplus = (Real)0.0;
477	>	}
478	>	if (fabs(fTminus) <= fEpsilon) {
479	>	fTminus = (Real)0.0;
480	>	}
481	>
482	>	if (fTplus >= (Real)0.0) {
483	>	Real fD = sqrt(fTplus);
484	>	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR+fD));
485	>	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fR-fD));
486	>	}
487	>	if (fTminus >= (Real)0.0) {
488	>	Real fE = sqrt(fTminus);
489	>	roots.push_back(-((Real)0.25)*fC3+((Real)0.5)*(fE-fR));
490	>	roots.push_back(-((Real)0.25)*fC3-((Real)0.5)*(fE+fR));
491	>	}
492	>	} else if (fDiscr < (Real)0.0) {
493	>	//roots.clear();
494	>	} else {
495	>	Real fT2 = fY*fY-((Real)4.0)*fC0;
496	>	if (fT2 >= -fEpsilon) {
497	>	if (fT2 < (Real)0.0) { // round to zero
498	>	fT2 = (Real)0.0;
499	>	}
500	>	fT2 = ((Real)2.0)*sqrt(fT2);
501	>	Real fT1 = ((Real)0.75)*fC3*fC3 - ((Real)2.0)*fC2;
502	>	if (fT1+fT2 >= fEpsilon) {
503	>	Real fD = sqrt(fT1+fT2);
504	>	roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fD);
505	>	roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fD);
506	>	}
507	>	if (fT1-fT2 >= fEpsilon) {
508	>	Real fE = sqrt(fT1-fT2);
509	>	roots.push_back( -((Real)0.25)*fC3+((Real)0.5)*fE);
510	>	roots.push_back( -((Real)0.25)*fC3-((Real)0.5)*fE);
511	>	}
512	>	}
513	>	}
514	>	}
515	>	return roots;
516	>	default: {
517	>	DynamicRectMatrix<Real> companion = CreateCompanion();
518	>	JAMA::Eigenvalue<Real> eig(companion);
519	>	DynamicVector<Real> reals, imags;
520	>	eig.getRealEigenvalues(reals);
521	>	eig.getImagEigenvalues(imags);
522	>
523	>	for (int i = 0; i < deg; i++) {
524	>	if (fabs(imags(i)) < fEpsilon)
525	>	roots.push_back(reals(i));
526	>	}
527	>	}
528	>	return roots;
529	>	}
530	>
531	>	return roots; // should be empty if you got here
532	>	}
533
534	+	private:
535	+
536	+	PolynomialPairMap polyPairMap_;
537	+	};
538	+
539	+
540	+	/**
541	+	* Generates and returns the product of two given Polynomials.
542	+	* @return A Polynomial containing the product of the two given Polynomial parameters
543	+	*/
544	+	template<typename Real>
545	+	Polynomial<Real> operator *(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
546	+	typename Polynomial<Real>::const_iterator i;
547	+	typename Polynomial<Real>::const_iterator j;
548	+	Polynomial<Real> p;
549	+
550		for (i = p1.begin(); i !=p1.end(); ++i) {
551	<	for (j = p2.begin(); j !=p2.end(); ++j) {
552	<	p.addCoefficient( i->first + j->first, i->second * j->second);
553	<	}
551	>	for (j = p2.begin(); j !=p2.end(); ++j) {
552	>	p.addCoefficient( i->first + j->first, i->second * j->second);
553	>	}
554		}
555
556		return p;
557	<	}
557	>	}
558
559	<	/**
560	<	* Generates and returns the sum of two given Polynomials.
561	<	* @param p1 the first polynomial
562	<	* @param p2 the second polynomial
563	<	*/
564	<	template<typename ElemType>
565	<	Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
566	<	Polynomial<ElemType> p(p1);
559	>	template<typename Real>
560	>	Polynomial<Real> operator *(const Polynomial<Real>& p, const Real v) {
561	>	typename Polynomial<Real>::const_iterator i;
562	>	Polynomial<Real> result;
563	>
564	>	for (i = p.begin(); i !=p.end(); ++i) {
565	>	result.setCoefficient( i->first , i->second * v);
566	>	}
567
568	<	typename Polynomial<ElemType>::const_iterator i;
568	>	return result;
569	>	}
570
571	+	template<typename Real>
572	+	Polynomial<Real> operator *( const Real v, const Polynomial<Real>& p) {
573	+	typename Polynomial<Real>::const_iterator i;
574	+	Polynomial<Real> result;
575	+
576	+	for (i = p.begin(); i !=p.end(); ++i) {
577	+	result.setCoefficient( i->first , i->second * v);
578	+	}
579	+
580	+	return result;
581	+	}
582	+
583	+	/**
584	+	* Generates and returns the sum of two given Polynomials.
585	+	* @param p1 the first polynomial
586	+	* @param p2 the second polynomial
587	+	*/
588	+	template<typename Real>
589	+	Polynomial<Real> operator +(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
590	+	Polynomial<Real> p(p1);
591	+
592	+	typename Polynomial<Real>::const_iterator i;
593	+
594		for (i =  p2.begin(); i  != p2.end(); ++i) {
595	<	p.addCoefficient(i->first, i->second);
595	>	p.addCoefficient(i->first, i->second);
596		}
597
598		return p;
599
600	<	}
600	>	}
601
602	<	/**
603	<	* Generates and returns the difference of two given Polynomials.
604	<	* @return
605	<	* @param p1 the first polynomial
606	<	* @param p2 the second polynomial
607	<	*/
608	<	template<typename ElemType>
609	<	Polynomial<ElemType> operator -(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
610	<	Polynomial<ElemType> p(p1);
602	>	/**
603	>	* Generates and returns the difference of two given Polynomials.
604	>	* @return
605	>	* @param p1 the first polynomial
606	>	* @param p2 the second polynomial
607	>	*/
608	>	template<typename Real>
609	>	Polynomial<Real> operator -(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
610	>	Polynomial<Real> p(p1);
611
612	<	typename Polynomial<ElemType>::const_iterator i;
612	>	typename Polynomial<Real>::const_iterator i;
613
614		for (i =  p2.begin(); i  != p2.end(); ++i) {
615	<	p.addCoefficient(i->first, -i->second);
615	>	p.addCoefficient(i->first, -i->second);
616		}
617
618		return p;
619
620	<	}
620	>	}
621
622	<	/**
623	<	* Tests if two polynomial have the same exponents
624	<	* @return true if these all of the exponents in these Polynomial are identical
625	<	* @param p1 the first polynomial
626	<	* @param p2 the second polynomial
627	<	* @note this function does not compare the coefficient
628	<	*/
629	<	template<typename ElemType>
630	<	bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
622	>	/**
623	>	* Returns the first derivative of this polynomial.
624	>	* @return the first derivative of this polynomial
625	>	*/
626	>	template<typename Real>
627	>	Polynomial<Real> getDerivative(const Polynomial<Real>& p1) {
628	>	Polynomial<Real> p;
629	>
630	>	typename Polynomial<Real>::const_iterator i;
631	>	int exponent;
632	>	Real coefficient;
633	>
634	>	for (i =  p1.begin(); i  != p1.end(); ++i) {
635	>	exponent = i->first;
636	>	coefficient = i->second;
637	>	p.setCoefficient(exponent-1, coefficient * exponent);
638	>	}
639	>
640	>	return p;
641	>	}
642
643	<	typename Polynomial<ElemType>::const_iterator i;
644	<	typename Polynomial<ElemType>::const_iterator j;
643	>	/**
644	>	* Tests if two polynomial have the same exponents
645	>	* @return true if all of the exponents in these Polynomial are identical
646	>	* @param p1 the first polynomial
647	>	* @param p2 the second polynomial
648	>	* @note this function does not compare the coefficient
649	>	*/
650	>	template<typename Real>
651	>	bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
652
653	+	typename Polynomial<Real>::const_iterator i;
654	+	typename Polynomial<Real>::const_iterator j;
655	+
656		if (p1.size() != p2.size() ) {
657	<	return false;
657	>	return false;
658		}
659
660		for (i =  p1.begin(), j = p2.begin(); i  != p1.end() && j != p2.end(); ++i, ++j) {
661	<	if (i->first != j->first) {
662	<	return false;
663	<	}
661	>	if (i->first != j->first) {
662	>	return false;
663	>	}
664		}
665
666		return true;
667	<	}
667	>	}
668
279	–	typedef Polynomial<double> DoublePolynomial;
669
670	<	} //end namespace oopse
670	>
671	>	typedef Polynomial<RealType> DoublePolynomial;
672	>
673	>	} //end namespace OpenMD
674		#endif //MATH_POLYNOMIAL_HPP

Comparing:
trunk/src/math/Polynomial.hpp (property svn:keywords), Revision 385 by tim, Tue Mar 1 20:10:14 2005 UTC vs.
branches/development/src/math/Polynomial.hpp (property svn:keywords), Revision 1827 by gezelter, Wed Jan 9 19:49:07 2013 UTC

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

Diff Legend

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