| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4JTPolynomialSolver
27 //
28 // Class description:
29 //
30 // G4JTPolynomialSolver implements the Jenkins-Traub algorithm
31 // for real polynomial root finding.
32 // The solver returns -1, if the leading coefficient is zero,
33 // the number of roots found, otherwise.
34 //
35 // ----------------------------- INPUT --------------------------------
36 //
37 // op - double precision vector of coefficients in order of
38 // decreasing powers
39 // degree - integer degree of polynomial
40 //
41 // ----------------------------- OUTPUT -------------------------------
42 //
43 // zeror,zeroi - double precision vectors of the
44 // real and imaginary parts of the zeros
45 //
46 // ---------------------------- EXAMPLE -------------------------------
47 //
48 // G4JTPolynomialSolver trapEq ;
49 // G4double coef[8] ;
50 // G4double zr[7] , zi[7] ;
51 // G4int num = trapEq.FindRoots(coef,7,zr,zi);
52 //
53 // Translated from original TOMS493 Fortran77 routine (ANSI C, by C.Bond).
54
55 // Author: Oliver Link, 15.02.2005
56 // Translated to C++ and adapted to use STL vectors.
57 // --------------------------------------------------------------------
58 #ifndef G4JTPOLYNOMIALSOLVER_HH
59 #define G4JTPOLYNOMIALSOLVER_HH 1
60
61 #include <cmath>
62 #include <vector>
63
64 #include "globals.hh"
65
66 class G4JTPolynomialSolver
67 {
68 public:
69 G4JTPolynomialSolver() = default;
70 ~G4JTPolynomialSolver() = default;
71
72 G4int FindRoots(G4double* op, G4int degree, G4double* zeror, G4double* zeroi);
73
74 private:
75 void Quadratic(G4double a, G4double b1, G4double c, G4double* sr,
76 G4double* si, G4double* lr, G4double* li);
77 void ComputeFixedShiftPolynomial(G4int l2, G4int* nz);
78 void QuadraticPolynomialIteration(G4double* uu, G4double* vv, G4int* nz);
79 void RealPolynomialIteration(G4double* sss, G4int* nz, G4int* iflag);
80 void ComputeScalarFactors(G4int* type);
81 void ComputeNextPolynomial(G4int* type);
82 void ComputeNewEstimate(G4int type, G4double* uu, G4double* vv);
83 void QuadraticSyntheticDivision(G4int n, G4double* u, G4double* v,
84 std::vector<G4double>& p,
85 std::vector<G4double>& q, G4double* a,
86 G4double* b);
87
88 private:
89 std::vector<G4double> p;
90 std::vector<G4double> qp;
91 std::vector<G4double> k;
92 std::vector<G4double> qk;
93 std::vector<G4double> svk;
94
95 G4double sr = 0.0;
96 G4double si = 0.0;
97 G4double u = 0.0, v = 0.0;
98 G4double a = 0.0, b = 0.0, c = 0.0, d = 0.0;
99 G4double a1 = 0.0, a3 = 0.0, a7 = 0.0;
100 G4double e = 0.0, f = 0.0, g = 0.0, h = 0.0;
101 G4double szr = 0.0, szi = 0.0;
102 G4double lzr = 0.0, lzi = 0.0;
103 G4int n = 0;
104
105 /* The following statements set machine constants */
106
107 static const G4double base;
108 static const G4double eta;
109 static const G4double infin;
110 static const G4double smalno;
111 static const G4double are;
112 static const G4double mre;
113 static const G4double lo;
114 };
115
116 #endif
117