completing generation of quad-double constants and math library struc…

…ture

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

include/universal/number/qd/math/complex.hpp → 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

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

include/universal/number/qd/math/logarithm.hpp → 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;
 	}
 
 	/// <summary>

include/universal/number/qd/math/trigonometry.hpp → 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 <universal/number/qd/math/sincos_tables.hpp>
+#include <universal/number/qd/math/constants/sincos_tables.hpp>
 
 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<int>(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<int>(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<int>(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)));

include/universal/number/qd/mathext/horners.hpp

@@ -20,8 +20,8 @@ namespace sw { namespace universal {
     /// <returns>polynomial at x</returns>
     inline qd polyeval(const std::vector<qd>& 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<unsigned>(n));
+        qd r{ coefficients[static_cast<unsigned>(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<qd> derivatives{ c };
         for (int i = 1; i <= n; ++i) {
-            double v = std::abs(double(c[i]));
+            double v = std::abs(double(c[static_cast<unsigned>(i)]));
             if (v > max_c) max_c = v;
-            derivatives[i - 1] = c[i] * static_cast<double>(i);
+            derivatives[i - 1] = c[static_cast<unsigned>(i)] * static_cast<double>(i);
         }
         threshold *= max_c;
 

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 <universal/number/qd/math/numerics.hpp>
+#include <universal/number/qd/math/functions/numerics.hpp>
 
-#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/fractional.hpp>
-#include <universal/number/qd/math/hyperbolic.hpp>
-#include <universal/number/qd/math/hypot.hpp>
-#include <universal/number/qd/math/logarithm.hpp>
-#include <universal/number/qd/math/minmax.hpp>
-#include <universal/number/qd/math/next.hpp>
-#include <universal/number/qd/math/pow.hpp>
-#include <universal/number/qd/math/sqrt.hpp>
-#include <universal/number/qd/math/trigonometry.hpp>
-#include <universal/number/qd/math/truncate.hpp>
+#include <universal/number/qd/math/functions/classify.hpp>
+#include <universal/number/qd/math/functions/error_and_gamma.hpp>
+#include <universal/number/qd/math/functions/exponent.hpp>
+#include <universal/number/qd/math/functions/fractional.hpp>
+#include <universal/number/qd/math/functions/hyperbolic.hpp>
+#include <universal/number/qd/math/functions/hypot.hpp>
+#include <universal/number/qd/math/functions/logarithm.hpp>
+#include <universal/number/qd/math/functions/minmax.hpp>
+#include <universal/number/qd/math/functions/next.hpp>
+#include <universal/number/qd/math/functions/pow.hpp>
+#include <universal/number/qd/math/functions/sqrt.hpp>
+#include <universal/number/qd/math/functions/trigonometry.hpp>
+#include <universal/number/qd/math/functions/truncate.hpp>

include/universal/number/qd/qd.hpp

@@ -79,6 +79,7 @@
 
 ///////////////////////////////////////////////////////////////////////////////////////
 /// elementary math functions library
+#include <universal/number/qd/math/constants/qd_constants.hpp>
 #include <universal/number/qd/mathlib.hpp>
 #include <universal/number/qd/mathext.hpp>
 

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