| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4JTPolynomialSolver class implementation. 26 // G4JTPolynomialSolver class implementation.
27 // 27 //
28 // Author: Oliver Link, 15.02.2005 28 // Author: Oliver Link, 15.02.2005
29 // Translated to C++ and adapted to us 29 // Translated to C++ and adapted to use STL vectors.
30 // ------------------------------------------- 30 // --------------------------------------------------------------------
31 31
32 #include "G4JTPolynomialSolver.hh" 32 #include "G4JTPolynomialSolver.hh"
33 #include "G4Pow.hh" 33 #include "G4Pow.hh"
34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
35 35
36 const G4double G4JTPolynomialSolver::base = 36 const G4double G4JTPolynomialSolver::base = 2;
37 const G4double G4JTPolynomialSolver::eta = 37 const G4double G4JTPolynomialSolver::eta = DBL_EPSILON;
38 const G4double G4JTPolynomialSolver::infin = 38 const G4double G4JTPolynomialSolver::infin = DBL_MAX;
39 const G4double G4JTPolynomialSolver::smalno = 39 const G4double G4JTPolynomialSolver::smalno = DBL_MIN;
40 const G4double G4JTPolynomialSolver::are = 40 const G4double G4JTPolynomialSolver::are = DBL_EPSILON;
41 const G4double G4JTPolynomialSolver::mre = 41 const G4double G4JTPolynomialSolver::mre = DBL_EPSILON;
42 const G4double G4JTPolynomialSolver::lo = 42 const G4double G4JTPolynomialSolver::lo = DBL_MIN / DBL_EPSILON;
43 43
>> 44 G4JTPolynomialSolver::G4JTPolynomialSolver() {}
>> 45
>> 46 G4JTPolynomialSolver::~G4JTPolynomialSolver() {}
>> 47
44 G4int G4JTPolynomialSolver::FindRoots(G4double 48 G4int G4JTPolynomialSolver::FindRoots(G4double* op, G4int degr, G4double* zeror,
45 G4double 49 G4double* zeroi)
46 { 50 {
47 G4double t = 0.0, aa = 0.0, bb = 0.0, cc = 0 51 G4double t = 0.0, aa = 0.0, bb = 0.0, cc = 0.0, factor = 1.0;
48 G4double max = 0.0, min = infin, xxx = 0.0, 52 G4double max = 0.0, min = infin, xxx = 0.0, x = 0.0, sc = 0.0, bnd = 0.0;
49 G4double xm = 0.0, ff = 0.0, df = 0.0, dx = 53 G4double xm = 0.0, ff = 0.0, df = 0.0, dx = 0.0;
50 G4int cnt = 0, nz = 0, i = 0, j = 0, jj = 0, 54 G4int cnt = 0, nz = 0, i = 0, j = 0, jj = 0, l = 0, nm1 = 0, zerok = 0;
51 G4Pow* power = G4Pow::GetInstance(); 55 G4Pow* power = G4Pow::GetInstance();
52 56
53 // Initialization of constants for shift rot 57 // Initialization of constants for shift rotation.
54 // 58 //
55 static const G4double xx = std::sqrt(0.5); 59 static const G4double xx = std::sqrt(0.5);
56 static const G4double rot = 94.0 * deg; 60 static const G4double rot = 94.0 * deg;
57 static const G4double cosr = std::cos(rot), 61 static const G4double cosr = std::cos(rot), sinr = std::sin(rot);
58 G4double xo = xx, yo = -xx; 62 G4double xo = xx, yo = -xx;
59 n = degr; 63 n = degr;
60 64
61 // Algorithm fails if the leading coefficie 65 // Algorithm fails if the leading coefficient is zero.
62 // 66 //
63 if(!(op[0] != 0.0)) 67 if(!(op[0] != 0.0))
64 { 68 {
65 return -1; 69 return -1;
66 } 70 }
67 71
68 // Remove the zeros at the origin, if any. 72 // Remove the zeros at the origin, if any.
69 // 73 //
70 while(!(op[n] != 0.0)) 74 while(!(op[n] != 0.0))
71 { 75 {
72 j = degr - n; 76 j = degr - n;
73 zeror[j] = 0.0; 77 zeror[j] = 0.0;
74 zeroi[j] = 0.0; 78 zeroi[j] = 0.0;
75 n--; 79 n--;
76 } 80 }
77 if(n < 1) 81 if(n < 1)
78 { 82 {
79 return -1; 83 return -1;
80 } 84 }
81 85
82 // Allocate buffers here 86 // Allocate buffers here
83 // 87 //
84 std::vector<G4double> temp(degr + 1); 88 std::vector<G4double> temp(degr + 1);
85 std::vector<G4double> pt(degr + 1); 89 std::vector<G4double> pt(degr + 1);
86 90
87 p.assign(degr + 1, 0); 91 p.assign(degr + 1, 0);
88 qp.assign(degr + 1, 0); 92 qp.assign(degr + 1, 0);
89 k.assign(degr + 1, 0); 93 k.assign(degr + 1, 0);
90 qk.assign(degr + 1, 0); 94 qk.assign(degr + 1, 0);
91 svk.assign(degr + 1, 0); 95 svk.assign(degr + 1, 0);
92 96
93 // Make a copy of the coefficients. 97 // Make a copy of the coefficients.
94 // 98 //
95 for(i = 0; i <= n; ++i) 99 for(i = 0; i <= n; ++i)
96 { 100 {
97 p[i] = op[i]; 101 p[i] = op[i];
98 } 102 }
99 103
100 do 104 do
101 { 105 {
102 if(n == 1) // Start the algorithm for one 106 if(n == 1) // Start the algorithm for one zero.
103 { 107 {
104 zeror[degr - 1] = -p[1] / p[0]; 108 zeror[degr - 1] = -p[1] / p[0];
105 zeroi[degr - 1] = 0.0; 109 zeroi[degr - 1] = 0.0;
106 n -= 1; 110 n -= 1;
107 return degr - n; 111 return degr - n;
108 } 112 }
109 if(n == 2) // Calculate the final zero or 113 if(n == 2) // Calculate the final zero or pair of zeros.
110 { 114 {
111 Quadratic(p[0], p[1], p[2], &zeror[degr 115 Quadratic(p[0], p[1], p[2], &zeror[degr - 2], &zeroi[degr - 2],
112 &zeror[degr - 1], &zeroi[degr 116 &zeror[degr - 1], &zeroi[degr - 1]);
113 n -= 2; 117 n -= 2;
114 return degr - n; 118 return degr - n;
115 } 119 }
116 120
117 // Find largest and smallest moduli of coe 121 // Find largest and smallest moduli of coefficients.
118 // 122 //
119 max = 0.0; 123 max = 0.0;
120 min = infin; 124 min = infin;
121 for(i = 0; i <= n; ++i) 125 for(i = 0; i <= n; ++i)
122 { 126 {
123 x = std::fabs(p[i]); 127 x = std::fabs(p[i]);
124 if(x > max) 128 if(x > max)
125 { 129 {
126 max = x; 130 max = x;
127 } 131 }
128 if(x != 0.0 && x < min) 132 if(x != 0.0 && x < min)
129 { 133 {
130 min = x; 134 min = x;
131 } 135 }
132 } 136 }
133 137
134 // Scale if there are large or very small 138 // Scale if there are large or very small coefficients.
135 // Computes a scale factor to multiply the 139 // Computes a scale factor to multiply the coefficients of the
136 // polynomial. The scaling is done to avoi 140 // polynomial. The scaling is done to avoid overflow and to
137 // avoid undetected underflow interfering 141 // avoid undetected underflow interfering with the convergence
138 // criterion. The factor is a power of the 142 // criterion. The factor is a power of the base.
139 // 143 //
140 sc = lo / min; 144 sc = lo / min;
141 145
142 if(((sc <= 1.0) && (max >= 10.0)) || ((sc 146 if(((sc <= 1.0) && (max >= 10.0)) || ((sc > 1.0) && (infin / sc >= max)) ||
143 ((infin / sc >= max) && (max >= 10))) 147 ((infin / sc >= max) && (max >= 10)))
144 { 148 {
145 if(!(sc != 0.0)) 149 if(!(sc != 0.0))
146 { 150 {
147 sc = smalno; 151 sc = smalno;
148 } 152 }
149 l = (G4int)(G4Log(sc) / G4Log(base) 153 l = (G4int)(G4Log(sc) / G4Log(base) + 0.5);
150 factor = power->powN(base, l); 154 factor = power->powN(base, l);
151 if(factor != 1.0) 155 if(factor != 1.0)
152 { 156 {
153 for(i = 0; i <= n; ++i) 157 for(i = 0; i <= n; ++i)
154 { 158 {
155 p[i] = factor * p[i]; 159 p[i] = factor * p[i];
156 } // Scale polynomial. 160 } // Scale polynomial.
157 } 161 }
158 } 162 }
159 163
160 // Compute lower bound on moduli of roots. 164 // Compute lower bound on moduli of roots.
161 // 165 //
162 for(i = 0; i <= n; ++i) 166 for(i = 0; i <= n; ++i)
163 { 167 {
164 pt[i] = (std::fabs(p[i])); 168 pt[i] = (std::fabs(p[i]));
165 } 169 }
166 pt[n] = -pt[n]; 170 pt[n] = -pt[n];
167 171
168 // Compute upper estimate of bound. 172 // Compute upper estimate of bound.
169 // 173 //
170 x = G4Exp((G4Log(-pt[n]) - G4Log(pt[0])) / 174 x = G4Exp((G4Log(-pt[n]) - G4Log(pt[0])) / (G4double) n);
171 175
172 // If Newton step at the origin is better, 176 // If Newton step at the origin is better, use it.
173 // 177 //
174 if(pt[n - 1] != 0.0) 178 if(pt[n - 1] != 0.0)
175 { 179 {
176 xm = -pt[n] / pt[n - 1]; 180 xm = -pt[n] / pt[n - 1];
177 if(xm < x) 181 if(xm < x)
178 { 182 {
179 x = xm; 183 x = xm;
180 } 184 }
181 } 185 }
182 186
183 // Chop the interval (0,x) until ff <= 0 187 // Chop the interval (0,x) until ff <= 0
184 // 188 //
185 while(true) << 189 while(1)
186 { 190 {
187 xm = x * 0.1; 191 xm = x * 0.1;
188 ff = pt[0]; 192 ff = pt[0];
189 for(i = 1; i <= n; ++i) 193 for(i = 1; i <= n; ++i)
190 { 194 {
191 ff = ff * xm + pt[i]; 195 ff = ff * xm + pt[i];
192 } 196 }
193 if(ff <= 0.0) 197 if(ff <= 0.0)
194 { 198 {
195 break; 199 break;
196 } 200 }
197 x = xm; 201 x = xm;
198 } 202 }
199 dx = x; 203 dx = x;
200 204
201 // Do Newton interation until x converges 205 // Do Newton interation until x converges to two decimal places.
202 // 206 //
203 while(std::fabs(dx / x) > 0.005) 207 while(std::fabs(dx / x) > 0.005)
204 { 208 {
205 ff = pt[0]; 209 ff = pt[0];
206 df = ff; 210 df = ff;
207 for(i = 1; i < n; ++i) 211 for(i = 1; i < n; ++i)
208 { 212 {
209 ff = ff * x + pt[i]; 213 ff = ff * x + pt[i];
210 df = df * x + ff; 214 df = df * x + ff;
211 } 215 }
212 ff = ff * x + pt[n]; 216 ff = ff * x + pt[n];
213 dx = ff / df; 217 dx = ff / df;
214 x -= dx; 218 x -= dx;
215 } 219 }
216 bnd = x; 220 bnd = x;
217 221
218 // Compute the derivative as the initial k 222 // Compute the derivative as the initial k polynomial
219 // and do 5 steps with no shift. 223 // and do 5 steps with no shift.
220 // 224 //
221 nm1 = n - 1; 225 nm1 = n - 1;
222 for(i = 1; i < n; ++i) 226 for(i = 1; i < n; ++i)
223 { 227 {
224 k[i] = (G4double)(n - i) * p[i] / (G4dou 228 k[i] = (G4double)(n - i) * p[i] / (G4double) n;
225 } 229 }
226 k[0] = p[0]; 230 k[0] = p[0];
227 aa = p[n]; 231 aa = p[n];
228 bb = p[n - 1]; 232 bb = p[n - 1];
229 zerok = static_cast<G4int>(k[n - 1] == 0); << 233 zerok = (k[n - 1] == 0);
230 for(jj = 0; jj < 5; ++jj) 234 for(jj = 0; jj < 5; ++jj)
231 { 235 {
232 cc = k[n - 1]; 236 cc = k[n - 1];
233 if(zerok == 0) // Use a scaled form of << 237 if(!zerok) // Use a scaled form of recurrence if k at 0 is nonzero.
234 { 238 {
235 // Use a scaled form of recurrence if 239 // Use a scaled form of recurrence if value of k at 0 is nonzero.
236 // 240 //
237 t = -aa / cc; 241 t = -aa / cc;
238 for(i = 0; i < nm1; ++i) 242 for(i = 0; i < nm1; ++i)
239 { 243 {
240 j = n - i - 1; 244 j = n - i - 1;
241 k[j] = t * k[j - 1] + p[j]; 245 k[j] = t * k[j - 1] + p[j];
242 } 246 }
243 k[0] = p[0]; 247 k[0] = p[0];
244 zerok = << 248 zerok = (std::fabs(k[n - 1]) <= std::fabs(bb) * eta * 10.0);
245 static_cast<G4int>(std::fabs(k[n - 1 <<
246 } 249 }
247 else // Use unscaled form of recurrence 250 else // Use unscaled form of recurrence.
248 { 251 {
249 for(i = 0; i < nm1; ++i) 252 for(i = 0; i < nm1; ++i)
250 { 253 {
251 j = n - i - 1; 254 j = n - i - 1;
252 k[j] = k[j - 1]; 255 k[j] = k[j - 1];
253 } 256 }
254 k[0] = 0.0; 257 k[0] = 0.0;
255 zerok = static_cast<G4int>(!(k[n - 1] << 258 zerok = (!(k[n - 1] != 0.0));
256 } 259 }
257 } 260 }
258 261
259 // Save k for restarts with new shifts. 262 // Save k for restarts with new shifts.
260 // 263 //
261 for(i = 0; i < n; ++i) 264 for(i = 0; i < n; ++i)
262 { 265 {
263 temp[i] = k[i]; 266 temp[i] = k[i];
264 } 267 }
265 268
266 // Loop to select the quadratic correspond 269 // Loop to select the quadratic corresponding to each new shift.
267 // 270 //
268 for(cnt = 0; cnt < 20; ++cnt) 271 for(cnt = 0; cnt < 20; ++cnt)
269 { 272 {
270 // Quadratic corresponds to a double shi 273 // Quadratic corresponds to a double shift to a
271 // non-real point and its complex conjug 274 // non-real point and its complex conjugate. The point
272 // has modulus bnd and amplitude rotated 275 // has modulus bnd and amplitude rotated by 94 degrees
273 // from the previous shift. 276 // from the previous shift.
274 // 277 //
275 xxx = cosr * xo - sinr * yo; 278 xxx = cosr * xo - sinr * yo;
276 yo = sinr * xo + cosr * yo; 279 yo = sinr * xo + cosr * yo;
277 xo = xxx; 280 xo = xxx;
278 sr = bnd * xo; 281 sr = bnd * xo;
279 si = bnd * yo; 282 si = bnd * yo;
280 u = -2.0 * sr; 283 u = -2.0 * sr;
281 v = bnd; 284 v = bnd;
282 ComputeFixedShiftPolynomial(20 * (cnt + 285 ComputeFixedShiftPolynomial(20 * (cnt + 1), &nz);
283 if(nz != 0) 286 if(nz != 0)
284 { 287 {
285 // The second stage jumps directly to 288 // The second stage jumps directly to one of the third
286 // stage iterations and returns here i 289 // stage iterations and returns here if successful.
287 // Deflate the polynomial, store the z 290 // Deflate the polynomial, store the zero or zeros and
288 // return to the main algorithm. 291 // return to the main algorithm.
289 // 292 //
290 j = degr - n; 293 j = degr - n;
291 zeror[j] = szr; 294 zeror[j] = szr;
292 zeroi[j] = szi; 295 zeroi[j] = szi;
293 n -= nz; 296 n -= nz;
294 for(i = 0; i <= n; ++i) 297 for(i = 0; i <= n; ++i)
295 { 298 {
296 p[i] = qp[i]; 299 p[i] = qp[i];
297 } 300 }
298 if(nz != 1) 301 if(nz != 1)
299 { 302 {
300 zeror[j + 1] = lzr; 303 zeror[j + 1] = lzr;
301 zeroi[j + 1] = lzi; 304 zeroi[j + 1] = lzi;
302 } 305 }
303 break; 306 break;
304 } 307 }
305 << 308 else
306 // If the iteration is unsuccessful anot <<
307 // is chosen after restoring k. <<
308 // <<
309 for(i = 0; i < n; ++i) <<
310 { 309 {
311 k[i] = temp[i]; << 310 // If the iteration is unsuccessful another quadratic
>> 311 // is chosen after restoring k.
>> 312 //
>> 313 for(i = 0; i < n; ++i)
>> 314 {
>> 315 k[i] = temp[i];
>> 316 }
312 } 317 }
313 <<
314 } 318 }
315 } while(nz != 0); // End of initial DO loop 319 } while(nz != 0); // End of initial DO loop
316 320
317 // Return with failure if no convergence wit 321 // Return with failure if no convergence with 20 shifts.
318 // 322 //
319 return degr - n; 323 return degr - n;
320 } 324 }
321 325
322 void G4JTPolynomialSolver::ComputeFixedShiftPo 326 void G4JTPolynomialSolver::ComputeFixedShiftPolynomial(G4int l2, G4int* nz)
323 { 327 {
324 // Computes up to L2 fixed shift k-polynomia 328 // Computes up to L2 fixed shift k-polynomials, testing for convergence
325 // in the linear or quadratic case. Initiate 329 // in the linear or quadratic case. Initiates one of the variable shift
326 // iterations and returns with the number of 330 // iterations and returns with the number of zeros found.
327 331
328 G4double svu = 0.0, svv = 0.0, ui = 0.0, vi 332 G4double svu = 0.0, svv = 0.0, ui = 0.0, vi = 0.0, xs = 0.0;
329 G4double betas = 0.25, betav = 0.25, oss = s 333 G4double betas = 0.25, betav = 0.25, oss = sr, ovv = v, ss = 0.0, vv = 0.0,
330 ts = 1.0, tv = 1.0; 334 ts = 1.0, tv = 1.0;
331 G4double ots = 0.0, otv = 0.0; 335 G4double ots = 0.0, otv = 0.0;
332 G4double tvv = 1.0, tss = 1.0; 336 G4double tvv = 1.0, tss = 1.0;
333 G4int type = 0, i = 0, j = 0, iflag = 0, vpa 337 G4int type = 0, i = 0, j = 0, iflag = 0, vpass = 0, spass = 0, vtry = 0,
334 stry = 0; 338 stry = 0;
335 339
336 *nz = 0; 340 *nz = 0;
337 341
338 // Evaluate polynomial by synthetic division 342 // Evaluate polynomial by synthetic division.
339 // 343 //
340 QuadraticSyntheticDivision(n, &u, &v, p, qp, 344 QuadraticSyntheticDivision(n, &u, &v, p, qp, &a, &b);
341 ComputeScalarFactors(&type); 345 ComputeScalarFactors(&type);
342 for(j = 0; j < l2; ++j) 346 for(j = 0; j < l2; ++j)
343 { 347 {
344 // Calculate next k polynomial and estimat 348 // Calculate next k polynomial and estimate v.
345 // 349 //
346 ComputeNextPolynomial(&type); 350 ComputeNextPolynomial(&type);
347 ComputeScalarFactors(&type); 351 ComputeScalarFactors(&type);
348 ComputeNewEstimate(type, &ui, &vi); 352 ComputeNewEstimate(type, &ui, &vi);
349 vv = vi; 353 vv = vi;
350 354
351 // Estimate xs. 355 // Estimate xs.
352 // 356 //
353 ss = 0.0; 357 ss = 0.0;
354 if(k[n - 1] != 0.0) 358 if(k[n - 1] != 0.0)
355 { 359 {
356 ss = -p[n] / k[n - 1]; 360 ss = -p[n] / k[n - 1];
357 } 361 }
358 tv = 1.0; 362 tv = 1.0;
359 ts = 1.0; 363 ts = 1.0;
360 if(j == 0 || type == 3) 364 if(j == 0 || type == 3)
361 { 365 {
362 ovv = vv; 366 ovv = vv;
363 oss = ss; 367 oss = ss;
364 otv = tv; 368 otv = tv;
365 ots = ts; 369 ots = ts;
366 continue; 370 continue;
367 } 371 }
368 372
369 // Compute relative measures of convergenc 373 // Compute relative measures of convergence of xs and v sequences.
370 // 374 //
371 if(vv != 0.0) 375 if(vv != 0.0)
372 { 376 {
373 tv = std::fabs((vv - ovv) / vv); 377 tv = std::fabs((vv - ovv) / vv);
374 } 378 }
375 if(ss != 0.0) 379 if(ss != 0.0)
376 { 380 {
377 ts = std::fabs((ss - oss) / ss); 381 ts = std::fabs((ss - oss) / ss);
378 } 382 }
379 383
380 // If decreasing, multiply two most recent 384 // If decreasing, multiply two most recent convergence measures.
381 tvv = 1.0; 385 tvv = 1.0;
382 if(tv < otv) 386 if(tv < otv)
383 { 387 {
384 tvv = tv * otv; 388 tvv = tv * otv;
385 } 389 }
386 tss = 1.0; 390 tss = 1.0;
387 if(ts < ots) 391 if(ts < ots)
388 { 392 {
389 tss = ts * ots; 393 tss = ts * ots;
390 } 394 }
391 395
392 // Compare with convergence criteria. 396 // Compare with convergence criteria.
393 vpass = static_cast<G4int>(tvv < betav); << 397 vpass = (tvv < betav);
394 spass = static_cast<G4int>(tss < betas); << 398 spass = (tss < betas);
395 if(!((spass != 0) || (vpass != 0))) << 399 if(!(spass || vpass))
396 { 400 {
397 ovv = vv; 401 ovv = vv;
398 oss = ss; 402 oss = ss;
399 otv = tv; 403 otv = tv;
400 ots = ts; 404 ots = ts;
401 continue; 405 continue;
402 } 406 }
403 407
404 // At least one sequence has passed the co 408 // At least one sequence has passed the convergence test.
405 // Store variables before iterating. 409 // Store variables before iterating.
406 // 410 //
407 svu = u; 411 svu = u;
408 svv = v; 412 svv = v;
409 for(i = 0; i < n; ++i) 413 for(i = 0; i < n; ++i)
410 { 414 {
411 svk[i] = k[i]; 415 svk[i] = k[i];
412 } 416 }
413 xs = ss; 417 xs = ss;
414 418
415 // Choose iteration according to the faste 419 // Choose iteration according to the fastest converging sequence.
416 // 420 //
417 vtry = 0; 421 vtry = 0;
418 stry = 0; 422 stry = 0;
419 if(((spass != 0) && (vpass == 0)) || (tss << 423 if((spass && (!vpass)) || (tss < tvv))
420 { 424 {
421 RealPolynomialIteration(&xs, nz, &iflag) 425 RealPolynomialIteration(&xs, nz, &iflag);
422 if(*nz > 0) 426 if(*nz > 0)
423 { 427 {
424 return; 428 return;
425 } 429 }
426 430
427 // Linear iteration has failed. Flag tha 431 // Linear iteration has failed. Flag that it has been
428 // tried and decrease the convergence cr 432 // tried and decrease the convergence criterion.
429 // 433 //
430 stry = 1; 434 stry = 1;
431 betas *= 0.25; 435 betas *= 0.25;
432 if(iflag == 0) 436 if(iflag == 0)
433 { 437 {
434 goto _restore_variables; 438 goto _restore_variables;
435 } 439 }
436 440
437 // If linear iteration signals an almost 441 // If linear iteration signals an almost double real
438 // zero attempt quadratic iteration. 442 // zero attempt quadratic iteration.
439 // 443 //
440 ui = -(xs + xs); 444 ui = -(xs + xs);
441 vi = xs * xs; 445 vi = xs * xs;
442 } 446 }
443 447
444 _quadratic_iteration: 448 _quadratic_iteration:
445 449
446 do 450 do
447 { 451 {
448 QuadraticPolynomialIteration(&ui, &vi, n 452 QuadraticPolynomialIteration(&ui, &vi, nz);
449 if(*nz > 0) 453 if(*nz > 0)
450 { 454 {
451 return; 455 return;
452 } 456 }
453 457
454 // Quadratic iteration has failed. Flag 458 // Quadratic iteration has failed. Flag that it has
455 // been tried and decrease the convergen 459 // been tried and decrease the convergence criterion.
456 // 460 //
457 vtry = 1; 461 vtry = 1;
458 betav *= 0.25; 462 betav *= 0.25;
459 463
460 // Try linear iteration if it has not be 464 // Try linear iteration if it has not been tried and
461 // the S sequence is converging. 465 // the S sequence is converging.
462 // 466 //
463 if((stry != 0) || (spass == 0)) << 467 if(stry || !spass)
464 { 468 {
465 break; 469 break;
466 } 470 }
467 for(i = 0; i < n; ++i) 471 for(i = 0; i < n; ++i)
468 { 472 {
469 k[i] = svk[i]; 473 k[i] = svk[i];
470 } 474 }
471 RealPolynomialIteration(&xs, nz, &iflag) 475 RealPolynomialIteration(&xs, nz, &iflag);
472 if(*nz > 0) 476 if(*nz > 0)
473 { 477 {
474 return; 478 return;
475 } 479 }
476 480
477 // Linear iteration has failed. Flag tha 481 // Linear iteration has failed. Flag that it has been
478 // tried and decrease the convergence cr 482 // tried and decrease the convergence criterion.
479 // 483 //
480 stry = 1; 484 stry = 1;
481 betas *= 0.25; 485 betas *= 0.25;
482 if(iflag == 0) 486 if(iflag == 0)
483 { 487 {
484 break; 488 break;
485 } 489 }
486 490
487 // If linear iteration signals an almost 491 // If linear iteration signals an almost double real
488 // zero attempt quadratic iteration. 492 // zero attempt quadratic iteration.
489 // 493 //
490 ui = -(xs + xs); 494 ui = -(xs + xs);
491 vi = xs * xs; 495 vi = xs * xs;
492 } while(iflag != 0); 496 } while(iflag != 0);
493 497
494 // Restore variables. 498 // Restore variables.
495 499
496 _restore_variables: 500 _restore_variables:
497 501
498 u = svu; 502 u = svu;
499 v = svv; 503 v = svv;
500 for(i = 0; i < n; ++i) 504 for(i = 0; i < n; ++i)
501 { 505 {
502 k[i] = svk[i]; 506 k[i] = svk[i];
503 } 507 }
504 508
505 // Try quadratic iteration if it has not b 509 // Try quadratic iteration if it has not been tried
506 // and the V sequence is converging. 510 // and the V sequence is converging.
507 // 511 //
508 if((vpass != 0) && (vtry == 0)) << 512 if(vpass && !vtry)
509 { 513 {
510 goto _quadratic_iteration; 514 goto _quadratic_iteration;
511 } 515 }
512 516
513 // Recompute QP and scalar values to conti 517 // Recompute QP and scalar values to continue the
514 // second stage. 518 // second stage.
515 // 519 //
516 QuadraticSyntheticDivision(n, &u, &v, p, q 520 QuadraticSyntheticDivision(n, &u, &v, p, qp, &a, &b);
517 ComputeScalarFactors(&type); 521 ComputeScalarFactors(&type);
518 522
519 ovv = vv; 523 ovv = vv;
520 oss = ss; 524 oss = ss;
521 otv = tv; 525 otv = tv;
522 ots = ts; 526 ots = ts;
523 } 527 }
524 } 528 }
525 529
526 void G4JTPolynomialSolver::QuadraticPolynomial 530 void G4JTPolynomialSolver::QuadraticPolynomialIteration(G4double* uu,
527 531 G4double* vv, G4int* nz)
528 { 532 {
529 // Variable-shift k-polynomial iteration for 533 // Variable-shift k-polynomial iteration for a
530 // quadratic factor converges only if the ze 534 // quadratic factor converges only if the zeros are
531 // equimodular or nearly so. 535 // equimodular or nearly so.
532 // uu, vv - coefficients of starting quadrat 536 // uu, vv - coefficients of starting quadratic.
533 // nz - number of zeros found. 537 // nz - number of zeros found.
534 // 538 //
535 G4double ui = 0.0, vi = 0.0; 539 G4double ui = 0.0, vi = 0.0;
536 G4double omp = 0.0; 540 G4double omp = 0.0;
537 G4double relstp = 0.0; 541 G4double relstp = 0.0;
538 G4double mp = 0.0, ee = 0.0, t = 0.0, zm = 0 542 G4double mp = 0.0, ee = 0.0, t = 0.0, zm = 0.0;
539 G4int type = 0, i = 1, j = 0, tried = 0; 543 G4int type = 0, i = 1, j = 0, tried = 0;
540 544
541 *nz = 0; 545 *nz = 0;
542 tried = 0; 546 tried = 0;
543 u = *uu; 547 u = *uu;
544 v = *vv; 548 v = *vv;
545 549
546 // Main loop. 550 // Main loop.
547 551
548 while(true) << 552 while(1)
549 { 553 {
550 Quadratic(1.0, u, v, &szr, &szi, &lzr, &lz 554 Quadratic(1.0, u, v, &szr, &szi, &lzr, &lzi);
551 555
552 // Return if roots of the quadratic are re 556 // Return if roots of the quadratic are real and not
553 // close to multiple or nearly equal and o 557 // close to multiple or nearly equal and of opposite
554 // sign. 558 // sign.
555 // 559 //
556 if(std::fabs(std::fabs(szr) - std::fabs(lz 560 if(std::fabs(std::fabs(szr) - std::fabs(lzr)) > 0.01 * std::fabs(lzr))
557 { 561 {
558 return; 562 return;
559 } 563 }
560 564
561 // Evaluate polynomial by quadratic synthe 565 // Evaluate polynomial by quadratic synthetic division.
562 // 566 //
563 QuadraticSyntheticDivision(n, &u, &v, p, q 567 QuadraticSyntheticDivision(n, &u, &v, p, qp, &a, &b);
564 mp = std::fabs(a - szr * b) + std::fabs(sz 568 mp = std::fabs(a - szr * b) + std::fabs(szi * b);
565 569
566 // Compute a rigorous bound on the roundin 570 // Compute a rigorous bound on the rounding error in evaluating p.
567 // 571 //
568 zm = std::sqrt(std::fabs(v)); 572 zm = std::sqrt(std::fabs(v));
569 ee = 2.0 * std::fabs(qp[0]); 573 ee = 2.0 * std::fabs(qp[0]);
570 t = -szr * b; 574 t = -szr * b;
571 for(i = 1; i < n; ++i) 575 for(i = 1; i < n; ++i)
572 { 576 {
573 ee = ee * zm + std::fabs(qp[i]); 577 ee = ee * zm + std::fabs(qp[i]);
574 } 578 }
575 ee = ee * zm + std::fabs(a + t); 579 ee = ee * zm + std::fabs(a + t);
576 ee *= (5.0 * mre + 4.0 * are); 580 ee *= (5.0 * mre + 4.0 * are);
577 ee = ee - (5.0 * mre + 2.0 * are) * (std:: 581 ee = ee - (5.0 * mre + 2.0 * are) * (std::fabs(a + t) + std::fabs(b) * zm) +
578 2.0 * are * std::fabs(t); 582 2.0 * are * std::fabs(t);
579 583
580 // Iteration has converged sufficiently if 584 // Iteration has converged sufficiently if the
581 // polynomial value is less than 20 times 585 // polynomial value is less than 20 times this bound.
582 // 586 //
583 if(mp <= 20.0 * ee) 587 if(mp <= 20.0 * ee)
584 { 588 {
585 *nz = 2; 589 *nz = 2;
586 return; 590 return;
587 } 591 }
588 j++; 592 j++;
589 593
590 // Stop iteration after 20 steps. 594 // Stop iteration after 20 steps.
591 // 595 //
592 if(j > 20) 596 if(j > 20)
593 { 597 {
594 return; 598 return;
595 } 599 }
596 if(j >= 2) 600 if(j >= 2)
597 { 601 {
598 if(!(relstp > 0.01 || mp < omp || (tried << 602 if(!(relstp > 0.01 || mp < omp || tried))
599 { 603 {
600 // A cluster appears to be stalling th 604 // A cluster appears to be stalling the convergence.
601 // Five fixed shift steps are taken wi 605 // Five fixed shift steps are taken with a u,v close to the cluster.
602 // 606 //
603 if(relstp < eta) 607 if(relstp < eta)
604 { 608 {
605 relstp = eta; 609 relstp = eta;
606 } 610 }
607 relstp = std::sqrt(relstp); 611 relstp = std::sqrt(relstp);
608 u = u - u * relstp; 612 u = u - u * relstp;
609 v = v + v * relstp; 613 v = v + v * relstp;
610 QuadraticSyntheticDivision(n, &u, &v, 614 QuadraticSyntheticDivision(n, &u, &v, p, qp, &a, &b);
611 for(i = 0; i < 5; ++i) 615 for(i = 0; i < 5; ++i)
612 { 616 {
613 ComputeScalarFactors(&type); 617 ComputeScalarFactors(&type);
614 ComputeNextPolynomial(&type); 618 ComputeNextPolynomial(&type);
615 } 619 }
616 tried = 1; 620 tried = 1;
617 j = 0; 621 j = 0;
618 } 622 }
619 } 623 }
620 omp = mp; 624 omp = mp;
621 625
622 // Calculate next k polynomial and new u a 626 // Calculate next k polynomial and new u and v.
623 // 627 //
624 ComputeScalarFactors(&type); 628 ComputeScalarFactors(&type);
625 ComputeNextPolynomial(&type); 629 ComputeNextPolynomial(&type);
626 ComputeScalarFactors(&type); 630 ComputeScalarFactors(&type);
627 ComputeNewEstimate(type, &ui, &vi); 631 ComputeNewEstimate(type, &ui, &vi);
628 632
629 // If vi is zero the iteration is not conv 633 // If vi is zero the iteration is not converging.
630 // 634 //
631 if(!(vi != 0.0)) 635 if(!(vi != 0.0))
632 { 636 {
633 return; 637 return;
634 } 638 }
635 relstp = std::fabs((vi - v) / vi); 639 relstp = std::fabs((vi - v) / vi);
636 u = ui; 640 u = ui;
637 v = vi; 641 v = vi;
638 } 642 }
639 } 643 }
640 644
641 void G4JTPolynomialSolver::RealPolynomialItera 645 void G4JTPolynomialSolver::RealPolynomialIteration(G4double* sss, G4int* nz,
642 646 G4int* iflag)
643 { 647 {
644 // Variable-shift H polynomial iteration for 648 // Variable-shift H polynomial iteration for a real zero.
645 // sss - starting iterate 649 // sss - starting iterate
646 // nz - number of zeros found 650 // nz - number of zeros found
647 // iflag - flag to indicate a pair of zeros 651 // iflag - flag to indicate a pair of zeros near real axis.
648 652
649 G4double t = 0.; 653 G4double t = 0.;
650 G4double omp = 0.; 654 G4double omp = 0.;
651 G4double pv = 0.0, kv = 0.0, xs = *sss; 655 G4double pv = 0.0, kv = 0.0, xs = *sss;
652 G4double mx = 0.0, mp = 0.0, ee = 0.0; 656 G4double mx = 0.0, mp = 0.0, ee = 0.0;
653 G4int i = 1, j = 0; 657 G4int i = 1, j = 0;
654 658
655 *nz = 0; 659 *nz = 0;
656 *iflag = 0; 660 *iflag = 0;
657 661
658 // Main loop 662 // Main loop
659 // 663 //
660 while(true) << 664 while(1)
661 { 665 {
662 pv = p[0]; 666 pv = p[0];
663 667
664 // Evaluate p at xs. 668 // Evaluate p at xs.
665 // 669 //
666 qp[0] = pv; 670 qp[0] = pv;
667 for(i = 1; i <= n; ++i) 671 for(i = 1; i <= n; ++i)
668 { 672 {
669 pv = pv * xs + p[i]; 673 pv = pv * xs + p[i];
670 qp[i] = pv; 674 qp[i] = pv;
671 } 675 }
672 mp = std::fabs(pv); 676 mp = std::fabs(pv);
673 677
674 // Compute a rigorous bound on the error i 678 // Compute a rigorous bound on the error in evaluating p.
675 // 679 //
676 mx = std::fabs(xs); 680 mx = std::fabs(xs);
677 ee = (mre / (are + mre)) * std::fabs(qp[0] 681 ee = (mre / (are + mre)) * std::fabs(qp[0]);
678 for(i = 1; i <= n; ++i) 682 for(i = 1; i <= n; ++i)
679 { 683 {
680 ee = ee * mx + std::fabs(qp[i]); 684 ee = ee * mx + std::fabs(qp[i]);
681 } 685 }
682 686
683 // Iteration has converged sufficiently if 687 // Iteration has converged sufficiently if the polynomial
684 // value is less than 20 times this bound. 688 // value is less than 20 times this bound.
685 // 689 //
686 if(mp <= 20.0 * ((are + mre) * ee - mre * 690 if(mp <= 20.0 * ((are + mre) * ee - mre * mp))
687 { 691 {
688 *nz = 1; 692 *nz = 1;
689 szr = xs; 693 szr = xs;
690 szi = 0.0; 694 szi = 0.0;
691 return; 695 return;
692 } 696 }
693 j++; 697 j++;
694 698
695 // Stop iteration after 10 steps. 699 // Stop iteration after 10 steps.
696 // 700 //
697 if(j > 10) 701 if(j > 10)
698 { 702 {
699 return; 703 return;
700 } 704 }
701 if(j >= 2) 705 if(j >= 2)
702 { 706 {
703 if(!(std::fabs(t) > 0.001 * std::fabs(xs 707 if(!(std::fabs(t) > 0.001 * std::fabs(xs - t) || mp < omp))
704 { 708 {
705 // A cluster of zeros near the real ax 709 // A cluster of zeros near the real axis has been encountered.
706 // Return with iflag set to initiate a 710 // Return with iflag set to initiate a quadratic iteration.
707 // 711 //
708 *iflag = 1; 712 *iflag = 1;
709 *sss = xs; 713 *sss = xs;
710 return; 714 return;
711 } // Return if the polynomial value has 715 } // Return if the polynomial value has increased significantly.
712 } 716 }
713 717
714 omp = mp; 718 omp = mp;
715 719
716 // Compute t, the next polynomial, and th 720 // Compute t, the next polynomial, and the new iterate.
717 // 721 //
718 kv = k[0]; 722 kv = k[0];
719 qk[0] = kv; 723 qk[0] = kv;
720 for(i = 1; i < n; ++i) 724 for(i = 1; i < n; ++i)
721 { 725 {
722 kv = kv * xs + k[i]; 726 kv = kv * xs + k[i];
723 qk[i] = kv; 727 qk[i] = kv;
724 } 728 }
725 if(std::fabs(kv) <= std::fabs(k[n - 1]) * 729 if(std::fabs(kv) <= std::fabs(k[n - 1]) * 10.0 * eta) // Use unscaled form.
726 { 730 {
727 k[0] = 0.0; 731 k[0] = 0.0;
728 for(i = 1; i < n; ++i) 732 for(i = 1; i < n; ++i)
729 { 733 {
730 k[i] = qk[i - 1]; 734 k[i] = qk[i - 1];
731 } 735 }
732 } 736 }
733 else // Use the scaled form of the recurr 737 else // Use the scaled form of the recurrence if k at xs is nonzero.
734 { 738 {
735 t = -pv / kv; 739 t = -pv / kv;
736 k[0] = qp[0]; 740 k[0] = qp[0];
737 for(i = 1; i < n; ++i) 741 for(i = 1; i < n; ++i)
738 { 742 {
739 k[i] = t * qk[i - 1] + qp[i]; 743 k[i] = t * qk[i - 1] + qp[i];
740 } 744 }
741 } 745 }
742 kv = k[0]; 746 kv = k[0];
743 for(i = 1; i < n; ++i) 747 for(i = 1; i < n; ++i)
744 { 748 {
745 kv = kv * xs + k[i]; 749 kv = kv * xs + k[i];
746 } 750 }
747 t = 0.0; 751 t = 0.0;
748 if(std::fabs(kv) > std::fabs(k[n - 1] * 10 752 if(std::fabs(kv) > std::fabs(k[n - 1] * 10.0 * eta))
749 { 753 {
750 t = -pv / kv; 754 t = -pv / kv;
751 } 755 }
752 xs += t; 756 xs += t;
753 } 757 }
754 } 758 }
755 759
756 void G4JTPolynomialSolver::ComputeScalarFactor 760 void G4JTPolynomialSolver::ComputeScalarFactors(G4int* type)
757 { 761 {
758 // This function calculates scalar quantitie 762 // This function calculates scalar quantities used to
759 // compute the next k polynomial and new est 763 // compute the next k polynomial and new estimates of
760 // the quadratic coefficients. 764 // the quadratic coefficients.
761 // type - integer variable set here indicati 765 // type - integer variable set here indicating how the
762 // calculations are normalized to avoid over 766 // calculations are normalized to avoid overflow.
763 767
764 // Synthetic division of k by the quadratic 768 // Synthetic division of k by the quadratic 1,u,v
765 // 769 //
766 QuadraticSyntheticDivision(n - 1, &u, &v, k, 770 QuadraticSyntheticDivision(n - 1, &u, &v, k, qk, &c, &d);
767 if(std::fabs(c) <= std::fabs(k[n - 1] * 100. 771 if(std::fabs(c) <= std::fabs(k[n - 1] * 100.0 * eta))
768 { 772 {
769 if(std::fabs(d) <= std::fabs(k[n - 2] * 10 773 if(std::fabs(d) <= std::fabs(k[n - 2] * 100.0 * eta))
770 { 774 {
771 *type = 3; // Type=3 indicates the quad 775 *type = 3; // Type=3 indicates the quadratic is almost a factor of k.
772 return; 776 return;
773 } 777 }
774 } 778 }
775 779
776 if(std::fabs(d) < std::fabs(c)) 780 if(std::fabs(d) < std::fabs(c))
777 { 781 {
778 *type = 1; // Type=1 indicates that all f 782 *type = 1; // Type=1 indicates that all formulas are divided by c.
779 e = a / c; 783 e = a / c;
780 f = d / c; 784 f = d / c;
781 g = u * e; 785 g = u * e;
782 h = v * b; 786 h = v * b;
783 a3 = a * e + (h / c + g) * b; 787 a3 = a * e + (h / c + g) * b;
784 a1 = b - a * (d / c); 788 a1 = b - a * (d / c);
785 a7 = a + g * d + h * f; 789 a7 = a + g * d + h * f;
786 return; 790 return;
787 } 791 }
788 *type = 2; // Type=2 indicates that all for 792 *type = 2; // Type=2 indicates that all formulas are divided by d.
789 e = a / d; 793 e = a / d;
790 f = c / d; 794 f = c / d;
791 g = u * b; 795 g = u * b;
792 h = v * b; 796 h = v * b;
793 a3 = (a + g) * e + h * (b / d); 797 a3 = (a + g) * e + h * (b / d);
794 a1 = b * f - a; 798 a1 = b * f - a;
795 a7 = (f + u) * a + h; 799 a7 = (f + u) * a + h;
796 } 800 }
797 801
798 void G4JTPolynomialSolver::ComputeNextPolynomi 802 void G4JTPolynomialSolver::ComputeNextPolynomial(G4int* type)
799 { 803 {
800 // Computes the next k polynomials using sca 804 // Computes the next k polynomials using scalars
801 // computed in ComputeScalarFactors. 805 // computed in ComputeScalarFactors.
802 806
803 G4int i = 2; 807 G4int i = 2;
804 808
805 if(*type == 3) // Use unscaled form of the 809 if(*type == 3) // Use unscaled form of the recurrence if type is 3.
806 { 810 {
807 k[0] = 0.0; 811 k[0] = 0.0;
808 k[1] = 0.0; 812 k[1] = 0.0;
809 for(i = 2; i < n; ++i) 813 for(i = 2; i < n; ++i)
810 { 814 {
811 k[i] = qk[i - 2]; 815 k[i] = qk[i - 2];
812 } 816 }
813 return; 817 return;
814 } 818 }
815 G4double temp = a; 819 G4double temp = a;
816 if(*type == 1) 820 if(*type == 1)
817 { 821 {
818 temp = b; 822 temp = b;
819 } 823 }
820 if(std::fabs(a1) <= std::fabs(temp) * eta * 824 if(std::fabs(a1) <= std::fabs(temp) * eta * 10.0)
821 { 825 {
822 // If a1 is nearly zero then use a special 826 // If a1 is nearly zero then use a special form of the recurrence.
823 // 827 //
824 k[0] = 0.0; 828 k[0] = 0.0;
825 k[1] = -a7 * qp[0]; 829 k[1] = -a7 * qp[0];
826 for(i = 2; i < n; ++i) 830 for(i = 2; i < n; ++i)
827 { 831 {
828 k[i] = a3 * qk[i - 2] - a7 * qp[i - 1]; 832 k[i] = a3 * qk[i - 2] - a7 * qp[i - 1];
829 } 833 }
830 return; 834 return;
831 } 835 }
832 836
833 // Use scaled form of the recurrence. 837 // Use scaled form of the recurrence.
834 // 838 //
835 a7 /= a1; 839 a7 /= a1;
836 a3 /= a1; 840 a3 /= a1;
837 k[0] = qp[0]; 841 k[0] = qp[0];
838 k[1] = qp[1] - a7 * qp[0]; 842 k[1] = qp[1] - a7 * qp[0];
839 for(i = 2; i < n; ++i) 843 for(i = 2; i < n; ++i)
840 { 844 {
841 k[i] = a3 * qk[i - 2] - a7 * qp[i - 1] + q 845 k[i] = a3 * qk[i - 2] - a7 * qp[i - 1] + qp[i];
842 } 846 }
843 } 847 }
844 848
845 void G4JTPolynomialSolver::ComputeNewEstimate( 849 void G4JTPolynomialSolver::ComputeNewEstimate(G4int type, G4double* uu,
846 850 G4double* vv)
847 { 851 {
848 // Compute new estimates of the quadratic co 852 // Compute new estimates of the quadratic coefficients
849 // using the scalars computed in calcsc. 853 // using the scalars computed in calcsc.
850 854
851 G4double a4 = 0.0, a5 = 0.0, b1 = 0.0, b2 = 855 G4double a4 = 0.0, a5 = 0.0, b1 = 0.0, b2 = 0.0, c1 = 0.0, c2 = 0.0, c3 = 0.0,
852 c4 = 0.0, temp = 0.0; 856 c4 = 0.0, temp = 0.0;
853 857
854 // Use formulas appropriate to setting of ty 858 // Use formulas appropriate to setting of type.
855 // 859 //
856 if(type == 3) // If type=3 the quadratic i 860 if(type == 3) // If type=3 the quadratic is zeroed.
857 { 861 {
858 *uu = 0.0; 862 *uu = 0.0;
859 *vv = 0.0; 863 *vv = 0.0;
860 return; 864 return;
861 } 865 }
862 if(type == 2) 866 if(type == 2)
863 { 867 {
864 a4 = (a + g) * f + h; 868 a4 = (a + g) * f + h;
865 a5 = (f + u) * c + v * d; 869 a5 = (f + u) * c + v * d;
866 } 870 }
867 else 871 else
868 { 872 {
869 a4 = a + u * b + h * f; 873 a4 = a + u * b + h * f;
870 a5 = c + (u + v * f) * d; 874 a5 = c + (u + v * f) * d;
871 } 875 }
872 876
873 // Evaluate new quadratic coefficients. 877 // Evaluate new quadratic coefficients.
874 // 878 //
875 b1 = -k[n - 1] / p[n]; 879 b1 = -k[n - 1] / p[n];
876 b2 = -(k[n - 2] + b1 * p[n - 1]) / p[n]; 880 b2 = -(k[n - 2] + b1 * p[n - 1]) / p[n];
877 c1 = v * b2 * a1; 881 c1 = v * b2 * a1;
878 c2 = b1 * a7; 882 c2 = b1 * a7;
879 c3 = b1 * b1 * a3; 883 c3 = b1 * b1 * a3;
880 c4 = c1 - c2 - c3; 884 c4 = c1 - c2 - c3;
881 temp = a5 + b1 * a4 - c4; 885 temp = a5 + b1 * a4 - c4;
882 if(!(temp != 0.0)) 886 if(!(temp != 0.0))
883 { 887 {
884 *uu = 0.0; 888 *uu = 0.0;
885 *vv = 0.0; 889 *vv = 0.0;
886 return; 890 return;
887 } 891 }
888 *uu = u - (u * (c3 + c2) + v * (b1 * a1 + b2 892 *uu = u - (u * (c3 + c2) + v * (b1 * a1 + b2 * a7)) / temp;
889 *vv = v * (1.0 + c4 / temp); 893 *vv = v * (1.0 + c4 / temp);
890 return; 894 return;
891 } 895 }
892 896
893 void G4JTPolynomialSolver::QuadraticSyntheticD 897 void G4JTPolynomialSolver::QuadraticSyntheticDivision(
894 G4int nn, G4double* uu, G4double* vv, std::v 898 G4int nn, G4double* uu, G4double* vv, std::vector<G4double>& pp,
895 std::vector<G4double>& qq, G4double* aa, G4d 899 std::vector<G4double>& qq, G4double* aa, G4double* bb)
896 { 900 {
897 // Divides pp by the quadratic 1,uu,vv placi 901 // Divides pp by the quadratic 1,uu,vv placing the quotient
898 // in qq and the remainder in aa,bb. 902 // in qq and the remainder in aa,bb.
899 903
900 G4double cc = 0.0; 904 G4double cc = 0.0;
901 *bb = pp[0]; 905 *bb = pp[0];
902 qq[0] = *bb; 906 qq[0] = *bb;
903 *aa = pp[1] - (*bb) * (*uu); 907 *aa = pp[1] - (*bb) * (*uu);
904 qq[1] = *aa; 908 qq[1] = *aa;
905 for(G4int i = 2; i <= nn; ++i) 909 for(G4int i = 2; i <= nn; ++i)
906 { 910 {
907 cc = pp[i] - (*aa) * (*uu) - (*bb) * (* 911 cc = pp[i] - (*aa) * (*uu) - (*bb) * (*vv);
908 qq[i] = cc; 912 qq[i] = cc;
909 *bb = *aa; 913 *bb = *aa;
910 *aa = cc; 914 *aa = cc;
911 } 915 }
912 } 916 }
913 917
914 void G4JTPolynomialSolver::Quadratic(G4double 918 void G4JTPolynomialSolver::Quadratic(G4double aa, G4double b1, G4double cc,
915 G4double* 919 G4double* ssr, G4double* ssi, G4double* lr,
916 G4double* 920 G4double* li)
917 { 921 {
918 // Calculate the zeros of the quadratic aa*z 922 // Calculate the zeros of the quadratic aa*z^2 + b1*z + cc.
919 // The quadratic formula, modified to avoid 923 // The quadratic formula, modified to avoid overflow, is used
920 // to find the larger zero if the zeros are 924 // to find the larger zero if the zeros are real and both
921 // are complex. The smaller real zero is fou 925 // are complex. The smaller real zero is found directly from
922 // the product of the zeros c/a. 926 // the product of the zeros c/a.
923 927
924 G4double bb = 0.0, dd = 0.0, ee = 0.0; 928 G4double bb = 0.0, dd = 0.0, ee = 0.0;
925 929
926 if(!(aa != 0.0)) // less than two roots 930 if(!(aa != 0.0)) // less than two roots
927 { 931 {
928 if(b1 != 0.0) 932 if(b1 != 0.0)
929 { 933 {
930 *ssr = -cc / b1; 934 *ssr = -cc / b1;
931 } 935 }
932 else 936 else
933 { 937 {
934 *ssr = 0.0; 938 *ssr = 0.0;
935 } 939 }
936 *lr = 0.0; 940 *lr = 0.0;
937 *ssi = 0.0; 941 *ssi = 0.0;
938 *li = 0.0; 942 *li = 0.0;
939 return; 943 return;
940 } 944 }
941 if(!(cc != 0.0)) // one real root, one zero 945 if(!(cc != 0.0)) // one real root, one zero root
942 { 946 {
943 *ssr = 0.0; 947 *ssr = 0.0;
944 *lr = -b1 / aa; 948 *lr = -b1 / aa;
945 *ssi = 0.0; 949 *ssi = 0.0;
946 *li = 0.0; 950 *li = 0.0;
947 return; 951 return;
948 } 952 }
949 953
950 // Compute discriminant avoiding overflow. 954 // Compute discriminant avoiding overflow.
951 // 955 //
952 bb = b1 / 2.0; 956 bb = b1 / 2.0;
953 if(std::fabs(bb) < std::fabs(cc)) 957 if(std::fabs(bb) < std::fabs(cc))
954 { 958 {
955 if(cc < 0.0) 959 if(cc < 0.0)
956 { 960 {
957 ee = -aa; 961 ee = -aa;
958 } 962 }
959 else 963 else
960 { 964 {
961 ee = aa; 965 ee = aa;
962 } 966 }
963 ee = bb * (bb / std::fabs(cc)) - ee; 967 ee = bb * (bb / std::fabs(cc)) - ee;
964 dd = std::sqrt(std::fabs(ee)) * std::sqrt( 968 dd = std::sqrt(std::fabs(ee)) * std::sqrt(std::fabs(cc));
965 } 969 }
966 else 970 else
967 { 971 {
968 ee = 1.0 - (aa / bb) * (cc / bb); 972 ee = 1.0 - (aa / bb) * (cc / bb);
969 dd = std::sqrt(std::fabs(ee)) * std::fabs( 973 dd = std::sqrt(std::fabs(ee)) * std::fabs(bb);
970 } 974 }
971 if(ee < 0.0) // complex conjugate zeros 975 if(ee < 0.0) // complex conjugate zeros
972 { 976 {
973 *ssr = -bb / aa; 977 *ssr = -bb / aa;
974 *lr = *ssr; 978 *lr = *ssr;
975 *ssi = std::fabs(dd / aa); 979 *ssi = std::fabs(dd / aa);
976 *li = -(*ssi); 980 *li = -(*ssi);
977 } 981 }
978 else 982 else
979 { 983 {
980 if(bb >= 0.0) // real zeros. 984 if(bb >= 0.0) // real zeros.
981 { 985 {
982 dd = -dd; 986 dd = -dd;
983 } 987 }
984 *lr = (-bb + dd) / aa; 988 *lr = (-bb + dd) / aa;
985 *ssr = 0.0; 989 *ssr = 0.0;
986 if(*lr != 0.0) 990 if(*lr != 0.0)
987 { 991 {
988 *ssr = (cc / *lr) / aa; 992 *ssr = (cc / *lr) / aa;
989 } 993 }
990 *ssi = 0.0; 994 *ssi = 0.0;
991 *li = 0.0; 995 *li = 0.0;
992 } 996 }
993 } 997 }
994 998