Skip to content

Commit

Permalink
adding exponent regression test for quad-double: exp algorithm is way…
Browse files Browse the repository at this point in the history
… off
  • Loading branch information
Ravenwater committed Aug 15, 2024
1 parent e97b3fd commit 8a4d57f
Show file tree
Hide file tree
Showing 5 changed files with 215 additions and 80 deletions.
2 changes: 1 addition & 1 deletion include/universal/number/dd/dd_impl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -965,7 +965,7 @@ inline dd div(double a, double b) {
return dd(s, e);
}

// double-double * double, where double is a power of 2. */
// double-double * double, where double is a power of 2
inline dd mul_pwr2(const dd& a, double b) {
return dd(a.high() * b, a.low() * b);
}
Expand Down
158 changes: 80 additions & 78 deletions include/universal/number/qd/math/exponent.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,93 +11,95 @@
namespace sw { namespace universal {

// fwd reference
qd ldexp(qd const&, int);

static const int n_inv_fact = 15;
static const double inv_fact[n_inv_fact][2] = {
{ 1.66666666666666657e-01, 9.25185853854297066e-18},
{ 4.16666666666666644e-02, 2.31296463463574266e-18},
{ 8.33333333333333322e-03, 1.15648231731787138e-19},
{ 1.38888888888888894e-03, -5.30054395437357706e-20},
{ 1.98412698412698413e-04, 1.72095582934207053e-22},
{ 2.48015873015873016e-05, 2.15119478667758816e-23},
{ 2.75573192239858925e-06, -1.85839327404647208e-22},
{ 2.75573192239858883e-07, 2.37677146222502973e-23},
{ 2.50521083854417202e-08, -1.44881407093591197e-24},
{ 2.08767569878681002e-09, -1.20734505911325997e-25},
{ 1.60590438368216133e-10, 1.25852945887520981e-26},
{ 1.14707455977297245e-11, 2.06555127528307454e-28},
{ 7.64716373181981641e-13, 7.03872877733453001e-30},
{ 4.77947733238738525e-14, 4.39920548583408126e-31},
{ 2.81145725434552060e-15, 1.65088427308614326e-31}
};

// Base-e exponential function
qd exp(const qd& a) {
/* Strategy: We first reduce the size of x by noting that
exp(kr + m * log(2)) = 2^m * exp(r)^k
where m and k are integers. By choosing m appropriately
we can make |kr| <= log(2) / 2 = 0.347. Then exp(r) is
evaluated using the familiar Taylor series. Reducing the
argument substantially speeds up the convergence. */

const double k = 512.0;
const double inv_k = 1.0 / k;

if (a.high() <= -709.0) return qd(0.0);

if (a.high() >= 709.0) return qd(SpecificValue::infpos);

if (a.iszero()) return qd(1.0);

if (a.isone()) return qd_e;

double m = std::floor(a.high() / qd_log2.high() + 0.5);
qd r = mul_pwr2(a - qd_log2 * m, inv_k);
qd s, t, p;

p = sqr(r);
s = r + mul_pwr2(p, 0.5);
p *= r;
t = p * qd(inv_fact[0][0], inv_fact[0][1]);
int i = 0;
do {
s += t;
p *= r;
++i;
t = p * qd(inv_fact[i][0], inv_fact[i][1]);
} while (std::abs(double(t)) > inv_k * d_eps && i < 5);

s += t;

s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s += 1.0;

return ldexp(s, static_cast<int>(m));
}
qd ldexp(const qd&, int);

constexpr unsigned n_inv_fact = 15;
static const qd inv_fact[n_inv_fact] = {
qd(1.66666666666666657e-01, 9.25185853854297066e-18, 5.13581318503262866e-34, 2.85094902409834186e-50),
qd(4.16666666666666644e-02, 2.31296463463574266e-18, 1.28395329625815716e-34, 7.12737256024585466e-51),
qd(8.33333333333333322e-03, 1.15648231731787138e-19, 1.60494162032269652e-36, 2.22730392507682967e-53),
qd(1.38888888888888894e-03, -5.30054395437357706e-20, -1.73868675534958776e-36, -1.63335621172300840e-52),
qd(1.98412698412698413e-04, 1.72095582934207053e-22, 1.49269123913941271e-40, 1.29470326746002471e-58),
qd(2.48015873015873016e-05, 2.15119478667758816e-23, 1.86586404892426588e-41, 1.61837908432503088e-59),
qd(2.75573192239858925e-06, -1.85839327404647208e-22, 8.49175460488199287e-39, -5.72661640789429621e-55),
qd(2.75573192239858883e-07, 2.37677146222502973e-23, -3.26318890334088294e-40, 1.61435111860404415e-56),
qd(2.50521083854417202e-08, -1.44881407093591197e-24, 2.04267351467144546e-41, -8.49632672007163175e-58),
qd(2.08767569878681002e-09, -1.20734505911325997e-25, 1.70222792889287100e-42, 1.41609532150396700e-58),
qd(1.60590438368216133e-10, 1.25852945887520981e-26, -5.31334602762985031e-43, 3.54021472597605528e-59),
qd(1.14707455977297245e-11, 2.06555127528307454e-28, 6.88907923246664603e-45, 5.72920002655109095e-61),
qd(7.64716373181981641e-13, 7.03872877733453001e-30, -7.82753927716258345e-48, 1.92138649443790242e-64),
qd(4.77947733238738525e-14, 4.39920548583408126e-31, -4.89221204822661465e-49, 1.20086655902368901e-65),
qd(2.81145725434552060e-15, 1.65088427308614326e-31, -2.87777179307447918e-50, 4.27110689256293549e-67)
};

qd exp(const qd& x) {
/* Strategy: We first reduce the size of x by noting that
exp(kr + m * log(2)) = 2^m * exp(r)^k
where m and k are integers. By choosing m appropriately
we can make |kr| <= log(2) / 2 = 0.347. Then exp(r) is
evaluated using the familiar Taylor series. Reducing the
argument substantially speeds up the convergence.
*/

constexpr double k = double(1ull << 16);
constexpr double inv_k = 1.0 / k;

if (x[0] <= -709.0) return 0.0;

if (x[0] >= 709.0) return qd(SpecificValue::infpos);

if (x.iszero()) return 1.0;

if (x.isone()) return qd_e;

double m = std::floor(x[0] / qd_log2[0] + 0.5);
qd r = mul_pwr2(x - qd_log2 * m, inv_k);
qd s, p, t;
double thresh = inv_k * qd_eps;

p = sqr(r);
s = r + mul_pwr2(p, 0.5);
int i = 0;
do {
p *= r;
t = p * inv_fact[i++];
s += t;
} while (std::abs(double(t)) > thresh && i < 9);

s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s = mul_pwr2(s, 2.0) + sqr(s);
s += 1.0;
return ldexp(s, static_cast<int>(m));
}

// Base-2 exponential function
qd exp2(qd x) {
qd exp2(const qd& x) {
return qd(std::exp2(double(x)));
}

// Base-10 exponential function
qd exp10(qd x) {
return qd(std::exp10(double(x)));
qd exp10(const qd& x) {
return qd(std::pow(10.0, double(x)));
}

// Base-e exponential function exp(x)-1
qd expm1(qd x) {
qd expm1(const qd& x) {
return qd(std::expm1(double(x)));
}

Expand Down
2 changes: 1 addition & 1 deletion include/universal/number/qd/mathlib.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
#include <universal/number/qd/math/classify.hpp>
//#include <universal/number/qd/math/complex.hpp>
#include <universal/number/qd/math/error_and_gamma.hpp>
//#include <universal/number/qd/math/exponent.hpp>
#include <universal/number/qd/math/exponent.hpp>
//#include <universal/number/qd/math/fractional.hpp>
//#include <universal/number/qd/math/hyperbolic.hpp>
//#include <universal/number/qd/math/hypot.hpp>
Expand Down
7 changes: 7 additions & 0 deletions include/universal/number/qd/qd_impl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -1094,6 +1094,8 @@ constexpr qd qd_inv_sqrt2(0.70710678118654757, -4.8336466567264567e-17, 2.069337

constexpr qd qd_max(1.79769313486231570815e+308, 9.97920154767359795037e+291);

constexpr double qd_eps = 4.93038065763132e-32; // 2^-104
constexpr double qd_min_normalized = 2.0041683600089728e-292; // = 2^(-1022 + 53)

//////////////////////// helper functions /////////////////////////////////

Expand Down Expand Up @@ -1246,6 +1248,11 @@ inline qd quick_nint(const qd& a) {
return r;
}

// quad-double * double, where double is a power of 2
inline qd mul_pwr2(const qd& a, double b) {
return qd(a[0] * b, a[1] * b, a[2] * b, a[3] * b);
}

/* quad-double ^ 2 = (x0 + x1 + x2 + x3) ^ 2
= x0 ^ 2 + 2 x0 * x1 + (2 x0 * x2 + x1 ^ 2)
+ (2 x0 * x3 + 2 x1 * x2) */
Expand Down
126 changes: 126 additions & 0 deletions static/qd/math/exponent.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@
// exponent.cpp: test suite runner for exponentiation function for quad-double (qd) floats
//
// Copyright (C) 2017 Stillwater Supercomputing, Inc.
// SPDX-License-Identifier: MIT
//
// This file is part of the universal numbers project, which is released under an MIT Open Source license.
#include <universal/utility/directives.hpp>
#include <universal/number/qd/qd.hpp>
#include <universal/verification/test_suite.hpp>

// generate specific test case
template<typename Ty>
void GenerateTestCase(Ty fa) {
unsigned precision = 25;
unsigned width = 30;
Ty fref;
sw::universal::qd a, ref, v;
a = fa;
fref = std::exp(fa);
ref = fref;
v = sw::universal::exp(a);
auto oldPrec = std::cout.precision();
std::cout << std::setprecision(precision);
std::cout << " -> exp(" << fa << ") = " << std::setw(width) << fref << std::endl;
std::cout << " -> exp( " << a << ") = " << v << '\n' << to_binary(v) << '\n';
std::cout << to_binary(ref) << "\n -> reference\n";
std::cout << (ref == v ? "PASS" : "FAIL") << std::endl << std::endl;
std::cout << std::setprecision(oldPrec);
}

// Regression testing guards: typically set by the cmake configuration, but MANUAL_TESTING is an override
#define MANUAL_TESTING 1
// REGRESSION_LEVEL_OVERRIDE is set by the cmake file to drive a specific regression intensity
// It is the responsibility of the regression test to organize the tests in a quartile progression.
//#undef REGRESSION_LEVEL_OVERRIDE
#ifndef REGRESSION_LEVEL_OVERRIDE
#undef REGRESSION_LEVEL_1
#undef REGRESSION_LEVEL_2
#undef REGRESSION_LEVEL_3
#undef REGRESSION_LEVEL_4
#define REGRESSION_LEVEL_1 1
#define REGRESSION_LEVEL_2 1
#define REGRESSION_LEVEL_3 1
#define REGRESSION_LEVEL_4 1
#endif

int main()
try {
using namespace sw::universal;

std::string test_suite = "quad-double mathlib exponentiation function validation";
std::string test_tag = "exp/exp2/exp10/expm1";
bool reportTestCases = false;
int nrOfFailedTestCases = 0;

ReportTestSuiteHeader(test_suite, reportTestCases);

#if MANUAL_TESTING
// generate individual testcases to hand trace/debug
GenerateTestCase(4.0);

auto oldPrec = std::cout.precision();
for (int i = 0; i < 30; ++i) {
std::string tag = "exp(" + std::to_string(i) + ")";
double exponentRef = std::exp(double(i));
qd exponent = exp(qd(i));
qd error = exponentRef - exponent;
std::cout << std::setw(20) << tag << " : " << std::setprecision(32) << exponentRef << " : " << exponent << " : " << std::setw(25) << error << '\n';
}

for (int i = 0; i < 30; ++i) {
std::string tag = "exp2(" + std::to_string(i) + ")";
double exponentRef = std::exp2(double(i));
qd exponent = exp2(qd(i));
qd error = exponentRef - exponent;
std::cout << std::setw(20) << tag << " : " << std::setprecision(32) << exponentRef << " : " << exponent << " : " << std::setw(25) << error << '\n';
}

for (int i = 0; i < 30; ++i) {
std::string tag = "exp10(" + std::to_string(i) + ")";
double exponentRef = std::pow(10.0, double(i));
qd exponent = exp10(qd(i));
qd error = exponentRef - exponent;
std::cout << std::setw(20) << tag << " : " << std::setprecision(32) << exponentRef << " : " << exponent << " : " << std::setw(25) << error << '\n';
}

for (int i = 0; i < 30; ++i) {
std::string tag = "expm1(" + std::to_string(i) + ")";
double exponentRef = std::expm1(double(i));
qd exponent = expm1(qd(i));
qd error = exponentRef - exponent;
std::cout << std::setw(20) << tag << " : " << std::setprecision(32) << exponentRef << " : " << exponent << " : " << std::setw(25) << error << '\n';
}

std::cout << std::setprecision(oldPrec);

ReportTestSuiteResults(test_suite, nrOfFailedTestCases);
return EXIT_SUCCESS; // ignore errors
#else


ReportTestSuiteResults(test_suite, nrOfFailedTestCases);
return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS);

#endif // MANUAL_TESTING
}
catch (char const* msg) {
std::cerr << "Caught ad-hoc exception: " << msg << std::endl;
return EXIT_FAILURE;
}
catch (const sw::universal::universal_arithmetic_exception& err) {
std::cerr << "Caught unexpected universal arithmetic exception : " << err.what() << std::endl;
return EXIT_FAILURE;
}
catch (const sw::universal::universal_internal_exception& err) {
std::cerr << "Caught unexpected universal internal exception: " << err.what() << std::endl;
return EXIT_FAILURE;
}
catch (const std::runtime_error& err) {
std::cerr << "Caught runtime exception: " << err.what() << std::endl;
return EXIT_FAILURE;
}
catch (...) {
std::cerr << "Caught unknown exception" << std::endl;
return EXIT_FAILURE;
}

0 comments on commit 8a4d57f

Please sign in to comment.