Changeset 7999 for trunk/MagicSoft/Mars/mbase

Timestamp:

10/01/06 22:19:55 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mbase

Files:

• MMath.cc (modified) (2 diffs)
• MMath.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -39 +39 @@
 #include <TArrayD.h>
 #endif
+
+#ifndef ROOT_TComplex
+#include <TComplex.h>
+#endif
+
 // --------------------------------------------------------------------------
 //
@@ -495 +500 @@
     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;
+}

trunk/MagicSoft/Mars/mbase/MMath.h

@@ -47 +47 @@
     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 +73 @@
     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; }
 }