00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id: G4JTPolynomialSolver.hh 69792 2013-05-15 12:59:26Z gcosmo $ 00028 // 00029 // Class description: 00030 // 00031 // G4JTPolynomialSolver implements the Jenkins-Traub algorithm 00032 // for real polynomial root finding. 00033 // The solver returns -1, if the leading coefficient is zero, 00034 // the number of roots found, otherwise. 00035 // 00036 // ----------------------------- INPUT -------------------------------- 00037 // 00038 // op - double precision vector of coefficients in order of 00039 // decreasing powers 00040 // degree - integer degree of polynomial 00041 // 00042 // ----------------------------- OUTPUT ------------------------------- 00043 // 00044 // zeror,zeroi - double precision vectors of the 00045 // real and imaginary parts of the zeros 00046 // 00047 // ---------------------------- EXAMPLE ------------------------------- 00048 // 00049 // G4JTPolynomialSolver trapEq ; 00050 // G4double coef[8] ; 00051 // G4double zr[7] , zi[7] ; 00052 // G4int num = trapEq.FindRoots(coef,7,zr,zi); 00053 00054 // ---------------------------- HISTORY ------------------------------- 00055 // 00056 // Translated from original TOMS493 Fortran77 routine (ANSI C, by C.Bond). 00057 // Translated to C++ and adapted to use STL vectors, 00058 // by Oliver Link (Oliver.Link@cern.ch) 00059 // 00060 // -------------------------------------------------------------------- 00061 00062 #ifndef G4JTPOLYNOMIALSOLVER_HH 00063 #define G4JTPOLYNOMIALSOLVER_HH 00064 00065 #include <cmath> 00066 #include <vector> 00067 00068 #include "globals.hh" 00069 00070 class G4JTPolynomialSolver 00071 { 00072 00073 public: 00074 00075 G4JTPolynomialSolver(); 00076 ~G4JTPolynomialSolver(); 00077 00078 G4int FindRoots(G4double *op, G4int degree, 00079 G4double *zeror, G4double *zeroi); 00080 00081 private: 00082 00083 std::vector<G4double> p; 00084 std::vector<G4double> qp; 00085 std::vector<G4double> k; 00086 std::vector<G4double> qk; 00087 std::vector<G4double> svk; 00088 00089 G4double sr; 00090 G4double si; 00091 G4double u,v; 00092 G4double a,b,c,d; 00093 G4double a1,a3,a7; 00094 G4double e,f,g,h; 00095 G4double szr,szi; 00096 G4double lzr,lzi; 00097 G4int n; 00098 00099 /* The following statements set machine constants */ 00100 00101 static const G4double base; 00102 static const G4double eta; 00103 static const G4double infin; 00104 static const G4double smalno; 00105 static const G4double are; 00106 static const G4double mre; 00107 static const G4double lo; 00108 00109 void Quadratic(G4double a,G4double b1,G4double c, 00110 G4double *sr,G4double *si, G4double *lr,G4double *li); 00111 void ComputeFixedShiftPolynomial(G4int l2, G4int *nz); 00112 void QuadraticPolynomialIteration(G4double *uu,G4double *vv,G4int *nz); 00113 void RealPolynomialIteration(G4double *sss, G4int *nz, G4int *iflag); 00114 void ComputeScalarFactors(G4int *type); 00115 void ComputeNextPolynomial(G4int *type); 00116 void ComputeNewEstimate(G4int type,G4double *uu,G4double *vv); 00117 void QuadraticSyntheticDivision(G4int n, G4double *u, G4double *v, 00118 std::vector<G4double> &p, 00119 std::vector<G4double> &q, 00120 G4double *a, G4double *b); 00121 }; 00122 00123 #endif