• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • Examples
  • File List
  • File Members

/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/math_polynomial_roots.h

Go to the documentation of this file.
00001 #ifndef OPENTISSUE_CORE_MATH_MATH_POLYNOMIAL_ROOTS_H
00002 #define OPENTISSUE_CORE_MATH_MATH_POLYNOMIAL_ROOTS_H
00003 //
00004 // OpenTissue Template Library
00005 // - A generic toolbox for physics-based modeling and simulation.
00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
00007 //
00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
00009 //
00010 #include <OpenTissue/configuration.h>
00011 
00012 #include <OpenTissue/core/math/math_precision.h>
00013 
00014 #include <OpenTissue/core/math/math_constants.h>
00015 
00016 #include <cmath>
00017 
00018 namespace OpenTissue
00019 {
00020 
00021   namespace math
00022   {
00023 
00038     template<typename real_type, typename real_array>
00039     inline bool compute_polynomial_roots(real_type const & c0, real_type const & c1, unsigned int & count, real_array & roots)
00040     {
00041       using std::fabs;
00042 
00043       static real_type const epsilon = math::working_precision<real_type>();
00044 
00045       if ( fabs(c1) >= epsilon )
00046       {
00047         roots[0] = -c0/c1;
00048         count = 1;
00049         return true;
00050       }
00051       count = 0;
00052       return false;
00053     }
00054 
00070     template<typename real_type, typename real_array>
00071     inline bool compute_polynomial_roots(real_type const & c0, real_type const & c1, real_type const & c2, unsigned int & count, real_array & roots)
00072     {
00073       using std::fabs;
00074       using std::sqrt;
00075 
00076       static real_type const epsilon = math::working_precision<real_type>();
00077       static real_type const four    = boost::numeric_cast<real_type>(4.0);
00078       static real_type const two     = detail::two<real_type>();
00079       static real_type const zero    = detail::zero<real_type>();
00080 
00081       if ( fabs(c2) <= epsilon )
00082         return compute_polynomial_roots(c0,c1,count,roots);
00083 
00084       real_type discr = c1*c1 - four*c0*c2;
00085 
00086       if ( fabs(discr) <= epsilon )
00087       {
00088         //discr = zero;
00089         roots[0] = -c1/(two*c2);
00090         count = 1;
00091         return true;
00092       }
00093 
00094       if ( discr < zero )
00095       {
00096         count = 0;
00097         return false;
00098       }
00099 
00100       assert( discr > zero || !"compute_polynomial_roots(): discriminant was expected to be positive");
00101       
00102       discr = sqrt(discr);
00103       roots[0] = (-c1 - discr)/(two*c2);
00104       roots[1] = (-c1 + discr)/(two*c2);
00105       count = 2;
00106       return true;
00107     }
00108 
00125     template<typename real_type, typename real_array>
00126     inline bool compute_polynomial_roots(real_type c0, real_type c1, real_type c2, real_type const & c3, unsigned int & count, real_array & roots)
00127     {
00128       using std::sqrt;
00129       using std::fabs;
00130       using std::pow;
00131       using std::atan2;
00132       using std::cos;
00133       using std::sin;
00134 
00135       static real_type const THIRD         = boost::numeric_cast<real_type>( 1.0/3.0   );
00136       static real_type const TWENTYSEVENTH = boost::numeric_cast<real_type>( 1.0/27.0  );
00137       static real_type const NINE          = boost::numeric_cast<real_type>( 9.0       );
00138       static real_type const TWO           = detail::two<real_type>();
00139       static real_type const ZERO          = detail::zero<real_type>();
00140       static real_type const SQRT3         = boost::numeric_cast<real_type>( sqrt(3.0) );
00141       static real_type const epsilon       = math::working_precision<real_type>();
00142 
00143       if ( fabs(c3) <= epsilon )
00144       {
00145         return compute_polynomial_roots(c0,c1,c2,count,roots);
00146       }
00147 
00148       // make polynomial monic, x^3+c2*x^2+c1*x+c0
00149       c0 /= c3;
00150       c1 /= c3;
00151       c2 /= c3;
00152 
00153       // convert to y^3 + a*y + b = 0 by x = y - c2/3
00154       real_type offset = THIRD*c2;
00155 
00156       real_type A      = c1 - c2*offset;
00157       real_type B      = c0 + c2*(TWO*c2*c2 - NINE*c1)*TWENTYSEVENTH;
00158       real_type half_B = B/TWO;
00159 
00160       real_type discr = half_B*half_B + A*A*A*TWENTYSEVENTH;
00161 
00162       if ( fabs(discr) <= epsilon )
00163       {
00164         //discr = ZERO;
00165 
00166         //if ( half_B >= ZERO )
00167         //  temp = -pow(half_B,THIRD);
00168         //else
00169         //  temp = pow(-half_B,THIRD);
00170         real_type temp = ( half_B >= ZERO )? -pow(half_B,THIRD) : pow(-half_B,THIRD);
00171 
00172         roots[0] = TWO*temp - offset;
00173         roots[1] = -temp - offset;
00174         roots[2] = roots[1];
00175         count = 3;
00176       }
00177       else if ( discr > ZERO )  // 1 real, 2 complex roots
00178       {
00179         discr = sqrt(discr);
00180         
00181         real_type temp = -half_B + discr;
00182         //if ( temp >= ZERO )
00183         //  roots[0] = pow(temp,THIRD);
00184         //else
00185         //  roots[0] = -pow(-temp,THIRD);
00186         roots[0] = ( temp >= ZERO ) ? pow(temp,THIRD) : -pow(-temp,THIRD);
00187         
00188         temp = - half_B - discr;
00189         
00190         //if ( temp >= ZERO )
00191         //  roots[0] += pow(temp,THIRD);
00192         //else
00193         //  roots[0] -= pow(-temp,THIRD);
00194         roots[0] += ( temp >= ZERO ) ? pow(temp,THIRD) : -pow(-temp,THIRD);
00195 
00196         roots[0] -= offset;
00197         count = 1;
00198       }
00199       else if ( discr < ZERO )
00200       {
00201         real_type dist  = sqrt(-THIRD*A);
00202         real_type angle = THIRD*atan2(sqrt(-discr), -half_B);
00203         real_type c     = cos(angle);
00204         real_type s     = sin(angle);
00205 
00206         roots[0] = TWO*dist*c-offset;
00207         roots[1] = -dist*(c+SQRT3*s)-offset;
00208         roots[2] = -dist*(c-SQRT3*s)-offset;
00209         count = 3;
00210       }
00211       return true;
00212     }
00213 
00231     template<typename real_type,typename real_array>
00232     inline bool compute_polynomial_roots(real_type c0, real_type c1, real_type c2, real_type c3, real_type const & c4, unsigned int & count, real_array & roots)
00233     {
00234       using std::fabs;
00235       using std::sqrt;
00236 
00237       static real_type const THREEQUATERS  = boost::numeric_cast<real_type>( 0.75 );
00238       static real_type const FOUR          = boost::numeric_cast<real_type>( 4.0 );
00239       static real_type const EIGHT         = boost::numeric_cast<real_type>( 8.0 );
00240       static real_type const ZERO          = detail::zero<real_type>();
00241       static real_type const TWO           = detail::two<real_type>();
00242       static real_type const ONE           = detail::one<real_type>();
00243 
00244       static real_type const epsilon = math::working_precision<real_type>();
00245 
00246       if ( fabs(c4) <= epsilon )
00247       {
00248         return compute_polynomial_roots(c0,c1,c2,c3,count,roots);
00249       }
00250 
00251       // make polynomial monic, x^4+c3*x^3+c2*x^2+c1*x+c0
00252       c0 /= c4;
00253       c1 /= c4;
00254       c2 /= c4;
00255       c3 /= c4;
00256 
00257       // reduction to resolvent cubic polynomial y^3+r2*y^2+r1*y+r0 = 0
00258       real_type r0 = -c3*c3*c0 + FOUR*c2*c0 - c1*c1;
00259       real_type r1 = c3*c1 - FOUR*c0;
00260       real_type r2 = -c2;
00261       compute_polynomial_roots(r0,r1,r2,ONE,count,roots);  // always produces at least one root
00262       real_type Y = roots[0];
00263 
00264       count = 0;
00265       real_type discr = (c3*c3)/FOUR - c2 + Y;
00266 
00267       if ( fabs(discr) <= epsilon )
00268         discr = ZERO;
00269 
00270       if ( discr < ZERO )
00271       {
00272         count = 0;
00273       }
00274       else if ( discr > ZERO )
00275       {
00276         real_type r       = sqrt(discr);
00277         real_type t1      = THREEQUATERS*c3*c3 - r*r - TWO*c2;
00278         real_type t2      = (FOUR*c3*c2 - EIGHT*c1 - c3*c3*c3) /(FOUR*r);
00279 
00280         real_type t_plus  = t1+t2;
00281         real_type t_minus = t1-t2;
00282 
00283         if ( fabs(t_plus) <= epsilon )
00284           t_plus = ZERO;
00285         if ( fabs(t_minus) <= epsilon )
00286           t_minus = ZERO;
00287 
00288         if ( t_plus >= ZERO )
00289         {
00290           real_type D = sqrt(t_plus);
00291           roots[0] = -c3/FOUR + (r+D)/TWO;
00292           roots[1] = -c3/FOUR + (r-D)/TWO;
00293           count += 2;
00294         }
00295         if ( t_minus >= ZERO )
00296         {
00297           real_type E = sqrt(t_minus);
00298           roots[count++] = -c3/FOUR + (E-r)/TWO;
00299           roots[count++] = -c3/FOUR - (E+r)/TWO;
00300         }
00301       }
00302       else
00303       {
00304         real_type t2 = Y*Y - FOUR*c0;
00305         if ( t2 >= - epsilon )
00306         {
00307           if ( t2 < ZERO ) // round to zero
00308             t2 = ZERO;
00309 
00310           t2 = TWO*sqrt(t2);
00311           real_type t1 = THREEQUATERS*c3*c3 - TWO*c2;
00312           if ( t1+t2 >= epsilon )
00313           {
00314             real_type D = sqrt(t1+t2);
00315             roots[0] = -c3/FOUR + D/TWO;
00316             roots[1] = -c3/FOUR - D/TWO;
00317             count += 2;
00318           }
00319           if ( t1-t2 >= epsilon )
00320           {
00321             real_type E = sqrt(t1-t2);
00322             roots[count++] = -c3/FOUR + E/TWO;
00323             roots[count++] = -c3/FOUR - E/TWO;
00324           }
00325         }
00326       }
00327 
00328       return count > 0;
00329     }
00330 
00331   } // namespace math
00332 
00333 } // namespace OpenTissue
00334 
00335 // OPENTISSUE_CORE_MATH_MATH_POLYNOMIAL_ROOTS_H
00336 #endif

Generated on Thu Dec 1 2011 12:52:08 for HUMIM Tracker by  doxygen 1.7.1