00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4JTPolynomialSolver.hh"
00037 #include "G4SystemOfUnits.hh"
00038
00039 const G4double G4JTPolynomialSolver::base = 2;
00040 const G4double G4JTPolynomialSolver::eta = DBL_EPSILON;
00041 const G4double G4JTPolynomialSolver::infin = DBL_MAX;
00042 const G4double G4JTPolynomialSolver::smalno = DBL_MIN;
00043 const G4double G4JTPolynomialSolver::are = DBL_EPSILON;
00044 const G4double G4JTPolynomialSolver::mre = DBL_EPSILON;
00045 const G4double G4JTPolynomialSolver::lo = DBL_MIN/DBL_EPSILON ;
00046
00047 G4JTPolynomialSolver::G4JTPolynomialSolver()
00048 : sr(0.), si(0.), u(0.),v(0.),
00049 a(0.), b(0.), c(0.), d(0.),
00050 a1(0.), a3(0.), a7(0.),
00051 e(0.), f(0.), g(0.), h(0.),
00052 szr(0.), szi(0.), lzr(0.), lzi(0.),
00053 n(0)
00054 {
00055 }
00056
00057 G4JTPolynomialSolver::~G4JTPolynomialSolver()
00058 {
00059 }
00060
00061 G4int G4JTPolynomialSolver::FindRoots(G4double *op, G4int degr,
00062 G4double *zeror, G4double *zeroi)
00063 {
00064 G4double t=0.0, aa=0.0, bb=0.0, cc=0.0, factor=1.0;
00065 G4double max=0.0, min=infin, xxx=0.0, x=0.0, sc=0.0, bnd=0.0;
00066 G4double xm=0.0, ff=0.0, df=0.0, dx=0.0;
00067 G4int cnt=0, nz=0, i=0, j=0, jj=0, l=0, nm1=0, zerok=0;
00068
00069
00070
00071 G4double xx = std::sqrt(0.5);
00072 G4double yy = -xx,
00073 rot = 94.0*deg;
00074 G4double cosr = std::cos(rot),
00075 sinr = std::sin(rot);
00076 n = degr;
00077
00078
00079
00080 if (!(op[0] != 0.0)) { return -1; }
00081
00082
00083
00084 while (!(op[n] != 0.0))
00085 {
00086 j = degr - n;
00087 zeror[j] = 0.0;
00088 zeroi[j] = 0.0;
00089 n--;
00090 }
00091 if (n < 1) { return -1; }
00092
00093
00094
00095 std::vector<G4double> temp(degr+1) ;
00096 std::vector<G4double> pt(degr+1) ;
00097
00098 p.assign(degr+1,0) ;
00099 qp.assign(degr+1,0) ;
00100 k.assign(degr+1,0) ;
00101 qk.assign(degr+1,0) ;
00102 svk.assign(degr+1,0) ;
00103
00104
00105
00106 for (i=0;i<=n;i++)
00107 { p[i] = op[i]; }
00108
00109 do
00110 {
00111 if (n == 1)
00112 {
00113 zeror[degr-1] = -p[1]/p[0];
00114 zeroi[degr-1] = 0.0;
00115 n -= 1;
00116 return degr - n ;
00117 }
00118 if (n == 2)
00119 {
00120 Quadratic(p[0],p[1],p[2],&zeror[degr-2],&zeroi[degr-2],
00121 &zeror[degr-1],&zeroi[degr-1]);
00122 n -= 2;
00123 return degr - n ;
00124 }
00125
00126
00127
00128 max = 0.0;
00129 min = infin;
00130 for (i=0;i<=n;i++)
00131 {
00132 x = std::fabs(p[i]);
00133 if (x > max) { max = x; }
00134 if (x != 0.0 && x < min) { min = x; }
00135 }
00136
00137
00138
00139
00140
00141
00142
00143 sc = lo/min;
00144
00145 if ( ((sc <= 1.0) && (max >= 10.0))
00146 || ((sc > 1.0) && (infin/sc >= max))
00147 || ((infin/sc >= max) && (max >= 10)) )
00148 {
00149 if (!( sc != 0.0 ))
00150 { sc = smalno ; }
00151 l = (G4int)(std::log(sc)/std::log(base) + 0.5);
00152 factor = std::pow(base*1.0,l);
00153 if (factor != 1.0)
00154 {
00155 for (i=0;i<=n;i++)
00156 { p[i] = factor*p[i]; }
00157 }
00158 }
00159
00160
00161
00162 for (i=0;i<=n;i++)
00163 {
00164 pt[i] = (std::fabs(p[i]));
00165 }
00166 pt[n] = - pt[n];
00167
00168
00169
00170 x = std::exp((std::log(-pt[n])-std::log(pt[0])) / (G4double)n);
00171
00172
00173
00174 if (pt[n-1] != 0.0)
00175 {
00176 xm = -pt[n]/pt[n-1];
00177 if (xm < x) { x = xm; }
00178 }
00179
00180
00181
00182 while (1)
00183 {
00184 xm = x*0.1;
00185 ff = pt[0];
00186 for (i=1;i<=n;i++)
00187 { ff = ff*xm + pt[i]; }
00188 if (ff <= 0.0) { break; }
00189 x = xm;
00190 }
00191 dx = x;
00192
00193
00194
00195 while (std::fabs(dx/x) > 0.005)
00196 {
00197 ff = pt[0];
00198 df = ff;
00199 for (i=1;i<n;i++)
00200 {
00201 ff = ff*x + pt[i];
00202 df = df*x + ff;
00203 }
00204 ff = ff*x + pt[n];
00205 dx = ff/df;
00206 x -= dx;
00207 }
00208 bnd = x;
00209
00210
00211
00212
00213 nm1 = n - 1;
00214 for (i=1;i<n;i++)
00215 { k[i] = (G4double)(n-i)*p[i]/(G4double)n; }
00216 k[0] = p[0];
00217 aa = p[n];
00218 bb = p[n-1];
00219 zerok = (k[n-1] == 0);
00220 for(jj=0;jj<5;jj++)
00221 {
00222 cc = k[n-1];
00223 if (!zerok)
00224 {
00225
00226
00227 t = -aa/cc;
00228 for (i=0;i<nm1;i++)
00229 {
00230 j = n-i-1;
00231 k[j] = t*k[j-1]+p[j];
00232 }
00233 k[0] = p[0];
00234 zerok = (std::fabs(k[n-1]) <= std::fabs(bb)*eta*10.0);
00235 }
00236 else
00237 {
00238 for (i=0;i<nm1;i++)
00239 {
00240 j = n-i-1;
00241 k[j] = k[j-1];
00242 }
00243 k[0] = 0.0;
00244 zerok = (!(k[n-1] != 0.0));
00245 }
00246 }
00247
00248
00249
00250 for (i=0;i<n;i++)
00251 { temp[i] = k[i]; }
00252
00253
00254
00255 for (cnt = 0;cnt < 20;cnt++)
00256 {
00257
00258
00259
00260
00261
00262 xxx = cosr*xx - sinr*yy;
00263 yy = sinr*xx + cosr*yy;
00264 xx = xxx;
00265 sr = bnd*xx;
00266 si = bnd*yy;
00267 u = -2.0 * sr;
00268 v = bnd;
00269 ComputeFixedShiftPolynomial(20*(cnt+1),&nz);
00270 if (nz != 0)
00271 {
00272
00273
00274
00275
00276
00277 j = degr - n;
00278 zeror[j] = szr;
00279 zeroi[j] = szi;
00280 n -= nz;
00281 for (i=0;i<=n;i++)
00282 { p[i] = qp[i]; }
00283 if (nz != 1)
00284 {
00285 zeror[j+1] = lzr;
00286 zeroi[j+1] = lzi;
00287 }
00288 break;
00289 }
00290 else
00291 {
00292
00293
00294
00295 for (i=0;i<n;i++)
00296 {
00297 k[i] = temp[i];
00298 }
00299 }
00300 }
00301 }
00302 while (nz != 0);
00303
00304
00305
00306 return degr - n;
00307 }
00308
00309 void G4JTPolynomialSolver::ComputeFixedShiftPolynomial(G4int l2, G4int *nz)
00310 {
00311
00312
00313
00314
00315 G4double svu=0.0, svv=0.0, ui=0.0, vi=0.0, xs=0.0;
00316 G4double betas=0.25, betav=0.25, oss=sr, ovv=v,
00317 ss=0.0, vv=0.0, ts=1.0, tv=1.0;
00318 G4double ots=0.0, otv=0.0;
00319 G4double tvv=1.0, tss=1.0;
00320 G4int type=0, i=0, j=0, iflag=0, vpass=0, spass=0, vtry=0, stry=0;
00321
00322 *nz = 0;
00323
00324
00325
00326 QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
00327 ComputeScalarFactors(&type);
00328 for (j=0;j<l2;j++)
00329 {
00330
00331
00332 ComputeNextPolynomial(&type);
00333 ComputeScalarFactors(&type);
00334 ComputeNewEstimate(type,&ui,&vi);
00335 vv = vi;
00336
00337
00338
00339 ss = 0.0;
00340 if (k[n-1] != 0.0) { ss = -p[n]/k[n-1]; }
00341 tv = 1.0;
00342 ts = 1.0;
00343 if (j == 0 || type == 3)
00344 {
00345 ovv = vv;
00346 oss = ss;
00347 otv = tv;
00348 ots = ts;
00349 continue;
00350 }
00351
00352
00353
00354 if (vv != 0.0) { tv = std::fabs((vv-ovv)/vv); }
00355 if (ss != 0.0) { ts = std::fabs((ss-oss)/ss); }
00356
00357
00358 tvv = 1.0;
00359 if (tv < otv) { tvv = tv*otv; }
00360 tss = 1.0;
00361 if (ts < ots) { tss = ts*ots; }
00362
00363
00364 vpass = (tvv < betav);
00365 spass = (tss < betas);
00366 if (!(spass || vpass))
00367 {
00368 ovv = vv;
00369 oss = ss;
00370 otv = tv;
00371 ots = ts;
00372 continue;
00373 }
00374
00375
00376
00377
00378 svu = u;
00379 svv = v;
00380 for (i=0;i<n;i++)
00381 {
00382 svk[i] = k[i];
00383 }
00384 xs = ss;
00385
00386
00387
00388 vtry = 0;
00389 stry = 0;
00390 if ((spass && (!vpass)) || (tss < tvv))
00391 {
00392 RealPolynomialIteration(&xs,nz,&iflag);
00393 if (*nz > 0) { return; }
00394
00395
00396
00397
00398 stry = 1;
00399 betas *=0.25;
00400 if (iflag == 0) { goto _restore_variables; }
00401
00402
00403
00404
00405 ui = -(xs+xs);
00406 vi = xs*xs;
00407 }
00408
00409 _quadratic_iteration:
00410
00411 do
00412 {
00413 QuadraticPolynomialIteration(&ui,&vi,nz);
00414 if (*nz > 0) { return; }
00415
00416
00417
00418
00419 vtry = 1;
00420 betav *= 0.25;
00421
00422
00423
00424
00425 if (stry || !spass) { break; }
00426 for (i=0;i<n;i++)
00427 {
00428 k[i] = svk[i];
00429 }
00430 RealPolynomialIteration(&xs,nz,&iflag);
00431 if (*nz > 0) { return; }
00432
00433
00434
00435
00436 stry = 1;
00437 betas *=0.25;
00438 if (iflag == 0) { break; }
00439
00440
00441
00442
00443 ui = -(xs+xs);
00444 vi = xs*xs;
00445 }
00446 while (iflag != 0);
00447
00448
00449
00450 _restore_variables:
00451
00452 u = svu;
00453 v = svv;
00454 for (i=0;i<n;i++)
00455 {
00456 k[i] = svk[i];
00457 }
00458
00459
00460
00461
00462 if (vpass && !vtry) { goto _quadratic_iteration; }
00463
00464
00465
00466
00467 QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
00468 ComputeScalarFactors(&type);
00469
00470 ovv = vv;
00471 oss = ss;
00472 otv = tv;
00473 ots = ts;
00474 }
00475 }
00476
00477 void G4JTPolynomialSolver::
00478 QuadraticPolynomialIteration(G4double *uu, G4double *vv, G4int *nz)
00479 {
00480
00481
00482
00483
00484
00485
00486 G4double ui=0.0, vi=0.0;
00487 G4double omp=0.0;
00488 G4double relstp=0.0;
00489 G4double mp=0.0, ee=0.0, t=0.0, zm=0.0;
00490 G4int type=0, i=1, j=0, tried=0;
00491
00492 *nz = 0;
00493 tried = 0;
00494 u = *uu;
00495 v = *vv;
00496
00497
00498
00499 while (1)
00500 {
00501 Quadratic(1.0,u,v,&szr,&szi,&lzr,&lzi);
00502
00503
00504
00505
00506
00507 if (std::fabs(std::fabs(szr)-std::fabs(lzr)) > 0.01 * std::fabs(lzr))
00508 { return; }
00509
00510
00511
00512 QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
00513 mp = std::fabs(a-szr*b) + std::fabs(szi*b);
00514
00515
00516
00517 zm = std::sqrt(std::fabs(v));
00518 ee = 2.0*std::fabs(qp[0]);
00519 t = -szr*b;
00520 for (i=1;i<n;i++)
00521 {
00522 ee = ee*zm + std::fabs(qp[i]);
00523 }
00524 ee = ee*zm + std::fabs(a+t);
00525 ee *= (5.0 *mre + 4.0*are);
00526 ee = ee - (5.0*mre+2.0*are)*(std::fabs(a+t)+std::fabs(b)*zm)
00527 + 2.0*are*std::fabs(t);
00528
00529
00530
00531
00532 if (mp <= 20.0*ee)
00533 {
00534 *nz = 2;
00535 return;
00536 }
00537 j++;
00538
00539
00540
00541 if (j > 20) { return; }
00542 if (j >= 2)
00543 {
00544 if (!(relstp > 0.01 || mp < omp || tried))
00545 {
00546
00547
00548
00549 if (relstp < eta) { relstp = eta; }
00550 relstp = std::sqrt(relstp);
00551 u = u - u*relstp;
00552 v = v + v*relstp;
00553 QuadraticSyntheticDivision(n,&u,&v,p,qp,&a,&b);
00554 for (i=0;i<5;i++)
00555 {
00556 ComputeScalarFactors(&type);
00557 ComputeNextPolynomial(&type);
00558 }
00559 tried = 1;
00560 j = 0;
00561 }
00562 }
00563 omp = mp;
00564
00565
00566
00567 ComputeScalarFactors(&type);
00568 ComputeNextPolynomial(&type);
00569 ComputeScalarFactors(&type);
00570 ComputeNewEstimate(type,&ui,&vi);
00571
00572
00573
00574 if (!(vi != 0.0)) { return; }
00575 relstp = std::fabs((vi-v)/vi);
00576 u = ui;
00577 v = vi;
00578 }
00579 }
00580
00581 void G4JTPolynomialSolver::
00582 RealPolynomialIteration(G4double *sss, G4int *nz, G4int *iflag)
00583 {
00584
00585
00586
00587
00588
00589 G4double t=0.;
00590 G4double omp=0.;
00591 G4double pv=0.0, kv=0.0, xs= *sss;
00592 G4double mx=0.0, mp=0.0, ee=0.0;
00593 G4int i=1, j=0;
00594
00595 *nz = 0;
00596 *iflag = 0;
00597
00598
00599
00600 while (1)
00601 {
00602 pv = p[0];
00603
00604
00605
00606 qp[0] = pv;
00607 for (i=1;i<=n;i++)
00608 {
00609 pv = pv*xs + p[i];
00610 qp[i] = pv;
00611 }
00612 mp = std::fabs(pv);
00613
00614
00615
00616 mx = std::fabs(xs);
00617 ee = (mre/(are+mre))*std::fabs(qp[0]);
00618 for (i=1;i<=n;i++)
00619 {
00620 ee = ee*mx + std::fabs(qp[i]);
00621 }
00622
00623
00624
00625
00626 if (mp <= 20.0*((are+mre)*ee-mre*mp))
00627 {
00628 *nz = 1;
00629 szr = xs;
00630 szi = 0.0;
00631 return;
00632 }
00633 j++;
00634
00635
00636
00637 if (j > 10) { return; }
00638 if (j >= 2)
00639 {
00640 if (!(std::fabs(t) > 0.001*std::fabs(xs-t) || mp < omp))
00641 {
00642
00643
00644
00645 *iflag = 1;
00646 *sss = xs;
00647 return;
00648 }
00649 }
00650
00651 omp = mp;
00652
00653
00654
00655 kv = k[0];
00656 qk[0] = kv;
00657 for (i=1;i<n;i++)
00658 {
00659 kv = kv*xs + k[i];
00660 qk[i] = kv;
00661 }
00662 if (std::fabs(kv) <= std::fabs(k[n-1])*10.0*eta)
00663 {
00664 k[0] = 0.0;
00665 for (i=1;i<n;i++)
00666 {
00667 k[i] = qk[i-1];
00668 }
00669 }
00670 else
00671 {
00672 t = -pv/kv;
00673 k[0] = qp[0];
00674 for (i=1;i<n;i++)
00675 {
00676 k[i] = t*qk[i-1] + qp[i];
00677 }
00678 }
00679 kv = k[0];
00680 for (i=1;i<n;i++)
00681 {
00682 kv = kv*xs + k[i];
00683 }
00684 t = 0.0;
00685 if (std::fabs(kv) > std::fabs(k[n-1]*10.0*eta)) { t = -pv/kv; }
00686 xs += t;
00687 }
00688 }
00689
00690 void G4JTPolynomialSolver::ComputeScalarFactors(G4int *type)
00691 {
00692
00693
00694
00695
00696
00697
00698
00699
00700 QuadraticSyntheticDivision(n-1,&u,&v,k,qk,&c,&d);
00701 if (std::fabs(c) <= std::fabs(k[n-1]*100.0*eta))
00702 {
00703 if (std::fabs(d) <= std::fabs(k[n-2]*100.0*eta))
00704 {
00705 *type = 3;
00706 return;
00707 }
00708 }
00709
00710 if (std::fabs(d) < std::fabs(c))
00711 {
00712 *type = 1;
00713 e = a/c;
00714 f = d/c;
00715 g = u*e;
00716 h = v*b;
00717 a3 = a*e + (h/c+g)*b;
00718 a1 = b - a*(d/c);
00719 a7 = a + g*d + h*f;
00720 return;
00721 }
00722 *type = 2;
00723 e = a/d;
00724 f = c/d;
00725 g = u*b;
00726 h = v*b;
00727 a3 = (a+g)*e + h*(b/d);
00728 a1 = b*f-a;
00729 a7 = (f+u)*a + h;
00730 }
00731
00732 void G4JTPolynomialSolver::ComputeNextPolynomial(G4int *type)
00733 {
00734
00735
00736
00737 G4int i=2;
00738
00739 if (*type == 3)
00740 {
00741 k[0] = 0.0;
00742 k[1] = 0.0;
00743 for (i=2;i<n;i++)
00744 {
00745 k[i] = qk[i-2];
00746 }
00747 return;
00748 }
00749 G4double temp = a;
00750 if (*type == 1) { temp = b; }
00751 if (std::fabs(a1) <= std::fabs(temp)*eta*10.0)
00752 {
00753
00754
00755 k[0] = 0.0;
00756 k[1] = -a7*qp[0];
00757 for(i=2;i<n;i++)
00758 {
00759 k[i] = a3*qk[i-2] - a7*qp[i-1];
00760 }
00761 return;
00762 }
00763
00764
00765
00766 a7 /= a1;
00767 a3 /= a1;
00768 k[0] = qp[0];
00769 k[1] = qp[1] - a7*qp[0];
00770 for (i=2;i<n;i++)
00771 {
00772 k[i] = a3*qk[i-2] - a7*qp[i-1] + qp[i];
00773 }
00774 }
00775
00776 void G4JTPolynomialSolver::
00777 ComputeNewEstimate(G4int type, G4double *uu, G4double *vv)
00778 {
00779
00780
00781
00782 G4double a4=0.0, a5=0.0, b1=0.0, b2=0.0,
00783 c1=0.0, c2=0.0, c3=0.0, c4=0.0, temp=0.0;
00784
00785
00786
00787 if (type == 3)
00788 {
00789 *uu = 0.0;
00790 *vv = 0.0;
00791 return;
00792 }
00793 if (type == 2)
00794 {
00795 a4 = (a+g)*f + h;
00796 a5 = (f+u)*c + v*d;
00797 }
00798 else
00799 {
00800 a4 = a + u*b +h*f;
00801 a5 = c + (u+v*f)*d;
00802 }
00803
00804
00805
00806 b1 = -k[n-1]/p[n];
00807 b2 = -(k[n-2]+b1*p[n-1])/p[n];
00808 c1 = v*b2*a1;
00809 c2 = b1*a7;
00810 c3 = b1*b1*a3;
00811 c4 = c1 - c2 - c3;
00812 temp = a5 + b1*a4 - c4;
00813 if (!(temp != 0.0))
00814 {
00815 *uu = 0.0;
00816 *vv = 0.0;
00817 return;
00818 }
00819 *uu = u - (u*(c3+c2)+v*(b1*a1+b2*a7))/temp;
00820 *vv = v*(1.0+c4/temp);
00821 return;
00822 }
00823
00824 void G4JTPolynomialSolver::
00825 QuadraticSyntheticDivision(G4int nn, G4double *uu, G4double *vv,
00826 std::vector<G4double> &pp, std::vector<G4double> &qq,
00827 G4double *aa, G4double *bb)
00828 {
00829
00830
00831
00832 G4double cc=0.0;
00833 *bb = pp[0];
00834 qq[0] = *bb;
00835 *aa = pp[1] - (*bb)*(*uu);
00836 qq[1] = *aa;
00837 for (G4int i=2;i<=nn;i++)
00838 {
00839 cc = pp[i] - (*aa)*(*uu) - (*bb)*(*vv);
00840 qq[i] = cc;
00841 *bb = *aa;
00842 *aa = cc;
00843 }
00844 }
00845
00846 void G4JTPolynomialSolver::Quadratic(G4double aa,G4double b1,
00847 G4double cc,G4double *ssr,G4double *ssi,
00848 G4double *lr,G4double *li)
00849 {
00850
00851
00852
00853
00854
00855
00856
00857 G4double bb=0.0, dd=0.0, ee=0.0;
00858
00859 if (!(aa != 0.0))
00860 {
00861 if (b1 != 0.0)
00862 { *ssr = -cc/b1; }
00863 else
00864 { *ssr = 0.0; }
00865 *lr = 0.0;
00866 *ssi = 0.0;
00867 *li = 0.0;
00868 return;
00869 }
00870 if (!(cc != 0.0))
00871 {
00872 *ssr = 0.0;
00873 *lr = -b1/aa;
00874 *ssi = 0.0;
00875 *li = 0.0;
00876 return;
00877 }
00878
00879
00880
00881 bb = b1/2.0;
00882 if (std::fabs(bb) < std::fabs(cc))
00883 {
00884 if (cc < 0.0)
00885 { ee = -aa; }
00886 else
00887 { ee = aa; }
00888 ee = bb*(bb/std::fabs(cc)) - ee;
00889 dd = std::sqrt(std::fabs(ee))*std::sqrt(std::fabs(cc));
00890 }
00891 else
00892 {
00893 ee = 1.0 - (aa/bb)*(cc/bb);
00894 dd = std::sqrt(std::fabs(ee))*std::fabs(bb);
00895 }
00896 if (ee < 0.0)
00897 {
00898 *ssr = -bb/aa;
00899 *lr = *ssr;
00900 *ssi = std::fabs(dd/aa);
00901 *li = -(*ssi);
00902 }
00903 else
00904 {
00905 if (bb >= 0.0)
00906 { dd = -dd; }
00907 *lr = (-bb+dd)/aa;
00908 *ssr = 0.0;
00909 if (*lr != 0.0)
00910 { *ssr = (cc/ *lr)/aa; }
00911 *ssi = 0.0;
00912 *li = 0.0;
00913 }
00914 }