lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Import Geant4 7.1.0 source tree

Browse Source
This commit is contained in:
Gabriele Cosmo
2016-06-09 12:11:21 +02:00
parent 516dbf1a58
commit d93e1e39a9
5384 changed files with 125662 additions and 82444 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/global/HEPNumerics/src/G4AnalyticalPolSolver.cc
+442
View File
@@ -0,0 +1,442 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4AnalyticalPolSolver.cc,v 1.5 2005/05/19 07:40:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "globals.hh"
#include <complex>
#include "G4AnalyticalPolSolver.hh"
//////////////////////////////////////////////////////////////////////////////
G4AnalyticalPolSolver::G4AnalyticalPolSolver() {;}
//////////////////////////////////////////////////////////////////////////////
G4AnalyticalPolSolver::~G4AnalyticalPolSolver() {;}
//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5] p[5]
// Roots of poly p[0] x^2 + p[1] x+p[2]=0
//
// x = r[1][k] + i r[2][k]; k = 1, 2
G4int G4AnalyticalPolSolver::QuadRoots( G4double p[5], G4double r[3][5] )
{
G4double b, c, d2, d;
b = -p[1]/p[0]/2.;
c = p[2]/p[0];
d2 = b*b - c;
if( d2 >= 0. )
{
d = std::sqrt(d2);
r[1][1] = b - d;
r[1][2] = b + d;
r[2][1] = 0.;
r[2][2] = 0.;
}
else
{
d = std::sqrt(-d2);
r[2][1] = d;
r[2][2] = -d;
r[1][1] = b;
r[1][2] = b;
}
return 2;
}
//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5] p[5]
// Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
// x=r[1][k] + i r[2][k] k=1,...,3
// Assumes 0<arctan(x)<pi/2 for x>0
G4int G4AnalyticalPolSolver::CubicRoots( G4double p[5], G4double r[3][5] )
{
G4double s,t,b,c,d;
G4int k;
if( p[0] != 1. )
{
for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
p[0] = 1.;
}
s = p[1]/3.0;
t = s*p[1];
b = 0.5*( s*( t/1.5 - p[2] ) + p[3] );
t = ( t - p[2] )/3.0;
c = t*t*t;
d = b*b - c;
if( d >= 0. )
{
d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
if( d != 0. )
{
if( b > 0. ) { b = -d; }
else { b = d; }
c = t/b;
}
d = std::sqrt(0.75)*(b - c);
r[2][2] = d;
b = b + c;
c = -0.5*b-s;
r[1][2] = c;
if( ( b > 0. && s <= 0. ) || ( b < 0. && s > 0. ) )
{
r[1][1] = c;
r[2][1] = -d;
r[1][3] = b - s;
r[2][3] = 0;
}
else
{
r[1][1] = b - s;
r[2][1] = 0.;
r[1][3] = c;
r[2][3] = -d;
}
} // end of 2 equal or complex roots
else
{
if( b == 0. ) { d = std::atan(1.0)/1.5; }
else { d = std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
if( b < 0. ) { b = std::sqrt(t)*2.0; }
else { b = -2.0*std::sqrt(t); }
c = std::cos(d)*b;
t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
d = -t - c - s;
c = c - s;
t = t - s;
if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
else
{
r[1][3] = t;
t = c;
}
if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
else
{
r[1][2] = t;
t = d;
}
r[1][1] = t;
for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
}
return 0;
}
//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5] p[5]
// Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
// x=r[1][k] + i r[2][k] k=1,...,4
G4int G4AnalyticalPolSolver::BiquadRoots( G4double p[5], G4double r[3][5] )
{
G4double a, b, c, d, e;
G4int i, k, j, noRoots;
if(p[0] != 1.0)
{
for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
p[0] = 1.;
}
e = 0.25*p[1];
b = 2*e;
c = b*b;
d = 0.75*c;
b = p[3] + b*( c - p[2] );
a = p[2] - d;
c = p[4] + e*( e*a - p[3] );
a = a - d;
p[1] = 0.5*a;
p[2] = (p[1]*p[1]-c)*0.25;
p[3] = b*b/(-64.0);
if( p[3] < 0. )
{
noRoots = CubicRoots(p,r);
for( k = 1; k < 4; k++ )
{
if( r[2][k] == 0. && r[1][k] > 0 )
{
d = r[1][k]*4;
a = a + d;
if ( a >= 0. && b >= 0.) { p[1] = std::sqrt(d); }
else if( a <= 0. && b <= 0.) { p[1] = std::sqrt(d); }
else { p[1] = -std::sqrt(d); }
b = 0.5*( a + b/p[1] );
p[2] = c/b;
noRoots = QuadRoots(p,r);
for( i = 1; i < 3; i++ )
{
for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
}
p[1] = -p[1];
p[2] = b;
noRoots = QuadRoots(p,r);
for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }
return 4;
}
}
}
if( p[2] < 0. )
{
b = std::sqrt(c);
d = b + b - a;
p[1] = 0.;
if( d > 0. ) { p[1] = std::sqrt(d); }
}
else
{
if( p[1] > 0.) { b = std::sqrt(p[2])*2.0 + p[1]; }
else { b = -std::sqrt(p[2])*2.0 + p[1]; }
if( b != 0.) { p[1] = 0; }
else
{
for(k = 1; k < 5; k++ )
{
r[1][k] = -e;
r[2][k] = 0;
}
return 0;
}
}
p[2] = c/b;
noRoots = QuadRoots(p,r);
for( k = 1; k < 3; k++ )
{
for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
}
p[1] = -p[1];
p[2] = b;
noRoots = QuadRoots(p,r);
for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }
return 4;
}
//////////////////////////////////////////////////////////////////////////////
G4int G4AnalyticalPolSolver::QuarticRoots( G4double p[5], G4double r[3][5])
{
G4double a0, a1, a2, a3, y1;
G4double R2, D2, E2, D, E, R = 0.;
G4double a, b, c, d, ds;
G4double reRoot[4];
G4int k, noRoots, noReRoots = 0;
for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }
if( p[0] != 1.0 )
{
for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
p[0] = 1.;
}
a3 = p[1];
a2 = p[2];
a1 = p[3];
a0 = p[4];
// resolvent cubic equation cofs:
p[1] = -a2;
p[2] = a1*a3 - 4*a0;
p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;
noRoots = CubicRoots(p,r);
for( k = 1; k < 4; k++ )
{
if( r[2][k] == 0. ) // find a real root
{
noReRoots++;
reRoot[k] = r[1][k];
}
else reRoot[k] = DBL_MAX; // kInfinity;
}
y1 = DBL_MAX; // kInfinity;
for( k = 1; k < 4; k++ )
{
if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
}
R2 = 0.25*a3*a3 - a2 + y1;
b = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
c = 0.75*a3*a3 - 2*a2;
a = c - R2;
d = 4*y1*y1 - 16*a0;
if( R2 > 0.)
{
R = std::sqrt(R2);
D2 = a + b/R;
E2 = a - b/R;
if( D2 >= 0. )
{
D = std::sqrt(D2);
r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
r[2][1] = 0.;
r[2][2] = 0.;
}
else
{
D = std::sqrt(-D2);
r[1][1] = -0.25*a3 + 0.5*R;
r[1][2] = -0.25*a3 + 0.5*R;
r[2][1] = 0.5*D;
r[2][2] = -0.5*D;
}
if( E2 >= 0. )
{
E = std::sqrt(E2);
r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
r[2][3] = 0.;
r[2][4] = 0.;
}
else
{
E = std::sqrt(-E2);
r[1][3] = -0.25*a3 - 0.5*R;
r[1][4] = -0.25*a3 - 0.5*R;
r[2][3] = 0.5*E;
r[2][4] = -0.5*E;
}
}
else if( R2 < 0.)
{
R = std::sqrt(-R2);
G4complex CD2(a,-b/R);
G4complex CD = std::sqrt(CD2);
r[1][1] = -0.25*a3 + 0.5*real(CD);
r[1][2] = -0.25*a3 - 0.5*real(CD);
r[2][1] = 0.5*R + 0.5*imag(CD);
r[2][2] = 0.5*R - 0.5*imag(CD);
G4complex CE2(a,b/R);
G4complex CE = std::sqrt(CE2);
r[1][3] = -0.25*a3 + 0.5*real(CE);
r[1][4] = -0.25*a3 - 0.5*real(CE);
r[2][3] = -0.5*R + 0.5*imag(CE);
r[2][4] = -0.5*R - 0.5*imag(CE);
}
else // R2=0 case
{
if(d >= 0.)
{
D2 = c + std::sqrt(d);
E2 = c - std::sqrt(d);
if( D2 >= 0. )
{
D = std::sqrt(D2);
r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
r[2][1] = 0.;
r[2][2] = 0.;
}
else
{
D = std::sqrt(-D2);
r[1][1] = -0.25*a3 + 0.5*R;
r[1][2] = -0.25*a3 + 0.5*R;
r[2][1] = 0.5*D;
r[2][2] = -0.5*D;
}
if( E2 >= 0. )
{
E = std::sqrt(E2);
r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
r[2][3] = 0.;
r[2][4] = 0.;
}
else
{
E = std::sqrt(-E2);
r[1][3] = -0.25*a3 - 0.5*R;
r[1][4] = -0.25*a3 - 0.5*R;
r[2][3] = 0.5*E;
r[2][4] = -0.5*E;
}
}
else
{
ds = std::sqrt(-d);
G4complex CD2(c,ds);
G4complex CD = std::sqrt(CD2);
r[1][1] = -0.25*a3 + 0.5*real(CD);
r[1][2] = -0.25*a3 - 0.5*real(CD);
r[2][1] = 0.5*R + 0.5*imag(CD);
r[2][2] = 0.5*R - 0.5*imag(CD);
G4complex CE2(c,-ds);
G4complex CE = std::sqrt(CE2);
r[1][3] = -0.25*a3 + 0.5*real(CE);
r[1][4] = -0.25*a3 - 0.5*real(CE);
r[2][3] = -0.5*R + 0.5*imag(CE);
r[2][4] = -0.5*R - 0.5*imag(CE);
}
}
return 4;
}
//
//
//////////////////////////////////////////////////////////////////////////////
source/global/HEPNumerics/src/G4ChebyshevApproximation.cc
+59 -79
View File
@@ -21,37 +21,31 @@
// ********************************************************************
//
//
// $Id: G4ChebyshevApproximation.cc,v 1.4 2004/11/12 17:38:32 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4ChebyshevApproximation.cc,v 1.5 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4ChebyshevApproximation.hh"
// Constructor for initialisation of the class data members. It creates the array
// fChebyshevCof[0,...,fNumber-1], fNumber = n ; which consists of Chebyshev
// coefficients describing the function pointed by pFunction. The values a and b
// fixe the interval of validity of Chebyshev approximation.
// Constructor for initialisation of the class data members.
// It creates the array fChebyshevCof[0,...,fNumber-1], fNumber = n ;
// which consists of Chebyshev coefficients describing the function
// pointed by pFunction. The values a and b fix the interval of validity
// of the Chebyshev approximation.
G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
G4int n,
G4double a,
G4double b )
G4double b )
: fFunction(pFunction), fNumber(n),
fChebyshevCof(new G4double[fNumber]),
fMean(0.5*(b+a)), fDiff(0.5*(b-a))
{
G4int i, j ;
G4double rootSum, cof, cofj, weight ;
fFunction = pFunction ;
fNumber = n ;
fDiff = 0.5*(b-a) ;
fMean = 0.5*(b+a) ;
fChebyshevCof = new G4double[fNumber] ;
G4int i=0, j=0 ;
G4double rootSum=0.0, cofj=0.0 ;
G4double* tempFunction = new G4double[fNumber] ;
weight = 2.0/fNumber ;
cof = 0.5*weight*pi ; // pi/n
G4double weight = 2.0/fNumber ;
G4double cof = 0.5*weight*pi ; // pi/n
for (i=0;i<fNumber;i++)
{
@@ -77,34 +71,26 @@ G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
// Constructor for creation of Chebyshev coefficients for m-derivative
// from pFunction. The value of m ! MUST BE ! < n , because the result
// array of fChebyshevCof will be of (n-m) size. The values a and b
// fixe the interval of validity of Chebyshev approximation.
// fix the interval of validity of the Chebyshev approximation.
G4ChebyshevApproximation::
G4ChebyshevApproximation( function pFunction,
G4int n,
G4int m,
G4double a,
G4double b )
G4int n, G4int m,
G4double a, G4double b )
: fFunction(pFunction), fNumber(n),
fChebyshevCof(new G4double[fNumber]),
fMean(0.5*(b+a)), fDiff(0.5*(b-a))
{
if(n <= m)
{
G4Exception
("Invalid arguments in G4ChebyshevApproximation::G4ChebyshevApproximation") ;
G4Exception("G4ChebyshevApproximation::G4ChebyshevApproximation()",
"InvalidCall", FatalException, "Invalid arguments !") ;
}
G4int i, j ;
G4double rootSum, cof, cofj, weight ;
fFunction = pFunction ;
fNumber = n ;
fDiff = 0.5*(b-a) ;
fMean = 0.5*(b+a) ;
fChebyshevCof = new G4double[fNumber] ;
G4int i=0, j=0 ;
G4double rootSum = 0.0, cofj=0.0;
G4double* tempFunction = new G4double[fNumber] ;
weight = 2.0/fNumber ;
cof = 0.5*weight*pi ; // pi/n
G4double weight = 2.0/fNumber ;
G4double cof = 0.5*weight*pi ; // pi/n
for (i=0;i<fNumber;i++)
{
@@ -130,7 +116,7 @@ G4ChebyshevApproximation( function pFunction,
fNumber-- ;
for(j=0;j<fNumber;j++)
{
fChebyshevCof[j] = tempFunction[j] ; // corresponds to (i)-derivative
fChebyshevCof[j] = tempFunction[j] ; // corresponds to (i)-derivative
}
}
delete[] tempFunction ; // delete of dynamically allocated tempFunction
@@ -140,25 +126,20 @@ G4ChebyshevApproximation( function pFunction,
//
// Constructor for creation of Chebyshev coefficients for integral
// from pFunction.
G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
G4double a,
G4double b,
G4double b,
G4int n )
: fFunction(pFunction), fNumber(n),
fChebyshevCof(new G4double[fNumber]),
fMean(0.5*(b+a)), fDiff(0.5*(b-a))
{
G4int i,j;
G4double rootSum, cof, cofj, weight ;
fFunction = pFunction ;
fNumber = n ;
fDiff = 0.5*(b-a) ;
fMean = 0.5*(b+a) ;
fChebyshevCof = new G4double[fNumber] ;
G4int i=0, j=0;
G4double rootSum=0.0, cofj=0.0;
G4double* tempFunction = new G4double[fNumber] ;
weight = 2.0/fNumber ;
cof = 0.5*weight*pi ; // pi/n
G4double weight = 2.0/fNumber;
G4double cof = 0.5*weight*pi ; // pi/n
for (i=0;i<fNumber;i++)
{
@@ -208,8 +189,8 @@ G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
{
if(number < 0 && number >= fNumber)
{
G4Exception
("Argument out of range in G4ChebyshevApproximation::GetChebyshevCof") ;
G4Exception("G4ChebyshevApproximation::GetChebyshevCof()",
"InvalidCall", FatalException, "Argument out of range !") ;
}
return fChebyshevCof[number] ;
}
@@ -222,22 +203,23 @@ G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
G4double
G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
{
G4int i;
G4double evaluate = 0.0, evaluate2 = 0.0, temp, xReduced, xReduced2 ;
G4double evaluate = 0.0, evaluate2 = 0.0, temp = 0.0,
xReduced = 0.0, xReduced2 = 0.0 ;
if ((x-fMean+fDiff)*(x-fMean-fDiff) > 0.0)
{
G4Exception("Invalid argument in G4ChebyshevApproximation::ChebyshevEvaluation");
}
xReduced = (x-fMean)/fDiff ;
xReduced2 = 2.0*xReduced ;
for (i=fNumber-1;i>=1;i--)
{
temp = evaluate ;
evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
evaluate2 = temp ;
}
return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
if ((x-fMean+fDiff)*(x-fMean-fDiff) > 0.0)
{
G4Exception("G4ChebyshevApproximation::ChebyshevEvaluation()",
"InvalidCall", FatalException, "Invalid argument !") ;
}
xReduced = (x-fMean)/fDiff ;
xReduced2 = 2.0*xReduced ;
for (G4int i=fNumber-1;i>=1;i--)
{
temp = evaluate ;
evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
evaluate2 = temp ;
}
return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
}
// ------------------------------------------------------------------
@@ -248,17 +230,16 @@ G4Exception("Invalid argument in G4ChebyshevApproximation::ChebyshevEvaluation")
void
G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
{
G4int i ;
G4double cof = 1.0/fDiff ;
derCof[fNumber-1] = 0.0 ;
derCof[fNumber-2] = 2*(fNumber-1)*fChebyshevCof[fNumber-1] ;
for(i=fNumber-3;i>=0;i--)
for(G4int i=fNumber-3;i>=0;i--)
{
derCof[i] = derCof[i+2] + 2*(i+1)*fChebyshevCof[i+1] ;
}
for(i=0;i<fNumber;i++)
for(G4int j=0;j<fNumber;j++)
{
derCof[i] *= cof ;
derCof[j] *= cof ;
}
}
@@ -274,9 +255,8 @@ G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
void
G4ChebyshevApproximation::IntegralChebyshevCof(G4double integralCof[]) const
{
G4int i ;
G4double cof = 0.5*fDiff, sum = 0.0, factor = 1.0 ;
for(i=1;i<fNumber-1;i++)
for(G4int i=1;i<fNumber-1;i++)
{
integralCof[i] = cof*(fChebyshevCof[i-1] - fChebyshevCof[i+1])/i ;
sum += factor*integralCof[i] ;
source/global/HEPNumerics/src/G4DataInterpolation.cc
+175 -184
View File
@@ -21,25 +21,25 @@
// ********************************************************************
//
//
// $Id: G4DataInterpolation.cc,v 1.6 2004/11/12 17:38:32 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4DataInterpolation.cc,v 1.7 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4DataInterpolation.hh"
//////////////////////////////////////////////////////////////////////////////
//
// Constructor for initializing of fArgument, fFunction and fNumber data members
// Constructor for initializing of fArgument, fFunction and fNumber
// data members
G4DataInterpolation::G4DataInterpolation( G4double pX[],
G4double pY[],
G4int number ):
fSecondDerivative (0)
G4double pY[],
G4int number )
: fArgument(new G4double[number]),
fFunction(new G4double[number]),
fSecondDerivative(0),
fNumber(number)
{
G4int i ;
fNumber = number ;
fArgument = new G4double[fNumber] ;
fFunction = new G4double[fNumber] ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
fArgument[i] = pX[i] ;
fFunction[i] = pY[i] ;
@@ -54,18 +54,18 @@ G4DataInterpolation::G4DataInterpolation( G4double pX[],
G4DataInterpolation::G4DataInterpolation( G4double pX[],
G4double pY[],
G4int number,
G4double pFirstDerStart,
G4double pFirstDerFinish )
G4double pY[],
G4int number,
G4double pFirstDerStart,
G4double pFirstDerFinish )
: fArgument(new G4double[number]),
fFunction(new G4double[number]),
fSecondDerivative(new G4double[number]),
fNumber(number)
{
G4int i, k ;
G4double p, qn, sig, un ;
G4int i=0 ;
G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ;
const G4double maxDerivative = 0.99e30 ;
fNumber = number ;
fArgument = new G4double[fNumber] ;
fFunction = new G4double[fNumber] ;
fSecondDerivative = new G4double[fNumber] ;
G4double* u = new G4double[fNumber - 1] ;
for(i=0;i<fNumber;i++)
@@ -81,13 +81,13 @@ G4DataInterpolation::G4DataInterpolation( G4double pX[],
else
{
fSecondDerivative[0] = -0.5 ;
u[0] = (3.0/(fArgument[1]-fArgument[0]))*
((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0]) -
pFirstDerStart) ;
u[0] = (3.0/(fArgument[1]-fArgument[0]))
* ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0])
- pFirstDerStart) ;
}
// Decomposition loop for tridiagonal algorithm. fSecondDerivative[i] and u[i]
// are used for temporary storage of the decomposed factors.
// Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
// and u[i] are used for temporary storage of the decomposed factors.
for(i=1;i<fNumber-1;i++)
{
@@ -95,7 +95,7 @@ G4DataInterpolation::G4DataInterpolation( G4double pX[],
p = sig*fSecondDerivative[i-1] + 2.0 ;
fSecondDerivative[i] = (sig - 1.0)/p ;
u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) -
(fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ;
(fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ;
u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ;
}
if(pFirstDerFinish > maxDerivative)
@@ -106,57 +106,58 @@ G4DataInterpolation::G4DataInterpolation( G4double pX[],
else
{
qn = 0.5 ;
un =(3.0/(fArgument[fNumber-1]-fArgument[fNumber-2]))*(pFirstDerFinish -
(fFunction[fNumber-1]-fFunction[fNumber-2])/
(fArgument[fNumber-1]-fArgument[fNumber-2])) ;
un = (3.0/(fArgument[fNumber-1]-fArgument[fNumber-2]))
* (pFirstDerFinish - (fFunction[fNumber-1]-fFunction[fNumber-2])
/ (fArgument[fNumber-1]-fArgument[fNumber-2])) ;
}
fSecondDerivative[fNumber-1] = (un - qn*u[fNumber-2])/
(qn*fSecondDerivative[fNumber-2] + 1.0) ;
// The backsubstitution loop for the triagonal algorithm of solving a linear
// system of equations.
// The backsubstitution loop for the triagonal algorithm of solving
// a linear system of equations.
for(k=fNumber-2;k>=0;k--)
for(G4int k=fNumber-2;k>=0;k--)
{
fSecondDerivative[k] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k] ;
fSecondDerivative[k] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k];
}
delete[] u ;
}
////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Destructor deletes dynamically created arrays for data members: fArgument,
// fFunction and fSecondDerivative, all have dimension of fNumber
G4DataInterpolation::~G4DataInterpolation()
{
delete[] fArgument ;
delete[] fFunction ;
if(fSecondDerivative) delete[] fSecondDerivative ;
delete [] fArgument ;
delete [] fFunction ;
if(fSecondDerivative) { delete [] fSecondDerivative; }
}
////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// This function returns the value P(pX), where P(x) is polynom of fNumber-1 degree
// such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1 . This is
// Lagrange's form of interpolation and it is based on Neville's algorithm
// This function returns the value P(pX), where P(x) is polynom of fNumber-1
// degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1.
// This is Lagrange's form of interpolation and it is based on Neville's
// algorithm
G4double
G4DataInterpolation::PolynomInterpolation(G4double pX,
G4double& deltaY ) const
G4double& deltaY ) const
{
G4int i, m, k = 0 ;
G4double mult, diff, difi, deltaLow, deltaUp, cd, y ;
G4int i=0, m=1, k=0 ;
G4double mult=0.0, difi=0.0, deltaLow=0.0, deltaUp=0.0, cd=0.0, y=0.0 ;
G4double* c = new G4double[fNumber] ;
G4double* d = new G4double[fNumber] ;
diff = std::fabs(pX-fArgument[0]) ;
G4double diff = std::fabs(pX-fArgument[0]) ;
for(i=0;i<fNumber;i++)
{
difi = std::fabs(pX-fArgument[i]) ;
if(difi <diff)
{
k = i ;
diff = difi ;
k = i ;
diff = difi ;
}
c[i] = fFunction[i] ;
d[i] = fFunction[i] ;
@@ -166,18 +167,18 @@ G4DataInterpolation::PolynomInterpolation(G4double pX,
{
for(i=0;i<fNumber-m;i++)
{
deltaLow = fArgument[i] - pX ;
deltaUp = fArgument[i+m] - pX ;
cd = c[i+1] - d[i] ;
mult = deltaLow - deltaUp ;
if(mult == 0.0)
{
G4Exception
("Coincident nodes in G4DataInterpolation::PolynomInterpolation") ;
}
mult = cd/mult ;
d[i] = deltaUp*mult ;
c[i] = deltaLow*mult ;
deltaLow = fArgument[i] - pX ;
deltaUp = fArgument[i+m] - pX ;
cd = c[i+1] - d[i] ;
mult = deltaLow - deltaUp ;
if (!(mult != 0.0))
{
G4Exception("G4DataInterpolation::PolynomInterpolation()",
"Error", FatalException, "Coincident nodes !") ;
}
mult = cd/mult ;
d[i] = deltaUp*mult ;
c[i] = deltaLow*mult ;
}
y += (deltaY = (2*k < (fNumber - m -1) ? c[k+1] : d[k--] )) ;
}
@@ -187,19 +188,19 @@ G4DataInterpolation::PolynomInterpolation(G4double pX,
return y ;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//
// Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1] , this
// Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this
// function calculates an array of coefficients. The coefficients don't provide
// usually (fNumber>10) better accuracy for polynom interpolation, as compared with
// PolynomInterpolation function. They could be used instead for derivate
// usually (fNumber>10) better accuracy for polynom interpolation, as compared
// with PolynomInterpolation function. They could be used instead for derivate
// calculations and some other applications.
void
G4DataInterpolation::PolIntCoefficient( G4double cof[]) const
{
G4int i, j ;
G4double factor, reducedY, mult ;
G4int i=0, j=0 ;
G4double factor=fNumber, reducedY=0.0, mult=1.0 ;
G4double* tempArgument = new G4double[fNumber] ;
for(i=0;i<fNumber;i++)
@@ -212,7 +213,7 @@ G4DataInterpolation::PolIntCoefficient( G4double cof[]) const
{
for(j=fNumber-1-i;j<fNumber-1;j++)
{
tempArgument[j] -= fArgument[i]*tempArgument[j+1] ;
tempArgument[j] -= fArgument[i]*tempArgument[j+1] ;
}
tempArgument[fNumber-1] -= fArgument[i] ;
}
@@ -221,54 +222,51 @@ G4DataInterpolation::PolIntCoefficient( G4double cof[]) const
factor = fNumber ;
for(j=fNumber-1;j>=1;j--)
{
factor = j*tempArgument[j] + factor*fArgument[i] ;
factor = j*tempArgument[j] + factor*fArgument[i] ;
}
reducedY = fFunction[i]/factor ;
mult = 1.0 ;
for(j=fNumber-1;j>=0;j--)
{
cof[j] += mult*reducedY ;
mult = tempArgument[j] + mult*fArgument[i] ;
cof[j] += mult*reducedY ;
mult = tempArgument[j] + mult*fArgument[i] ;
}
}
delete[] tempArgument ;
}
///////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// The function returns diagonal rational function (Bulirsch and Stoer algorithm
// of Neville type) Pn(x)/Qm(x) where P and Q are polynoms.
// Tests showed the method is not stable and hasn't advantage if compared with
// polynomial interpolation ?!
// The function returns diagonal rational function (Bulirsch and Stoer
// algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms.
// Tests showed the method is not stable and hasn't advantage if compared
// with polynomial interpolation ?!
G4double
G4DataInterpolation::RationalPolInterpolation(G4double pX,
G4double& deltaY ) const
G4double& deltaY ) const
{
G4int i, m, k = 0 ;
G4int i=0, m=1, k=0 ;
const G4double tolerance = 1.6e-24 ;
G4double mult, difi, diff, cd, y, cof ;
G4double mult=0.0, difi=0.0, cd=0.0, y=0.0, cof=0.0 ;
G4double* c = new G4double[fNumber] ;
G4double* d = new G4double[fNumber] ;
diff = std::fabs(pX-fArgument[0]) ;
G4double diff = std::fabs(pX-fArgument[0]) ;
for(i=0;i<fNumber;i++)
{
difi = std::fabs(pX-fArgument[i]) ;
if(difi == 0.0)
if (!(difi != 0.0))
{
y = fFunction[i] ;
deltaY = 0.0 ;
y = fFunction[i] ;
deltaY = 0.0 ;
delete[] c ;
delete[] d ;
return y ;
return y ;
}
else if(difi < diff)
{
k = i ;
diff = difi ;
k = i ;
diff = difi ;
}
c[i] = fFunction[i] ;
d[i] = fFunction[i] + tolerance ; // to prevent rare zero/zero cases
@@ -278,17 +276,18 @@ G4DataInterpolation::RationalPolInterpolation(G4double pX,
{
for(i=0;i<fNumber-m;i++)
{
cd = c[i+1] - d[i] ;
difi = fArgument[i+m] - pX ;
cof = (fArgument[i] - pX)*d[i]/difi ;
mult = cof - c[i+1] ;
if(mult == 0.0) // function to be interpolated has pole at pX
{
G4Exception("Error in G4DataInterpolation::RationalPolInterpolation") ;
}
mult = cd/mult ;
d[i] = c[i+1]*mult ;
c[i] = cof*mult ;
cd = c[i+1] - d[i] ;
difi = fArgument[i+m] - pX ;
cof = (fArgument[i] - pX)*d[i]/difi ;
mult = cof - c[i+1] ;
if (!(mult != 0.0)) // function to be interpolated has pole at pX
{
G4Exception("G4DataInterpolation::RationalPolInterpolation()",
"Error", FatalException, "Coincident nodes !") ;
}
mult = cd/mult ;
d[i] = c[i+1]*mult ;
c[i] = cof*mult ;
}
y += (deltaY = (2*k < (fNumber - m - 1) ? c[k+1] : d[k--] )) ;
}
@@ -298,44 +297,41 @@ G4DataInterpolation::RationalPolInterpolation(G4double pX,
return y ;
}
//////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Cubic spline interpolation in point pX for function given by the table:
// fArgument, fFunction. The constructor, which creates fSecondDerivative, must be
// called before. The function works optimal, if sequential calls are in random
// values of pX.
// fArgument, fFunction. The constructor, which creates fSecondDerivative,
// must be called before. The function works optimal, if sequential calls
// are in random values of pX.
G4double
G4DataInterpolation::CubicSplineInterpolation(G4double pX) const
{
G4int kLow, kHigh, k ;
G4double deltaHL, a, b ;
G4int kLow=0, kHigh=fNumber-1, k=0 ;
// Searching in the table by means of bisection method.
// fArgument must be monotonic, either increasing or decreasing
kLow = 0 ;
kHigh = fNumber - 1 ;
while((kHigh - kLow) > 1)
{
k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
if(fArgument[k] > pX)
{
kHigh = k ;
kHigh = k ;
}
else
{
kLow = k ;
kLow = k ;
}
} // kLow and kHigh now bracket the input value of pX
deltaHL = fArgument[kHigh] - fArgument[kLow] ;
if(deltaHL == 0.0)
G4double deltaHL = fArgument[kHigh] - fArgument[kLow] ;
if (!(deltaHL != 0.0))
{
G4Exception(
"Bad fArgument input in G4DataInterpolation::CubicSplineInterpolation") ;
G4Exception("G4DataInterpolation::CubicSplineInterpolation()",
"Error", FatalException, "Bad fArgument input !") ;
}
a = (fArgument[kHigh] - pX)/deltaHL ;
b = (pX - fArgument[kLow])/deltaHL ;
G4double a = (fArgument[kHigh] - pX)/deltaHL ;
G4double b = (pX - fArgument[kLow])/deltaHL ;
// Final evaluation of cubic spline polynomial for return
@@ -347,22 +343,21 @@ G4DataInterpolation::CubicSplineInterpolation(G4double pX) const
///////////////////////////////////////////////////////////////////////////
//
// Return cubic spline interpolation in the point pX which is located between
// fArgument[index] and fArgument[index+1]. It is usually called in sequence of
// known from external analysis values of index.
// fArgument[index] and fArgument[index+1]. It is usually called in sequence
// of known from external analysis values of index.
G4double
G4DataInterpolation::FastCubicSpline(G4double pX,
G4int index) const
G4int index) const
{
G4double delta, a, b ;
delta = fArgument[index+1] - fArgument[index] ;
if(delta == 0.0)
G4double delta = fArgument[index+1] - fArgument[index] ;
if (!(delta != 0.0))
{
G4Exception("Bad fArgument input in G4DataInterpolation::FastCubicSpline") ;
G4Exception("G4DataInterpolation::FastCubicSpline()",
"Error", FatalException, "Bad fArgument input !") ;
}
a = (fArgument[index+1] - pX)/delta ;
b = (pX - fArgument[index])/delta ;
G4double a = (fArgument[index+1] - pX)/delta ;
G4double b = (pX - fArgument[index])/delta ;
// Final evaluation of cubic spline polynomial for return
@@ -371,56 +366,52 @@ G4DataInterpolation::FastCubicSpline(G4double pX,
(b*b*b - b)*fSecondDerivative[index+1])*delta*delta/6.0 ;
}
///////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//
// Given argument pX, returns index k, so that pX bracketed by fArgument[k] and
// fArgument[k+1]
// Given argument pX, returns index k, so that pX bracketed by fArgument[k]
// and fArgument[k+1]
G4int
G4DataInterpolation::LocateArgument(G4double pX) const
{
G4int kLow, kHigh, k ;
G4bool ascend ;
kLow = -1 ;
kHigh = fNumber ;
ascend = (fArgument[fNumber-1] >= fArgument[0]) ;
G4int kLow=-1, kHigh=fNumber, k=0 ;
G4bool ascend=(fArgument[fNumber-1] >= fArgument[0]) ;
while((kHigh - kLow) > 1)
{
k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
if(pX >= fArgument[k] == ascend)
{
kLow = k ;
kLow = k ;
}
else
{
kHigh = k ;
kHigh = k ;
}
}
if(pX == fArgument[0])
if (!(pX != fArgument[0]))
{
return 1 ;
}
else if(pX == fArgument[fNumber-1])
else if (!(pX != fArgument[fNumber-1]))
{
return fNumber - 2 ;
}
else return kLow ;
}
/////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Given a value pX, returns a value 'index' such that pX is between fArgument[index]
// and fArgument[index+1]. fArgument MUST BE MONOTONIC, either increasing or
// decreasing. If index = -1 or fNumber, this indicates that pX is out of range.
// The value index on input is taken as the initial approximation for index on
// output.
// Given a value pX, returns a value 'index' such that pX is between
// fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC,
// either increasing or decreasing. If index = -1 or fNumber, this indicates
// that pX is out of range. The value index on input is taken as the initial
// approximation for index on output.
void
G4DataInterpolation::CorrelatedSearch( G4double pX,
G4int& index ) const
G4int& index ) const
{
G4int kHigh, k, Increment ;
G4int kHigh=0, k=0, Increment=0 ;
// ascend = true for ascending order of table, false otherwise
G4bool ascend = (fArgument[fNumber-1] >= fArgument[0]) ;
if(index < 0 || index > fNumber-1)
@@ -433,45 +424,45 @@ G4DataInterpolation::CorrelatedSearch( G4double pX,
Increment = 1 ; // What value would be the best ?
if((pX >= fArgument[index]) == ascend)
{
if(index == fNumber -1)
{
index = fNumber ;
return ;
}
kHigh = index + 1 ;
while((pX >= fArgument[kHigh]) == ascend)
{
index = kHigh ;
Increment += Increment ; // double the Increment
kHigh = index + Increment ;
if(kHigh > (fNumber - 1))
{
kHigh = fNumber ;
break ;
}
}
if(index == fNumber -1)
{
index = fNumber ;
return ;
}
kHigh = index + 1 ;
while((pX >= fArgument[kHigh]) == ascend)
{
index = kHigh ;
Increment += Increment ; // double the Increment
kHigh = index + Increment ;
if(kHigh > (fNumber - 1))
{
kHigh = fNumber ;
break ;
}
}
}
else
{
if(index == 0)
{
index = -1 ;
return ;
}
if(index == 0)
{
index = -1 ;
return ;
}
kHigh = index-- ;
while((pX < fArgument[index]) == ascend)
{
kHigh = index ;
Increment <<= 1 ; // double the Increment
if(Increment >= kHigh)
{
index = -1 ;
break ;
}
else
{
index = kHigh - Increment ;
}
kHigh = index ;
Increment <<= 1 ; // double the Increment
if(Increment >= kHigh)
{
index = -1 ;
break ;
}
else
{
index = kHigh - Increment ;
}
}
} // Value bracketed
}
@@ -482,18 +473,18 @@ G4DataInterpolation::CorrelatedSearch( G4double pX,
k = (kHigh + index) >> 1 ;
if((pX >= fArgument[k]) == ascend)
{
index = k ;
index = k ;
}
else
{
kHigh = k ;
kHigh = k ;
}
}
if(pX == fArgument[fNumber-1])
if (!(pX != fArgument[fNumber-1]))
{
index = fNumber - 2 ;
}
if(pX == fArgument[0])
if (!(pX != fArgument[0]))
{
index = 0 ;
}
@@ -502,4 +493,4 @@ G4DataInterpolation::CorrelatedSearch( G4double pX,
//
//
//////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
source/global/HEPNumerics/src/G4GaussChebyshevQ.cc
+8 -11
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4GaussChebyshevQ.cc,v 1.4 2004/11/12 17:38:32 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4GaussChebyshevQ.cc,v 1.5 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4GaussChebyshevQ.hh"
@@ -31,15 +31,14 @@
// Constructor for Gauss-Chebyshev quadrature method
G4GaussChebyshevQ::G4GaussChebyshevQ( function pFunction ,
G4int nChebyshev )
G4int nChebyshev )
: G4VGaussianQuadrature(pFunction)
{
G4int i ;
fNumber = nChebyshev ; // Try to reduce fNumber twice ??
G4double cof = pi/fNumber ;
fAbscissa = new G4double[fNumber] ;
fWeight = new G4double[fNumber] ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
@@ -51,7 +50,6 @@ G4GaussChebyshevQ::G4GaussChebyshevQ( function pFunction ,
G4GaussChebyshevQ::~G4GaussChebyshevQ()
{
;
}
// -------------------------------------------------------------------------------
@@ -62,12 +60,11 @@ G4GaussChebyshevQ::~G4GaussChebyshevQ()
G4double
G4GaussChebyshevQ::Integral(G4double a, G4double b) const
{
G4int i ;
G4double xDiff, xMean, dx, integral = 0.0 ;
G4double xDiff=0.5*(b - a),
xMean=0.5*(a + b),
dx=0.0, integral=0.0 ;
xMean = 0.5*(a + b) ;
xDiff = 0.5*(b - a) ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
dx = xDiff*fAbscissa[i] ;
integral += fWeight[i]*fFunction(xMean + dx) ;
source/global/HEPNumerics/src/G4GaussHermiteQ.cc
+37 -35
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4GaussHermiteQ.cc,v 1.5 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4GaussHermiteQ.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4GaussHermiteQ.hh"
@@ -31,17 +31,17 @@
//
// Constructor for Gauss-Hermite
G4GaussHermiteQ::G4GaussHermiteQ( function pFunction,
G4int nHermite )
G4GaussHermiteQ::G4GaussHermiteQ ( function pFunction,
G4int nHermite )
: G4VGaussianQuadrature(pFunction)
{
const G4double tolerance = 1.0e-12 ;
const G4int maxNumber = 12 ;
G4int i, j, k ;
G4int i=1, j=1, k=1 ;
G4double newton=0.;
G4double newton1, temp1, temp2, temp3, temp ;
G4double piInMinusQ = std::pow(pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
G4double newton1=0.0, temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0 ;
G4double piInMinusQ = std::pow(pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
fNumber = (nHermite +1)/2 ;
fAbscissa = new G4double[fNumber] ;
@@ -51,46 +51,49 @@ G4GaussHermiteQ::G4GaussHermiteQ( function pFunction,
{
if(i == 1)
{
newton = std::sqrt((G4double)(2*nHermite + 1)) -
1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
newton = std::sqrt((G4double)(2*nHermite + 1)) -
1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
}
else if(i == 2)
{
newton -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton ;
newton -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton ;
}
else if(i == 3)
{
newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
newton = 1.86002*newton - 0.86002*fAbscissa[0] ;
}
else if(i == 4)
{
newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
newton = 1.91001*newton - 0.91001*fAbscissa[1] ;
}
else
{
newton = 2.0*newton - fAbscissa[i - 3] ;
newton = 2.0*newton - fAbscissa[i - 3] ;
}
for(k=1;k<=maxNumber;k++)
{
temp1 = piInMinusQ ;
temp2 = 0.0 ;
for(j=1;j<=nHermite;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = newton*std::sqrt(2.0/j)*temp2 - std::sqrt(((G4double)(j - 1))/j)*temp3 ;
}
temp = std::sqrt((G4double)2*nHermite)*temp2 ;
newton1 = newton ;
newton = newton1 - temp1/temp ;
temp1 = piInMinusQ ;
temp2 = 0.0 ;
for(j=1;j<=nHermite;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = newton*std::sqrt(2.0/j)*temp2
- std::sqrt(((G4double)(j - 1))/j)*temp3 ;
}
temp = std::sqrt((G4double)2*nHermite)*temp2 ;
newton1 = newton ;
newton = newton1 - temp1/temp ;
if(std::fabs(newton - newton1) <= tolerance)
{
break ;
}
{
break ;
}
}
if(k > maxNumber)
{
G4Exception("Too many iterations in Gauss-Hermite constructor") ;
G4Exception("G4GaussHermiteQ::G4GaussHermiteQ()",
"OutOfRange", FatalException,
"Too many iterations in Gauss-Hermite constructor.") ;
}
fAbscissa[i-1] = newton ;
fWeight[i-1] = 2.0/(temp*temp) ;
@@ -100,17 +103,16 @@ G4GaussHermiteQ::G4GaussHermiteQ( function pFunction,
// ----------------------------------------------------------
//
// Gauss-Hermite method for integration of std::exp(-x*x)*nFunction(x) from minus infinity
// to plus infinity .
// Gauss-Hermite method for integration of std::exp(-x*x)*nFunction(x)
// from minus infinity to plus infinity .
G4double
G4GaussHermiteQ::Integral() const
G4double G4GaussHermiteQ::Integral() const
{
G4int i ;
G4double integral = 0.0 ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
integral += fWeight[i]*(fFunction(fAbscissa[i]) + fFunction(-fAbscissa[i])) ;
integral += fWeight[i]*(fFunction(fAbscissa[i])
+ fFunction(-fAbscissa[i])) ;
}
return integral ;
}
source/global/HEPNumerics/src/G4GaussJacobiQ.cc
+66 -61
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4GaussJacobiQ.cc,v 1.5 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4GaussJacobiQ.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4GaussJacobiQ.hh"
@@ -33,18 +33,21 @@
//
G4GaussJacobiQ::G4GaussJacobiQ( function pFunction,
G4double alpha,
G4double alpha,
G4double beta,
G4int nJacobi )
G4int nJacobi )
: G4VGaussianQuadrature(pFunction)
{
const G4double tolerance = 1.0e-12 ;
const G4double maxNumber = 12 ;
G4int i, k, j ;
G4int i=1, k=1 ;
G4double root=0.;
G4double alphaBeta, alphaReduced, betaReduced, root1, root2, root3 ;
G4double a, b, c, newton1, newton2, newton3, newton, temp, rootTemp ;
G4double alphaBeta=0.0, alphaReduced=0.0, betaReduced=0.0,
root1=0.0, root2=0.0, root3=0.0 ;
G4double a=0.0, b=0.0, c=0.0,
newton1=0.0, newton2=0.0, newton3=0.0, newton=0.0,
temp=0.0, rootTemp=0.0 ;
fNumber = nJacobi ;
fAbscissa = new G4double[fNumber] ;
@@ -54,98 +57,100 @@ G4GaussJacobiQ::G4GaussJacobiQ( function pFunction,
{
if (i == 1)
{
alphaReduced = alpha/nJacobi ;
betaReduced = beta/nJacobi ;
root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
0.767999*alphaReduced/nJacobi) ;
root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
0.451998*alphaReduced*alphaReduced+0.83001*alphaReduced*betaReduced ;
root = 1.0-root1/root2 ;
alphaReduced = alpha/nJacobi ;
betaReduced = beta/nJacobi ;
root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
0.767999*alphaReduced/nJacobi) ;
root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced
+ 0.451998*alphaReduced*alphaReduced
+ 0.83001*alphaReduced*betaReduced ;
root = 1.0-root1/root2 ;
}
else if (i == 2)
{
root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
root -= (1.0-root)*root1*root2*root3 ;
root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
root -= (1.0-root)*root1*root2*root3 ;
}
else if (i == 3)
{
root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
root -= (fAbscissa[0]-root)*root1*root2*root3 ;
root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
root -= (fAbscissa[0]-root)*root1*root2*root3 ;
}
else if (i == nJacobi-1)
{
root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
}
else if (i == nJacobi)
{
root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
}
else
{
root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
}
alphaBeta = alpha + beta ;
for (k=1;k<=maxNumber;k++)
{
temp = 2.0 + alphaBeta ;
newton1 = (alpha-beta+temp*root)/2.0 ;
newton2 = 1.0 ;
for (j=2;j<=nJacobi;j++)
{
newton3 = newton2 ;
newton2 = newton1 ;
temp = 2*j+alphaBeta ;
a = 2*j*(j+alphaBeta)*(temp-2.0) ;
b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
newton1 = (b*newton2-c*newton3)/a ;
}
newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
(temp*(1.0 - root*root)) ;
rootTemp = root ;
root = rootTemp - newton1/newton ;
if (std::fabs(root-rootTemp) <= tolerance)
{
break ;
}
temp = 2.0 + alphaBeta ;
newton1 = (alpha-beta+temp*root)/2.0 ;
newton2 = 1.0 ;
for (G4int j=2;j<=nJacobi;j++)
{
newton3 = newton2 ;
newton2 = newton1 ;
temp = 2*j+alphaBeta ;
a = 2*j*(j+alphaBeta)*(temp-2.0) ;
b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
newton1 = (b*newton2-c*newton3)/a ;
}
newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
(temp*(1.0 - root*root)) ;
rootTemp = root ;
root = rootTemp - newton1/newton ;
if (std::fabs(root-rootTemp) <= tolerance)
{
break ;
}
}
if (k > maxNumber)
{
G4Exception("Too many iterations in G4GaussJacobiQ::G4GaussJacobiQ") ;
G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange",
FatalException, "Too many iterations in constructor.") ;
}
fAbscissa[i-1] = root ;
fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
GammaLogarithm((G4double)(beta+nJacobi)) -
GammaLogarithm((G4double)(nJacobi+1.0)) -
GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
*temp*std::pow(2.0,alphaBeta)/(newton*newton2) ;
GammaLogarithm((G4double)(beta+nJacobi)) -
GammaLogarithm((G4double)(nJacobi+1.0)) -
GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
*temp*std::pow(2.0,alphaBeta)/(newton*newton2) ;
}
}
// ----------------------------------------------------------
//
// Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*pFunction(x)
// Gauss-Jacobi method for integration of
// ((1-x)^alpha)*((1+x)^beta)*pFunction(x)
// from minus unit to plus unit .
G4double
G4GaussJacobiQ::Integral() const
{
G4int i ;
G4double integral = 0.0 ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
integral += fWeight[i]*fFunction(fAbscissa[i]) ;
}
source/global/HEPNumerics/src/G4GaussLaguerreQ.cc
+42 -36
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4GaussLaguerreQ.cc,v 1.5 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4GaussLaguerreQ.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4GaussLaguerreQ.hh"
@@ -31,21 +31,22 @@
// ------------------------------------------------------------
//
// Constructor for Gauss-Laguerre quadrature method: integral from zero to
// infinity of std::pow(x,alpha)*std::exp(-x)*f(x). The value of nLaguerre sets the accuracy.
// infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
// The value of nLaguerre sets the accuracy.
// The constructor creates arrays fAbscissa[0,..,nLaguerre-1] and
// fWeight[0,..,nLaguerre-1] .
//
G4GaussLaguerreQ::G4GaussLaguerreQ( function pFunction,
G4double alpha,
G4int nLaguerre )
G4double alpha,
G4int nLaguerre )
: G4VGaussianQuadrature(pFunction)
{
const G4double tolerance = 1.0e-10 ;
const G4int maxNumber = 12 ;
G4int i, j, k ;
G4double newton=0.;
G4double newton1, temp1, temp2, temp3, temp, cofi ;
G4int i=1, k=1 ;
G4double newton=0.0, newton1=0.0,
temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0, cofi=0.0 ;
fNumber = nLaguerre ;
fAbscissa = new G4double[fNumber] ;
@@ -55,60 +56,65 @@ G4GaussLaguerreQ::G4GaussLaguerreQ( function pFunction,
{
if(i == 1)
{
newton = (1.0 + alpha)*(3.0 + 0.92*alpha)/(1.0 + 2.4*fNumber + 1.8*alpha) ;
newton = (1.0 + alpha)*(3.0 + 0.92*alpha)
/ (1.0 + 2.4*fNumber + 1.8*alpha) ;
}
else if(i == 2)
{
newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
newton += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
}
else
{
cofi = i - 2 ;
newton += ((1.0+2.55*cofi)/(1.9*cofi) + 1.26*cofi*alpha/(1.0+3.5*cofi))*
(newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
cofi = i - 2 ;
newton += ((1.0+2.55*cofi)/(1.9*cofi)
+ 1.26*cofi*alpha/(1.0+3.5*cofi))
* (newton - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
}
for(k=1;k<=maxNumber;k++)
{
temp1 = 1.0 ;
temp2 = 0.0 ;
for(j=1;j<=fNumber;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = ((2*j - 1 + alpha - newton)*temp2 - (j - 1 + alpha)*temp3)/j ;
}
temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
newton1 = newton ;
newton = newton1 - temp1/temp ;
temp1 = 1.0 ;
temp2 = 0.0 ;
for(G4int j=1;j<=fNumber;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = ((2*j - 1 + alpha - newton)*temp2
- (j - 1 + alpha)*temp3)/j ;
}
temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton ;
newton1 = newton ;
newton = newton1 - temp1/temp ;
if(std::fabs(newton - newton1) <= tolerance)
{
break ;
}
{
break ;
}
}
if(k > maxNumber)
{
G4Exception("Too many iterations in Gauss-Laguerre constructor") ;
G4Exception("G4GaussLaguerreQ::G4GaussLaguerreQ()",
"OutOfRange", FatalException,
"Too many iterations in Gauss-Laguerre constructor") ;
}
fAbscissa[i-1] = newton ;
fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) -
GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber)
- GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
}
}
// -----------------------------------------------------------------
//
// Gauss-Laguerre method for integration of std::pow(x,alpha)*std::exp(-x)*pFunction(x)
// from zero up to infinity. pFunction is evaluated in fNumber points for which
// fAbscissa[i] and fWeight[i] arrays were created in
// Gauss-Laguerre method for integration of
// std::pow(x,alpha)*std::exp(-x)*pFunction(x)
// from zero up to infinity. pFunction is evaluated in fNumber points
// for which fAbscissa[i] and fWeight[i] arrays were created in
// G4VGaussianQuadrature(double,int) constructor
G4double
G4GaussLaguerreQ::Integral() const
{
G4int i ;
G4double integral = 0.0 ;
for(i=0;i<fNumber;i++)
for(G4int i=0;i<fNumber;i++)
{
integral += fWeight[i]*fFunction(fAbscissa[i]) ;
}
source/global/HEPNumerics/src/G4GaussLegendreQ.cc
+120 -131
View File
@@ -21,88 +21,83 @@
// ********************************************************************
//
//
// $Id: G4GaussLegendreQ.cc,v 1.4 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4GaussLegendreQ.cc,v 1.6 2005/03/16 13:45:16 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
#include "G4GaussLegendreQ.hh"
G4GaussLegendreQ::G4GaussLegendreQ( function pFunction )
: G4VGaussianQuadrature(pFunction)
{
;
}
// ----------------------------------------------------------------------------
// --------------------------------------------------------------------------
//
// Constructor for GaussLegendre quadrature method. The value nLegendre set the
// accuracy required, i.e the number of points where the function pFunction will
// be evaluated during integration. The constructor creates the arrays for
// abscissas and weights that used in Gauss-Legendre quadrature method.
// Constructor for GaussLegendre quadrature method. The value nLegendre sets
// the accuracy required, i.e the number of points where the function pFunction
// will be evaluated during integration. The constructor creates the arrays for
// abscissas and weights that are used in Gauss-Legendre quadrature method.
// The values a and b are the limits of integration of the pFunction.
// nLegendre MUST BE EVEN !!!
G4GaussLegendreQ::G4GaussLegendreQ( function pFunction,
G4int nLegendre )
G4int nLegendre )
: G4VGaussianQuadrature(pFunction)
{
const G4double tolerance = 1.6e-10 ;
G4int i, j, k = nLegendre ;
G4int k = nLegendre ;
fNumber = (nLegendre + 1)/2 ;
if(2*fNumber != k)
{
G4Exception("Invalid nLegendre in G4GaussLegendreQ::G4GaussLegendreQ") ;
G4Exception("G4GaussLegendreQ::G4GaussLegendreQ()", "InvalidCall",
FatalException, "Invalid nLegendre argument !") ;
}
G4double newton, newton1, temp1, temp2, temp3, temp ;
G4double newton=0.0, newton1=0.0,
temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0 ;
fAbscissa = new G4double[fNumber] ;
fWeight = new G4double[fNumber] ;
for(i=1;i<=fNumber;i++) // Loop over the desired roots
for(G4int i=1;i<=fNumber;i++) // Loop over the desired roots
{
newton = std::cos(pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation
do
newton = std::cos(pi*(i - 0.25)/(k + 0.5)) ; // Initial root
do // approximation
{ // loop of Newton's method
temp1 = 1.0 ;
temp2 = 0.0 ;
for(j=1;j<=k;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
}
temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
newton1 = newton ;
newton = newton1 - temp1/temp ; // Newton's method
temp1 = 1.0 ;
temp2 = 0.0 ;
for(G4int j=1;j<=k;j++)
{
temp3 = temp2 ;
temp2 = temp1 ;
temp1 = ((2.0*j - 1.0)*newton*temp2 - (j - 1.0)*temp3)/j ;
}
temp = k*(newton*temp1 - temp2)/(newton*newton - 1.0) ;
newton1 = newton ;
newton = newton1 - temp1/temp ; // Newton's method
}
while(std::fabs(newton - newton1) > tolerance) ;
fAbscissa[fNumber-i] = newton ;
fWeight[fNumber-i] = 2.0/((1.0 - newton*newton)*temp*temp) ;
}
}
// -------------------------------------------------------------------------------
// --------------------------------------------------------------------------
//
// Returns the integral of the function to be pointed by fFunction between a and b,
// by 2*fNumber point Gauss-Legendre integration: the function is evaluated exactly
// 2*fNumber Times at interior points in the range of integration. Since the weights
// and abscissas are, in this case, symmetric around the midpoint of the range of
// integration, there are actually only fNumber distinct values of each.
// Returns the integral of the function to be pointed by fFunction between a
// and b, by 2*fNumber point Gauss-Legendre integration: the function is
// evaluated exactly 2*fNumber times at interior points in the range of
// integration. Since the weights and abscissas are, in this case, symmetric
// around the midpoint of the range of integration, there are actually only
// fNumber distinct values of each.
G4double
G4GaussLegendreQ::Integral(G4double a, G4double b) const
{
G4int i ;
G4double xDiff, xMean, dx, integral ;
xMean = 0.5*(a + b) ;
xDiff = 0.5*(b - a) ;
integral = 0.0 ;
for(i=0;i<fNumber;i++)
G4double xMean = 0.5*(a + b),
xDiff = 0.5*(b - a),
integral = 0.0, dx = 0.0 ;
for(G4int i=0;i<fNumber;i++)
{
dx = xDiff*fAbscissa[i] ;
integral += fWeight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
@@ -110,33 +105,31 @@ G4GaussLegendreQ::Integral(G4double a, G4double b) const
return integral *= xDiff ;
}
// -------------------------------------------------------------------------------
// --------------------------------------------------------------------------
//
// Returns the integral of the function to be pointed by fFunction between a and b,
// by ten point Gauss-Legendre integration: the function is evaluated exactly
// ten Times at interior points in the range of integration. Since the weights
// and abscissas are, in this case, symmetric around the midpoint of the range of
// integration, there are actually only five distinct values of each
// Returns the integral of the function to be pointed by fFunction between a
// and b, by ten point Gauss-Legendre integration: the function is evaluated
// exactly ten times at interior points in the range of integration. Since the
// weights and abscissas are, in this case, symmetric around the midpoint of
// the range of integration, there are actually only five distinct values of
// each.
G4double
G4GaussLegendreQ::QuickIntegral(G4double a, G4double b) const
{
G4int i ;
G4double xDiff, xMean, dx, integral ;
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
0.679409568299024, 0.865063366688985,
0.973906528517172 } ;
0.973906528517172 } ;
static G4double weight[] = { 0.295524224714753, 0.269266719309996,
0.219086362515982, 0.149451349150581,
0.066671344308688 } ;
xMean = 0.5*(a + b) ;
xDiff = 0.5*(b - a) ;
integral = 0.0 ;
for(i=0;i<5;i++)
0.066671344308688 } ;
G4double xMean = 0.5*(a + b),
xDiff = 0.5*(b - a),
integral = 0.0, dx = 0.0 ;
for(G4int i=0;i<5;i++)
{
dx = xDiff*abscissa[i] ;
integral += weight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
@@ -144,21 +137,18 @@ G4double
return integral *= xDiff ;
}
// -------------------------------------------------------------------------
//
// Returns the integral of the function to be pointed by fFunction between a and b,
// by 96 point Gauss-Legendre integration: the function is evaluated exactly
// ten Times at interior points in the range of integration. Since the weights
// and abscissas are, in this case, symmetric around the midpoint of the range of
// integration, there are actually only five distinct values of each
// Returns the integral of the function to be pointed by fFunction between a
// and b, by 96 point Gauss-Legendre integration: the function is evaluated
// exactly ten times at interior points in the range of integration. Since the
// weights and abscissas are, in this case, symmetric around the midpoint of
// the range of integration, there are actually only five distinct values of
// each.
G4double
G4GaussLegendreQ::AccurateIntegral(G4double a, G4double b) const
{
G4int i ;
G4double xDiff, xMean, dx, integral ;
{
// From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
static
@@ -167,77 +157,76 @@ G4double
0.081297495464425558994, 0.113695850110665920911,
0.145973714654896941989, 0.178096882367618602759, // 6
0.210031310460567203603, 0.241743156163840012328,
0.273198812591049141487, 0.304364944354496353024,
0.335208522892625422616, 0.365696861472313635031, // 12
0.395797649828908603285, 0.425478988407300545365,
0.454709422167743008636, 0.483457973920596359768,
0.511694177154667673586, 0.539388108324357436227, // 18
0.566510418561397168404, 0.593032364777572080684,
0.618925840125468570386, 0.644163403784967106798,
0.668718310043916153953, 0.692564536642171561344, // 24
0.715676812348967626225, 0.738030643744400132851,
0.759602341176647498703, 0.780369043867433217604,
0.800308744139140817229, 0.819400310737931675539, // 30
0.837623511228187121494, 0.854959033434601455463,
0.871388505909296502874, 0.886894517402420416057,
0.901460635315852341319, 0.915071423120898074206, // 36
0.927712456722308690965, 0.939370339752755216932,
0.950032717784437635756, 0.959688291448742539300,
0.968326828463264212174, 0.975939174585136466453, // 42
0.982517263563014677447, 0.988054126329623799481,
0.992543900323762624572, 0.995981842987209290650,
0.998364375863181677724, 0.999689503883230766828 // 48
0.210031310460567203603, 0.241743156163840012328,
0.273198812591049141487, 0.304364944354496353024,
0.335208522892625422616, 0.365696861472313635031, // 12
0.395797649828908603285, 0.425478988407300545365,
0.454709422167743008636, 0.483457973920596359768,
0.511694177154667673586, 0.539388108324357436227, // 18
0.566510418561397168404, 0.593032364777572080684,
0.618925840125468570386, 0.644163403784967106798,
0.668718310043916153953, 0.692564536642171561344, // 24
0.715676812348967626225, 0.738030643744400132851,
0.759602341176647498703, 0.780369043867433217604,
0.800308744139140817229, 0.819400310737931675539, // 30
0.837623511228187121494, 0.854959033434601455463,
0.871388505909296502874, 0.886894517402420416057,
0.901460635315852341319, 0.915071423120898074206, // 36
0.927712456722308690965, 0.939370339752755216932,
0.950032717784437635756, 0.959688291448742539300,
0.968326828463264212174, 0.975939174585136466453, // 42
0.982517263563014677447, 0.988054126329623799481,
0.992543900323762624572, 0.995981842987209290650,
0.998364375863181677724, 0.999689503883230766828 // 48
} ;
static
G4double weight[] = {
0.032550614492363166242, 0.032516118713868835987,
0.032447163714064269364, 0.032343822568575928429,
0.032206204794030250669, 0.032034456231992663218, // 6
0.031828758894411006535, 0.031589330770727168558,
0.031316425596862355813, 0.031010332586313837423,
0.030671376123669149014, 0.030299915420827593794, // 12
0.029896344136328385984, 0.029461089958167905970,
0.028994614150555236543, 0.028497411065085385646,
0.027970007616848334440, 0.027412962726029242823, // 18
0.026826866725591762198, 0.026212340735672413913,
0.025570036005349361499, 0.024900633222483610288,
0.024204841792364691282, 0.023483399085926219842, // 24
0.022737069658329374001, 0.021966644438744349195,
0.021172939892191298988, 0.020356797154333324595,
0.019519081140145022410, 0.018660679627411467385, // 30
0.017782502316045260838, 0.016885479864245172450,
0.015970562902562291381, 0.015038721026994938006,
0.014090941772314860916, 0.013128229566961572637, // 36
0.012151604671088319635, 0.011162102099838498591,
0.010160770535008415758, 0.009148671230783386633,
0.008126876925698759217, 0.007096470791153865269, // 42
0.006058545504235961683, 0.005014202742927517693,
0.003964554338444686674, 0.002910731817934946408,
0.001853960788946921732, 0.000796792065552012429 // 48
0.032206204794030250669, 0.032034456231992663218, // 6
0.031828758894411006535, 0.031589330770727168558,
0.031316425596862355813, 0.031010332586313837423,
0.030671376123669149014, 0.030299915420827593794, // 12
0.029896344136328385984, 0.029461089958167905970,
0.028994614150555236543, 0.028497411065085385646,
0.027970007616848334440, 0.027412962726029242823, // 18
0.026826866725591762198, 0.026212340735672413913,
0.025570036005349361499, 0.024900633222483610288,
0.024204841792364691282, 0.023483399085926219842, // 24
0.022737069658329374001, 0.021966644438744349195,
0.021172939892191298988, 0.020356797154333324595,
0.019519081140145022410, 0.018660679627411467385, // 30
0.017782502316045260838, 0.016885479864245172450,
0.015970562902562291381, 0.015038721026994938006,
0.014090941772314860916, 0.013128229566961572637, // 36
0.012151604671088319635, 0.011162102099838498591,
0.010160770535008415758, 0.009148671230783386633,
0.008126876925698759217, 0.007096470791153865269, // 42
0.006058545504235961683, 0.005014202742927517693,
0.003964554338444686674, 0.002910731817934946408,
0.001853960788946921732, 0.000796792065552012429 // 48
} ;
xMean = 0.5*(a + b) ;
xDiff = 0.5*(b - a) ;
integral = 0.0 ;
for(i=0;i<48;i++)
G4double xMean = 0.5*(a + b),
xDiff = 0.5*(b - a),
integral = 0.0, dx = 0.0 ;
for(G4int i=0;i<48;i++)
{
dx = xDiff*abscissa[i] ;
integral += weight[i]*(fFunction(xMean + dx) + fFunction(xMean - dx)) ;
}
return integral *= xDiff ;
}
source/global/HEPNumerics/src/G4JTPolynomialSolver.cc
+905
View File
@@ -0,0 +1,905 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
// $Id: G4JTPolynomialSolver.cc,v 1.4 2005/03/21 18:27:40 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
// --------------------------------------------------------------------
// GEANT 4 class source file
//
// G4JTPolynomialSolver
//
// Implementation based on Jenkins-Traub algorithm.
// --------------------------------------------------------------------
#include "G4JTPolynomialSolver.hh"
const G4double G4JTPolynomialSolver::base = 2;
const G4double G4JTPolynomialSolver::eta = DBL_EPSILON;
const G4double G4JTPolynomialSolver::infin = DBL_MAX;
const G4double G4JTPolynomialSolver::smalno = DBL_MIN;
const G4double G4JTPolynomialSolver::are = DBL_EPSILON;
const G4double G4JTPolynomialSolver::mre = DBL_EPSILON;
const G4double G4JTPolynomialSolver::lo = DBL_MIN/DBL_EPSILON ;
G4JTPolynomialSolver::G4JTPolynomialSolver()
{
}
G4JTPolynomialSolver::~G4JTPolynomialSolver()
{
}
G4int G4JTPolynomialSolver::FindRoots(G4double *op, G4int degree,
G4double *zeror, G4double *zeroi)
{
G4double t=0.0, aa=0.0, bb=0.0, cc=0.0, factor=1.0;
G4double max=0.0, min=infin, xxx=0.0, x=0.0, sc=0.0, bnd=0.0;
G4double xm=0.0, ff=0.0, df=0.0, dx=0.0;
G4int cnt=0, nz=0, i=0, j=0, jj=0, l=0, nm1=0, zerok=0;
// Initialization of constants for shift rotation.
//
G4double xx = std::sqrt(0.5);
G4double yy = -xx,
rot = 94.0*deg;
G4double cosr = std::cos(rot),
sinr = std::sin(rot);
n = degree;
// Algorithm fails if the leading coefficient is zero.
//
if (!(op[0] != 0.0)) { return -1; }
// Remove the zeros at the origin, if any.
//
while (!(op[n] != 0.0))
{
j = degree - n;
zeror[j] = 0.0;
zeroi[j] = 0.0;
n--;
}
if (n < 1) { return -1; }
// Allocate buffers here
//
std::vector<G4double> temp(degree+1) ;
std::vector<G4double> pt(degree+1) ;
p.assign(degree+1,0) ;
qp.assign(degree+1,0) ;
k.assign(degree+1,0) ;
qk.assign(degree+1,0) ;
svk.assign(degree+1,0) ;
// Make a copy of the coefficients.
//
for (i=0;i<=n;i++)
{ p[i] = op[i]; }
do
{
if (n == 1) // Start the algorithm for one zero.
{
zeror[degree-1] = -p[1]/p[0];
zeroi[degree-1] = 0.0;
n -= 1;
return degree - n ;
}
if (n == 2) // Calculate the final zero or pair of zeros.
{
Quadratic(p[0],p[1],p[2],&zeror[degree-2],&zeroi[degree-2],
&zeror[degree-1],&zeroi[degree-1]);
n -= 2;
return degree - n ;
}
// Find largest and smallest moduli of coefficients.
//
max = 0.0;
min = infin;
for (i=0;i<=n;i++)
{
x = std::fabs(p[i]);
if (x > max) { max = x; }
if (x != 0.0 && x < min) { min = x; }
}
// Scale if there are large or very small coefficients.
// Computes a scale factor to multiply the coefficients of the
// polynomial. The scaling is done to avoid overflow and to
// avoid undetected underflow interfering with the convergence
// criterion. The factor is a power of the base.
//
sc = lo/min;
if ( sc <= 1.0 && max >= 10.0
|| sc > 1.0 && infin/sc >= max
|| infin/sc >= max && max >= 10 )
{
if (!( sc != 0.0 ))
{ sc = smalno ; }
l = (G4int)(std::log(sc)/std::log(base) + 0.5);
factor = std::pow(base*1.0,l);
if (factor != 1.0)
{
for (i=0;i<=n;i++)
{ p[i] = factor*p[i]; } // Scale polynomial.
}
}
// Compute lower bound on moduli of roots.
//
for (i=0;i<=n;i++)
{
pt[i] = (std::fabs(p[i]));
}
pt[n] = - pt[n];
// Compute upper estimate of bound.
//
x = std::exp((std::log(-pt[n])-std::log(pt[0])) / (G4double)n);
// If Newton step at the origin is better, use it.
//
if (pt[n-1] != 0.0)
{
xm = -pt[n]/pt[n-1];
if (xm < x) { x = xm; }
}
// Chop the interval (0,x) until ff <= 0
//
while (1)
{
xm = x*0.1;
ff = pt[0];
for (i=1;i<=n;i++)
{ ff = ff*xm + pt[i]; }
if (ff <= 0.0) { break; }
x = xm;
}
dx = x;
// Do Newton interation until x converges to two decimal places.
//
while (std::fabs(dx/x) > 0.005)
{
ff = pt[0];
df = ff;
for (i=1;i<n;i++)
{
ff = ff*x + pt[i];
df = df*x + ff;
}
ff = ff*x + pt[n];
dx = ff/df;
x -= dx;
}
bnd = x;
// Compute the derivative as the initial k polynomial
// and do 5 steps with no shift.
//
nm1 = n - 1;
for (i=1;i<n;i++)
{ k[i] = (G4double)(n-i)*p[i]/(G4double)n; }
k[0] = p[0];
aa = p[n];
bb = p[n-1];
zerok = (k[n-1] == 0);
for(jj=0;jj<5;jj++)
{
cc = k[n-1];
if (!zerok) // Use a scaled form of recurrence if k at 0 is nonzero.
{
// Use a scaled form of recurrence if value of k at 0 is nonzero.
//
t = -aa/cc;
for (i=0;i<nm1;i++)
{
j = n-i-1;
k[j] = t*k[j-1]+p[j];
}
k[0] = p[0];
zerok = (std::fabs(k[n-1]) <= std::fabs(bb)*eta*10.0);
}
else // Use unscaled form of recurrence.
{
for (i=0;i<nm1;i++)
{
j = n-i-1;
k[j] = k[j-1];
}
k[0] = 0.0;
zerok = (!(k[n-1] != 0.0));
}
}
// Save k for restarts with new shifts.
//
for (i=0;i<n;i++)
{ temp[i] = k[i]; }
// Loop to select the quadratic corresponding to each new shift.
//
for (cnt = 0;cnt < 20;cnt++)
{
// Quadratic corresponds to a double shift to a
// non-real point and its complex conjugate. The point
// has modulus bnd and amplitude rotated by 94 degrees
// from the previous shift.
//
xxx = cosr*xx - sinr*yy;
yy = sinr*xx + cosr*yy;
xx = xxx;
sr = bnd*xx;
si = bnd*yy;
u = -2.0 * sr;
v = bnd;
ComputeFixedShiftPolynomial(20*(cnt+1),&nz);
if (nz != 0)
{
// The second stage jumps directly to one of the third
// stage iterations and returns here if successful.
// Deflate the polynomial, store the zero or zeros and
// return to the main algorithm.
//
j = degree - n;
zeror[j] = szr;
zeroi[j] = szi;
n -= nz;
for (i=0;i<=n;i++)
{ p[i] = qp[i]; }
if (nz != 1)
{
zeror[j+1] = lzr;
zeroi[j+1] = lzi;
}
break;
}
else
{
// If the iteration is unsuccessful another quadratic
// is chosen after restoring k.
//
for (i=0;i<n;i++)
{
k[i] = temp[i];
}
}
}
}
while (nz != 0); // End of initial DO loop
// Return with failure if no convergence with 20 shifts.
//
return degree - n;
}
void G4JTPolynomialSolver::ComputeFixedShiftPolynomial(G4int l2, G4int *nz)
{
// Computes up to L2 fixed shift k-polynomials, testing for convergence
// in the linear or quadratic case. Initiates one of the variable shift
// iterations and returns with the number of zeros found.
G4double svu=0.0, svv=0.0, ui=0.0, vi=0.0, s=0.0;
G4double betas=0.25, betav=0.25, oss=sr, ovv=v,
ss=0.0, vv=0.0, ts=1.0, tv=1.0;
G4double ots=0.0, otv=0.0;
G4double tvv=1.0, tss=1.0;
G4int type=0, i=0, j=0, iflag=0, vpass=0, spass=0, vtry=0, stry=0;
*nz = 0;
// Evaluate polynomial by synthetic division.
//
QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
ComputeScalarFactors(&type);
for (j=0;j<l2;j++)
{
// Calculate next k polynomial and estimate v.
//
ComputeNextPolynomial(&type);
ComputeScalarFactors(&type);
ComputeNewEstimate(type,&ui,&vi);
vv = vi;
// Estimate s.
//
ss = 0.0;
if (k[n-1] != 0.0) { ss = -p[n]/k[n-1]; }
tv = 1.0;
ts = 1.0;
if (j == 0 || type == 3)
{
ovv = vv;
oss = ss;
otv = tv;
ots = ts;
continue;
}
// Compute relative measures of convergence of s and v sequences.
//
if (vv != 0.0) { tv = std::fabs((vv-ovv)/vv); }
if (ss != 0.0) { ts = std::fabs((ss-oss)/ss); }
// If decreasing, multiply two most recent convergence measures.
tvv = 1.0;
if (tv < otv) { tvv = tv*otv; }
tss = 1.0;
if (ts < ots) { tss = ts*ots; }
// Compare with convergence criteria.
vpass = (tvv < betav);
spass = (tss < betas);
if (!(spass || vpass))
{
ovv = vv;
oss = ss;
otv = tv;
ots = ts;
continue;
}
// At least one sequence has passed the convergence test.
// Store variables before iterating.
//
svu = u;
svv = v;
for (i=0;i<n;i++)
{
svk[i] = k[i];
}
s = ss;
// Choose iteration according to the fastest converging sequence.
//
vtry = 0;
stry = 0;
if (spass && (!vpass) || tss < tvv)
{
RealPolynomialIteration(&s,nz,&iflag);
if (*nz > 0) { return; }
// Linear iteration has failed. Flag that it has been
// tried and decrease the convergence criterion.
//
stry = 1;
betas *=0.25;
if (iflag == 0) { goto _restore_variables; }
// If linear iteration signals an almost double real
// zero attempt quadratic iteration.
//
ui = -(s+s);
vi = s*s;
}
_quadratic_iteration:
do
{
QuadraticPolynomialIteration(&ui,&vi,nz);
if (*nz > 0) { return; }
// Quadratic iteration has failed. Flag that it has
// been tried and decrease the convergence criterion.
//
vtry = 1;
betav *= 0.25;
// Try linear iteration if it has not been tried and
// the S sequence is converging.
//
if (stry || !spass) { break; }
for (i=0;i<n;i++)
{
k[i] = svk[i];
}
RealPolynomialIteration(&s,nz,&iflag);
if (*nz > 0) { return; }
// Linear iteration has failed. Flag that it has been
// tried and decrease the convergence criterion.
//
stry = 1;
betas *=0.25;
if (iflag == 0) { break; }
// If linear iteration signals an almost double real
// zero attempt quadratic iteration.
//
ui = -(s+s);
vi = s*s;
}
while (iflag != 0);
// Restore variables.
_restore_variables:
u = svu;
v = svv;
for (i=0;i<n;i++)
{
k[i] = svk[i];
}
// Try quadratic iteration if it has not been tried
// and the V sequence is converging.
//
if (vpass && !vtry) { goto _quadratic_iteration; }
// Recompute QP and scalar values to continue the
// second stage.
//
QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
ComputeScalarFactors(&type);
ovv = vv;
oss = ss;
otv = tv;
ots = ts;
}
}
void G4JTPolynomialSolver::
QuadraticPolynomialIteration(G4double *uu, G4double *vv, G4int *nz)
{
// Variable-shift k-polynomial iteration for a
// quadratic factor converges only if the zeros are
// equimodular or nearly so.
// uu, vv - coefficients of starting quadratic.
// nz - number of zeros found.
//
G4double ui=0.0, vi=0.0;
G4double omp=0.0;
G4double relstp=0.0;
G4double mp=0.0, ee=0.0, t=0.0, zm=0.0;
G4int type=0, i=1, j=0, tried=0;
*nz = 0;
tried = 0;
u = *uu;
v = *vv;
// Main loop.
while (1)
{
Quadratic(1.0,u,v,&szr,&szi,&lzr,&lzi);
// Return if roots of the quadratic are real and not
// close to multiple or nearly equal and of opposite
// sign.
//
if (std::fabs(std::fabs(szr)-std::fabs(lzr)) > 0.01 * std::fabs(lzr))
{ return; }
// Evaluate polynomial by quadratic synthetic division.
//
QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
mp = std::fabs(a-szr*b) + std::fabs(szi*b);
// Compute a rigorous bound on the rounding error in evaluating p.
//
zm = std::sqrt(std::fabs(v));
ee = 2.0*std::fabs(qp[0]);
t = -szr*b;
for (i=1;i<n;i++)
{
ee = ee*zm + std::fabs(qp[i]);
}
ee = ee*zm + std::fabs(a+t);
ee *= (5.0 *mre + 4.0*are);
ee = ee - (5.0*mre+2.0*are)*(std::fabs(a+t)+std::fabs(b)*zm)
+ 2.0*are*std::fabs(t);
// Iteration has converged sufficiently if the
// polynomial value is less than 20 times this bound.
//
if (mp <= 20.0*ee)
{
*nz = 2;
return;
}
j++;
// Stop iteration after 20 steps.
//
if (j > 20) { return; }
if (j >= 2)
{
if (!(relstp > 0.01 || mp < omp || tried))
{
// A cluster appears to be stalling the convergence.
// Five fixed shift steps are taken with a u,v close to the cluster.
//
if (relstp < eta) { relstp = eta; }
relstp = std::sqrt(relstp);
u = u - u*relstp;
v = v + v*relstp;
QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
for (i=0;i<5;i++)
{
ComputeScalarFactors(&type);
ComputeNextPolynomial(&type);
}
tried = 1;
j = 0;
}
}
omp = mp;
// Calculate next k polynomial and new u and v.
//
ComputeScalarFactors(&type);
ComputeNextPolynomial(&type);
ComputeScalarFactors(&type);
ComputeNewEstimate(type,&ui,&vi);
// If vi is zero the iteration is not converging.
//
if (!(vi != 0.0)) { return; }
relstp = std::fabs((vi-v)/vi);
u = ui;
v = vi;
}
}
void G4JTPolynomialSolver::
RealPolynomialIteration(G4double *sss, G4int *nz, G4int *iflag)
{
// Variable-shift H polynomial iteration for a real zero.
// sss - starting iterate
// nz - number of zeros found
// iflag - flag to indicate a pair of zeros near real axis.
G4double t=0.;
G4double omp=0.;
G4double pv=0.0, kv=0.0, s= *sss;
G4double ms=0.0, mp=0.0, ee=0.0;
G4int i=1, j=0;
*nz = 0;
*iflag = 0;
// Main loop
//
while (1)
{
pv = p[0];
// Evaluate p at s.
//
qp[0] = pv;
for (i=1;i<=n;i++)
{
pv = pv*s + p[i];
qp[i] = pv;
}
mp = std::fabs(pv);
// Compute a rigorous bound on the error in evaluating p.
//
ms = std::fabs(s);
ee = (mre/(are+mre))*std::fabs(qp[0]);
for (i=1;i<=n;i++)
{
ee = ee*ms + std::fabs(qp[i]);
}
// Iteration has converged sufficiently if the polynomial
// value is less than 20 times this bound.
//
if (mp <= 20.0*((are+mre)*ee-mre*mp))
{
*nz = 1;
szr = s;
szi = 0.0;
return;
}
j++;
// Stop iteration after 10 steps.
//
if (j > 10) { return; }
if (j >= 2)
{
if (!(std::fabs(t) > 0.001*std::fabs(s-t) || mp < omp))
{
// A cluster of zeros near the real axis has been encountered.
// Return with iflag set to initiate a quadratic iteration.
//
*iflag = 1;
*sss = s;
return;
} // Return if the polynomial value has increased significantly.
}
omp = mp;
// Compute t, the next polynomial, and the new iterate.
//
kv = k[0];
qk[0] = kv;
for (i=1;i<n;i++)
{
kv = kv*s + k[i];
qk[i] = kv;
}
if (std::fabs(kv) <= std::fabs(k[n-1])*10.0*eta) // Use unscaled form.
{
k[0] = 0.0;
for (i=1;i<n;i++)
{
k[i] = qk[i-1];
}
}
else // Use the scaled form of the recurrence if k at s is nonzero.
{
t = -pv/kv;
k[0] = qp[0];
for (i=1;i<n;i++)
{
k[i] = t*qk[i-1] + qp[i];
}
}
kv = k[0];
for (i=1;i<n;i++)
{
kv = kv*s + k[i];
}
t = 0.0;
if (std::fabs(kv) > std::fabs(k[n-1]*10.0*eta)) { t = -pv/kv; }
s += t;
}
}
void G4JTPolynomialSolver::ComputeScalarFactors(G4int *type)
{
// This function calculates scalar quantities used to
// compute the next k polynomial and new estimates of
// the quadratic coefficients.
// type - integer variable set here indicating how the
// calculations are normalized to avoid overflow.
// Synthetic division of k by the quadratic 1,u,v
//
QuadraticSyntheticDivision(n-1,&u,&v,k,qk,&c,&d);
if (std::fabs(c) <= std::fabs(k[n-1]*100.0*eta))
{
if (std::fabs(d) <= std::fabs(k[n-2]*100.0*eta))
{
*type = 3; // Type=3 indicates the quadratic is almost a factor of k.
return;
}
}
if (std::fabs(d) < std::fabs(c))
{
*type = 1; // Type=1 indicates that all formulas are divided by c.
e = a/c;
f = d/c;
g = u*e;
h = v*b;
a3 = a*e + (h/c+g)*b;
a1 = b - a*(d/c);
a7 = a + g*d + h*f;
return;
}
*type = 2; // Type=2 indicates that all formulas are divided by d.
e = a/d;
f = c/d;
g = u*b;
h = v*b;
a3 = (a+g)*e + h*(b/d);
a1 = b*f-a;
a7 = (f+u)*a + h;
}
void G4JTPolynomialSolver::ComputeNextPolynomial(G4int *type)
{
// Computes the next k polynomials using scalars
// computed in ComputeScalarFactors.
G4int i=2;
if (*type == 3) // Use unscaled form of the recurrence if type is 3.
{
k[0] = 0.0;
k[1] = 0.0;
for (i=2;i<n;i++)
{
k[i] = qk[i-2];
}
return;
}
G4double temp = a;
if (*type == 1) { temp = b; }
if (std::fabs(a1) <= std::fabs(temp)*eta*10.0)
{
// If a1 is nearly zero then use a special form of the recurrence.
//
k[0] = 0.0;
k[1] = -a7*qp[0];
for(i=2;i<n;i++)
{
k[i] = a3*qk[i-2] - a7*qp[i-1];
}
return;
}
// Use scaled form of the recurrence.
//
a7 /= a1;
a3 /= a1;
k[0] = qp[0];
k[1] = qp[1] - a7*qp[0];
for (i=2;i<n;i++)
{
k[i] = a3*qk[i-2] - a7*qp[i-1] + qp[i];
}
}
void G4JTPolynomialSolver::
ComputeNewEstimate(G4int type, G4double *uu, G4double *vv)
{
// Compute new estimates of the quadratic coefficients
// using the scalars computed in calcsc.
G4double a4=0.0, a5=0.0, b1=0.0, b2=0.0,
c1=0.0, c2=0.0, c3=0.0, c4=0.0, temp=0.0;
// Use formulas appropriate to setting of type.
//
if (type == 3) // If type=3 the quadratic is zeroed.
{
*uu = 0.0;
*vv = 0.0;
return;
}
if (type == 2)
{
a4 = (a+g)*f + h;
a5 = (f+u)*c + v*d;
}
else
{
a4 = a + u*b +h*f;
a5 = c + (u+v*f)*d;
}
// Evaluate new quadratic coefficients.
//
b1 = -k[n-1]/p[n];
b2 = -(k[n-2]+b1*p[n-1])/p[n];
c1 = v*b2*a1;
c2 = b1*a7;
c3 = b1*b1*a3;
c4 = c1 - c2 - c3;
temp = a5 + b1*a4 - c4;
if (!(temp != 0.0))
{
*uu = 0.0;
*vv = 0.0;
return;
}
*uu = u - (u*(c3+c2)+v*(b1*a1+b2*a7))/temp;
*vv = v*(1.0+c4/temp);
return;
}
void G4JTPolynomialSolver::
QuadraticSyntheticDivision(G4int nn, G4double *uu, G4double *vv,
std::vector<G4double> &pp, std::vector<G4double> &qq,
G4double *aa, G4double *bb)
{
// Divides pp by the quadratic 1,uu,vv placing the quotient
// in qq and the remainder in aa,bb.
G4double cc=0.0;
*bb = pp[0];
qq[0] = *bb;
*aa = pp[1] - (*bb)*(*uu);
qq[1] = *aa;
for (G4int i=2;i<=nn;i++)
{
cc = pp[i] - (*aa)*(*uu) - (*bb)*(*vv);
qq[i] = cc;
*bb = *aa;
*aa = cc;
}
}
void G4JTPolynomialSolver::Quadratic(G4double aa,G4double b1,
G4double cc,G4double *ssr,G4double *ssi,
G4double *lr,G4double *li)
{
// Calculate the zeros of the quadratic aa*z^2 + b1*z + cc.
// The quadratic formula, modified to avoid overflow, is used
// to find the larger zero if the zeros are real and both
// are complex. The smaller real zero is found directly from
// the product of the zeros c/a.
G4double bb=0.0, dd=0.0, ee=0.0;
if (!(aa != 0.0)) // less than two roots
{
if (b1 != 0.0)
{ *ssr = -cc/b1; }
else
{ *ssr = 0.0; }
*lr = 0.0;
*ssi = 0.0;
*li = 0.0;
return;
}
if (!(cc != 0.0)) // one real root, one zero root
{
*ssr = 0.0;
*lr = -b1/aa;
*ssi = 0.0;
*li = 0.0;
return;
}
// Compute discriminant avoiding overflow.
//
bb = b1/2.0;
if (std::fabs(bb) < std::fabs(cc))
{
if (cc < 0.0)
{ ee = -aa; }
else
{ ee = aa; }
ee = bb*(bb/std::fabs(cc)) - ee;
dd = std::sqrt(std::fabs(ee))*std::sqrt(std::fabs(cc));
}
else
{
ee = 1.0 - (aa/bb)*(cc/bb);
dd = std::sqrt(std::fabs(ee))*std::fabs(bb);
}
if (ee < 0.0) // complex conjugate zeros
{
*ssr = -bb/aa;
*lr = *ssr;
*ssi = std::fabs(dd/aa);
*li = -(*ssi);
}
else
{
if (bb >= 0.0) // real zeros.
{ dd = -dd; }
*lr = (-bb+dd)/aa;
*ssr = 0.0;
if (*lr != 0.0)
{ *ssr = (cc/ *lr)/aa; }
*ssi = 0.0;
*li = 0.0;
}
}
source/global/HEPNumerics/src/G4SimpleIntegration.cc
+21 -24
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4SimpleIntegration.cc,v 1.5 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4SimpleIntegration.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
// Implementation file for simple integration methods
//
@@ -35,21 +35,21 @@ G4int G4SimpleIntegration::fMaxDepth = 100 ;
G4SimpleIntegration::G4SimpleIntegration( function pFunction )
: fFunction(pFunction),
fTolerance(.0001)
{
fFunction = pFunction ;
}
G4SimpleIntegration::G4SimpleIntegration( function pFunction,
G4double pTolerance)
G4double pTolerance)
: fFunction(pFunction),
fTolerance(pTolerance)
{
fFunction = pFunction ;
fTolerance = pTolerance ;
}
G4SimpleIntegration::~G4SimpleIntegration()
{
;
}
// Simple integration methods
@@ -57,13 +57,12 @@ G4SimpleIntegration::~G4SimpleIntegration()
G4double
G4SimpleIntegration::Trapezoidal(G4double xInitial,
G4double xFinal,
G4int iterationNumber )
G4int iterationNumber )
{
G4int i ;
G4double Step = (xFinal - xInitial)/iterationNumber ;
G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
G4double x = xInitial ;
for(i=1;i<iterationNumber;i++)
for(G4int i=1;i<iterationNumber;i++)
{
x += Step ;
mean += fFunction(x) ;
@@ -74,13 +73,12 @@ G4SimpleIntegration::Trapezoidal(G4double xInitial,
G4double
G4SimpleIntegration::MidPoint(G4double xInitial,
G4double xFinal,
G4int iterationNumber )
G4int iterationNumber )
{
G4int i ;
G4double Step = (xFinal - xInitial)/iterationNumber ;
G4double x = xInitial + 0.5*Step;
G4double mean = fFunction(x) ;
for(i=1;i<iterationNumber;i++)
for(G4int i=1;i<iterationNumber;i++)
{
x += Step ;
mean += fFunction(x) ;
@@ -91,15 +89,14 @@ G4SimpleIntegration::MidPoint(G4double xInitial,
G4double
G4SimpleIntegration::Gauss(G4double xInitial,
G4double xFinal,
G4int iterationNumber )
G4int iterationNumber )
{
G4int i ;
G4double x ;
G4double x=0.;
static G4double root = 1.0/std::sqrt(3.0) ;
G4double Step = (xFinal - xInitial)/(2.0*iterationNumber) ;
G4double delta = Step*root ;
G4double mean = 0.0 ;
for(i=0;i<iterationNumber;i++)
for(G4int i=0;i<iterationNumber;i++)
{
x = (2*i + 1)*Step ;
mean += (fFunction(x+delta) + fFunction(x-delta)) ;
@@ -110,15 +107,14 @@ G4SimpleIntegration::Gauss(G4double xInitial,
G4double
G4SimpleIntegration::Simpson(G4double xInitial,
G4double xFinal,
G4int iterationNumber )
G4int iterationNumber )
{
G4int i ;
G4double Step = (xFinal - xInitial)/iterationNumber ;
G4double x = xInitial ;
G4double xPlus = xInitial + 0.5*Step ;
G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
G4double sum = fFunction(xPlus) ;
for(i=1;i<iterationNumber;i++)
for(G4int i=1;i<iterationNumber;i++)
{
x += Step ;
xPlus += Step ;
@@ -142,7 +138,7 @@ G4SimpleIntegration::AdaptGaussIntegration( G4double xInitial,
AdaptGauss(xInitial,xFinal,sum,depth) ;
return sum ;
}
G4double
G4SimpleIntegration::Gauss( G4double xInitial,
@@ -162,12 +158,13 @@ G4SimpleIntegration::Gauss( G4double xInitial,
void
G4SimpleIntegration::AdaptGauss( G4double xInitial,
G4double xFinal,
G4double& sum,
G4int& depth )
G4double& sum,
G4int& depth )
{
if(depth >fMaxDepth)
{
G4Exception("Function varies too rapidly in G4SimpleIntegration::AdaptGauss") ;
G4Exception("G4SimpleIntegration::AdaptGauss()", "Error",
FatalException, "Function varies too rapidly !") ;
}
G4double xMean = (xInitial + xFinal)/2.0 ;
G4double leftHalf = Gauss(xInitial,xMean) ;
source/global/HEPNumerics/src/G4VGaussianQuadrature.cc
+8 -14
View File
@@ -21,8 +21,8 @@
// ********************************************************************
//
//
// $Id: G4VGaussianQuadrature.cc,v 1.5 2004/11/12 17:38:33 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
// $Id: G4VGaussianQuadrature.cc,v 1.6 2005/03/15 19:11:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-01 $
//
// Implementation file for G4VGaussianQuadrature virtual base class
//
@@ -31,14 +31,9 @@
#include "globals.hh"
#include "G4VGaussianQuadrature.hh"
G4VGaussianQuadrature::G4VGaussianQuadrature( function pFunction )
: fFunction(pFunction), fAbscissa(0), fWeight(0), fNumber(0)
{
fFunction = pFunction ;
fAbscissa = 0;
fWeight = 0;
}
// -------------------------------------------------------------------
@@ -54,7 +49,6 @@ G4VGaussianQuadrature::~G4VGaussianQuadrature()
// -------------------------- Access functions ----------------------------------
G4double
G4VGaussianQuadrature::GetAbscissa(G4int index) const
{
@@ -67,6 +61,10 @@ G4VGaussianQuadrature::GetWeight(G4int index) const
return fWeight[index] ;
}
G4int G4VGaussianQuadrature::GetNumber() const
{
return fNumber ;
}
// ----------------------------------------------------------------------------
//
@@ -84,19 +82,15 @@ G4VGaussianQuadrature::GammaLogarithm(G4double xx)
static G4double cof[6] = { 76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2, -0.5395239384953e-5 } ;
register G4int j;
G4double x = xx - 1.0;
G4double tmp = x + 5.5;
tmp -= (x + 0.5) * std::log(tmp);
G4double ser = 1.000000000190015;
for ( j = 0; j <= 5; j++ )
for ( size_t j = 0; j <= 5; j++ )
{
x += 1.0;
ser += cof[j]/x;
}
return -tmp + std::log(2.5066282746310005*ser);
}