synaps/usolve/solve_by_radical.h

00001 // Copyright (c) 2003  INRIA Sophia-Antipolis (France) and
00002 //                     Department of Informatics and Telecommunications
00003 //                     University of Athens (Greece).
00004 //
00005 // Authors : Elias P. TSIGARIDAS <et@di.uoa.gr>
00006 
00007 // Partially supported by INRIA's project "CALAMATA", a
00008 // bilateral collaboration with National Kapodistrian University of
00009 // Athens.
00010 // Partially supported by the IST Programme of the EU as a Shared-cost
00011 // RTD (FET Open) Project under Contract No  IST-2000-26473
00012 // (ECG - Effective Computational Geometry for Curves and Surfaces)
00013 
00014 
00015 // This code is a modification of the code from
00016 // Polynomial root finder (polynoms up to degree 4)
00017 // AUTHOR: Jochen SCHARZE
00018 //           (See 'Cubic & Quartic Roots' in 'Graphics Gems 1', AP)
00019 // CUBIC: Modification from the GNU gsl library
00020 // QUARTIC: We use the gnu gsl-solver
00021 
00022 // 27/11/2003: Modification mourrain@sophia.inria.fr for generic coefficients.
00023 // 05/12/2003: Modification et@di.uoa.gr multiple roots and vectors as input
00024 
00025 
00026 #ifndef SYNAPS_ARITHM_ALGEBRAIC_SOLVE_BY_RADICAL_H
00027 #define SYNAPS_ARITHM_ALGEBRAIC_SOLVE_BY_RADICAL_H
00028 
00029 #include <synaps/init.h>
00030 #include <synaps/usolve/Traits.h>
00031 
00032 __BEGIN_NAMESPACE_SYNAPS
00033 
00034 #define EQN_EPS 1e-5
00035 
00036 static int isZero(double x)
00037 {
00038         return x > -EQN_EPS && x < EQN_EPS;
00039 }
00040 
00041 
00042 template<class V>
00043 int solve_quadratic(std::vector<double>& s, const V& c)
00044 {
00045         typedef typename V::value_type                RT;
00046         typedef typename Integer_type_traits<RT>::FT  FT;
00047         typedef Rational_traits<FT>                   Rat_traits;
00048 
00049         RT D = c[1] * c[1] - RT(4) *  c[2] * c[0];
00050         if (D == RT(0)) {
00051                 // one double root
00052                 FT r = Rat_traits::make_rational(RT(-c[1]),  RT(c[2] * 2));
00053                 s.push_back(to_double(r));
00054                 s.push_back(to_double(r));
00055                 return 2;
00056         }
00057 
00058         if (D < RT(0)) return 0;
00059 
00060         // two real roots
00061         double sqrt_D = std::sqrt(to_double(D));
00062 
00063 
00064 //      FT r1 = Rat_traits::make_rational( RT(-c[1] - sqrt_D),   RT(c[2] * 2) );
00065 //      FT r2 = Rat_traits::make_rational( RT(-c[1] + sqrt_D),   RT(c[2] * 2) );
00066 
00067         FT r1 = Rat_traits::make_rational( RT(-c[1] ),   RT(c[2] * 2) );
00068         FT r2 = Rat_traits::make_rational( RT(-c[1] ),   RT(c[2] * 2) );
00069 
00070         r1 += FT(sqrt_D) /  RT(c[2] * 2);
00071         r1 -= FT(sqrt_D) /  RT(c[2] * 2);
00072 
00073         s.push_back(to_double(r1));
00074         s.push_back(to_double(r2));
00075         return 2;
00076 }
00077 
00078         template<class NT, class VV>
00079 int solve_cubic(VV& s, const VV& coeff)
00080 {
00081         double a = coeff[2] / coeff[3];
00082         double b = coeff[1] / coeff[3];
00083         double c = coeff[0] / coeff[3];
00084 
00085         double q = (a * a - 3 * b);
00086         double r = (2 * a * a * a - 9 * a * b + 27 * c);
00087 
00088         double Q = q / 9;
00089         double R = r / 54;
00090 
00091         double Q3 = Q * Q * Q;
00092         double R2 = R * R;
00093 
00094         double CR2 = 729 * r * r;
00095         double CQ3 = 2916 * q * q * q;
00096 
00097         if (R == 0 && Q == 0) {
00098                 s.push_back( - a / 3 );
00099                 s.push_back( - a / 3 );
00100                 s.push_back( - a / 3 );
00101                 return 3 ;
00102         } else if (CR2 == CQ3) {
00103                 /* this test is actually R2 == Q3, written in a form suitable
00104                    for exact computation with integers */
00105 
00106                 /* Due to finite precision some double roots may be missed, and
00107                    considered to be a pair of complex roots z = x +/- epsilon i
00108                    close to the real axis. */
00109                 double sqrtQ = sqrt (Q);
00110                 s.push_back( -2 * sqrtQ  - a / 3 );
00111                 s.push_back( sqrtQ - a / 3 );
00112                 s.push_back( sqrtQ - a / 3 );
00113                 sort(s.begin(), s.end());
00114                 return 3 ;
00115         } else if (CR2 < CQ3)  { /* equivalent to R2 < Q3 */
00116                 double sqrtQ = sqrt (Q);
00117                 double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
00118                 double theta = acos (R / sqrtQ3);
00119                 double norm = -2 * sqrtQ;
00120                 s.push_back( norm * cos (theta / 3) - a / 3 );
00121                 s.push_back( norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3 );
00122                 s.push_back( norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3 );
00123 
00124                 sort(s.begin(), s.end());
00125                 return 3;
00126         } else {
00127                 double sgnR = (R >= 0 ? 1 : -1);
00128                 double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
00129                 double B = Q / A ;
00130                 s.push_back( A + B - a / 3 );
00131                 return 1;
00132         }
00133 }
00134 
00135 __END_NAMESPACE_SYNAPS
00136 
00137 #endif // SYNAPS_ARITHM_ALGEBRAIC_SOLVE_BY_RADICAL_H
00138 

SYNAPS DOCUMENTATION
logo