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