G4JTPolynomialSolver Class Reference

#include <G4JTPolynomialSolver.hh>

Public Member Functions
	G4JTPolynomialSolver ()=default
	~G4JTPolynomialSolver ()=default
G4int	FindRoots (G4double *op, G4int degree, G4double *zeror, G4double *zeroi)

Detailed Description

Definition at line 66 of file G4JTPolynomialSolver.hh.

Constructor & Destructor Documentation

G4JTPolynomialSolver()

G4JTPolynomialSolver::G4JTPolynomialSolver ( )

default

~G4JTPolynomialSolver()

G4JTPolynomialSolver::~G4JTPolynomialSolver ( )

default

Member Function Documentation

FindRoots()

G4int G4JTPolynomialSolver::FindRoots	(	G4double *	op,
		G4int	degree,
		G4double *	zeror,
		G4double *	zeroi )

Definition at line 44 of file G4JTPolynomialSolver.cc.

{
  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(true)
    {
      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 = static_cast<G4int>(k[n - 1] == 0);
    for(jj = 0; jj < 5; ++jj)
    {
      cc = k[n - 1];
      if(zerok == 0)  // 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 =
          static_cast<G4int>(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 = static_cast<G4int>(!(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;
      }
      
      // 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;
}

Referenced by G4TwistBoxSide::DistanceToSurface(), G4TwistTrapAlphaSide::DistanceToSurface(), and G4TwistTrapParallelSide::DistanceToSurface().

---

The documentation for this class was generated from the following files:

• G4JTPolynomialSolver.hh
• G4JTPolynomialSolver.cc