ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/math/Polynomial.hpp
(Generate patch)

Comparing trunk/src/math/Polynomial.hpp (file contents):
Revision 385 by tim, Tue Mar 1 20:10:14 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 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]  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;
66 >      result *= x;
67      }
68  
69      return result;
70 < }
70 >  }
71  
72 < /**
73 < * @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
74 < * A generic Polynomial class
75 < */
76 < template<typename ElemType>
77 < class Polynomial {
72 >  /**
73 >   * @class Polynomial Polynomial.hpp "math/Polynomial.hpp"
74 >   * A generic Polynomial class
75 >   */
76 >  template<typename Real>
77 >  class Polynomial {
78  
79 <    public:
80 <        
81 <        typedef int ExponentType;
82 <        typedef ElemType CoefficientType;
83 <        typedef std::map<ExponentType, CoefficientType> PolynomialPairMap;
84 <        typedef typename PolynomialPairMap::iterator iterator;
85 <        typedef typename PolynomialPairMap::const_iterator const_iterator;
86 <        /**
87 <         * Calculates the value of this Polynomial evaluated at the given x value.
88 <         * @return The value of this Polynomial evaluates at the given x value
89 <         * @param x the value of the independent variable for this Polynomial function
90 <         */
91 <        ElemType evaluate(const ElemType& x) {
92 <            ElemType result = ElemType();
93 <            ExponentType exponent;
94 <            CoefficientType coefficient;
79 >  public:
80 >    typedef Polynomial<Real> PolynomialType;    
81 >    typedef int ExponentType;
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(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
93 >     * Polynomial function
94 >     */
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;
104 <            }
100 >      for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
101 >        exponent = i->first;
102 >        coefficient = i->second;
103 >        result  += fastpow(x, exponent) * coefficient;
104 >      }
105  
106 <            return result;
107 <        }
106 >      return result;
107 >    }
108  
109 <        /**
110 <         * Returns the first derivative of this polynomial.
111 <         * @return the first derivative of this polynomial
112 <         * @param x
113 <         */
114 <        ElemType evaluateDerivative(const ElemType& x) {
115 <            ElemType result = ElemType();
116 <            ExponentType exponent;
117 <            CoefficientType coefficient;
109 >    /**
110 >     * Returns the first derivative of this polynomial.
111 >     * @return the first derivative of this polynomial
112 >     * @param x
113 >     */
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;
123 <            }
119 >      for (iterator i = polyPairMap_.begin(); i != polyPairMap_.end(); ++i) {
120 >        exponent = i->first;
121 >        coefficient = i->second;
122 >        result  += fastpow(x, exponent - 1) * coefficient * exponent;
123 >      }
124  
125 <            return result;
126 <        }
125 >      return result;
126 >    }
127  
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 –        }
128  
129 <        /**
130 <         * Set the coefficent of the specified exponent. If the coefficient is already there,  just add the
131 <         * new coefficient to the old one, otherwise,  just call setCoefficent
132 <         * @param exponent exponent of a term in this Polynomial
133 <         * @param coefficient multiplier of a term in this Polynomial
134 <         */
135 <        
136 <        void addCoefficient(int exponent, const ElemType& coefficient) {
137 <            iterator i = polyPairMap_.find(exponent);
129 >    /**
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 >    void setCoefficient(int exponent, const Real& coefficient) {
136 >      polyPairMap_[exponent] = coefficient;
137 >    }
138 >    
139 >    /**
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 >    void addCoefficient(int exponent, const Real& coefficient) {
147 >      iterator i = polyPairMap_.find(exponent);
148  
149 <            if (i != end()) {
150 <                i->second += coefficient;
151 <            } else {
152 <                setCoefficient(exponent, coefficient);
153 <            }
154 <        }
149 >      if (i != end()) {
150 >        i->second += coefficient;
151 >      } else {
152 >        setCoefficient(exponent, coefficient);
153 >      }
154 >    }
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
160 <         * @exponent exponent of any term in this Polynomial
161 <         */
162 <        ElemType getCoefficient(ExponentType exponent) {
163 <            iterator i = polyPairMap_.find(exponent);
164 <
158 <            if (i != end()) {
159 <                return i->second;
160 <            } else {
161 <                return ElemType(0);
162 <            }
163 <        }
156 >    /**
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 >    Real getCoefficient(ExponentType exponent) {
164 >      iterator i = polyPairMap_.find(exponent);
165  
166 <        iterator begin() {
167 <            return polyPairMap_.begin();
168 <        }
166 >      if (i != end()) {
167 >        return i->second;
168 >      } else {
169 >        return Real(0);
170 >      }
171 >    }
172  
173 <        const_iterator begin() const{
174 <            return polyPairMap_.begin();
175 <        }
173 >    iterator begin() {
174 >      return polyPairMap_.begin();
175 >    }
176 >
177 >    const_iterator begin() const{
178 >      return polyPairMap_.begin();
179 >    }
180          
181 <        iterator end() {
182 <            return polyPairMap_.end();
183 <        }
181 >    iterator end() {
182 >      return polyPairMap_.end();
183 >    }
184  
185 <        const_iterator end() const{
186 <            return polyPairMap_.end();
185 >    const_iterator end() const{
186 >      return polyPairMap_.end();
187 >    }
188 >
189 >    iterator find(ExponentType exponent) {
190 >      return polyPairMap_.find(exponent);
191 >    }
192 >
193 >    size_t size() {
194 >      return polyPairMap_.size();
195 >    }
196 >
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 >    PolynomialType& operator = (const PolynomialType& p) {
207 >
208 >      if (this != &p)  // protect against invalid self-assignment
209 >        {
210 >          typename Polynomial<Real>::const_iterator i;
211 >
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 <        iterator find(ExponentType exponent) {
223 <            return polyPairMap_.find(exponent);
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 >    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 <        size_t size() {
255 <            return polyPairMap_.size();
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 <    private:
306 <        
307 <        PolynomialPairMap polyPairMap_;
308 < };
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 < /**
326 < * Generates and returns the product of two given Polynomials.
327 < * @return A Polynomial containing the product of the two given Polynomial parameters
328 < */
329 < template<typename ElemType>
330 < Polynomial<ElemType> operator *(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
331 <    typename Polynomial<ElemType>::const_iterator i;
332 <    typename Polynomial<ElemType>::const_iterator j;
333 <    Polynomial<ElemType> p;
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 >  
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 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) {
557 <            p.addCoefficient( i->first + j->first, i->second * j->second);
558 <        }
556 >      for (j = p2.begin(); j !=p2.end(); ++j) {
557 >        p.addCoefficient( i->first + j->first, i->second * j->second);
558 >      }
559      }
560  
561      return p;
562 < }
562 >  }
563  
564 < /**
565 < * Generates and returns the sum of two given Polynomials.
566 < * @param p1 the first polynomial
567 < * @param p2 the second polynomial
568 < */
569 < template<typename ElemType>
570 < Polynomial<ElemType> operator +(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
571 <    Polynomial<ElemType> p(p1);
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.setCoefficient( i->first , i->second * v);
571 >    }
572  
573 <    typename Polynomial<ElemType>::const_iterator i;
573 >    return result;
574 >  }
575  
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.setCoefficient( i->first , i->second * v);
583 +    }
584 +
585 +    return result;
586 +  }
587 +  
588 +  /**
589 +   * Generates and returns the sum of two given Polynomials.
590 +   * @param p1 the first polynomial
591 +   * @param p2 the second polynomial
592 +   */
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<Real>::const_iterator i;
598 +
599      for (i =  p2.begin(); i  != p2.end(); ++i) {
600 <        p.addCoefficient(i->first, i->second);
600 >      p.addCoefficient(i->first, i->second);
601      }
602  
603      return p;
604  
605 < }
605 >  }
606  
607 < /**
608 < * Generates and returns the difference of two given Polynomials.
609 < * @return
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);
607 >  /**
608 >   * Generates and returns the difference of two given Polynomials.
609 >   * @return
610 >   * @param p1 the first polynomial
611 >   * @param p2 the second polynomial
612 >   */
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);
620 >      p.addCoefficient(i->first, -i->second);
621      }
622  
623      return p;
624  
625 < }
625 >  }
626  
627 < /**
628 < * Tests if two polynomial have the same exponents
629 < * @return true if these all of the exponents in these Polynomial are identical
630 < * @param p1 the first polynomial
631 < * @param p2 the second polynomial
632 < * @note this function does not compare the coefficient
633 < */
634 < template<typename ElemType>
635 < bool equal(const Polynomial<ElemType>& p1, const Polynomial<ElemType>& p2) {
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 <    typename Polynomial<ElemType>::const_iterator i;
649 <    typename Polynomial<ElemType>::const_iterator j;
648 >  /**
649 >   * Tests if two polynomial have the same exponents
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 Real>
656 >  bool equal(const Polynomial<Real>& p1, const Polynomial<Real>& p2) {
657  
658 +    typename Polynomial<Real>::const_iterator i;
659 +    typename Polynomial<Real>::const_iterator j;
660 +
661      if (p1.size() != p2.size() ) {
662 <        return false;
662 >      return false;
663      }
664      
665      for (i =  p1.begin(), j = p2.begin(); i  != p1.end() && j != p2.end(); ++i, ++j) {
666 <        if (i->first != j->first) {
667 <            return false;
668 <        }
666 >      if (i->first != j->first) {
667 >        return false;
668 >      }
669      }
670  
671      return true;
672 < }
672 >  }
673  
279 – 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

Diff Legend

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