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

Comparing trunk/src/math/Polynomial.hpp (file contents):
Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
Revision 1358 by gezelter, Thu Aug 13 21:21:51 2009 UTC

# 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 {
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 <        
131 <    void setCoefficient(int exponent, const ElemType& 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 150 | Line 154 | namespace oopse {
154      }
155  
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 188 | Line 194 | namespace oopse {
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<ElemType>::const_iterator i;
209 >        {
210 >          typename Polynomial<Real>::const_iterator i;
211  
212 <        polyPairMap_.clear();  // clear out the old map
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);
214 >          for (i =  p.begin(); i != p.end(); ++i) {
215 >            this->setCoefficient(i->first, i->second);
216 >          }
217          }
202 –      }
218        // by convention, always return *this
219        return *this;
220      }
221  
222      PolynomialType& operator += (const PolynomialType& p) {
223 <        typename Polynomial<ElemType>::const_iterator i;
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 <        }
225 >      for (i =  p.begin(); i  != p.end(); ++i) {
226 >        this->addCoefficient(i->first, i->second);
227 >      }
228  
229 <        return *this;        
229 >      return *this;        
230      }
231  
232      PolynomialType& operator -= (const PolynomialType& p) {
233 <        typename Polynomial<ElemType>::const_iterator i;
234 <        for (i =  p.begin(); i  != p.end(); ++i) {
235 <          this->addCoefficient(i->first, -i->second);
236 <        }        
237 <        return *this;
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 <
239 >    
240      PolynomialType& operator *= (const PolynomialType& p) {
241 <    typename Polynomial<ElemType>::const_iterator i;
242 <    typename Polynomial<ElemType>::const_iterator j;
243 <    Polynomial<ElemType> 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);
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 <    return *this;
265 >
266 >    PolynomialType& operator += (const Real v) {    
267 >      this->addCoefficient( 0, v);
268 >      return *this;
269      }
270  
271 <    //PolynomialType& operator *= (const ElemType v)
272 <    PolynomialType& operator *= (const ElemType v) {
273 <    typename Polynomial<ElemType>::const_iterator i;
274 <    //Polynomial<ElemType> result;
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 <    for (i = this->begin(); i != this->end(); ++i) {
292 <        this->setCoefficient( i->first, i->second*v);
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(reals(i), imags(i)));
320 +      }
321  
322 <    return *this;
322 >      return roots;
323      }
324  
325 <    PolynomialType& operator += (const ElemType v) {    
326 <    this->addCoefficient( 0, v);
327 <    return *this;
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:
# Line 258 | Line 541 | namespace oopse {
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 278 | 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.setCoefficient( i->first , i->second * v);
# Line 290 | Line 573 | namespace oopse {
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.setCoefficient( i->first , i->second * v);
# Line 307 | 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 327 | 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 342 | 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 +    ExponentType exponent;
637 +    CoefficientType 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 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 367 | Line 671 | namespace oopse {
671      return true;
672    }
673  
674 +
675 +
676    typedef Polynomial<RealType> DoublePolynomial;
677  
678   } //end namespace oopse

Diff Legend

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