Gabriele Cosmo
916 lines
22 KiB
C++
916 lines
22 KiB
C++
//
|
|
// ********************************************************************
|
|
// * License and Disclaimer *
|
|
// * *
|
|
// * The Geant4 software is copyright of the Copyright Holders of *
|
|
// * the Geant4 Collaboration. It is provided under the terms and *
|
|
// * conditions of the Geant4 Software License, included in the file *
|
|
// * LICENSE and available at http://cern.ch/geant4/license . These *
|
|
// * include a list of copyright holders. *
|
|
// * *
|
|
// * 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. Please see the license in the file LICENSE and URL above *
|
|
// * for the full disclaimer and the limitation of liability. *
|
|
// * *
|
|
// * This code implementation is the result of the scientific and *
|
|
// * technical work of the GEANT4 collaboration. *
|
|
// * By using, copying, modifying or distributing the software (or *
|
|
// * any work based on the software) you agree to acknowledge its *
|
|
// * use in resulting scientific publications, and indicate your *
|
|
// * acceptance of all terms of the Geant4 Software license. *
|
|
// ********************************************************************
|
|
//
|
|
//
|
|
// --------------------------------------------------------------------
|
|
// GEANT 4 class source file
|
|
//
|
|
// G4JTPolynomialSolver
|
|
//
|
|
// Implementation based on Jenkins-Traub algorithm.
|
|
// --------------------------------------------------------------------
|
|
|
|
#include "G4JTPolynomialSolver.hh"
|
|
#include "G4SystemOfUnits.hh"
|
|
#include "G4Pow.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()
|
|
: sr(0.), si(0.), u(0.),v(0.),
|
|
a(0.), b(0.), c(0.), d(0.),
|
|
a1(0.), a3(0.), a7(0.),
|
|
e(0.), f(0.), g(0.), h(0.),
|
|
szr(0.), szi(0.), lzr(0.), lzi(0.),
|
|
n(0)
|
|
{
|
|
}
|
|
|
|
G4JTPolynomialSolver::~G4JTPolynomialSolver()
|
|
{
|
|
}
|
|
|
|
G4int G4JTPolynomialSolver::FindRoots(G4double *op, G4int degr,
|
|
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;
|
|
G4Pow* power = G4Pow::GetInstance();
|
|
|
|
// Initialization of constants for shift rotation.
|
|
//
|
|
static const G4double xx = std::sqrt(0.5);
|
|
static const G4double rot = 94.0*deg;
|
|
static const G4double cosr = std::cos(rot),
|
|
sinr = std::sin(rot);
|
|
G4double xo = xx, yo = -xx;
|
|
n = degr;
|
|
|
|
// 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 = degr - n;
|
|
zeror[j] = 0.0;
|
|
zeroi[j] = 0.0;
|
|
n--;
|
|
}
|
|
if (n < 1) { return -1; }
|
|
|
|
// Allocate buffers here
|
|
//
|
|
std::vector<G4double> temp(degr+1) ;
|
|
std::vector<G4double> pt(degr+1) ;
|
|
|
|
p.assign(degr+1,0) ;
|
|
qp.assign(degr+1,0) ;
|
|
k.assign(degr+1,0) ;
|
|
qk.assign(degr+1,0) ;
|
|
svk.assign(degr+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[degr-1] = -p[1]/p[0];
|
|
zeroi[degr-1] = 0.0;
|
|
n -= 1;
|
|
return degr - n ;
|
|
}
|
|
if (n == 2) // Calculate the final zero or pair of zeros.
|
|
{
|
|
Quadratic(p[0],p[1],p[2],&zeror[degr-2],&zeroi[degr-2],
|
|
&zeror[degr-1],&zeroi[degr-1]);
|
|
n -= 2;
|
|
return degr - 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)(G4Log(sc)/G4Log(base) + 0.5);
|
|
factor = power->powN(base,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 = G4Exp((G4Log(-pt[n])-G4Log(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*xo - sinr*yo;
|
|
yo = sinr*xo + cosr*yo;
|
|
xo = xxx;
|
|
sr = bnd*xo;
|
|
si = bnd*yo;
|
|
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 = degr - 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 degr - 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, xs=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 xs.
|
|
//
|
|
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 xs 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];
|
|
}
|
|
xs = ss;
|
|
|
|
// Choose iteration according to the fastest converging sequence.
|
|
//
|
|
vtry = 0;
|
|
stry = 0;
|
|
if ((spass && (!vpass)) || (tss < tvv))
|
|
{
|
|
RealPolynomialIteration(&xs,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 = -(xs+xs);
|
|
vi = xs*xs;
|
|
}
|
|
|
|
_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(&xs,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 = -(xs+xs);
|
|
vi = xs*xs;
|
|
}
|
|
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, xs= *sss;
|
|
G4double mx=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 xs.
|
|
//
|
|
qp[0] = pv;
|
|
for (i=1;i<=n;i++)
|
|
{
|
|
pv = pv*xs + p[i];
|
|
qp[i] = pv;
|
|
}
|
|
mp = std::fabs(pv);
|
|
|
|
// Compute a rigorous bound on the error in evaluating p.
|
|
//
|
|
mx = std::fabs(xs);
|
|
ee = (mre/(are+mre))*std::fabs(qp[0]);
|
|
for (i=1;i<=n;i++)
|
|
{
|
|
ee = ee*mx + 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 = xs;
|
|
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(xs-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 = xs;
|
|
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*xs + 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 xs 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*xs + k[i];
|
|
}
|
|
t = 0.0;
|
|
if (std::fabs(kv) > std::fabs(k[n-1]*10.0*eta)) { t = -pv/kv; }
|
|
xs += 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;
|
|
}
|
|
}
|