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