#include <G4JTPolynomialSolver.hh>
Public Member Functions | |
G4JTPolynomialSolver () | |
~G4JTPolynomialSolver () | |
G4int | FindRoots (G4double *op, G4int degree, G4double *zeror, G4double *zeroi) |
Definition at line 70 of file G4JTPolynomialSolver.hh.
G4JTPolynomialSolver::G4JTPolynomialSolver | ( | ) |
Definition at line 47 of file G4JTPolynomialSolver.cc.
00048 : sr(0.), si(0.), u(0.),v(0.), 00049 a(0.), b(0.), c(0.), d(0.), 00050 a1(0.), a3(0.), a7(0.), 00051 e(0.), f(0.), g(0.), h(0.), 00052 szr(0.), szi(0.), lzr(0.), lzi(0.), 00053 n(0) 00054 { 00055 }
G4JTPolynomialSolver::~G4JTPolynomialSolver | ( | ) |
G4int G4JTPolynomialSolver::FindRoots | ( | G4double * | op, | |
G4int | degree, | |||
G4double * | zeror, | |||
G4double * | zeroi | |||
) |
Definition at line 61 of file G4JTPolynomialSolver.cc.
Referenced by G4TwistTrapParallelSide::DistanceToSurface(), G4TwistTrapAlphaSide::DistanceToSurface(), and G4TwistBoxSide::DistanceToSurface().
00063 { 00064 G4double t=0.0, aa=0.0, bb=0.0, cc=0.0, factor=1.0; 00065 G4double max=0.0, min=infin, xxx=0.0, x=0.0, sc=0.0, bnd=0.0; 00066 G4double xm=0.0, ff=0.0, df=0.0, dx=0.0; 00067 G4int cnt=0, nz=0, i=0, j=0, jj=0, l=0, nm1=0, zerok=0; 00068 00069 // Initialization of constants for shift rotation. 00070 // 00071 G4double xx = std::sqrt(0.5); 00072 G4double yy = -xx, 00073 rot = 94.0*deg; 00074 G4double cosr = std::cos(rot), 00075 sinr = std::sin(rot); 00076 n = degr; 00077 00078 // Algorithm fails if the leading coefficient is zero. 00079 // 00080 if (!(op[0] != 0.0)) { return -1; } 00081 00082 // Remove the zeros at the origin, if any. 00083 // 00084 while (!(op[n] != 0.0)) 00085 { 00086 j = degr - n; 00087 zeror[j] = 0.0; 00088 zeroi[j] = 0.0; 00089 n--; 00090 } 00091 if (n < 1) { return -1; } 00092 00093 // Allocate buffers here 00094 // 00095 std::vector<G4double> temp(degr+1) ; 00096 std::vector<G4double> pt(degr+1) ; 00097 00098 p.assign(degr+1,0) ; 00099 qp.assign(degr+1,0) ; 00100 k.assign(degr+1,0) ; 00101 qk.assign(degr+1,0) ; 00102 svk.assign(degr+1,0) ; 00103 00104 // Make a copy of the coefficients. 00105 // 00106 for (i=0;i<=n;i++) 00107 { p[i] = op[i]; } 00108 00109 do 00110 { 00111 if (n == 1) // Start the algorithm for one zero. 00112 { 00113 zeror[degr-1] = -p[1]/p[0]; 00114 zeroi[degr-1] = 0.0; 00115 n -= 1; 00116 return degr - n ; 00117 } 00118 if (n == 2) // Calculate the final zero or pair of zeros. 00119 { 00120 Quadratic(p[0],p[1],p[2],&zeror[degr-2],&zeroi[degr-2], 00121 &zeror[degr-1],&zeroi[degr-1]); 00122 n -= 2; 00123 return degr - n ; 00124 } 00125 00126 // Find largest and smallest moduli of coefficients. 00127 // 00128 max = 0.0; 00129 min = infin; 00130 for (i=0;i<=n;i++) 00131 { 00132 x = std::fabs(p[i]); 00133 if (x > max) { max = x; } 00134 if (x != 0.0 && x < min) { min = x; } 00135 } 00136 00137 // Scale if there are large or very small coefficients. 00138 // Computes a scale factor to multiply the coefficients of the 00139 // polynomial. The scaling is done to avoid overflow and to 00140 // avoid undetected underflow interfering with the convergence 00141 // criterion. The factor is a power of the base. 00142 // 00143 sc = lo/min; 00144 00145 if ( ((sc <= 1.0) && (max >= 10.0)) 00146 || ((sc > 1.0) && (infin/sc >= max)) 00147 || ((infin/sc >= max) && (max >= 10)) ) 00148 { 00149 if (!( sc != 0.0 )) 00150 { sc = smalno ; } 00151 l = (G4int)(std::log(sc)/std::log(base) + 0.5); 00152 factor = std::pow(base*1.0,l); 00153 if (factor != 1.0) 00154 { 00155 for (i=0;i<=n;i++) 00156 { p[i] = factor*p[i]; } // Scale polynomial. 00157 } 00158 } 00159 00160 // Compute lower bound on moduli of roots. 00161 // 00162 for (i=0;i<=n;i++) 00163 { 00164 pt[i] = (std::fabs(p[i])); 00165 } 00166 pt[n] = - pt[n]; 00167 00168 // Compute upper estimate of bound. 00169 // 00170 x = std::exp((std::log(-pt[n])-std::log(pt[0])) / (G4double)n); 00171 00172 // If Newton step at the origin is better, use it. 00173 // 00174 if (pt[n-1] != 0.0) 00175 { 00176 xm = -pt[n]/pt[n-1]; 00177 if (xm < x) { x = xm; } 00178 } 00179 00180 // Chop the interval (0,x) until ff <= 0 00181 // 00182 while (1) 00183 { 00184 xm = x*0.1; 00185 ff = pt[0]; 00186 for (i=1;i<=n;i++) 00187 { ff = ff*xm + pt[i]; } 00188 if (ff <= 0.0) { break; } 00189 x = xm; 00190 } 00191 dx = x; 00192 00193 // Do Newton interation until x converges to two decimal places. 00194 // 00195 while (std::fabs(dx/x) > 0.005) 00196 { 00197 ff = pt[0]; 00198 df = ff; 00199 for (i=1;i<n;i++) 00200 { 00201 ff = ff*x + pt[i]; 00202 df = df*x + ff; 00203 } 00204 ff = ff*x + pt[n]; 00205 dx = ff/df; 00206 x -= dx; 00207 } 00208 bnd = x; 00209 00210 // Compute the derivative as the initial k polynomial 00211 // and do 5 steps with no shift. 00212 // 00213 nm1 = n - 1; 00214 for (i=1;i<n;i++) 00215 { k[i] = (G4double)(n-i)*p[i]/(G4double)n; } 00216 k[0] = p[0]; 00217 aa = p[n]; 00218 bb = p[n-1]; 00219 zerok = (k[n-1] == 0); 00220 for(jj=0;jj<5;jj++) 00221 { 00222 cc = k[n-1]; 00223 if (!zerok) // Use a scaled form of recurrence if k at 0 is nonzero. 00224 { 00225 // Use a scaled form of recurrence if value of k at 0 is nonzero. 00226 // 00227 t = -aa/cc; 00228 for (i=0;i<nm1;i++) 00229 { 00230 j = n-i-1; 00231 k[j] = t*k[j-1]+p[j]; 00232 } 00233 k[0] = p[0]; 00234 zerok = (std::fabs(k[n-1]) <= std::fabs(bb)*eta*10.0); 00235 } 00236 else // Use unscaled form of recurrence. 00237 { 00238 for (i=0;i<nm1;i++) 00239 { 00240 j = n-i-1; 00241 k[j] = k[j-1]; 00242 } 00243 k[0] = 0.0; 00244 zerok = (!(k[n-1] != 0.0)); 00245 } 00246 } 00247 00248 // Save k for restarts with new shifts. 00249 // 00250 for (i=0;i<n;i++) 00251 { temp[i] = k[i]; } 00252 00253 // Loop to select the quadratic corresponding to each new shift. 00254 // 00255 for (cnt = 0;cnt < 20;cnt++) 00256 { 00257 // Quadratic corresponds to a double shift to a 00258 // non-real point and its complex conjugate. The point 00259 // has modulus bnd and amplitude rotated by 94 degrees 00260 // from the previous shift. 00261 // 00262 xxx = cosr*xx - sinr*yy; 00263 yy = sinr*xx + cosr*yy; 00264 xx = xxx; 00265 sr = bnd*xx; 00266 si = bnd*yy; 00267 u = -2.0 * sr; 00268 v = bnd; 00269 ComputeFixedShiftPolynomial(20*(cnt+1),&nz); 00270 if (nz != 0) 00271 { 00272 // The second stage jumps directly to one of the third 00273 // stage iterations and returns here if successful. 00274 // Deflate the polynomial, store the zero or zeros and 00275 // return to the main algorithm. 00276 // 00277 j = degr - n; 00278 zeror[j] = szr; 00279 zeroi[j] = szi; 00280 n -= nz; 00281 for (i=0;i<=n;i++) 00282 { p[i] = qp[i]; } 00283 if (nz != 1) 00284 { 00285 zeror[j+1] = lzr; 00286 zeroi[j+1] = lzi; 00287 } 00288 break; 00289 } 00290 else 00291 { 00292 // If the iteration is unsuccessful another quadratic 00293 // is chosen after restoring k. 00294 // 00295 for (i=0;i<n;i++) 00296 { 00297 k[i] = temp[i]; 00298 } 00299 } 00300 } 00301 } 00302 while (nz != 0); // End of initial DO loop 00303 00304 // Return with failure if no convergence with 20 shifts. 00305 // 00306 return degr - n; 00307 }