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