/******************************************************************************* * This file is part of SWIFT. * Copyright (C) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ #include <config.h> /* Local includes. */ #include "swift.h" /* Standard includes */ #include <fenv.h> #include <math.h> /** * Compute erfcf(u) using eq. 7.1.26 of * Abramowitz & Stegun, 1972. * * This has a *relative* error of less than 3.4e-3 over * the range of interest (0 < x < 5) */ float optimized_erfcf(const float x) { const float x2 = x * x; const float exp_x2 = expf(-x2); const float t = 1.f / (1.f + 0.3275911f * x); const float a1 = 0.254829592f; const float a2 = -0.284496736f; const float a3 = 1.421413741f; const float a4 = -1.453152027; const float a5 = 1.061405429f; /* a1 * t + a2 * t^2 + a3 * t^3 + a4 * t^4 + a5 * t^5 */ float a = a5 * t + a4; a = a * t + a3; a = a * t + a2; a = a * t + a1; a = a * t; return a * exp_x2; } /** * @brief Check that a and b are consistent (up to some relative error) * * @param a First value * @param b Second value * @param rel_tol Relative tolerance * @param abs_tol Absolute tolerance * @param x Value for which we tested the function (for error messages) */ void check_value(const double a, const double b, const double rel_tol, const double abs_tol, const double x) { if (fabs(a - b) / fabs(a + b) > rel_tol) error("Values are inconsistent: %12.15e %12.15e rel=%e (for x=%e).", a, b, fabs(a - b) / fabs(a + b), x); if (fabs(a - b) > abs_tol) error("Values are inconsistent: %12.15e %12.15e abs=%e (for x=%e).", a, b, fabs(a - b), x); } int main(int argc, char* argv[]) { /* Initialize CPU frequency, this also starts time. */ unsigned long long cpufreq = 0; clocks_set_cpufreq(cpufreq); /* Choke on FPEs */ #ifdef HAVE_FE_ENABLE_EXCEPT feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); #endif /* Loop over some values */ for (float x = 0.f; x < 5.f; x += 0.000001f) { const double exact = erfc(x); const double swift_erfcf = optimized_erfcf(x); check_value(exact, swift_erfcf, 3.358e-3, 6.1e-7, x); } return 0; }