From 377ee77486a2631d1d3f006d81ee0bebf34eaab4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 20 Dec 2021 09:28:21 -0500 Subject: [PATCH] Use quadric solver from ACM TOMS 46(2), pp1-28 (algorithm 1010) --- CMakeLists.txt | 4 + include/openmc/external/quartic_solver.h | 10 + src/external/LICENSE | 32 ++ src/external/quartic_solver.c | 621 +++++++++++++++++++++++ src/surface.cpp | 308 +++-------- 5 files changed, 732 insertions(+), 243 deletions(-) create mode 100644 include/openmc/external/quartic_solver.h create mode 100644 src/external/LICENSE create mode 100644 src/external/quartic_solver.c diff --git a/CMakeLists.txt b/CMakeLists.txt index d0c2569c34..9cfb55b7a3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -364,6 +364,10 @@ list(APPEND libopenmc_SOURCES src/xml_interface.cpp src/xsdata.cpp) +# Add bundled external dependencies +list(APPEND libopenmc_SOURCES + src/external/quartic_solver.c) + # For Visual Studio compilers if(MSVC) # Use static library (otherwise explicit symbol portings are needed) diff --git a/include/openmc/external/quartic_solver.h b/include/openmc/external/quartic_solver.h new file mode 100644 index 0000000000..98159afb07 --- /dev/null +++ b/include/openmc/external/quartic_solver.h @@ -0,0 +1,10 @@ +#ifndef OPENMC_EXTERNAL_QUARTIC_SOLVER_H +#define OPENMC_EXTERNAL_QUARTIC_SOLVER_H + +#include + +extern "C" { +void oqs_quartic_solver(double coeff[5], std::complex roots[4]); +} + +#endif // OPENMC_EXTERNAL_QUARTIC_SOLVER_H diff --git a/src/external/LICENSE b/src/external/LICENSE new file mode 100644 index 0000000000..ac09c7fcb7 --- /dev/null +++ b/src/external/LICENSE @@ -0,0 +1,32 @@ +The quartic solver was obtained from the paper: Alberto Giacomo Orellana and +Cristiano De Michele, "Algorithm 1010: Boosting Efficiency in Solving Quartic +Equations with No Compromise in Accuracy," ACM Transactions on Mathematical +Software, 46 (2), pp. 1-28. https://doi.org/10.1145/3386241 + +OpenMC developers contacted the authors, who have agreed to license their +software under the simplified BSD license, reproduced below: + +------------------------------------------------------------------------------- +Copyright (c) 2020 Alberto Giacomo Orellana and Cristiano De Michele +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE PYNE DEVELOPMENT TEAM ``AS IS'' AND ANY EXPRESS +OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/src/external/quartic_solver.c b/src/external/quartic_solver.c new file mode 100644 index 0000000000..71b24b8a4a --- /dev/null +++ b/src/external/quartic_solver.c @@ -0,0 +1,621 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define Sqr(x) ((x) * (x)) +#ifndef CMPLX +#define CMPLX(x, y) (x) + (y)*I +#endif +const double cubic_rescal_fact = + 3.488062113727083E+102; //= pow(DBL_MAX,1.0/3.0)/1.618034; +const double quart_rescal_fact = + 7.156344627944542E+76; // = pow(DBL_MAX,1.0/4.0)/1.618034; +const double macheps = 2.2204460492503131E-16; // DBL_EPSILON +double oqs_max2(double a, double b) +{ + if (a >= b) + return a; + else + return b; +} +double oqs_max3(double a, double b, double c) +{ + double t; + t = oqs_max2(a, b); + return oqs_max2(t, c); +} + +void oqs_solve_cubic_analytic_depressed_handle_inf( + double b, double c, double* sol) +{ + /* find analytically the dominant root of a depressed cubic x^3+b*x+c + * where coefficients b and c are large (see sec. 2.2 in the manuscript) */ + double Q, R, theta, A, B, QR, QRSQ, KK, sqrtQ, RQ; + ; + const double PI2 = M_PI / 2.0, TWOPI = 2.0 * M_PI; + Q = -b / 3.0; + R = 0.5 * c; + if (R == 0) { + if (b <= 0) { + *sol = sqrt(-b); + } else { + *sol = 0; + } + return; + } + + if (fabs(Q) < fabs(R)) { + QR = Q / R; + QRSQ = QR * QR; + KK = 1.0 - Q * QRSQ; + } else { + RQ = R / Q; + KK = copysign(1.0, Q) * (RQ * RQ / Q - 1.0); + } + + if (KK < 0.0) { + sqrtQ = sqrt(Q); + theta = acos((R / fabs(Q)) / sqrtQ); + if (theta < PI2) + *sol = -2.0 * sqrtQ * cos(theta / 3.0); + else + *sol = -2.0 * sqrtQ * cos((theta + TWOPI) / 3.0); + } else { + if (fabs(Q) < fabs(R)) + A = -copysign(1.0, R) * cbrt(fabs(R) * (1.0 + sqrt(KK))); + else { + A = + -copysign(1.0, R) * cbrt(fabs(R) + sqrt(fabs(Q)) * fabs(Q) * sqrt(KK)); + } + if (A == 0.0) + B = 0.0; + else + B = Q / A; + *sol = A + B; + } +} +void oqs_solve_cubic_analytic_depressed(double b, double c, double* sol) +{ + /* find analytically the dominant root of a depressed cubic x^3+b*x+c + * (see sec. 2.2 in the manuscript) */ + double Q, R, theta, Q3, R2, A, B, sqrtQ; + Q = -b / 3.0; + R = 0.5 * c; + if (fabs(Q) > 1E102 || fabs(R) > 1E154) { + oqs_solve_cubic_analytic_depressed_handle_inf(b, c, sol); + return; + } + Q3 = Sqr(Q) * Q; + R2 = Sqr(R); + if (R2 < Q3) { + theta = acos(R / sqrt(Q3)); + sqrtQ = -2.0 * sqrt(Q); + if (theta < M_PI / 2) + *sol = sqrtQ * cos(theta / 3.0); + else + *sol = sqrtQ * cos((theta + 2.0 * M_PI) / 3.0); + } else { + A = -copysign(1.0, R) * pow(fabs(R) + sqrt(R2 - Q3), 1.0 / 3.0); + if (A == 0.0) + B = 0.0; + else + B = Q / A; + *sol = A + B; /* this is always largest root even if A=B */ + } +} +void oqs_calc_phi0( + double a, double b, double c, double d, double* phi0, int scaled) +{ + /* find phi0 as the dominant root of the depressed and shifted cubic + * in eq. (79) (see also the discussion in sec. 2.2 of the manuscript) */ + double rmax, g, h, gg, hh, aq, bq, cq, dq, s, diskr; + double maxtt, xxx, gx, x, xold, f, fold, df, xsq; + double ggss, hhss, dqss, aqs, bqs, cqs, rfact, rfactsq; + int iter; + diskr = 9 * a * a - 24 * b; + /* eq. (87) */ + if (diskr > 0.0) { + diskr = sqrt(diskr); + if (a > 0.0) + s = -2 * b / (3 * a + diskr); + else + s = -2 * b / (3 * a - diskr); + } else { + s = -a / 4; + } + /* eqs. (83) */ + aq = a + 4 * s; + bq = b + 3 * s * (a + 2 * s); + cq = c + s * (2 * b + s * (3 * a + 4 * s)); + dq = d + s * (c + s * (b + s * (a + s))); + gg = bq * bq / 9; + hh = aq * cq; + + g = hh - 4 * dq - 3 * gg; /* eq. (85) */ + h = (8 * dq + hh - 2 * gg) * bq / 3 - cq * cq - dq * aq * aq; /* eq. (86) */ + oqs_solve_cubic_analytic_depressed(g, h, &rmax); + if (isnan(rmax) || isinf(rmax)) { + oqs_solve_cubic_analytic_depressed_handle_inf(g, h, &rmax); + if ((isnan(rmax) || isinf(rmax)) && scaled) { + // try harder: rescale also the depressed cubic if quartic has been + // already rescaled + rfact = cubic_rescal_fact; + rfactsq = rfact * rfact; + ggss = gg / rfactsq; + hhss = hh / rfactsq; + dqss = dq / rfactsq; + aqs = aq / rfact; + bqs = bq / rfact; + cqs = cq / rfact; + ggss = bqs * bqs / 9.0; + hhss = aqs * cqs; + g = hhss - 4.0 * dqss - 3.0 * ggss; + h = (8.0 * dqss + hhss - 2.0 * ggss) * bqs / 3 - cqs * (cqs / rfact) - + (dq / rfact) * aqs * aqs; + oqs_solve_cubic_analytic_depressed(g, h, &rmax); + if (isnan(rmax) || isinf(rmax)) { + oqs_solve_cubic_analytic_depressed_handle_inf(g, h, &rmax); + } + rmax *= rfact; + } + } + /* Newton-Raphson used to refine phi0 (see end of sec. 2.2 in the manuscript) + */ + x = rmax; + xsq = x * x; + xxx = x * xsq; + gx = g * x; + f = x * (xsq + g) + h; + if (fabs(xxx) > fabs(gx)) + maxtt = fabs(xxx); + else + maxtt = fabs(gx); + if (fabs(h) > maxtt) + maxtt = fabs(h); + + if (fabs(f) > macheps * maxtt) { + for (iter = 0; iter < 8; iter++) { + df = 3.0 * xsq + g; + if (df == 0) { + break; + } + xold = x; + x += -f / df; + fold = f; + xsq = x * x; + f = x * (xsq + g) + h; + if (f == 0) { + break; + } + + if (fabs(f) >= fabs(fold)) { + x = xold; + break; + } + } + } + *phi0 = x; +} +double oqs_calc_err_ldlt( + double b, double c, double d, double d2, double l1, double l2, double l3) +{ + /* Eqs. (29) and (30) in the manuscript */ + double sum; + sum = (b == 0) ? fabs(d2 + l1 * l1 + 2.0 * l3) + : fabs(((d2 + l1 * l1 + 2.0 * l3) - b) / b); + sum += (c == 0) ? fabs(2.0 * d2 * l2 + 2.0 * l1 * l3) + : fabs(((2.0 * d2 * l2 + 2.0 * l1 * l3) - c) / c); + sum += (d == 0) ? fabs(d2 * l2 * l2 + l3 * l3) + : fabs(((d2 * l2 * l2 + l3 * l3) - d) / d); + return sum; +} +double oqs_calc_err_abcd_cmplx(double a, double b, double c, double d, + complex double aq, complex double bq, complex double cq, complex double dq) +{ + /* Eqs. (68) and (69) in the manuscript for complex alpha1 (aq), beta1 (bq), + * alpha2 (cq) and beta2 (dq) */ + double sum; + sum = (d == 0) ? cabs(bq * dq) : cabs((bq * dq - d) / d); + sum += + (c == 0) ? cabs(bq * cq + aq * dq) : cabs(((bq * cq + aq * dq) - c) / c); + sum += + (b == 0) ? cabs(bq + aq * cq + dq) : cabs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? cabs(aq + cq) : cabs(((aq + cq) - a) / a); + return sum; +} +double oqs_calc_err_abcd(double a, double b, double c, double d, double aq, + double bq, double cq, double dq) +{ + /* Eqs. (68) and (69) in the manuscript for real alpha1 (aq), beta1 (bq), + * alpha2 (cq) and beta2 (dq)*/ + double sum; + sum = (d == 0) ? fabs(bq * dq) : fabs((bq * dq - d) / d); + sum += + (c == 0) ? fabs(bq * cq + aq * dq) : fabs(((bq * cq + aq * dq) - c) / c); + sum += + (b == 0) ? fabs(bq + aq * cq + dq) : fabs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? fabs(aq + cq) : fabs(((aq + cq) - a) / a); + return sum; +} +double oqs_calc_err_abc( + double a, double b, double c, double aq, double bq, double cq, double dq) +{ + /* Eqs. (48)-(51) in the manuscript */ + double sum; + sum = + (c == 0) ? fabs(bq * cq + aq * dq) : fabs(((bq * cq + aq * dq) - c) / c); + sum += + (b == 0) ? fabs(bq + aq * cq + dq) : fabs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? fabs(aq + cq) : fabs(((aq + cq) - a) / a); + return sum; +} +void oqs_NRabcd(double a, double b, double c, double d, double* AQ, double* BQ, + double* CQ, double* DQ) +{ + /* Newton-Raphson described in sec. 2.3 of the manuscript for complex + * coefficients a,b,c,d */ + int iter, k1, k2; + double x02, errf, errfold, xold[4], x[4], dx[4], det, Jinv[4][4], fvec[4], + vr[4]; + x[0] = *AQ; + x[1] = *BQ; + x[2] = *CQ; + x[3] = *DQ; + vr[0] = d; + vr[1] = c; + vr[2] = b; + vr[3] = a; + fvec[0] = x[1] * x[3] - d; + fvec[1] = x[1] * x[2] + x[0] * x[3] - c; + fvec[2] = x[1] + x[0] * x[2] + x[3] - b; + fvec[3] = x[0] + x[2] - a; + errf = 0; + for (k1 = 0; k1 < 4; k1++) { + errf += (vr[k1] == 0) ? fabs(fvec[k1]) : fabs(fvec[k1] / vr[k1]); + } + for (iter = 0; iter < 8; iter++) { + x02 = x[0] - x[2]; + det = x[1] * x[1] + x[1] * (-x[2] * x02 - 2.0 * x[3]) + + x[3] * (x[0] * x02 + x[3]); + if (det == 0.0) + break; + Jinv[0][0] = x02; + Jinv[0][1] = x[3] - x[1]; + Jinv[0][2] = x[1] * x[2] - x[0] * x[3]; + Jinv[0][3] = -x[1] * Jinv[0][1] - x[0] * Jinv[0][2]; + Jinv[1][0] = x[0] * Jinv[0][0] + Jinv[0][1]; + Jinv[1][1] = -x[1] * Jinv[0][0]; + Jinv[1][2] = -x[1] * Jinv[0][1]; + Jinv[1][3] = -x[1] * Jinv[0][2]; + Jinv[2][0] = -Jinv[0][0]; + Jinv[2][1] = -Jinv[0][1]; + Jinv[2][2] = -Jinv[0][2]; + Jinv[2][3] = Jinv[0][2] * x[2] + Jinv[0][1] * x[3]; + Jinv[3][0] = -x[2] * Jinv[0][0] - Jinv[0][1]; + Jinv[3][1] = Jinv[0][0] * x[3]; + Jinv[3][2] = x[3] * Jinv[0][1]; + Jinv[3][3] = x[3] * Jinv[0][2]; + for (k1 = 0; k1 < 4; k1++) { + dx[k1] = 0; + for (k2 = 0; k2 < 4; k2++) + dx[k1] += Jinv[k1][k2] * fvec[k2]; + } + for (k1 = 0; k1 < 4; k1++) + xold[k1] = x[k1]; + + for (k1 = 0; k1 < 4; k1++) { + x[k1] += -dx[k1] / det; + } + fvec[0] = x[1] * x[3] - d; + fvec[1] = x[1] * x[2] + x[0] * x[3] - c; + fvec[2] = x[1] + x[0] * x[2] + x[3] - b; + fvec[3] = x[0] + x[2] - a; + errfold = errf; + errf = 0; + for (k1 = 0; k1 < 4; k1++) { + errf += (vr[k1] == 0) ? fabs(fvec[k1]) : fabs(fvec[k1] / vr[k1]); + } + if (errf == 0) + break; + if (errf >= errfold) { + for (k1 = 0; k1 < 4; k1++) + x[k1] = xold[k1]; + break; + } + } + *AQ = x[0]; + *BQ = x[1]; + *CQ = x[2]; + *DQ = x[3]; +} +void oqs_solve_quadratic(double a, double b, complex double roots[2]) +{ + double div, sqrtd, diskr, zmax, zmin; + diskr = a * a - 4 * b; + if (diskr >= 0.0) { + if (a >= 0.0) + div = -a - sqrt(diskr); + else + div = -a + sqrt(diskr); + + zmax = div / 2; + + if (zmax == 0.0) + zmin = 0.0; + else + zmin = b / zmax; + + roots[0] = CMPLX(zmax, 0.0); + roots[1] = CMPLX(zmin, 0.0); + } else { + sqrtd = sqrt(-diskr); + roots[0] = CMPLX(-a / 2, sqrtd / 2); + roots[1] = CMPLX(-a / 2, -sqrtd / 2); + } +} +void oqs_quartic_solver(double coeff[5], complex double roots[4]) +{ + /* USAGE: + * + * This routine calculates the roots of the quartic equation + * + * coeff[4]*x^4 + coeff[3]*x^3 + coeff[2]*x^2 + coeff[1]*x + coeff[0] = 0 + * + * if coeff[4] != 0 + * + * the four roots will be stored in the complex array roots[] + * + * */ + complex double acx1, bcx1, ccx1, dcx1, acx, bcx, ccx, dcx, cdiskr, zx1, zx2, + zxmax, zxmin, qroots[2]; + double l2m[12], d2m[12], res[12], resmin, bl311, dml3l3, err0 = 0, err1 = 0, + aq1, bq1, cq1, dq1; + double a, b, c, d, phi0, aq, bq, cq, dq, d2, d3, l1, l2, l3, errmin, errv[3], + aqv[3], cqv[3], gamma, del2; + int realcase[2], whichcase, k1, k, kmin, nsol; + double rfactsq, rfact = 1.0; + + if (coeff[4] == 0.0) { + printf("That's not a quartic!\n"); + return; + } + a = coeff[3] / coeff[4]; + b = coeff[2] / coeff[4]; + c = coeff[1] / coeff[4]; + d = coeff[0] / coeff[4]; + oqs_calc_phi0(a, b, c, d, &phi0, 0); + + // simple polynomial rescaling + if (isnan(phi0) || isinf(phi0)) { + rfact = quart_rescal_fact; + a /= rfact; + rfactsq = rfact * rfact; + b /= rfactsq; + c /= rfactsq * rfact; + d /= rfactsq * rfactsq; + oqs_calc_phi0(a, b, c, d, &phi0, 1); + } + l1 = a / 2; /* eq. (16) */ + l3 = b / 6 + phi0 / 2; /* eq. (18) */ + del2 = c - a * l3; /* defined just after eq. (27) */ + nsol = 0; + bl311 = 2. * b / 3. - phi0 - l1 * l1; /* This is d2 as defined in eq. (20)*/ + dml3l3 = d - l3 * l3; /* dml3l3 is d3 as defined in eq. (15) with d2=0 */ + + /* Three possible solutions for d2 and l2 (see eqs. (28) and discussion which + * follows) */ + if (bl311 != 0.0) { + d2m[nsol] = bl311; + l2m[nsol] = del2 / (2.0 * d2m[nsol]); + res[nsol] = oqs_calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + if (del2 != 0) { + l2m[nsol] = 2 * dml3l3 / del2; + if (l2m[nsol] != 0) { + d2m[nsol] = del2 / (2 * l2m[nsol]); + res[nsol] = oqs_calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + + d2m[nsol] = bl311; + l2m[nsol] = 2.0 * dml3l3 / del2; + res[nsol] = oqs_calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + + if (nsol == 0) { + l2 = d2 = 0.0; + } else { + /* we select the (d2,l2) pair which minimizes errors */ + for (k1 = 0; k1 < nsol; k1++) { + if (k1 == 0 || res[k1] < resmin) { + resmin = res[k1]; + kmin = k1; + } + } + d2 = d2m[kmin]; + l2 = l2m[kmin]; + } + whichcase = 0; + if (d2 < 0.0) { + /* Case I eqs. (37)-(40) */ + gamma = sqrt(-d2); + aq = l1 + gamma; + bq = l3 + gamma * l2; + + cq = l1 - gamma; + dq = l3 - gamma * l2; + if (fabs(dq) < fabs(bq)) + dq = d / bq; + else if (fabs(dq) > fabs(bq)) + bq = d / dq; + if (fabs(aq) < fabs(cq)) { + nsol = 0; + if (dq != 0) { + aqv[nsol] = (c - bq * cq) / dq; /* see eqs. (47) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + } + if (cq != 0) { + aqv[nsol] = (b - dq - bq) / cq; /* see eqs. (47) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + } + aqv[nsol] = a - cq; /* see eqs. (47) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + /* we select the value of aq (i.e. alpha1 in the manuscript) which + * minimizes errors */ + for (k = 0; k < nsol; k++) { + if (k == 0 || errv[k] < errmin) { + kmin = k; + errmin = errv[k]; + } + } + aq = aqv[kmin]; + } else { + nsol = 0; + if (bq != 0) { + cqv[nsol] = (c - aq * dq) / bq; /* see eqs. (53) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + } + if (aq != 0) { + cqv[nsol] = (b - bq - dq) / aq; /* see eqs. (53) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + } + cqv[nsol] = a - aq; /* see eqs. (53) */ + errv[nsol] = oqs_calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + /* we select the value of cq (i.e. alpha2 in the manuscript) which + * minimizes errors */ + for (k = 0; k < nsol; k++) { + if (k == 0 || errv[k] < errmin) { + kmin = k; + errmin = errv[k]; + } + } + cq = cqv[kmin]; + } + realcase[0] = 1; + } else if (d2 > 0) { + /* Case II eqs. (53)-(56) */ + gamma = sqrt(d2); + acx = CMPLX(l1, gamma); + bcx = CMPLX(l3, gamma * l2); + ccx = conj(acx); + dcx = conj(bcx); + realcase[0] = 0; + } else + realcase[0] = -1; // d2=0 + /* Case III: d2 is 0 or approximately 0 (in this case check which solution is + * better) */ + if (realcase[0] == -1 || (fabs(d2) <= macheps * oqs_max3(fabs(2. * b / 3.), + fabs(phi0), l1 * l1))) { + d3 = d - l3 * l3; + if (realcase[0] == 1) + err0 = oqs_calc_err_abcd(a, b, c, d, aq, bq, cq, dq); + else if (realcase[0] == 0) + err0 = oqs_calc_err_abcd_cmplx(a, b, c, d, acx, bcx, ccx, dcx); + if (d3 <= 0) { + realcase[1] = 1; + aq1 = l1; + bq1 = l3 + sqrt(-d3); + cq1 = l1; + dq1 = l3 - sqrt(-d3); + if (fabs(dq1) < fabs(bq1)) + dq1 = d / bq1; + else if (fabs(dq1) > fabs(bq1)) + bq1 = d / dq1; + err1 = oqs_calc_err_abcd(a, b, c, d, aq1, bq1, cq1, dq1); /* eq. (68) */ + } else /* complex */ + { + realcase[1] = 0; + acx1 = l1; + bcx1 = l3 + I * sqrt(d3); + ccx1 = l1; + dcx1 = conj(bcx1); + err1 = oqs_calc_err_abcd_cmplx(a, b, c, d, acx1, bcx1, ccx1, dcx1); + } + if (realcase[0] == -1 || err1 < err0) { + whichcase = 1; // d2 = 0 + if (realcase[1] == 1) { + aq = aq1; + bq = bq1; + cq = cq1; + dq = dq1; + } else { + acx = acx1; + bcx = bcx1; + ccx = ccx1; + dcx = dcx1; + } + } + } + if (realcase[whichcase] == 1) { + /* if alpha1, beta1, alpha2 and beta2 are real first refine + * the coefficient through a Newton-Raphson */ + oqs_NRabcd(a, b, c, d, &aq, &bq, &cq, &dq); + /* finally calculate the roots as roots of p1(x) and p2(x) (see end of + * sec. 2.1) */ + oqs_solve_quadratic(aq, bq, qroots); + roots[0] = qroots[0]; + roots[1] = qroots[1]; + oqs_solve_quadratic(cq, dq, qroots); + roots[2] = qroots[0]; + roots[3] = qroots[1]; + } else { + /* complex coefficients of p1 and p2 */ + if (whichcase == 0) // d2!=0 + { + cdiskr = acx * acx / 4 - bcx; + /* calculate the roots as roots of p1(x) and p2(x) (see end of sec. 2.1) + */ + zx1 = -acx / 2 + csqrt(cdiskr); + zx2 = -acx / 2 - csqrt(cdiskr); + if (cabs(zx1) > cabs(zx2)) + zxmax = zx1; + else + zxmax = zx2; + zxmin = bcx / zxmax; + roots[0] = zxmin; + roots[1] = conj(zxmin); + roots[2] = zxmax; + roots[3] = conj(zxmax); + } else // d2 ~ 0 + { + /* never gets here! */ + cdiskr = csqrt(acx * acx - 4.0 * bcx); + zx1 = -0.5 * (acx + cdiskr); + zx2 = -0.5 * (acx - cdiskr); + if (cabs(zx1) > cabs(zx2)) + zxmax = zx1; + else + zxmax = zx2; + zxmin = bcx / zxmax; + roots[0] = zxmax; + roots[1] = zxmin; + cdiskr = csqrt(ccx * ccx - 4.0 * dcx); + zx1 = -0.5 * (ccx + cdiskr); + zx2 = -0.5 * (ccx - cdiskr); + if (cabs(zx1) > cabs(zx2)) + zxmax = zx1; + else + zxmax = zx2; + zxmin = dcx / zxmax; + roots[2] = zxmax; + roots[3] = zxmin; + } + } + if (rfact != 1.0) { + for (k = 0; k < 4; k++) + roots[k] *= rfact; + } +} diff --git a/src/surface.cpp b/src/surface.cpp index c63b54b9ad..38a70fcdb1 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1,6 +1,7 @@ #include "openmc/surface.h" #include +#include #include #include @@ -10,6 +11,7 @@ #include "openmc/array.h" #include "openmc/container_util.h" #include "openmc/error.h" +#include "openmc/external/quartic_solver.h" #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" #include "openmc/random_lcg.h" @@ -1004,169 +1006,6 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -//============================================================================== -//============================================================================== -// Generic functions for quadratic, cubic & quartic solver -//============================================================================== - -int quadratic_solve(double a, double b, double c, std::array& x) -{ - - double func = (b * b) - 4 * a * c; - - if (func < 0) { - // this would be imaginary - } else { - x[0] = -b / (2. * a) - std::sqrt(func) / (2. * a); - x[1] = -b / (2. * a) + std::sqrt(func) / (2. * a); - } - - return 0; -} - -const double M_2PI = 2 * PI; -const double eps = FP_COINCIDENT; - -// typedef std::complex DComplex; - -//--------------------------------------------------------------------------- -// solve cubic equation x^3 + a*x^2 + b*x + c -// x - array of size 3 -// In case 3 real roots: => x[0], x[1], x[2], return 3 -// 2 real roots: x[0], x[1], return 2 -// 1 real root : x[0], x[1] ± i*x[2], return 1 -unsigned int solve_cubic( - const double a, const double b, const double c, std::array& x) -{ - double a2 = a * a; - double q = (a2 - 3 * b) / 9; - double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54; - double r2 = r * r; - double q3 = std::pow(q, 3); - double A, B; - double a_prime = 0.; - if (r2 < q3) // 3 roots - { - double t = r / sqrt(q3); - if (t < -1) - t = -1; - if (t > 1) - t = 1; - t = std::acos(t); - a_prime = a / 3.; - q = -2 * std::sqrt(q); - x[0] = q * std::cos(t / 3) - a_prime; - x[1] = q * std::cos((t + M_2PI) / 3) - a_prime; - x[2] = q * std::cos((t - M_2PI) / 3) - a_prime; - return 3; - } else { - A = -std::pow(std::fabs(r) + std::sqrt(r2 - q3), 1. / 3); - if (r < 0) - A = -A; - B = (0 == A ? 0 : q / A); - - a_prime = a / 3; - x[0] = (A + B) - a_prime; - x[1] = -0.5 * (A + B) - a_prime; - x[2] = 0.5 * std::sqrt(3.) * (A - B); - // 2 real roots - if (std::fabs(x[2]) < eps) { - x[2] = x[1]; - return 2; - } - // one real root - return 1; - } -} -//--------------------------------------------------------------------------- -// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d -// Attention - this function returns dynamically allocated array. It has to be -// released afterwards. -void quartic_solve( - double a, double b, double c, double d, std::array& real_roots) -{ - double a3 = -b; - double b3 = a * c - 4. * d; - double c3 = -a * a * d - c * c + 4. * b * d; - - // cubic resolvent - // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0 - - std::array cube_roots; - unsigned int num_roots = solve_cubic(a3, b3, c3, cube_roots); - - double q1, q2, p1, p2, D, sqD, y; - - y = cube_roots[0]; - // The essence - choosing Y with maximal absolute value. - if (num_roots != 1) { - if (std::fabs(cube_roots[1]) > std::fabs(y)) - y = cube_roots[1]; - if (std::fabs(cube_roots[2]) > std::fabs(y)) - y = cube_roots[2]; - } - - // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q) - - D = y * y - 4 * d; - if (std::fabs(D) < eps) // in other words - D==0 - { - q1 = q2 = y * 0.5; - // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g) - D = a * a - 4 * (b - y); - if (std::fabs(D) < eps) { - p1 = p2 = a * 0.5; - } else { - sqD = std::sqrt(D); - p1 = (a + sqD) * 0.5; - p2 = (a - sqD) * 0.5; - } - } else { - sqD = std::sqrt(D); - q1 = (y + sqD) * 0.5; - q2 = (y - sqD) * 0.5; - p1 = (a * q1 - c) / (q1 - q2); - p2 = (c - a * q2) / (q1 - q2); - } - - std::array, 4> roots; // the roots to return - - // solving quadratic eq. - x^2 + p1*x + q1 = 0 - D = p1 * p1 - 4 * q1; - if (D < 0.0) { - roots[0].real(-p1 * 0.5); - roots[0].imag(std::sqrt(-D) * 0.5); - roots[1] = std::conj(roots[0]); - } else { - sqD = std::sqrt(D); - roots[0].real((-p1 + sqD) * 0.5); - roots[1].real((-p1 - sqD) * 0.5); - } - - // solving quadratic eq. - x^2 + p2*x + q2 = 0 - D = p2 * p2 - 4 * q2; - if (D < 0.0) { - roots[2].real(-p2 * 0.5); - roots[2].imag(std::sqrt(-D) * 0.5); - roots[3] = std::conj(roots[2]); - } else { - sqD = std::sqrt(D); - roots[2].real((-p2 + sqD) * 0.5); - roots[3].real((-p2 - sqD) * 0.5); - } - - for (int i = 0; i < 4; i++) { - if (roots[i].imag() == 0.) - real_roots[i] = roots[i].real(); - else - real_roots[i] = 0.; - } - - std::sort(real_roots.begin(), real_roots.end()); - - return; -} - //============================================================================== SurfaceXTorus::SurfaceXTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) @@ -1212,37 +1051,31 @@ double SurfaceXTorus::distance(Position r, Direction ang, bool coincident) const double d0 = four_A2 * (y * y + z * z); // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0 - double a = c2 * c2; - double b = 2 * c1 * c2; - double c = c1 * c1 + 2 * c0 * c2 - d2; - double d = 2 * c0 * c1 - d1; - double e = c0 * c0 - d0; + double coeff[5]; + coeff[0] = c0 * c0 - d0; + coeff[1] = 2 * c0 * c1 - d1; + coeff[2] = c1 * c1 + 2 * c0 * c2 - d2; + coeff[3] = 2 * c1 * c2; + coeff[4] = c2 * c2; - std::array roots; - quartic_solve(b, c, d, e, roots); + std::complex roots[4]; + oqs_quartic_solver(coeff, roots); - if (coincident) - r += ang * TINY_BIT; - - // special degerenate case two sets of repated - // roots - if (b == 0.0 && d == 0) { - if (roots[1] - roots[0] < 1e-5) - return INFTY; - if (roots[3] - roots[2] < 1e-5) - return INFTY; - } - - for (int i = 0; i < 4; i++) { - // need something better than just raw tolerance - // use fastQS to get back lost precision - if (roots[i] > 1e-6) { - return roots[i]; + // Find smallest positive, real root. In the case where the particle is + // coincident with the surface, we are sure to have one root very close to + // zero but possibly small and positive. A tolerance is set to discard that + // zero. + double distance = INFTY; + double cutoff = coincident ? 1e-9 : 0.0; + for (int i = 0; i < 4; ++i) { + if (roots[i].imag() == 0) { + double root = roots[i].real(); + if (root > cutoff && root < distance) { + distance = root; + } } } - - // otherwise no hit - return INFTY; + return distance; } Direction SurfaceXTorus::normal(Position r) const @@ -1306,36 +1139,31 @@ double SurfaceYTorus::distance(Position r, Direction ang, bool coincident) const double d0 = four_A2 * (x * x + z * z); // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0 - double a = c2 * c2; - double b = 2 * c1 * c2; - double c = c1 * c1 + 2 * c0 * c2 - d2; - double d = 2 * c0 * c1 - d1; - double e = c0 * c0 - d0; + double coeff[5]; + coeff[0] = c0 * c0 - d0; + coeff[1] = 2 * c0 * c1 - d1; + coeff[2] = c1 * c1 + 2 * c0 * c2 - d2; + coeff[3] = 2 * c1 * c2; + coeff[4] = c2 * c2; - std::array roots; - quartic_solve(b, c, d, e, roots); + std::complex roots[4]; + oqs_quartic_solver(coeff, roots); - if (coincident) - r += ang * TINY_BIT; - - // special degerenate case two sets of repated - if (b == 0.0 && d == 0) { - if (roots[1] - roots[0] < 1e-5) - return INFTY; - if (roots[3] - roots[2] < 1e-5) - return INFTY; - } - - for (int i = 0; i < 4; i++) { - // need something better than just raw tolerance - // use fastQS to get back lost precision - if (roots[i] > 1e-6) { - return roots[i]; + // Find smallest positive, real root. In the case where the particle is + // coincident with the surface, we are sure to have one root very close to + // zero but possibly small and positive. A tolerance is set to discard that + // zero. + double distance = INFTY; + double cutoff = coincident ? 1e-9 : 0.0; + for (int i = 0; i < 4; ++i) { + if (roots[i].imag() == 0) { + double root = roots[i].real(); + if (root > cutoff && root < distance) { + distance = root; + } } } - - // otherwise no hit - return INFTY; + return distance; } Direction SurfaceYTorus::normal(Position r) const @@ -1399,37 +1227,31 @@ double SurfaceZTorus::distance(Position r, Direction ang, bool coincident) const double d0 = four_A2 * (x * x + y * y); // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0 - double a = c2 * c2; - double b = 2 * c1 * c2; - double c = c1 * c1 + 2 * c0 * c2 - d2; - double d = 2 * c0 * c1 - d1; - double e = c0 * c0 - d0; + double coeff[5]; + coeff[0] = c0 * c0 - d0; + coeff[1] = 2 * c0 * c1 - d1; + coeff[2] = c1 * c1 + 2 * c0 * c2 - d2; + coeff[3] = 2 * c1 * c2; + coeff[4] = c2 * c2; - std::array roots; - quartic_solve(b, c, d, e, roots); + std::complex roots[4]; + oqs_quartic_solver(coeff, roots); - if (coincident) - r += ang * TINY_BIT; - - // special degerenate case two sets of repated - // roots - if (b == 0.0 && d == 0) { - if (roots[1] - roots[0] < 1e-5) - return INFTY; - if (roots[3] - roots[2] < 1e-5) - return INFTY; - } - - for (int i = 0; i < 4; i++) { - // need something better than just raw tolerance - // use fastQS to get back lost precision - if (roots[i] > 1e-6) { - return roots[i]; + // Find smallest positive, real root. In the case where the particle is + // coincident with the surface, we are sure to have one root very close to + // zero but possibly small and positive. A tolerance is set to discard that + // zero. + double distance = INFTY; + double cutoff = coincident ? 1e-9 : 0.0; + for (int i = 0; i < 4; ++i) { + if (roots[i].imag() == 0) { + double root = roots[i].real(); + if (root > cutoff && root < distance) { + distance = root; + } } } - - // otherwise no hit - return INFTY; + return distance; } Direction SurfaceZTorus::normal(Position r) const