#include <G4JTPolynomialSolver.hh>

Public Member Functions
	G4JTPolynomialSolver ()

	~G4JTPolynomialSolver ()

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

Detailed Description

Definition at line 70 of file G4JTPolynomialSolver.hh.

Constructor & Destructor Documentation

◆ G4JTPolynomialSolver()

G4JTPolynomialSolver::G4JTPolynomialSolver ( )

Definition at line 47 of file G4JTPolynomialSolver.cc.

  : sr(0.), si(0.), u(0.),v(0.),
    a(0.), b(0.), c(0.), d(0.),
    a1(0.), a2(0.), a3(0.), a6(0.), a7(0.),
    e(0.), f(0.), g(0.), h(0.),
    szr(0.), szi(0.), lzr(0.), lzi(0.),
    n(0), nmi(0)
{
}

◆ ~G4JTPolynomialSolver()

G4JTPolynomialSolver::~G4JTPolynomialSolver ( )

Definition at line 57 of file G4JTPolynomialSolver.cc.

58{

59}

Member Function Documentation

◆ FindRoots()

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

Definition at line 61 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;
 
  // 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 = 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)(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 = 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;
}

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