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