crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Use quadric solver from ACM TOMS 46(2), pp1-28 (algorithm 1010)

Browse source
This commit is contained in:
Paul Romano 2021-12-20 09:28:21 -05:00
parent b557e62b5d
commit 377ee77486
5 changed files with 732 additions and 243 deletions
Download patch file Download diff file Expand all files Collapse all files

4
CMakeLists.txt
View file

@ -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)

10
include/openmc/external/quartic_solver.h vendored Normal file
View file

@ -0,0 +1,10 @@
#ifndef OPENMC_EXTERNAL_QUARTIC_SOLVER_H
#define OPENMC_EXTERNAL_QUARTIC_SOLVER_H
#include <complex>
extern "C" {
void oqs_quartic_solver(double coeff[5], std::complex<double> roots[4]);
}
#endif // OPENMC_EXTERNAL_QUARTIC_SOLVER_H

32
src/external/LICENSE vendored Normal file
View file

@ -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.

621
src/external/quartic_solver.c vendored Normal file
View file

@ -0,0 +1,621 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#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;
}
}

308
src/surface.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/surface.h"
#include <cmath>
#include <complex>
#include <set>
#include <utility>
@ -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<double, 2>& 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<double> 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<double, 3>& 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<double, 4>& 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<double, 3> 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<std::complex<double>, 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<double, 4> roots;
quartic_solve(b, c, d, e, roots);
std::complex<double> 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<double, 4> roots;
quartic_solve(b, c, d, e, roots);
std::complex<double> 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<double, 4> roots;
quartic_solve(b, c, d, e, roots);
std::complex<double> 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