Index: trunk/MagicSoft/Mars/mbase/MMath.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MMath.cc (revision 7971) +++ trunk/MagicSoft/Mars/mbase/MMath.cc (revision 7999) @@ -39,4 +39,9 @@ #include #endif + +#ifndef ROOT_TComplex +#include +#endif + // -------------------------------------------------------------------------- // @@ -495,2 +500,163 @@ return rc; } + +Int_t MMath::SolvePol2(Double_t a, Double_t b, Double_t &x1, Double_t &x2) +{ + const Double_t r = a*a - 4*b; + if (r<0) + return 0; + + if (r==0) + { + x1 = -a/2; + return 1; + } + + const Double_t s = TMath::Sqrt(r); + + x1 = (-a+s)/2; + x2 = (-a-s)/2; + + return 2; +} + +// -------------------------------------------------------------------------- +// +// This is a helper function making the execution of SolverPol3 a bit faster +// +static inline Double_t ReMul(const TComplex &c1, const TComplex &th) +{ + const TComplex c2 = TComplex::Cos(th/3.); + return c1.Re() * c2.Re() - c1.Im() * c2.Im(); +} + +// -------------------------------------------------------------------------- +// +// Solves: x^3 + ax^2 + bx + c = 0; +// Return number of the real solutions returns as z1, z2, z3 +// +// Algorithm adapted from http://home.att.net/~srschmitt/cubizen.heml +// Which is based on the solution given in +// http://mathworld.wolfram.com/CubicEquation.html +// +// ------------------------------------------------------------------------- +// +// Exact solutions of cubic polynomial equations +// by Stephen R. Schmitt Algorithm +// +// An exact solution of the cubic polynomial equation: +// +// x^3 + a*x^2 + b*x + c = 0 +// +// was first published by Gerolamo Cardano (1501-1576) in his treatise, +// Ars Magna. He did not discoverer of the solution; a professor of +// mathematics at the University of Bologna named Scipione del Ferro (ca. +// 1465-1526) is credited as the first to find an exact solution. In the +// years since, several improvements to the original solution have been +// discovered. Zeno source code +// +// % compute real or complex roots of cubic polynomial +// function cubic( var z1, z2, z3 : real, a, b, c : real ) : real +// +// var Q, R, D, S, T : real +// var im, th : real +// +// Q := (3*b - a^2)/9 +// R := (9*b*a - 27*c - 2*a^3)/54 +// D := Q^3 + R^2 % polynomial discriminant +// +// if (D >= 0) then % complex or duplicate roots +// +// S := sgn(R + sqrt(D))*abs(R + sqrt(D))^(1/3) +// T := sgn(R - sqrt(D))*abs(R - sqrt(D))^(1/3) +// +// z1 := -a/3 + (S + T) % real root +// z2 := -a/3 - (S + T)/2 % real part of complex root +// z3 := -a/3 - (S + T)/2 % real part of complex root +// im := abs(sqrt(3)*(S - T)/2) % complex part of root pair +// +// else % distinct real roots +// +// th := arccos(R/sqrt( -Q^3)) +// +// z1 := 2*sqrt(-Q)*cos(th/3) - a/3 +// z2 := 2*sqrt(-Q)*cos((th + 2*pi)/3) - a/3 +// z3 := 2*sqrt(-Q)*cos((th + 4*pi)/3) - a/3 +// im := 0 +// +// end if +// +// return im % imaginary part +// +// end function +// +Int_t MMath::SolvePol3(Double_t a, Double_t b, Double_t c, + Double_t &x1, Double_t &x2, Double_t &x3) +{ + const Double_t Q = (a*a - 3*b)/9; + const Double_t R = (9*b*a - 27*c - 2*a*a*a)/54; + const Double_t D = R*R - Q*Q*Q; // polynomial discriminant + + // ----- The single-real / duplicate-roots solution ----- + + if (D==0) + { + const Double_t r = MMath::Sqrt3(R); + + x1 = 2*r - a/3.; // real root + x2 = r - a/3.; // real root + + return 2; + } + + if (D>0) // complex or duplicate roots + { + const Double_t sqrtd = TMath::Sqrt(D); + + const Double_t S = MMath::Sqrt3(R + sqrtd); + const Double_t T = MMath::Sqrt3(R - sqrtd); + + x1 = (S+T) - a/3.; // real root + + return 1; + + //z2 = (S + T)/2 - a/3.; // real part of complex root + //z3 = (S + T)/2 - a/3.; // real part of complex root + //im = fabs(sqrt(3)*(S - T)/2) // complex part of root pair + } + + // ----- The general solution with three roots --- + + if (Q==0) + return 0; + + if (Q>0) // This is here for speed reasons + { + const Double_t sqrtq = TMath::Sqrt(Q); + const Double_t rq = R/TMath::Abs(Q); + + const Double_t th1 = TMath::ACos(rq/sqrtq); + const Double_t th2 = th1 + TMath::TwoPi(); + const Double_t th3 = th2 + TMath::TwoPi(); + + x1 = 2.*sqrtq * TMath::Cos(th1/3.) - a/3.; + x2 = 2.*sqrtq * TMath::Cos(th2/3.) - a/3.; + x3 = 2.*sqrtq * TMath::Cos(th3/3.) - a/3.; + + return 3; + } + + const TComplex sqrtq = TComplex::Sqrt(Q); + const Double_t rq = R/TMath::Abs(Q); + + const TComplex th1 = TComplex::ACos(rq/sqrtq); + const TComplex th2 = th1 + TMath::TwoPi(); + const TComplex th3 = th2 + TMath::TwoPi(); + + // For ReMul, see bove + x1 = ReMul(2.*sqrtq, th1) - a/3.; + x2 = ReMul(2.*sqrtq, th2) - a/3.; + x3 = ReMul(2.*sqrtq, th3) - a/3.; + + return 3; +} Index: trunk/MagicSoft/Mars/mbase/MMath.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MMath.h (revision 7971) +++ trunk/MagicSoft/Mars/mbase/MMath.h (revision 7999) @@ -47,4 +47,23 @@ Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x); + inline Int_t SolvePol1(Double_t c, Double_t d, Double_t &x1) + { + if (c==0) + return 0; + + x1 = -d/c; + return 1; + } + Int_t SolvePol2(Double_t c, Double_t d, Double_t &x1, Double_t &x2); + inline Int_t SolvePol2(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2) + { + return b==0 ? SolvePol1(c, d, x1) : SolvePol2(c/b, d/b, x1, x2); + } + Int_t SolvePol3(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2, Double_t &x3); + inline Int_t SolvePol3(Double_t a, Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2, Double_t &x3) + { + return a==0 ? SolvePol2(b, c, d, x1, x2) : SolvePol3(b/a, c/a, d/a, x1, x2, x3); + } + TArrayD LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y); TArrayD LeastSqFitExp(Int_t n, Double_t *x, Double_t *y); @@ -54,4 +73,6 @@ inline Int_t ModF(Double_t dbl, Double_t &frac) { Double_t rc; frac = modf(dbl, &rc); return TMath::Nint(rc); } + inline Double_t Sqrt3(Double_t x) { return TMath::Sign(TMath::Power(TMath::Abs(x), 1./3), x); } + inline Double_t Sgn(Double_t d) { return d<0 ? -1 : 1; } }