G4JTPolynomialSolver.hh

Go to the documentation of this file.
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

Generated on Mon May 27 17:48:42 2013 for Geant4 by  doxygen 1.4.7