diff --git a/include/universal/math/constants/ieee754_constants.hpp b/include/universal/math/constants/ieee754_constants.hpp new file mode 100644 index 000000000..7874c2aab --- /dev/null +++ b/include/universal/math/constants/ieee754_constants.hpp @@ -0,0 +1,115 @@ +#pragma once +// ieee754_constants.hpp: definition of math constants in different ieee-754 floating-piont precisions +// +// 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. + +namespace sw { namespace universal { + +// best practice for C++11 +// inject mathematical constants in our namespace + +// a_b reads a over b, as in 1_pi being 1 over pi + + /////////////////////////////////////////////////////////////////////////////////////////////////////// +/// float values +constexpr double f_pi_4 = 0.785398163f; // pi/4 +constexpr double f_pi_3 = 1.047197551f; // pi/3 +constexpr double f_pi_2 = 1.570796327f; // pi/2 +constexpr double f_pi = 3.141592653f; // pi +constexpr double f_3pi_4 = 4.712388980f; // 3*pi/4 +constexpr double f_2pi = 6.283185307f; // 2*pi + +constexpr double f_phi = 1.618033989f; // phi == golden ratio == most irrational number + +constexpr double f_e = 2.718281828f; // e +constexpr double f_log2e = 1.442695041f; // log2(e) +constexpr double f_log10e = 0.434294482f; // log10(e) + +constexpr double f_ln2 = 0.693147181f; // ln(2) +constexpr double f_ln10 = 2.302585093f; // ln(10) + +constexpr double f_sqrt2 = 1.4142135624f; // sqrt(2) +constexpr double f_sqrt3 = 1.7320508076f; // sqrt(3) +constexpr double f_sqrt5 = 2.2360679775f; // sqrt(5) + +constexpr double f_1_phi = 0.6180339887f; // 1/phi + +constexpr double f_1_e = 0.3678794411f; // 1/e + +constexpr double f_1_pi = 0.318309886f; // 1/pi +constexpr double f_2_pi = 0.636619772f; // 2/pi + +constexpr double f_1_sqrt2 = 0.707106781f; // 1/sqrt(2) +constexpr double f_2_sqrtpi = 1.128379167f; // 2/sqrt(pi) + + +/////////////////////////////////////////////////////////////////////////////////////////////////////// +/// double values +constexpr double d_pi_4 = 0.785398163397448309616; // pi/4 +constexpr double d_pi_3 = 1.04719755119659774615; // pi/3 +constexpr double d_pi_2 = 1.57079632679489661923; // pi/2 +constexpr double d_pi = 3.14159265358979323846; // pi +constexpr double d_3pi_4 = 4.71238898038468985769; // 3*pi/4 +constexpr double d_2pi = 6.28318530717958647693; // 2*pi + +constexpr double d_phi = 1.61803398874989484820; // phi == golden ratio == most irrational number + +constexpr double d_e = 2.71828182845904523536; // e +constexpr double d_log2e = 1.44269504088896340736; // log2(e) +constexpr double d_log10e = 0.434294481903251827651; // log10(e) + +constexpr double d_ln2 = 0.693147180559945309417; // ln(2) +constexpr double d_ln10 = 2.30258509299404568402; // ln(10) + +constexpr double d_sqrt2 = 1.414213562373095048801; // sqrt(2) +constexpr double d_sqrt3 = 1.732050807568877293527; // sqrt(3) +constexpr double d_sqrt5 = 2.236067977499789696409; // sqrt(5) + +constexpr double d_1_phi = 0.618033988749894791503; // 1/phi + +constexpr double d_1_e = 0.367879441171442321596; // 1/e + +constexpr double d_1_pi = 0.318309886183790671538; // 1/pi +constexpr double d_2_pi = 0.636619772367581343076; // 2/pi + +constexpr double d_1_sqrt2 = 0.707106781186547524401; // 1/sqrt(2) +constexpr double d_2_sqrtpi = 1.12837916709551257390; // 2/sqrt(pi) == 1 / sqrt(pi / 4) + + +/////////////////////////////////////////////////////////////////////////////////////////////////////// +/// long double values +constexpr long double ld_pi_4 = 0.78539816339744830961566084581988l; // pi/4 +constexpr long double ld_pi_3 = 1.0471975511965977461542144610932l; // pi/3 +constexpr long double ld_pi_2 = 1.5707963267948966192313216916398l; // pi/2 +constexpr long double ld_pi = 3.1415926535897932384626433832795l; // pi +constexpr long double ld_3pi_4 = 4.7123889803846898576939650749193l; // 3*pi/4 +constexpr long double ld_2pi = 6.283185307179586476925286766559l; // 2*pi + +// TODO: generate long double constants for the following constants +constexpr long double ld_phi = 1.61803398874989484820; // phi == golden ratio == most irrational number + +constexpr long double ld_e = 2.71828182845904523536l; // e +constexpr long double ld_log2e = 1.44269504088896340736l; // log2(e) +constexpr long double ld_log10e = 0.434294481903251827651l; // log10(e) + +constexpr long double ld_ln2 = 0.693147180559945309417l; // ln(2) +constexpr long double ld_ln10 = 2.30258509299404568402l; // ln(10) + +constexpr long double ld_sqrt2 = 1.41421356237309504880l; // sqrt(2) +constexpr long double ld_sqrt3 = 1.732050807568877293527l; // sqrt(2) +constexpr long double ld_sqrt5 = 2.236067977499789696409l; // sqrt(2) + +constexpr long double ld_1_phi = 0.618033988749894791503l; // 1/phi + +constexpr long double ld_1_e = 0.367879441171442321596l; // 1/e + +constexpr long double ld_1_pi = 0.318309886183790671538l; // 1/pi +constexpr long double ld_2_pi = 0.636619772367581343076l; // 2/pi + +constexpr long double ld_1_sqrt2 = 0.707106781186547524401l; // 1/sqrt(2) +constexpr long double ld_2_sqrtpi = 1.12837916709551257390l; // 2/sqrt(pi) == 1 / sqrt(pi / 4) + +}} // namespace sw::universal diff --git a/include/universal/number/qd/math/complex.hpp b/include/universal/number/qd/math/complex/complex.hpp similarity index 81% rename from include/universal/number/qd/math/complex.hpp rename to include/universal/number/qd/math/complex/complex.hpp index 79adf6c46..4c7c4e9cc 100644 --- a/include/universal/number/qd/math/complex.hpp +++ b/include/universal/number/qd/math/complex/complex.hpp @@ -1,5 +1,5 @@ #pragma once -// complex.hpp: complex support for double-double floating-point +// complex.hpp: complex support for quad-double floating-point // // Copyright (C) 2017 Stillwater Supercomputing, Inc. // SPDX-License-Identifier: MIT diff --git a/include/universal/number/qd/math/constants/qd_constants.hpp b/include/universal/number/qd/math/constants/qd_constants.hpp new file mode 100644 index 000000000..1217717b8 --- /dev/null +++ b/include/universal/number/qd/math/constants/qd_constants.hpp @@ -0,0 +1,61 @@ +#pragma once +// qd_constants.hpp: definition of math constants in quad-double precision +// +// 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. + +namespace sw { namespace universal { + +// best practice for C++11 +// inject mathematical constants in our namespace + +// a_b reads a over b, as in 1_pi being 1 over pi + +// precomputed quad-double constants courtesy Scibuilders, Jack Poulson +// and Stillwater Supercomputing, Theodore Omtzigt +// +// pi multiples and fractions +constexpr qd qd_pi_4 (0.78539816339744828, 3.061616997868383e-17, -7.4869245242958492e-34, 2.7811355521584142e-50); // pi/4 +constexpr qd qd_pi_3 (1.570796326794896558, 6.123233995736766036e-17); // pi/3 +constexpr qd qd_pi_2 (1.5707963267948966, 6.123233995736766e-17, -1.4973849048591698e-33, 5.5622711043168283e-50); // pi/2 +constexpr qd qd_pi (3.1415926535897931, 1.2246467991473532e-16, -2.9947698097183397e-33, 1.1124542208633657e-49); // pi +constexpr qd qd_3pi_4 (2.3561944901923448, 9.1848509936051484e-17, 3.9168984647504003e-33, -2.586798163270486e-49); // 3*pi/4 +constexpr qd qd_2pi (6.2831853071795862, 2.4492935982947064e-16, -5.9895396194366793e-33, 2.2249084417267313e-49); // 2*pi + +// Golden ratio PHI +constexpr qd qd_phi (1.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975020988e-50); // phi == golden ratio == most irrational number + +// Euler's number e +constexpr qd qd_e (2.7182818284590451, 1.4456468917292502e-16, -2.1277171080381768e-33, 1.515630159841219e-49); // e + +// natural logarithm (base = e) +constexpr qd qd_ln2 (0.69314718055994529, 2.3190468138462996e-17, 5.7077084384162121e-34, -3.5824322106018105e-50); // ln(2) = elog(2) +constexpr qd qd_ln10 (2.3025850929940459, -2.1707562233822494e-16, -9.9842624544657766e-33, -4.0233574544502071e-49); // ln(10) = elog(10) + +// binary logarithm (base = 2) +constexpr qd qd_lge (1.4426950408889634, 2.0355273740931033e-17, -1.0614659956117258e-33, -1.3836716780181395e-50); // 2log(e) +constexpr qd qd_lg10 (3.3219280948873622, 1.6616175169735920e-16, +1.2215512178458181e-32, +5.9551189702782481e-49); // 2log(10) +// common logarithm (base = 10) +constexpr qd qd_log2 (0.3010299956639812, -2.8037281277851704e-18, 5.4719484023146385e-35, 5.1051389831070996e-51); // 10log(2) +constexpr qd qd_loge (0.43429448190325182, 1.0983196502167651e-17, 3.7171812331109590e-34, 7.7344843465042927e-51); // 10log(3) + +// square roots +constexpr qd qd_sqrt2 (1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683509e-50); // sqrt(2) +constexpr qd qd_sqrt3 (1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33, 5.3061475632961675e-50); // sqrt(3) +constexpr qd qd_sqrt5 (2.2360679774997898, -1.0864230407365012e-16, 5.3086504167631309e-33, -6.6099839950042175e-50); // sqrt(5) + +constexpr qd qd_1_phi (0.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975021111e-50); // 1/phi + +constexpr qd qd_1_e (0.36787944117144233, -1.2428753672788363e-17, -5.830044851072742e-34, -2.8267977849017436e-50); // 1/e + +constexpr qd qd_1_pi (0.31830988618379069, -1.9678676675182486e-17, -1.0721436282893004e-33, 8.053563926594112e-50); // 1/pi +constexpr qd qd_2_pi (0.63661977236758138, -3.9357353350364972e-17, -2.1442872565786008e-33, 1.6107127853188224e-49);// 2/pi == 1 / (pi/2) + +//constexpr qd qd_1_sqrt2 (0.70710678118654757, -4.8336466312625432e-17, -3.0392667356269840e-34, -1.350504809842679e-50); // 1 / sqrt(2) +//constexpr qd qd_1_sqrt2 (0.70710678118654757, -4.8336466567264567e-17, +2.0693376543497068e-33, +2.4677734957341745e-50); // 1 / sqrt(2) +constexpr qd qd_1_sqrt2 (0.70710678118654757, -4.8336466567264567e-17, +2.0693376543497068e-33, +2.4677734957341755e-50); // 1 / sqrt(2) +constexpr qd qd_2_sqrtpi (1.1283791670955126, 1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50);// 2 / sqrt(pi) = 1 / (sqrt(pi/4) + +}} // namespace sw::universal diff --git a/include/universal/number/qd/math/sincos_tables.hpp b/include/universal/number/qd/math/constants/sincos_tables.hpp similarity index 100% rename from include/universal/number/qd/math/sincos_tables.hpp rename to include/universal/number/qd/math/constants/sincos_tables.hpp diff --git a/include/universal/number/qd/math/classify.hpp b/include/universal/number/qd/math/functions/classify.hpp similarity index 100% rename from include/universal/number/qd/math/classify.hpp rename to include/universal/number/qd/math/functions/classify.hpp diff --git a/include/universal/number/qd/math/error_and_gamma.hpp b/include/universal/number/qd/math/functions/error_and_gamma.hpp similarity index 100% rename from include/universal/number/qd/math/error_and_gamma.hpp rename to include/universal/number/qd/math/functions/error_and_gamma.hpp diff --git a/include/universal/number/qd/math/exponent.hpp b/include/universal/number/qd/math/functions/exponent.hpp similarity index 100% rename from include/universal/number/qd/math/exponent.hpp rename to include/universal/number/qd/math/functions/exponent.hpp diff --git a/include/universal/number/qd/math/fractional.hpp b/include/universal/number/qd/math/functions/fractional.hpp similarity index 100% rename from include/universal/number/qd/math/fractional.hpp rename to include/universal/number/qd/math/functions/fractional.hpp diff --git a/include/universal/number/qd/math/hyperbolic.hpp b/include/universal/number/qd/math/functions/hyperbolic.hpp similarity index 100% rename from include/universal/number/qd/math/hyperbolic.hpp rename to include/universal/number/qd/math/functions/hyperbolic.hpp diff --git a/include/universal/number/qd/math/hypot.hpp b/include/universal/number/qd/math/functions/hypot.hpp similarity index 100% rename from include/universal/number/qd/math/hypot.hpp rename to include/universal/number/qd/math/functions/hypot.hpp diff --git a/include/universal/number/qd/math/logarithm.hpp b/include/universal/number/qd/math/functions/logarithm.hpp similarity index 99% rename from include/universal/number/qd/math/logarithm.hpp rename to include/universal/number/qd/math/functions/logarithm.hpp index 53f2de81c..5d45ed722 100644 --- a/include/universal/number/qd/math/logarithm.hpp +++ b/include/universal/number/qd/math/functions/logarithm.hpp @@ -88,7 +88,7 @@ namespace sw { namespace universal { } - return log(a) / qd_log10; + return log(a) / qd_ln10; } /// diff --git a/include/universal/number/qd/math/minmax.hpp b/include/universal/number/qd/math/functions/minmax.hpp similarity index 100% rename from include/universal/number/qd/math/minmax.hpp rename to include/universal/number/qd/math/functions/minmax.hpp diff --git a/include/universal/number/qd/math/next.hpp b/include/universal/number/qd/math/functions/next.hpp similarity index 100% rename from include/universal/number/qd/math/next.hpp rename to include/universal/number/qd/math/functions/next.hpp diff --git a/include/universal/number/qd/math/numerics.hpp b/include/universal/number/qd/math/functions/numerics.hpp similarity index 100% rename from include/universal/number/qd/math/numerics.hpp rename to include/universal/number/qd/math/functions/numerics.hpp diff --git a/include/universal/number/qd/math/pow.hpp b/include/universal/number/qd/math/functions/pow.hpp similarity index 100% rename from include/universal/number/qd/math/pow.hpp rename to include/universal/number/qd/math/functions/pow.hpp diff --git a/include/universal/number/qd/math/sqrt.hpp b/include/universal/number/qd/math/functions/sqrt.hpp similarity index 100% rename from include/universal/number/qd/math/sqrt.hpp rename to include/universal/number/qd/math/functions/sqrt.hpp diff --git a/include/universal/number/qd/math/trigonometry.hpp b/include/universal/number/qd/math/functions/trigonometry.hpp similarity index 96% rename from include/universal/number/qd/math/trigonometry.hpp rename to include/universal/number/qd/math/functions/trigonometry.hpp index 712f0e6f4..fb032f8c2 100644 --- a/include/universal/number/qd/math/trigonometry.hpp +++ b/include/universal/number/qd/math/functions/trigonometry.hpp @@ -7,7 +7,7 @@ // SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. -#include +#include namespace sw { namespace universal { @@ -160,8 +160,8 @@ inline qd sin(const qd& a) { qd r = a - qd_2pi * z; // approximately reduce modulo pi/2 and then modulo pi/1024 - double q = std::floor(r[0] / qd_pi2[0] + 0.5); - qd t = r - qd_pi2 * q; + double q = std::floor(r[0] / qd_pi_2[0] + 0.5); + qd t = r - qd_pi_2 * q; int j = static_cast(q); q = std::floor(t[0] / qd_pi1024[0] + 0.5); t -= qd_pi1024 * q; @@ -240,8 +240,8 @@ inline qd cos(const qd& a) { qd r = a - qd_2pi * z; // approximately reduce modulo pi/2 and then modulo pi/1024 - double q = std::floor(r[0] / qd_pi2[0] + 0.5); - qd t = r - qd_pi2 * q; + double q = std::floor(r[0] / qd_pi_2[0] + 0.5); + qd t = r - qd_pi_2 * q; int j = static_cast(q); q = std::floor(t[0] / qd_pi1024[0] + 0.5); t -= qd_pi1024 * q; @@ -326,8 +326,8 @@ inline void sincos(const qd& a, qd& sin_a, qd& cos_a) { qd t = a - qd_2pi * z; // approximately reduce by pi/2 and then by pi/1024. - double q = std::floor(t[0] / qd_pi2[0] + 0.5); - t -= qd_pi2 * q; + double q = std::floor(t[0] / qd_pi_2[0] + 0.5); + t -= qd_pi_2 * q; int j = static_cast(q); q = std::floor(t[0] / qd_pi1024[0] + 0.5); t -= qd_pi1024 * q; @@ -440,18 +440,18 @@ inline qd atan2(const qd& y, const qd& x) { return qd(SpecificValue::snan); } - return (y.ispos()) ? qd_pi2 : -qd_pi2; + return (y.ispos()) ? qd_pi_2 : -qd_pi_2; } else if (y.iszero()) { return (x.ispos()) ? qd(0.0) : qd_pi; } if (x == y) { - return (y.ispos()) ? qd_pi4 : -qd_3pi4; + return (y.ispos()) ? qd_pi_4 : -qd_3pi_4; } if (x == -y) { - return (y.ispos()) ? qd_3pi4 : -qd_pi4; + return (y.ispos()) ? qd_3pi_4 : -qd_pi_4; } qd r = sqrt(sqr(x) + sqr(y)); @@ -503,7 +503,7 @@ inline qd asin(const qd& a) { } if (abs_a.isone()) { - return (a.ispos()) ? qd_pi2 : -qd_pi2; + return (a.ispos()) ? qd_pi_2 : -qd_pi_2; } return atan2(a, sqrt(1.0 - sqr(a))); diff --git a/include/universal/number/qd/math/truncate.hpp b/include/universal/number/qd/math/functions/truncate.hpp similarity index 100% rename from include/universal/number/qd/math/truncate.hpp rename to include/universal/number/qd/math/functions/truncate.hpp diff --git a/include/universal/number/qd/mathext/horners.hpp b/include/universal/number/qd/mathext/horners.hpp index 3857d2ebf..d89f753e0 100644 --- a/include/universal/number/qd/mathext/horners.hpp +++ b/include/universal/number/qd/mathext/horners.hpp @@ -20,8 +20,8 @@ namespace sw { namespace universal { /// polynomial at x inline qd polyeval(const std::vector& coefficients, int n, const qd& x) { // Horner's method of polynomial evaluation - assert(coefficients.size() > n); - qd r{ coefficients[n] }; + assert(coefficients.size() > static_cast(n)); + qd r{ coefficients[static_cast(n)] }; for (int i = n - 1; i >= 0; --i) { r *= x; @@ -48,9 +48,9 @@ namespace sw { namespace universal { // Compute the coefficients of the derivatives std::vector derivatives{ c }; for (int i = 1; i <= n; ++i) { - double v = std::abs(double(c[i])); + double v = std::abs(double(c[static_cast(i)])); if (v > max_c) max_c = v; - derivatives[i - 1] = c[i] * static_cast(i); + derivatives[i - 1] = c[static_cast(i)] * static_cast(i); } threshold *= max_c; diff --git a/include/universal/number/qd/mathlib.hpp b/include/universal/number/qd/mathlib.hpp index cb923ea57..90f9134a3 100644 --- a/include/universal/number/qd/mathlib.hpp +++ b/include/universal/number/qd/mathlib.hpp @@ -6,19 +6,18 @@ // // This file is part of the universal numbers project, which is released under an MIT Open Source license. -#include +#include -#include -//#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/include/universal/number/qd/qd.hpp b/include/universal/number/qd/qd.hpp index aa1654e93..12d4fa6bc 100644 --- a/include/universal/number/qd/qd.hpp +++ b/include/universal/number/qd/qd.hpp @@ -79,6 +79,7 @@ /////////////////////////////////////////////////////////////////////////////////////// /// elementary math functions library +#include #include #include diff --git a/include/universal/number/qd/qd_impl.hpp b/include/universal/number/qd/qd_impl.hpp index 42c690526..79cfa2ae2 100644 --- a/include/universal/number/qd/qd_impl.hpp +++ b/include/universal/number/qd/qd_impl.hpp @@ -1133,37 +1133,16 @@ class qd { // precomputed quad-double constants courtesy of constants example program -// Golden ratio PHI -constexpr qd qd_phi (1.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975020988e-50); -constexpr qd qd_inv_phi(0.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975021111e-50); -// Euler's number e -constexpr qd qd_e (2.7182818284590451, 1.4456468917292502e-16, -2.1277171080381768e-33, 1.515630159841219e-49); -constexpr qd qd_inv_e (0.36787944117144233, -1.2428753672788363e-17, -5.830044851072742e-34, -2.8267977849017436e-50); - -// pi multiples and fractions -constexpr qd qd_2pi (6.2831853071795862, 2.4492935982947064e-16, -5.9895396194366793e-33, 2.2249084417267313e-49); -constexpr qd qd_pi (3.1415926535897931, 1.2246467991473532e-16, -2.9947698097183397e-33, 1.1124542208633657e-49); -constexpr qd qd_pi2 (1.5707963267948966, 6.123233995736766e-17, -1.4973849048591698e-33, 5.5622711043168283e-50); -constexpr qd qd_pi4 (0.78539816339744828, 3.061616997868383e-17, -7.4869245242958492e-34, 2.7811355521584142e-50); -constexpr qd qd_3pi4 (2.3561944901923448, 9.1848509936051484e-17, 3.9168984647504003e-33, -2.586798163270486e-49); -constexpr qd qd_inv_pi (0.31830988618379069, -1.9678676675182486e-17, -1.0721436282893004e-33, 8.053563926594112e-50); -constexpr qd qd_inv_pi2(0.63661977236758138, -3.9357353350364972e-17, -2.1442872565786008e-33, 1.6107127853188224e-49); - -// natural logarithm (base = e) -constexpr qd qd_ln2 (0.69314718055994529, 2.3190468138462996e-17, 5.7077084384162121e-34, -3.5824322106018105e-50); -constexpr qd qd_lne (1.0, 0.0, 0.0, 0.0); -constexpr qd qd_ln10 (2.3025850929940459, -2.1707562233822494e-16, -9.9842624544657766e-33, -4.0233574544502071e-49); -// binary logarithm (base = 2) -constexpr qd qd_lg2 (1.0, 0.0, 0.0, 0.0); -constexpr qd qd_lge (1.4426950408889634, 2.0355273740931033e-17, -1.0614659956117258e-33, -1.3836716780181395e-50); -constexpr qd qd_lg10 (3.3219280948873622, 1.661617516973592e-16, 1.2215512178458181e-32, 5.9551189702782481e-49); -// common logarithm (base = 10) -constexpr qd qd_log2 (0.3010299956639812, -2.8037281277851704e-18, 5.4719484023146385e-35, 5.1051389831070996e-51); -constexpr qd qd_loge (0.43429448190325182, 1.0983196502167651e-17, 3.717181233110959e-34, 7.7344843465042927e-51); -constexpr qd qd_log10 (1.0, 0.0, 0.0, 0.0); - -constexpr qd qd_sqrt2 (1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683538e-50); -constexpr qd qd_inv_sqrt2(0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341745e-50); + + + + + + + + + + constexpr qd qd_max(1.79769313486231570815e+308, 9.97920154767359795037e+291); diff --git a/static/qd/api/constants.cpp b/static/qd/api/constants.cpp index 74f801401..4d82f7616 100644 --- a/static/qd/api/constants.cpp +++ b/static/qd/api/constants.cpp @@ -284,53 +284,53 @@ try { qd value; } constant_symbol_table[] = { { "qd_phi" , "1.6180339887498948482045868343656381177203091798057628621354486228", qd_phi }, - { "qd_inv_phi", "0.6180339887498948482045868343656381177203091798057628621354486227", qd_inv_phi }, + { "qd_1_phi", "0.6180339887498948482045868343656381177203091798057628621354486227", qd_1_phi }, { "qd_e" , "2.7182818284590452353602874713526624977572470936999595749669676277", qd_e }, - { "qd_inv_e" , "0.3678794411714423215955237701614608674458111310317678345078368017", qd_inv_e }, + { "qd_1_e" , "0.3678794411714423215955237701614608674458111310317678345078368017", qd_1_e }, { "qd_2pi" , "6.2831853071795864769252867665590057683943387987502116419498891847", qd_2pi }, { "qd_pi" , "3.1415926535897932384626433832795028841971693993751058209749445923", qd_pi }, - { "qd_pi2" , "1.5707963267948966192313216916397514420985846996875529104874722962", qd_pi2 }, - { "qd_pi4" , "0.7853981633974483096156608458198757210492923498437764552437361481", qd_pi4 }, - { "qd_3pi4" , "2.3561944901923449288469825374596271631478770495313293657312084443", qd_3pi4 }, + { "qd_pi2" , "1.5707963267948966192313216916397514420985846996875529104874722962", qd_pi_2 }, + { "qd_pi4" , "0.7853981633974483096156608458198757210492923498437764552437361481", qd_pi_4 }, + { "qd_3pi4" , "2.3561944901923449288469825374596271631478770495313293657312084443", qd_3pi_4 }, - { "qd_inv_pi" , "0.31830988618379067153776752674502872406891929148091289749533468812", qd_inv_pi }, - { "qd_inv_pi2", "0.63661977236758134307553505349005744813783858296182579499066937624", qd_inv_pi2 }, + { "qd_1_pi" , "0.31830988618379067153776752674502872406891929148091289749533468812", qd_1_pi }, + { "qd_2_pi" , "0.63661977236758134307553505349005744813783858296182579499066937624", qd_2_pi }, { "qd_ln2" , "0.69314718055994530941723212145817656807550013436025525412068000950", qd_ln2 }, - { "qd_lne" , "1.00000000000000000000000000000000000000000000000000000000000000000", qd_lne }, + { "qd_lne" , "1.00000000000000000000000000000000000000000000000000000000000000000", qd(1.0)}, { "qd_ln10" , "2.30258509299404568401799145468436420760110148862877297603332790097", qd_ln10 }, - { "qd_lg2" , "1.0000000000000000000000000000000000000000000000000000000000000000", qd_lg2 }, + { "qd_lg2" , "1.0000000000000000000000000000000000000000000000000000000000000000", qd(1.0)}, { "qd_lge" , "1.4426950408889634073599246810018921374266459541529859341354494069", qd_lge }, { "qd_lg10" , "3.3219280948873623478703194294893901758648313930245806120547563956", qd_lg10 }, { "qd_log2" , "3.0102999566398119521373889472449302676818988146210854131042746113e-01", qd_log2 }, { "qd_loge" , "4.3429448190325182765112891891660508229439700580366656611445378316e-01", qd_loge }, - { "qd_log10" , "1.00000000000000000000000000000000000000000000000000000000000000000", qd_log10 }, + { "qd_log10" , "1.0000000000000000000000000000000000000000000000000000000000000000", qd(1.0)}, { "qd_sqrt2" , "1.4142135623730950488016887242096980785696718753769480731766797380", qd_sqrt2 }, - { "qd_inv_sqrt2", "7.0710678118654752440084436210484903928483593768847403658833986899e-01", qd_inv_sqrt2 }, + { "qd_1_sqrt2", "7.0710678118654752440084436210484903928483593768847403658833986899e-01", qd_1_sqrt2 }, }; /* * - * ETLO August 14, 2024 + * ETLO August 31, 2024 * Need to verify if these are the most accurate quad-double approximations available. - * + * This is Debug, Release cuts the precision in half verifying constants qd_phi : 1.61803398874989484820458683436564e+00 vs 1.61803398874989484820458683436564e+00 : ( 1.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975020983e-50) : 4.74778387287989937373662113478098e-66 -qd_inv_phi : 6.18033988749894848204586834365638e-01 vs 6.18033988749894848204586834365638e-01 : ( 0.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975021083e-50) : 2.84867032372793962424197268086859e-65 +qd_1_phi : 6.18033988749894848204586834365638e-01 vs 6.18033988749894848204586834365638e-01 : ( 0.6180339887498949, -5.4321152036825061e-17, 2.6543252083815655e-33, -3.3049919975021083e-50) : 2.84867032372793962424197268086859e-65 qd_e : 2.71828182845904523536028747135266e+00 vs 2.71828182845904523536028747135266e+00 : ( 2.7182818284590451, 1.4456468917292502e-16, -2.1277171080381768e-33, 1.5156301598412188e-49) : -1.89911354915195974949464845391239e-65 -qd_inv_e : 3.67879441171442321595523770161461e-01 vs 3.67879441171442321595523770161461e-01 : ( 0.36787944117144233, -1.2428753672788363e-17, -5.830044851072742e-34, -2.8267977849017436e-50) : 0.00000000000000000000000000000000e+00 +qd_1_e : 3.67879441171442321595523770161461e-01 vs 3.67879441171442321595523770161461e-01 : ( 0.36787944117144233, -1.2428753672788363e-17, -5.830044851072742e-34, -2.8267977849017436e-50) : 0.00000000000000000000000000000000e+00 qd_2pi : 6.28318530717958647692528676655901e+00 vs 6.28318530717958647692528676655901e+00 : ( 6.2831853071795862, 2.4492935982947064e-16, -5.9895396194366793e-33, 2.2249084417267317e-49) : 3.79822709830391949898929690782478e-65 qd_pi : 3.14159265358979323846264338327950e+00 vs 3.14159265358979323846264338327950e+00 : ( 3.1415926535897931, 1.2246467991473532e-16, -2.9947698097183397e-33, 1.1124542208633653e-49) : -3.79822709830391949898929690782478e-65 qd_pi2 : 1.57079632679489661923132169163975e+00 vs 1.57079632679489661923132169163975e+00 : ( 1.5707963267948966, 6.123233995736766e-17, -1.4973849048591698e-33, 5.5622711043168312e-50) : 2.84867032372793962424197268086859e-65 qd_pi4 : 7.85398163397448309615660845819876e-01 vs 7.85398163397448309615660845819876e-01 : ( 0.78539816339744828, 3.061616997868383e-17, -7.4869245242958492e-34, 2.7811355521584156e-50) : 1.42433516186396981212098634043429e-65 qd_3pi4 : 2.35619449019234492884698253745963e+00 vs 2.35619449019234492884698253745963e+00 : ( 2.3561944901923448, 9.1848509936051484e-17, 3.9168984647504003e-33, -2.5867981632704857e-49) : 3.79822709830391949898929690782478e-65 -qd_inv_pi : 3.18309886183790671537767526745029e-01 vs 3.18309886183790671537767526745029e-01 : ( 0.31830988618379069, -1.9678676675182486e-17, -1.0721436282893004e-33, 8.053563926594112e-50) : 0.00000000000000000000000000000000e+00 -qd_inv_pi2 : 6.36619772367581343075535053490057e-01 vs 6.36619772367581343075535053490057e-01 : ( 0.63661977236758138, -3.9357353350364972e-17, -2.1442872565786008e-33, 1.6107127853188224e-49) : 0.00000000000000000000000000000000e+00 +qd_1_pi : 3.18309886183790671537767526745029e-01 vs 3.18309886183790671537767526745029e-01 : ( 0.31830988618379069, -1.9678676675182486e-17, -1.0721436282893004e-33, 8.053563926594112e-50) : 0.00000000000000000000000000000000e+00 +qd_2_pi : 6.36619772367581343075535053490057e-01 vs 6.36619772367581343075535053490057e-01 : ( 0.63661977236758138, -3.9357353350364972e-17, -2.1442872565786008e-33, 1.6107127853188224e-49) : 0.00000000000000000000000000000000e+00 qd_ln2 : 6.93147180559945309417232121458177e-01 vs 6.93147180559945309417232121458177e-01 : ( 0.69314718055994529, 2.3190468138462996e-17, 5.7077084384162121e-34, -3.5824322106018109e-50) : -4.74778387287989937373662113478098e-66 qd_lne : 1.00000000000000000000000000000000e+00 vs 1.00000000000000000000000000000000e+00 : ( 1, 0, 0, 0) : 0.00000000000000000000000000000000e+00 qd_ln10 : 2.30258509299404568401799145468436e+00 vs 2.30258509299404568401799145468436e+00 : ( 2.3025850929940459, -2.1707562233822494e-16, -9.9842624544657766e-33, -4.0233574544502064e-49) : 7.59645419660783899797859381564957e-65 @@ -340,8 +340,8 @@ qd_lg10 : 3.32192809488736234787031942948939e+00 vs 3.321928094887362347 qd_log2 : 3.01029995663981195213738894724493e-01 vs 3.01029995663981195213738894724493e-01 : ( 0.3010299956639812, -2.8037281277851704e-18, 5.4719484023146385e-35, 5.1051389831070954e-51) : -4.15431088876991195201954349293336e-66 qd_loge : 4.34294481903251827651128918916605e-01 vs 4.34294481903251827651128918916605e-01 : ( 0.43429448190325182, 1.0983196502167651e-17, 3.717181233110959e-34, 7.7344843465042927e-51) : 0.00000000000000000000000000000000e+00 qd_log10 : 1.00000000000000000000000000000000e+00 vs 1.00000000000000000000000000000000e+00 : ( 1, 0, 0, 0) : 0.00000000000000000000000000000000e+00 -qd_sqrt2 : 1.41421356237309504880168872420970e+00 vs 1.41421356237309504880168872420970e+00 : ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683519e-50) : -1.89911354915195974949464845391239e-65 -qd_inv_sqrt2 : 7.07106781186547524400844362104849e-01 vs 7.07106781186547524400844362104849e-01 : ( 0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341759e-50) : 1.42433516186396981212098634043429e-65 +qd_sqrt2 : 1.41421356237309504880168872420970e+00 vs 1.41421356237309504880168872420970e+00 : ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683519e-50) : 9.49556774575979874747324226956196e-66 +qd_1_sqrt2 : 7.07106781186547524400844362104849e-01 vs 7.07106781186547524400844362104849e-01 : ( 0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341759e-50) : 4.74778387287989937373662113478098e-66 */ auto oldPrec = std::cout.precision(); std::cout << std::setprecision(32); @@ -350,6 +350,22 @@ qd_inv_sqrt2 : 7.07106781186547524400844362104849e-01 vs 7.071067811865475244 qd error = (c - e.value); std::cout << std::left << std::setw(15) << e.name << " : " << c << " vs " << e.value << " : " << to_quad(c) << " : " << error << '\n'; } + + { + + qd a{ 2.0 }; + std::cout << "sqrt(2.0) " << to_quad(sqrt(a)) << '\n'; + std::cout << "1/sqrt(2.0) " << to_quad(reciprocal(sqrt(a))) << '\n'; + a = 3.0; + std::cout << "sqrt(3.0) " << to_quad(sqrt(a)) << '\n'; + a = 5.0; + std::cout << "sqrt(5.0) " << to_quad(sqrt(a)) << '\n'; + + a = sqrt(qd_pi); + std::cout << "2 / sqrtpi " << to_quad(2.0 / a) << '\n'; + } + + std::cout << std::setprecision(oldPrec); ReportTestSuiteResults(test_suite, nrOfFailedTestCases); diff --git a/static/qd/math/hyperbolic.cpp b/static/qd/math/hyperbolic.cpp index cbdd4640f..2bbacaf5d 100644 --- a/static/qd/math/hyperbolic.cpp +++ b/static/qd/math/hyperbolic.cpp @@ -55,7 +55,7 @@ try { { std::cout << "quad-double reference\n"; - qd x = qd_pi4; + qd x = qd_pi_4; std::cout << "sinh( " << x << " ) = " << sinh(x) << '\n'; std::cout << "cosh( " << x << " ) = " << cosh(x) << '\n'; std::cout << "tanh( " << x << " ) = " << tanh(x) << '\n'; diff --git a/static/qd/math/trigonometry.cpp b/static/qd/math/trigonometry.cpp index 440c9c1ff..8608ffad6 100644 --- a/static/qd/math/trigonometry.cpp +++ b/static/qd/math/trigonometry.cpp @@ -9,9 +9,8 @@ #include #include -#include -#include -#include +// for the constant d_2pi used in the verification functions +#include template int VerifySinFunction(bool reportTestCases) { @@ -19,16 +18,16 @@ int VerifySinFunction(bool reportTestCases) { constexpr bool bTraceError{ false }; int nrOfFailedTestCases{ 0 }; - const double d2pi = 6.283185307179586476925286766559; - //const double piOver4 = 0.78539816339744830961566084581988; - //const double piOver8 = 0.39269908169872415480783042290994; - //const double piOver16 = 0.19634954084936207740391521145497; - const double piOver32 = 0.01227184630308512983774470071594; + //const double d2pi = 6.283185307179586476925286766559; + //const double pi_4 = 0.78539816339744830961566084581988; + //const double pi_8 = 0.39269908169872415480783042290994; + //const double pi_16 = 0.19634954084936207740391521145497; + const double pi_32 = 0.01227184630308512983774470071594; // walk the unit circle in steps of pi/32 - double dinc{ piOver32 }; - unsigned samples{ static_cast(d2pi / dinc) }; - Real increment{ piOver32 }; + double dinc{ pi_32 }; + unsigned samples{ static_cast(sw::universal::d_2pi / dinc) }; + Real increment{ pi_32 }; for (unsigned i = 0; i < samples; ++i) { Real angle = Real(i) * increment; double dangle = double(i) * dinc; @@ -53,16 +52,16 @@ int VerifyCosFunction(bool reportTestCases) { constexpr bool bTraceError{ false }; int nrOfFailedTestCases{ 0 }; - const double d2pi = 6.283185307179586476925286766559; - //const double piOver4 = 0.78539816339744830961566084581988; - //const double piOver8 = 0.39269908169872415480783042290994; - //const double piOver16 = 0.19634954084936207740391521145497; - const double piOver32 = 0.01227184630308512983774470071594; + //const double d2pi = 6.283185307179586476925286766559; + //const double pi_4 = 0.78539816339744830961566084581988; + //const double pi_8 = 0.39269908169872415480783042290994; + //const double pi_16 = 0.19634954084936207740391521145497; + const double pi_32 = 0.01227184630308512983774470071594; // walk the unit circle in steps of pi/32 - double dinc{ piOver32 }; - unsigned samples{ static_cast(d2pi / dinc) }; - Real increment{ piOver32 }; + double dinc{ pi_32 }; + unsigned samples{ static_cast(sw::universal::d_2pi / dinc) }; + Real increment{ pi_32 }; for (unsigned i = 0; i < samples; ++i) { Real angle = Real(i) * increment; double dangle = double(i) * dinc; @@ -87,17 +86,17 @@ int VerifyTanFunction(bool reportTestCases) { constexpr bool bTraceError{ false }; int nrOfFailedTestCases{ 0 }; - const double d2pi = 6.283185307179586476925286766559; - //const double piOver2 = 1.5707963267948966192313216916398; - //const double piOver4 = 0.78539816339744830961566084581988; - //const double piOver8 = 0.39269908169872415480783042290994; - //const double piOver16 = 0.19634954084936207740391521145497; - const double piOver32 = 0.01227184630308512983774470071594; + //const double d2pi = 6.283185307179586476925286766559; + //const double pi_2 = 1.5707963267948966192313216916398; + //const double pi_4 = 0.78539816339744830961566084581988; + //const double pi_8 = 0.39269908169872415480783042290994; + //const double pi_16 = 0.19634954084936207740391521145497; + const double pi_32 = 0.01227184630308512983774470071594; // walk the unit circle in steps of pi/32 - double dinc{ piOver32 }; - unsigned samples{ static_cast(d2pi / dinc) }; - Real increment{ piOver32 }; + double dinc{ pi_32 }; + unsigned samples{ static_cast(sw::universal::d_2pi / dinc) }; + Real increment{ pi_32 }; // tan(x) is inf at pi/2 and 3pi/4 // they are at 1/4 and 3/4s of the sample sequence for (unsigned i = 0; i < samples; ++i) { @@ -190,13 +189,13 @@ int VerifyArctanFunction(bool reportTestCases) { // walk the domain of arctan = [ -inf, inf ] to the range of [ -pi/2, pi/2 ] // we are going to use tan(x) to generate the values to inverse - const double d2pi = 6.283185307179586476925286766559; - const double piOver32 = 0.01227184630308512983774470071594; + //const double d2pi = 6.283185307179586476925286766559; + const double pi_32 = 0.01227184630308512983774470071594; // walk the unit circle in steps of pi/32 - double dinc{ piOver32 }; - unsigned samples{ static_cast(d2pi / dinc) }; - Real increment{ piOver32 }; + double dinc{ pi_32 }; + unsigned samples{ static_cast(sw::universal::d_2pi / dinc) }; + Real increment{ pi_32 }; // tan(x) is inf at pi/2 and 3pi/4 // they are at 1/4 and 3/4s of the sample sequence for (unsigned i = 0; i < samples; ++i) { @@ -223,7 +222,7 @@ int VerifyArctanFunction(bool reportTestCases) { } // Regression testing guards: typically set by the cmake configuration, but MANUAL_TESTING is an override -#define MANUAL_TESTING 1 +#define MANUAL_TESTING 0 // 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 @@ -270,17 +269,17 @@ try { VerifySinFunction(reportTestCases); - qd piOver4("0.78539816339744830961566084581988"); - qd piOver8("0.39269908169872415480783042290994"); - qd piOver16("0.19634954084936207740391521145497"); - qd piOver32("0.01227184630308512983774470071594"); + qd pi_4("0.78539816339744830961566084581988"); + qd pi_8("0.39269908169872415480783042290994"); + qd pi_16("0.19634954084936207740391521145497"); + qd pi_32("0.01227184630308512983774470071594"); - qd a = sin(piOver4); + qd a = sin(pi_4); - std::cout << "pi/4 : " << std::setprecision(32) << piOver4 << '\n'; - std::cout << "pi/8 : " << std::setprecision(32) << piOver8 << '\n'; - std::cout << "pi/16 : " << std::setprecision(32) << piOver16 << '\n'; - std::cout << "pi/32 : " << std::setprecision(32) << piOver32 << '\n'; + std::cout << "pi/4 : " << std::setprecision(32) << pi_4 << '\n'; + std::cout << "pi/8 : " << std::setprecision(32) << pi_8 << '\n'; + std::cout << "pi/16 : " << std::setprecision(32) << pi_16 << '\n'; + std::cout << "pi/32 : " << std::setprecision(32) << pi_32 << '\n'; qd b{}; b = asin(qd(0));