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