00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
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
00061 double sqrt_D = std::sqrt(to_double(D));
00062
00063
00064
00065
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
00104
00105
00106
00107
00108
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) {
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