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