Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4JTPolynomialSolver class implementation. 27 // 28 // Author: Oliver Link, 15.02.2005 29 // Translated to C++ and adapted to us 30 // ------------------------------------------- 31 32 #include "G4JTPolynomialSolver.hh" 33 #include "G4Pow.hh" 34 #include "G4SystemOfUnits.hh" 35 36 const G4double G4JTPolynomialSolver::base = 37 const G4double G4JTPolynomialSolver::eta = 38 const G4double G4JTPolynomialSolver::infin = 39 const G4double G4JTPolynomialSolver::smalno = 40 const G4double G4JTPolynomialSolver::are = 41 const G4double G4JTPolynomialSolver::mre = 42 const G4double G4JTPolynomialSolver::lo = 43 44 G4int G4JTPolynomialSolver::FindRoots(G4double 45 G4double 46 { 47 G4double t = 0.0, aa = 0.0, bb = 0.0, cc = 0 48 G4double max = 0.0, min = infin, xxx = 0.0, 49 G4double xm = 0.0, ff = 0.0, df = 0.0, dx = 50 G4int cnt = 0, nz = 0, i = 0, j = 0, jj = 0, 51 G4Pow* power = G4Pow::GetInstance(); 52 53 // Initialization of constants for shift rot 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), 58 G4double xo = xx, yo = -xx; 59 n = degr; 60 61 // Algorithm fails if the leading coefficie 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 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 110 { 111 Quadratic(p[0], p[1], p[2], &zeror[degr 112 &zeror[degr - 1], &zeroi[degr 113 n -= 2; 114 return degr - n; 115 } 116 117 // Find largest and smallest moduli of coe 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 135 // Computes a scale factor to multiply the 136 // polynomial. The scaling is done to avoi 137 // avoid undetected underflow interfering 138 // criterion. The factor is a power of the 139 // 140 sc = lo / min; 141 142 if(((sc <= 1.0) && (max >= 10.0)) || ((sc 143 ((infin / sc >= max) && (max >= 10))) 144 { 145 if(!(sc != 0.0)) 146 { 147 sc = smalno; 148 } 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 } 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])) / 171 172 // If Newton step at the origin is better, 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 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 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] / (G4dou 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 234 { 235 // Use a scaled form of recurrence if 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 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] 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 correspond 267 // 268 for(cnt = 0; cnt < 20; ++cnt) 269 { 270 // Quadratic corresponds to a double shi 271 // non-real point and its complex conjug 272 // has modulus bnd and amplitude rotated 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 + 283 if(nz != 0) 284 { 285 // The second stage jumps directly to 286 // stage iterations and returns here i 287 // Deflate the polynomial, store the z 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 anot 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 wit 318 // 319 return degr - n; 320 } 321 322 void G4JTPolynomialSolver::ComputeFixedShiftPo 323 { 324 // Computes up to L2 fixed shift k-polynomia 325 // in the linear or quadratic case. Initiate 326 // iterations and returns with the number of 327 328 G4double svu = 0.0, svv = 0.0, ui = 0.0, vi 329 G4double betas = 0.25, betav = 0.25, oss = s 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, vpa 334 stry = 0; 335 336 *nz = 0; 337 338 // Evaluate polynomial by synthetic division 339 // 340 QuadraticSyntheticDivision(n, &u, &v, p, qp, 341 ComputeScalarFactors(&type); 342 for(j = 0; j < l2; ++j) 343 { 344 // Calculate next k polynomial and estimat 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 convergenc 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 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 co 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 faste 416 // 417 vtry = 0; 418 stry = 0; 419 if(((spass != 0) && (vpass == 0)) || (tss 420 { 421 RealPolynomialIteration(&xs, nz, &iflag) 422 if(*nz > 0) 423 { 424 return; 425 } 426 427 // Linear iteration has failed. Flag tha 428 // tried and decrease the convergence cr 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 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, n 449 if(*nz > 0) 450 { 451 return; 452 } 453 454 // Quadratic iteration has failed. Flag 455 // been tried and decrease the convergen 456 // 457 vtry = 1; 458 betav *= 0.25; 459 460 // Try linear iteration if it has not be 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 tha 478 // tried and decrease the convergence cr 479 // 480 stry = 1; 481 betas *= 0.25; 482 if(iflag == 0) 483 { 484 break; 485 } 486 487 // If linear iteration signals an almost 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 b 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 conti 514 // second stage. 515 // 516 QuadraticSyntheticDivision(n, &u, &v, p, q 517 ComputeScalarFactors(&type); 518 519 ovv = vv; 520 oss = ss; 521 otv = tv; 522 ots = ts; 523 } 524 } 525 526 void G4JTPolynomialSolver::QuadraticPolynomial 527 528 { 529 // Variable-shift k-polynomial iteration for 530 // quadratic factor converges only if the ze 531 // equimodular or nearly so. 532 // uu, vv - coefficients of starting quadrat 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 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, &lz 551 552 // Return if roots of the quadratic are re 553 // close to multiple or nearly equal and o 554 // sign. 555 // 556 if(std::fabs(std::fabs(szr) - std::fabs(lz 557 { 558 return; 559 } 560 561 // Evaluate polynomial by quadratic synthe 562 // 563 QuadraticSyntheticDivision(n, &u, &v, p, q 564 mp = std::fabs(a - szr * b) + std::fabs(sz 565 566 // Compute a rigorous bound on the roundin 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:: 578 2.0 * are * std::fabs(t); 579 580 // Iteration has converged sufficiently if 581 // polynomial value is less than 20 times 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 599 { 600 // A cluster appears to be stalling th 601 // Five fixed shift steps are taken wi 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, 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 a 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 conv 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::RealPolynomialItera 642 643 { 644 // Variable-shift H polynomial iteration for 645 // sss - starting iterate 646 // nz - number of zeros found 647 // iflag - flag to indicate a pair of zeros 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 i 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 684 // value is less than 20 times this bound. 685 // 686 if(mp <= 20.0 * ((are + mre) * ee - mre * 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 704 { 705 // A cluster of zeros near the real ax 706 // Return with iflag set to initiate a 707 // 708 *iflag = 1; 709 *sss = xs; 710 return; 711 } // Return if the polynomial value has 712 } 713 714 omp = mp; 715 716 // Compute t, the next polynomial, and th 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]) * 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 recurr 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 749 { 750 t = -pv / kv; 751 } 752 xs += t; 753 } 754 } 755 756 void G4JTPolynomialSolver::ComputeScalarFactor 757 { 758 // This function calculates scalar quantitie 759 // compute the next k polynomial and new est 760 // the quadratic coefficients. 761 // type - integer variable set here indicati 762 // calculations are normalized to avoid over 763 764 // Synthetic division of k by the quadratic 765 // 766 QuadraticSyntheticDivision(n - 1, &u, &v, k, 767 if(std::fabs(c) <= std::fabs(k[n - 1] * 100. 768 { 769 if(std::fabs(d) <= std::fabs(k[n - 2] * 10 770 { 771 *type = 3; // Type=3 indicates the quad 772 return; 773 } 774 } 775 776 if(std::fabs(d) < std::fabs(c)) 777 { 778 *type = 1; // Type=1 indicates that all f 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 for 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::ComputeNextPolynomi 799 { 800 // Computes the next k polynomials using sca 801 // computed in ComputeScalarFactors. 802 803 G4int i = 2; 804 805 if(*type == 3) // Use unscaled form of the 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 * 821 { 822 // If a1 is nearly zero then use a special 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] + q 842 } 843 } 844 845 void G4JTPolynomialSolver::ComputeNewEstimate( 846 847 { 848 // Compute new estimates of the quadratic co 849 // using the scalars computed in calcsc. 850 851 G4double a4 = 0.0, a5 = 0.0, b1 = 0.0, b2 = 852 c4 = 0.0, temp = 0.0; 853 854 // Use formulas appropriate to setting of ty 855 // 856 if(type == 3) // If type=3 the quadratic i 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 889 *vv = v * (1.0 + c4 / temp); 890 return; 891 } 892 893 void G4JTPolynomialSolver::QuadraticSyntheticD 894 G4int nn, G4double* uu, G4double* vv, std::v 895 std::vector<G4double>& qq, G4double* aa, G4d 896 { 897 // Divides pp by the quadratic 1,uu,vv placi 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) * (* 908 qq[i] = cc; 909 *bb = *aa; 910 *aa = cc; 911 } 912 } 913 914 void G4JTPolynomialSolver::Quadratic(G4double 915 G4double* 916 G4double* 917 { 918 // Calculate the zeros of the quadratic aa*z 919 // The quadratic formula, modified to avoid 920 // to find the larger zero if the zeros are 921 // are complex. The smaller real zero is fou 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 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( 965 } 966 else 967 { 968 ee = 1.0 - (aa / bb) * (cc / bb); 969 dd = std::sqrt(std::fabs(ee)) * std::fabs( 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