From 94481da233a7856acd7feb9cf8608e29267181f4 Mon Sep 17 00:00:00 2001 From: Ravenwater Date: Wed, 28 Aug 2024 16:46:25 -0400 Subject: [PATCH] resolving all the missing inline keywords on the dd and qd math functions --- CMakeLists.txt | 15 +- include/universal/number/dd/attributes.hpp | 8 +- include/universal/number/dd/dd_impl.hpp | 32 ++-- include/universal/number/dd/manipulators.hpp | 4 +- include/universal/number/dd/math/classify.hpp | 2 +- .../number/dd/math/error_and_gamma.hpp | 4 +- include/universal/number/dd/math/exponent.hpp | 50 +++--- .../universal/number/dd/math/fractional.hpp | 4 +- .../universal/number/dd/math/hyperbolic.hpp | 12 +- include/universal/number/dd/math/hypot.hpp | 2 +- .../universal/number/dd/math/logarithm.hpp | 167 +++++++++--------- include/universal/number/dd/math/minmax.hpp | 4 +- include/universal/number/dd/math/next.hpp | 82 ++++----- include/universal/number/dd/math/numerics.hpp | 6 +- include/universal/number/dd/math/pow.hpp | 11 +- include/universal/number/dd/math/sqrt.hpp | 4 +- .../universal/number/dd/math/trigonometry.hpp | 8 +- include/universal/number/dd/math/truncate.hpp | 4 +- include/universal/number/qd/attributes.hpp | 8 +- include/universal/number/qd/math/classify.hpp | 84 ++++----- .../number/qd/math/error_and_gamma.hpp | 4 +- include/universal/number/qd/math/exponent.hpp | 62 +++---- .../universal/number/qd/math/fractional.hpp | 4 +- .../universal/number/qd/math/hyperbolic.hpp | 12 +- include/universal/number/qd/math/hypot.hpp | 2 +- .../universal/number/qd/math/logarithm.hpp | 16 +- include/universal/number/qd/math/minmax.hpp | 4 +- include/universal/number/qd/math/next.hpp | 70 ++------ include/universal/number/qd/math/numerics.hpp | 6 +- include/universal/number/qd/math/pow.hpp | 10 +- include/universal/number/qd/math/sqrt.hpp | 2 +- .../universal/number/qd/math/trigonometry.hpp | 35 ++-- include/universal/number/qd/math/truncate.hpp | 4 +- include/universal/number/qd/qd_impl.hpp | 6 +- static/appenv/multifile/areals.cpp | 3 +- static/appenv/multifile/cfloats.cpp | 3 +- static/appenv/multifile/dbns.cpp | 3 +- static/appenv/multifile/dd.cpp | 25 +++ static/appenv/multifile/fixpnts.cpp | 3 +- static/appenv/multifile/integers.cpp | 5 +- static/appenv/multifile/logs.cpp | 3 +- static/appenv/multifile/main.cpp | 44 +++-- static/appenv/multifile/posits.cpp | 3 +- static/appenv/multifile/qd.cpp | 25 +++ static/appenv/multifile/unums.cpp | 3 +- 45 files changed, 457 insertions(+), 411 deletions(-) create mode 100644 static/appenv/multifile/dd.cpp create mode 100644 static/appenv/multifile/qd.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ea1082a90..a29c117e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -622,13 +622,16 @@ if(BUILD_ALL) set(BUILD_VALIDATION_HW ON) endif(BUILD_ALL) -# set the grouped components to build (will trigger builds when tested) +# set the grouped components to build Continuous Integration regression suite +# we are disabling COMPLEX because we need a solution to the +# disappearing support of complex<> for user-defined types if(BUILD_CI) set(BUILD_NUMBERS ON) set(BUILD_DEMONSTRATION OFF) set(BUILD_NUMERICS OFF) set(BUILD_MIXEDPRECISION_SDK OFF) set(BUILD_APP_ENVIRONMENT ON) + set(BUILD_COMPLEX OFF) endif(BUILD_CI) # core demonstration example applications that use the library @@ -639,8 +642,16 @@ if(BUILD_DEMONSTRATION) set(BUILD_PLAYGROUND ON) endif(BUILD_DEMONSTRATION) +# enable complex environment components if(BUILD_COMPLEX) - +# right now, the complex arithmetic and math functions (imag, real, conj) +# are maintained in the individual number system regression suites. +# Their respective build systems will use BUILD_COMPLEX to add/remove +# the complex<> regression suites. + +# I am leaving this segment here as a pattern that might drive the development +# of a complex<> replacement that might need its own regression suite +# independent of user-defined types. endif(BUILD_COMPLEX) if(BUILD_NUMBERS) diff --git a/include/universal/number/dd/attributes.hpp b/include/universal/number/dd/attributes.hpp index 4973332d6..384a578d4 100644 --- a/include/universal/number/dd/attributes.hpp +++ b/include/universal/number/dd/attributes.hpp @@ -18,7 +18,7 @@ namespace sw { namespace universal { } // generate the maxneg through maxpos value range of a double-double configuration - std::string dd_range() { + inline std::string dd_range() { dd v; std::stringstream s; s << std::setw(80) << type_tag(v) << " : [ " @@ -47,17 +47,17 @@ namespace sw { namespace universal { } */ - int minpos_scale(const dd& b) { + inline int minpos_scale(const dd& b) { dd c(b); return c.minpos().scale(); } - int maxpos_scale(const dd& b) { + inline int maxpos_scale(const dd& b) { dd c(b); return c.maxpos().scale(); } - int max_negative_scale(const dd& b) { + inline int max_negative_scale(const dd& b) { dd c(b); return c.maxneg().scale(); } diff --git a/include/universal/number/dd/dd_impl.hpp b/include/universal/number/dd/dd_impl.hpp index 3ebdc16d6..4c1f91641 100644 --- a/include/universal/number/dd/dd_impl.hpp +++ b/include/universal/number/dd/dd_impl.hpp @@ -34,7 +34,7 @@ namespace sw { namespace universal { // this is debug infrastructure: TODO: remove when decimal conversion is solved reliably constexpr bool bTraceDecimalConversion = false; constexpr bool bTraceDecimalRounding = false; - std::ostream& operator<<(std::ostream& ostr, const std::vector& s) { + inline std::ostream& operator<<(std::ostream& ostr, const std::vector& s) { for (auto c : s) { ostr << c; } @@ -42,14 +42,14 @@ namespace sw { namespace universal { } // fwd references to free functions -dd operator-(const dd&, const dd&); -dd operator*(const dd&, const dd&); -dd operator/(const dd&, const dd&); -std::ostream& operator<<(std::ostream&, const dd&); -bool signbit(const dd&); -dd pown(const dd&, int); -dd frexp(const dd&, int*); -dd ldexp(const dd&, int); + inline dd operator-(const dd&, const dd&); +inline dd operator*(const dd&, const dd&); +inline dd operator/(const dd&, const dd&); +inline std::ostream& operator<<(std::ostream&, const dd&); +inline bool signbit(const dd&); +inline dd pown(const dd&, int); +inline dd frexp(const dd&, int*); +inline dd ldexp(const dd&, int); // dd is an unevaluated pair of IEEE-754 doubles that provides a (1,11,106) floating-point triple class dd { @@ -1184,7 +1184,7 @@ inline dd mul_pwr2(const dd& a, double b) { // quad-double operators // quad-double + double-double -void qd_add(double const a[4], const dd& b, double s[4]) { +inline void qd_add(double const a[4], const dd& b, double s[4]) { double t[5]; s[0] = two_sum(a[0], b.high(), t[0]); // s0 - O( 1 ); t0 - O( e ) s[1] = two_sum(a[1], b.low(), t[1]); // s1 - O( e ); t1 - O( e^2 ) @@ -1201,7 +1201,7 @@ void qd_add(double const a[4], const dd& b, double s[4]) { } // quad-double = double-double * double-double -void qd_mul(const dd& a, const dd& b, double p[4]) { +inline void qd_mul(const dd& a, const dd& b, double p[4]) { double p4, p5, p6, p7; // powers of e - 0, 1, 1, 1, 2, 2, 2, 3 @@ -1278,7 +1278,13 @@ inline dd reciprocal(const dd& a) { ///////////////////////////////////////////////////////////////////////////// // power functions -dd cbrt(const dd& a) { +/// +/// cbrt is cube root function, that is, x^1/3 +/// +/// input +/// cube root of a +inline dd cbrt(const dd& a) { + using std::pow; if (!a.isfinite() || a.iszero()) return a; // NaN returns NaN; +/-Inf returns +/-Inf, +/-0.0 returns +/-0.0 @@ -1377,7 +1383,7 @@ inline std::istream& operator>>(std::istream& istr, dd& v) { ////////////////// string operators // parse a decimal ASCII floating-point format and make a doubledouble (dd) out of it -bool parse(const std::string& number, dd& value) { +inline bool parse(const std::string& number, dd& value) { char const* p = number.c_str(); // Skip any leading spaces diff --git a/include/universal/number/dd/manipulators.hpp b/include/universal/number/dd/manipulators.hpp index 4d200c190..75cbeac1f 100644 --- a/include/universal/number/dd/manipulators.hpp +++ b/include/universal/number/dd/manipulators.hpp @@ -14,12 +14,12 @@ namespace sw { namespace universal { // Generate a type tag for a doubledouble - std::string type_tag(const dd& = {}) { + inline std::string type_tag(const dd& = {}) { return std::string("double-double"); } // generate a binary, color-coded representation of the doubledouble - std::string color_print(const dd& r, bool nibbleMarker = false) { + inline std::string color_print(const dd& r, bool nibbleMarker = false) { std::stringstream s; double high = r.high(); double low = r.low(); diff --git a/include/universal/number/dd/math/classify.hpp b/include/universal/number/dd/math/classify.hpp index 3ce43fdcf..94f7336fb 100644 --- a/include/universal/number/dd/math/classify.hpp +++ b/include/universal/number/dd/math/classify.hpp @@ -9,7 +9,7 @@ namespace sw { namespace universal { // STD LIB function for IEEE floats: Categorizes floating point value arg into the following categories: zero, subnormal, normal, infinite, NAN, or implementation-defined category. -int fpclassify(const dd& a) { +inline int fpclassify(const dd& a) { return (std::fpclassify(a.high())); } diff --git a/include/universal/number/dd/math/error_and_gamma.hpp b/include/universal/number/dd/math/error_and_gamma.hpp index 6260e53c4..714b9937b 100644 --- a/include/universal/number/dd/math/error_and_gamma.hpp +++ b/include/universal/number/dd/math/error_and_gamma.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // Compute the error function erf(x) = 2 over sqrt(PI) times Integral from 0 to x of e ^ (-t)^2 dt - dd erf(dd x) { + inline dd erf(dd x) { return dd(std::erf(double(x.high()))); } // Compute the complementary error function: 1 - erf(x) - dd erfc(dd x) { + inline dd erfc(dd x) { return dd(std::erfc(double(x.high()))); } diff --git a/include/universal/number/dd/math/exponent.hpp b/include/universal/number/dd/math/exponent.hpp index db614b6af..9ac4207f4 100644 --- a/include/universal/number/dd/math/exponent.hpp +++ b/include/universal/number/dd/math/exponent.hpp @@ -11,29 +11,29 @@ namespace sw { namespace universal { // fwd reference - dd ldexp(dd const&, int); - -static const int n_inv_fact = 15; -static const double inv_fact[n_inv_fact][2] = { - { 1.66666666666666657e-01, 9.25185853854297066e-18}, // 1/3! - { 4.16666666666666644e-02, 2.31296463463574266e-18}, // 1/4! - { 8.33333333333333322e-03, 1.15648231731787138e-19}, // 1/5! - { 1.38888888888888894e-03, -5.30054395437357706e-20}, // 1/6! - { 1.98412698412698413e-04, 1.72095582934207053e-22}, // 1/7! - { 2.48015873015873016e-05, 2.15119478667758816e-23}, // 1/8! - { 2.75573192239858925e-06, -1.85839327404647208e-22}, // 1/9! - { 2.75573192239858883e-07, 2.37677146222502973e-23}, // 1/10! - { 2.50521083854417202e-08, -1.44881407093591197e-24}, // 1/11! - { 2.08767569878681002e-09, -1.20734505911325997e-25}, // 1/12! - { 1.60590438368216133e-10, 1.25852945887520981e-26}, // 1/13! - { 1.14707455977297245e-11, 2.06555127528307454e-28}, // 1/14! - { 7.64716373181981641e-13, 7.03872877733453001e-30}, // 1/15! - { 4.77947733238738525e-14, 4.39920548583408126e-31}, // 1/16! - { 2.81145725434552060e-15, 1.65088427308614326e-31} // 1/17! + inline dd ldexp(const dd&, int); + +static const int dd_inverse_factorial_table_size = 15; +static const dd dd_inverse_factorial[dd_inverse_factorial_table_size] = { + dd(1.66666666666666657e-01, 9.25185853854297066e-18), // 1/3! + dd(4.16666666666666644e-02, 2.31296463463574266e-18), // 1/4! + dd(8.33333333333333322e-03, 1.15648231731787138e-19), // 1/5! + dd(1.38888888888888894e-03, -5.30054395437357706e-20), // 1/6! + dd(1.98412698412698413e-04, 1.72095582934207053e-22), // 1/7! + dd(2.48015873015873016e-05, 2.15119478667758816e-23), // 1/8! + dd(2.75573192239858925e-06, -1.85839327404647208e-22), // 1/9! + dd(2.75573192239858883e-07, 2.37677146222502973e-23), // 1/10! + dd(2.50521083854417202e-08, -1.44881407093591197e-24), // 1/11! + dd(2.08767569878681002e-09, -1.20734505911325997e-25), // 1/12! + dd(1.60590438368216133e-10, 1.25852945887520981e-26), // 1/13! + dd(1.14707455977297245e-11, 2.06555127528307454e-28), // 1/14! + dd(7.64716373181981641e-13, 7.03872877733453001e-30), // 1/15! + dd(4.77947733238738525e-14, 4.39920548583408126e-31), // 1/16! + dd(2.81145725434552060e-15, 1.65088427308614326e-31) // 1/17! }; // Base-e exponential function -dd exp(const dd& a) { +inline dd exp(const dd& a) { /* Strategy: We first reduce the size of x by noting that exp(kr + m * log(2)) = 2^m * exp(r)^k @@ -61,13 +61,13 @@ dd exp(const dd& a) { p = sqr(r); s = r + mul_pwr2(p, 0.5); p *= r; - t = p * dd(inv_fact[0][0], inv_fact[0][1]); + t = p * dd_inverse_factorial[0]; int i = 0; do { s += t; p *= r; ++i; - t = p * dd(inv_fact[i][0], inv_fact[i][1]); + t = p * dd_inverse_factorial[i]; } while (std::abs(double(t)) > inv_k * dd_eps && i < 5); s += t; @@ -87,17 +87,17 @@ dd exp(const dd& a) { } // Base-2 exponential function -dd exp2(dd x) { +inline dd exp2(dd x) { return dd(std::exp2(double(x))); } // Base-10 exponential function -dd exp10(dd x) { +inline dd exp10(dd x) { return dd(std::pow(10.0, double(x))); } // Base-e exponential function exp(x)-1 -dd expm1(dd x) { +inline dd expm1(dd x) { return dd(std::expm1(double(x))); } diff --git a/include/universal/number/dd/math/fractional.hpp b/include/universal/number/dd/math/fractional.hpp index 15148342e..09bf11224 100644 --- a/include/universal/number/dd/math/fractional.hpp +++ b/include/universal/number/dd/math/fractional.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // fmod retuns x - n*y where n = x/y with the fractional part truncated - dd fmod(dd x, dd y) { + inline dd fmod(dd x, dd y) { return dd(std::fmod(double(x), double(y))); } // shim to stdlib - dd remainder(dd x, dd y) { + inline dd remainder(dd x, dd y) { return dd(std::remainder(double(x), double(y))); } diff --git a/include/universal/number/dd/math/hyperbolic.hpp b/include/universal/number/dd/math/hyperbolic.hpp index 1793be766..c922fc112 100644 --- a/include/universal/number/dd/math/hyperbolic.hpp +++ b/include/universal/number/dd/math/hyperbolic.hpp @@ -9,32 +9,32 @@ namespace sw { namespace universal { // hyperbolic sine of an angle of x radians - dd sinh(dd x) { + inline dd sinh(dd x) { return dd(std::sinh(double(x))); } // hyperbolic cosine of an angle of x radians - dd cosh(dd x) { + inline dd cosh(dd x) { return dd(std::cosh(double(x))); } // hyperbolic tangent of an angle of x radians - dd tanh(dd x) { + inline dd tanh(dd x) { return dd(std::tanh(double(x))); } // hyperbolic cotangent of an angle of x radians - dd atanh(dd x) { + inline dd atanh(dd x) { return dd(std::atanh(double(x))); } // hyperbolic cosecant of an angle of x radians - dd acosh(dd x) { + inline dd acosh(dd x) { return dd(std::acosh(double(x))); } // hyperbolic secant of an angle of x radians - dd asinh(dd x) { + inline dd asinh(dd x) { return dd(std::asinh(double(x))); } diff --git a/include/universal/number/dd/math/hypot.hpp b/include/universal/number/dd/math/hypot.hpp index 99bcc4636..14f254290 100644 --- a/include/universal/number/dd/math/hypot.hpp +++ b/include/universal/number/dd/math/hypot.hpp @@ -8,7 +8,7 @@ namespace sw { namespace universal { - dd hypot(dd x, dd y) { + inline dd hypot(dd x, dd y) { return dd(std::hypot(double(x), double(y))); } diff --git a/include/universal/number/dd/math/logarithm.hpp b/include/universal/number/dd/math/logarithm.hpp index 8905ce872..aac8cde20 100644 --- a/include/universal/number/dd/math/logarithm.hpp +++ b/include/universal/number/dd/math/logarithm.hpp @@ -11,113 +11,112 @@ namespace sw { namespace universal { -/// -/// Natural logarithm (base = e) -/// -/// input -/// natural logarithm of a -dd log(const dd& a) { - if (a.isnan() || a.isinf()) return a; - - if (a.iszero()) return dd(SpecificValue::infneg); - - if (a.isone()) return 0.0; - - if (a.sign()) { - std::cerr << "log: non-positive argument\n"; - errno = EDOM; - return dd(SpecificValue::qnan); - } - - /* Strategy. The Taylor series for log converges much more - slowly than that of exp, due to the lack of the factorial - term in the denominator. Hence this routine instead tries - to determine the root of the function + /// + /// Natural logarithm (base = e) + /// + /// input + /// natural logarithm of a + inline dd log(const dd& a) { + if (a.isnan() || a.isinf()) return a; - f(x) = exp(x) - a + if (a.iszero()) return dd(SpecificValue::infneg); - using Newton iteration. The iteration is given by + if (a.isone()) return 0.0; - x' = x - f(x)/f'(x) - = x - (1 - a * exp(-x)) - = x + a * exp(-x) - 1. + if (a.sign()) { + std::cerr << "log: non-positive argument\n"; + errno = EDOM; + return dd(SpecificValue::qnan); + } - Only one iteration is needed, since Newton's iteration - approximately doubles the number of digits per iteration. - */ + /* Strategy. The Taylor series for log converges much more + slowly than that of exp, due to the lack of the factorial + term in the denominator. Hence this routine instead tries + to determine the root of the function - dd x = std::log(a.high()); // Initial approximation - x = x + a * exp(-x) - 1.0; - return x; -} + f(x) = exp(x) - a -/// -/// binary logarithm (base = 2) -/// -/// input -/// binary logarithm of a -dd log2(const dd& a) { - if (a.isnan() || a.isinf()) return a; + using Newton iteration. The iteration is given by - if (a.iszero()) return dd(SpecificValue::infneg); + x' = x - f(x)/f'(x) + = x - (1 - a * exp(-x)) + = x + a * exp(-x) - 1. - if (a.isone()) return 0.0; + Only one iteration is needed, since Newton's iteration + approximately doubles the number of digits per iteration. + */ - if (a.sign()) { - std::cerr << "log2: non-positive argument\n"; - errno = EDOM; - return dd(SpecificValue::qnan); + dd x = std::log(a.high()); // Initial approximation + x = x + a * exp(-x) - 1.0; + return x; } - return log(a) * dd_lge; -} + /// + /// binary logarithm (base = 2) + /// + /// input + /// binary logarithm of a + inline dd log2(const dd& a) { + if (a.isnan() || a.isinf()) return a; -/// -/// decimal logarithm (base = 10) -/// -/// input -/// binary logarithm of a -dd log10(const dd& a) { - if (a.isnan() || a.isinf()) return a; + if (a.iszero()) return dd(SpecificValue::infneg); - if (a.iszero()) return dd(SpecificValue::infneg); + if (a.isone()) return 0.0; - if (a.isone()) return 0.0; + if (a.sign()) { + std::cerr << "log2: non-positive argument\n"; + errno = EDOM; + return dd(SpecificValue::qnan); + } - if (a.sign()) { - std::cerr << "log10: non-positive argument\n"; - errno = EDOM; - return dd(SpecificValue::qnan); + return log(a) * dd_lge; } - return log(a) / dd_log10; -} + /// + /// decimal logarithm (base = 10) + /// + /// input + /// binary logarithm of a + inline dd log10(const dd& a) { + if (a.isnan() || a.isinf()) return a; -/// -/// Natural logarithm of 1+x -/// -/// input -/// binary logarithm of a -dd log1p(const dd& a) -{ - if (a.isnan() || a.isinf()) return a; + if (a.iszero()) return dd(SpecificValue::infneg); - if (a.iszero()) return dd(0.0); + if (a.isone()) return 0.0; - if (a == -1.0) return dd(SpecificValue::infneg); + if (a.sign()) { + std::cerr << "log10: non-positive argument\n"; + errno = EDOM; + return dd(SpecificValue::qnan); + } - if (a < -1.0) { - std::cerr << "log1p: non-positive argument\n"; - errno = EDOM; - return dd(SpecificValue::qnan); + return log(a) / dd_log10; } - if ((a >= 2.0) || (a <= -0.5)) // a >= 2.0 - no loss of significant bits - use log() - return log(1.0 + a); + /// + /// Natural logarithm of 1+x + /// + /// input + /// binary logarithm of a + inline dd log1p(const dd& a) { + if (a.isnan() || a.isinf()) return a; + + if (a.iszero()) return dd(0.0); + + if (a == -1.0) return dd(SpecificValue::infneg); - // at this point, -1.0 < a < 2.0 - // return _log1p(a); - return log(1.0 + a); // TODO: evaluate loss of precision -} + if (a < -1.0) { + std::cerr << "log1p: non-positive argument\n"; + errno = EDOM; + return dd(SpecificValue::qnan); + } + + if ((a >= 2.0) || (a <= -0.5)) // a >= 2.0 - no loss of significant bits - use log() + return log(1.0 + a); + + // at this point, -1.0 < a < 2.0 + // return _log1p(a); + return log(1.0 + a); // TODO: evaluate loss of precision + } }} // namespace sw::universal diff --git a/include/universal/number/dd/math/minmax.hpp b/include/universal/number/dd/math/minmax.hpp index a4771bde8..db55829ea 100644 --- a/include/universal/number/dd/math/minmax.hpp +++ b/include/universal/number/dd/math/minmax.hpp @@ -8,11 +8,11 @@ namespace sw { namespace universal { - dd min(dd x, dd y) { + inline dd min(dd x, dd y) { return dd(std::min(double(x), double(y))); } - dd max(dd x, dd y) { + inline dd max(dd x, dd y) { return dd(std::max(double(x), double(y))); } diff --git a/include/universal/number/dd/math/next.hpp b/include/universal/number/dd/math/next.hpp index 399194719..eb9595b3f 100644 --- a/include/universal/number/dd/math/next.hpp +++ b/include/universal/number/dd/math/next.hpp @@ -8,57 +8,57 @@ namespace sw { namespace universal { -/* -Parameters - x Base value. - t Value toward which the return value is approximated. -If both parameters compare equal, the function returns t. + /* + Parameters + x Base value. + t Value toward which the return value is approximated. + If both parameters compare equal, the function returns t. -Return Value - The next representable value after x in the direction of t. + Return Value + The next representable value after x in the direction of t. - If x is the largest finite value representable in the type, - and the result is infinite or not representable, an overflow range error occurs. + If x is the largest finite value representable in the type, + and the result is infinite or not representable, an overflow range error occurs. - If an overflow range error occurs: - - And math_errhandling has MATH_ERRNO set: the global variable errno is set to ERANGE. - - And math_errhandling has MATH_ERREXCEPT set: FE_OVERFLOW is raised. - */ -dd nextafter(const dd& x, const dd& target) { - double hi{ x.high() }; - double lo; - if (x < target) { - lo = std::nextafter(x.low(), +INFINITY); - } - else { - lo = std::nextafter(x.low(), -INFINITY); - } - return dd(hi, lo); -} - -/* TODO -dd nexttoward(dd x, dd target) { - double _x(x); - if (_x == target) return x; - if (target.isnan()) { - if (_x.isneg()) { - --_x; + If an overflow range error occurs: + - And math_errhandling has MATH_ERRNO set: the global variable errno is set to ERANGE. + - And math_errhandling has MATH_ERREXCEPT set: FE_OVERFLOW is raised. + */ + inline dd nextafter(const dd& x, const dd& target) { + double hi{ x.high() }; + double lo; + if (x < target) { + lo = std::nextafter(x.low(), +INFINITY); } else { - ++_x; + lo = std::nextafter(x.low(), -INFINITY); } + return dd(hi, lo); } - else { - if (_x > target) { - --_x; + + /* TODO + inline dd nexttoward(dd x, dd target) { + double _x(x); + if (_x == target) return x; + if (target.isnan()) { + if (_x.isneg()) { + --_x; + } + else { + ++_x; + } } else { - ++_x; + if (_x > target) { + --_x; + } + else { + ++_x; + } } + x = _x; + return x; } - x = _x; - return x; -} -*/ + */ }} // namespace sw::universal diff --git a/include/universal/number/dd/math/numerics.hpp b/include/universal/number/dd/math/numerics.hpp index c3a473630..c4901d550 100644 --- a/include/universal/number/dd/math/numerics.hpp +++ b/include/universal/number/dd/math/numerics.hpp @@ -13,7 +13,7 @@ namespace sw { namespace universal { // copysign returns a value with the magnitude of a, and the sign of b - dd copysign(dd const& a, dd const& b) { + inline dd copysign(dd const& a, dd const& b) { auto signA = std::copysign(1.0, a.high()); auto signB = std::copysign(1.0, b.high()); @@ -21,14 +21,14 @@ namespace sw { namespace universal { } // decompose doubledouble into a fraction and an exponent - dd frexp(dd const& a, int* pexp) { + inline dd frexp(dd const& a, int* pexp) { double hi = std::frexp(a.high(), pexp); double lo = std::ldexp(a.low(), -*pexp); return dd(hi, lo); } // recompose doubledouble from a fraction and an exponent - dd ldexp(dd const& a, int exp) { + inline dd ldexp(dd const& a, int exp) { static_assert(std::numeric_limits< dd >::radix == 2, "CONFIGURATION: dd radix must be 2!"); static_assert(std::numeric_limits< double >::radix == 2, "CONFIGURATION: double radix must be 2!"); diff --git a/include/universal/number/dd/math/pow.hpp b/include/universal/number/dd/math/pow.hpp index 285a8263e..94ce5e237 100644 --- a/include/universal/number/dd/math/pow.hpp +++ b/include/universal/number/dd/math/pow.hpp @@ -12,21 +12,21 @@ namespace sw { namespace universal { // fwd reference - dd exp(const dd&); + inline dd exp(const dd&); // power function - dd pow(const dd& a, const dd& b) { + inline dd pow(const dd& a, const dd& b) { return exp(b * log(a)); } // power function of a dd to double - dd pow(dd x, double y) { + inline dd pow(dd x, double y) { return pow(x, dd(y)); } // Computes the n-th power of a double-double number. // NOTE: 0^0 causes an error. - dd npwr(const dd& a, int n) { + inline dd npwr(const dd& a, int n) { if (n == 0) { #if DOUBLEDOUBLE_THROW_ARITHMETIC_EXCEPTION if (a.iszero()) throw dd_invalid_argument(); @@ -61,9 +61,8 @@ namespace sw { namespace universal { return s; } - dd pow(const dd& a, int n) { + inline dd pow(const dd& a, int n) { return npwr(a, n); } - }} // namespace sw::universal diff --git a/include/universal/number/dd/math/sqrt.hpp b/include/universal/number/dd/math/sqrt.hpp index 2c819cbb7..3600ffa90 100644 --- a/include/universal/number/dd/math/sqrt.hpp +++ b/include/universal/number/dd/math/sqrt.hpp @@ -57,7 +57,7 @@ namespace sw { namespace universal { #endif // ! DOUBLEDOUBLE_NATIVE_SQRT // Computes the square root of a double in double-double precision. - dd sqrt(double d) { + inline dd sqrt(double d) { return sw::universal::sqrt(dd(d)); } @@ -71,7 +71,7 @@ namespace sw { namespace universal { /* Computes the n-th root of the double-double number a. NOTE: n must be a positive integer. NOTE: If n is even, then a must not be negative. */ - dd nroot(const dd& a, int n) { + inline dd nroot(const dd& a, int n) { /* Strategy: Use Newton iteration for the function f(x) = x^(-n) - a diff --git a/include/universal/number/dd/math/trigonometry.hpp b/include/universal/number/dd/math/trigonometry.hpp index 66505eae8..77402e9f0 100644 --- a/include/universal/number/dd/math/trigonometry.hpp +++ b/include/universal/number/dd/math/trigonometry.hpp @@ -42,10 +42,10 @@ namespace sw { namespace universal { int i = 0; do { r *= x; - t = r * dd(inv_fact[i][0], inv_fact[i][1]); + t = r * dd_inverse_factorial[i]; s += t; i += 2; - } while (i < n_inv_fact && std::abs(double(t)) > thresh); + } while (i < dd_inverse_factorial_table_size && std::abs(double(t)) > thresh); return s; } @@ -62,10 +62,10 @@ namespace sw { namespace universal { int i = 1; do { r *= x; - t = r * dd(inv_fact[i][0], inv_fact[i][1]); + t = r * dd_inverse_factorial[i]; s += t; i += 2; - } while (i < n_inv_fact && std::abs(double(t)) > thresh); + } while (i < dd_inverse_factorial_table_size && std::abs(double(t)) > thresh); return s; } diff --git a/include/universal/number/dd/math/truncate.hpp b/include/universal/number/dd/math/truncate.hpp index 85fed9e30..e714da859 100644 --- a/include/universal/number/dd/math/truncate.hpp +++ b/include/universal/number/dd/math/truncate.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // Truncate value by rounding toward zero, returning the nearest integral value that is not larger in magnitude than x - dd trunc(dd x) { + inline dd trunc(dd x) { return dd(std::trunc(double(x))); } // Round to nearest: returns the integral value that is nearest to x, with halfway cases rounded away from zero - dd round(dd x) { + inline dd round(dd x) { return dd(std::round(double(x))); } diff --git a/include/universal/number/qd/attributes.hpp b/include/universal/number/qd/attributes.hpp index 7dbd3ca18..0252289bf 100644 --- a/include/universal/number/qd/attributes.hpp +++ b/include/universal/number/qd/attributes.hpp @@ -18,7 +18,7 @@ namespace sw { namespace universal { } // generate the maxneg through maxpos value range of a quad-double configuration - std::string qd_range() { + inline std::string qd_range() { qd v; std::stringstream s; s << std::setw(80) << type_tag(v) << " : [ " @@ -47,17 +47,17 @@ namespace sw { namespace universal { } */ - int minpos_scale(const qd& b) { + inline int minpos_scale(const qd& b) { qd c(b); return c.minpos().scale(); } - int maxpos_scale(const qd& b) { + inline int maxpos_scale(const qd& b) { qd c(b); return c.maxpos().scale(); } - int max_negative_scale(const qd& b) { + inline int max_negative_scale(const qd& b) { qd c(b); return c.maxneg().scale(); } diff --git a/include/universal/number/qd/math/classify.hpp b/include/universal/number/qd/math/classify.hpp index 9f02eaef4..701e45821 100644 --- a/include/universal/number/qd/math/classify.hpp +++ b/include/universal/number/qd/math/classify.hpp @@ -8,48 +8,48 @@ namespace sw { namespace universal { -// STD LIB function for IEEE floats: Categorizes floating point value arg into the following categories: zero, subnormal, normal, infinite, NAN, or implementation-defined category. -int fpclassify(const qd& a) { - return (std::fpclassify(a[0])); -} + // Categorizes floating point value arg into the following categories: zero, subnormal, normal, infinite, NAN, or implementation-defined category. + inline int fpclassify(const qd& a) { + return (std::fpclassify(a[0])); + } -// STD LIB function for IEEE floats: Determines if the given floating point number arg is a cfloative or negative infinity. -// specialized for quad-double (qd) -inline bool isinf(const qd& a) { - return (std::fpclassify(a[0]) == FP_INFINITE); -} - -// STD LIB function for IEEE floats: Determines if the given floating point number arg is a not-a-number (NaN) value. -// specialized for quad-double (qd) -inline bool isnan(const qd& a) { - return (std::fpclassify(a[0]) == FP_NAN); -} - -// STD LIB function for IEEE floats: Determines if the given floating point number arg has finite value i.e. it is normal, subnormal or zero, but not infinite or NaN. -// specialized for quad-double (qd) -inline bool isfinite(const qd& a) { - return (std::fpclassify(a[0]) != FP_INFINITE) && (std::fpclassify(a[0]) != FP_NAN); -} - -// STD LIB function for IEEE floats: Determines if the given floating point number arg is normal, i.e. is neither zero, subnormal, infinite, nor NaN. -// specialized for quad-double (qd) -inline bool isnormal(const qd& a) { - return (std::fpclassify(a[0]) == FP_NORMAL); -} - -// STD LIB function for IEEE floats: Determines if the given floating point number arg is denormal, i.e. is neither zero, normal, infinite, nor NaN. -// specialized for quad-double (qd) -inline bool isdenorm(const qd& a) { - return (std::fpclassify(a[0]) == FP_SUBNORMAL); -} - -inline bool iszero(const qd& a) { - return (std::fpclassify(a[0]) == FP_ZERO); -} - -bool signbit(qd const& a) { - auto signA = std::copysign(1.0, a[0]); - return signA < 0.0; -} + // Determines if the given floating point number arg is a cfloative or negative infinity. + // specialized for quad-double (qd) + inline bool isinf(const qd& a) { + return (std::fpclassify(a[0]) == FP_INFINITE); + } + + // Determines if the given floating point number arg is a not-a-number (NaN) value. + // specialized for quad-double (qd) + inline bool isnan(const qd& a) { + return (std::fpclassify(a[0]) == FP_NAN); + } + + // Determines if the given floating point number arg has finite value i.e. it is normal, subnormal or zero, but not infinite or NaN. + // specialized for quad-double (qd) + inline bool isfinite(const qd& a) { + return (std::fpclassify(a[0]) != FP_INFINITE) && (std::fpclassify(a[0]) != FP_NAN); + } + + // Determines if the given floating point number arg is normal, i.e. is neither zero, subnormal, infinite, nor NaN. + // specialized for quad-double (qd) + inline bool isnormal(const qd& a) { + return (std::fpclassify(a[0]) == FP_NORMAL); + } + + // Determines if the given floating point number arg is denormal, i.e. is neither zero, normal, infinite, nor NaN. + // specialized for quad-double (qd) + inline bool isdenorm(const qd& a) { + return (std::fpclassify(a[0]) == FP_SUBNORMAL); + } + + inline bool iszero(const qd& a) { + return (std::fpclassify(a[0]) == FP_ZERO); + } + + inline bool signbit(qd const& a) { + auto signA = std::copysign(1.0, a[0]); + return signA < 0.0; + } }} // namespace sw::universal diff --git a/include/universal/number/qd/math/error_and_gamma.hpp b/include/universal/number/qd/math/error_and_gamma.hpp index c3f78063a..9652d836d 100644 --- a/include/universal/number/qd/math/error_and_gamma.hpp +++ b/include/universal/number/qd/math/error_and_gamma.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // Compute the error function erf(x) = 2 over sqrt(PI) times Integral from 0 to x of e ^ (-t)^2 dt - qd erf(qd x) { + inline qd erf(qd x) { return qd(std::erf(x[0])); } // Compute the complementary error function: 1 - erf(x) - qd erfc(qd x) { + inline qd erfc(qd x) { return qd(std::erfc(x[0])); } diff --git a/include/universal/number/qd/math/exponent.hpp b/include/universal/number/qd/math/exponent.hpp index 5c42e8c98..fb6349b32 100644 --- a/include/universal/number/qd/math/exponent.hpp +++ b/include/universal/number/qd/math/exponent.hpp @@ -13,26 +13,26 @@ namespace sw { namespace universal { // fwd reference 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) + constexpr unsigned qd_inverse_factorial_table_size = 15; + static const qd qd_inverse_factorial[qd_inverse_factorial_table_size] = { + qd(1.66666666666666657e-01, 9.25185853854297066e-18, 5.13581318503262866e-34, 2.85094902409834186e-50), // 1/3! + qd(4.16666666666666644e-02, 2.31296463463574266e-18, 1.28395329625815716e-34, 7.12737256024585466e-51), // 1/4! + qd(8.33333333333333322e-03, 1.15648231731787138e-19, 1.60494162032269652e-36, 2.22730392507682967e-53), // 1/5! + qd(1.38888888888888894e-03, -5.30054395437357706e-20, -1.73868675534958776e-36, -1.63335621172300840e-52), // 1/6! + qd(1.98412698412698413e-04, 1.72095582934207053e-22, 1.49269123913941271e-40, 1.29470326746002471e-58), // 1/7! + qd(2.48015873015873016e-05, 2.15119478667758816e-23, 1.86586404892426588e-41, 1.61837908432503088e-59), // 1/8! + qd(2.75573192239858925e-06, -1.85839327404647208e-22, 8.49175460488199287e-39, -5.72661640789429621e-55), // 1/9! + qd(2.75573192239858883e-07, 2.37677146222502973e-23, -3.26318890334088294e-40, 1.61435111860404415e-56), // 1/10! + qd(2.50521083854417202e-08, -1.44881407093591197e-24, 2.04267351467144546e-41, -8.49632672007163175e-58), // 1/11! + qd(2.08767569878681002e-09, -1.20734505911325997e-25, 1.70222792889287100e-42, 1.41609532150396700e-58), // 1/12! + qd(1.60590438368216133e-10, 1.25852945887520981e-26, -5.31334602762985031e-43, 3.54021472597605528e-59), // 1/13! + qd(1.14707455977297245e-11, 2.06555127528307454e-28, 6.88907923246664603e-45, 5.72920002655109095e-61), // 1/14! + qd(7.64716373181981641e-13, 7.03872877733453001e-30, -7.82753927716258345e-48, 1.92138649443790242e-64), // 1/15! + qd(4.77947733238738525e-14, 4.39920548583408126e-31, -4.89221204822661465e-49, 1.20086655902368901e-65), // 1/16! + qd(2.81145725434552060e-15, 1.65088427308614326e-31, -2.87777179307447918e-50, 4.27110689256293549e-67) // 1/17! }; - qd exp(const qd& x) { + inline 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 @@ -64,7 +64,7 @@ namespace sw { namespace universal { int i = 0; do { p *= r; - t = p * inv_fact[i++]; + t = p * qd_inverse_factorial[i++]; s += t; } while (std::abs(double(t)) > thresh && i < 9); @@ -88,20 +88,20 @@ namespace sw { namespace universal { return ldexp(s, static_cast(m)); } -// Base-2 exponential function -qd exp2(const qd& x) { - return qd(std::exp2(double(x))); -} + // Base-2 exponential function + inline qd exp2(const qd& x) { + return qd(std::exp2(double(x))); + } -// Base-10 exponential function -qd exp10(const qd& x) { - return qd(std::pow(10.0, double(x))); -} + // Base-10 exponential function + inline qd exp10(const qd& x) { + return qd(std::pow(10.0, double(x))); + } -// Base-e exponential function exp(x)-1 -qd expm1(const qd& x) { - return qd(std::expm1(double(x))); -} + // Base-e exponential function exp(x)-1 + inline qd expm1(const qd& x) { + return qd(std::expm1(double(x))); + } }} // namespace sw::universal diff --git a/include/universal/number/qd/math/fractional.hpp b/include/universal/number/qd/math/fractional.hpp index 637f03d73..0911eb335 100644 --- a/include/universal/number/qd/math/fractional.hpp +++ b/include/universal/number/qd/math/fractional.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // fmod retuns x - n*y where n = x/y with the fractional part truncated - qd fmod(const qd& x, const qd& y) { + inline qd fmod(const qd& x, const qd& y) { return qd(std::fmod(double(x), double(y))); } // shim to stdlib - qd remainder(const qd& x, const qd& y) { + inline qd remainder(const qd& x, const qd& y) { return qd(std::remainder(double(x), double(y))); } diff --git a/include/universal/number/qd/math/hyperbolic.hpp b/include/universal/number/qd/math/hyperbolic.hpp index 33f7b0a49..2b43c015b 100644 --- a/include/universal/number/qd/math/hyperbolic.hpp +++ b/include/universal/number/qd/math/hyperbolic.hpp @@ -9,32 +9,32 @@ namespace sw { namespace universal { // hyperbolic sine of an angle of x radians - qd sinh(qd x) { + inline qd sinh(qd x) { return qd(std::sinh(double(x))); } // hyperbolic cosine of an angle of x radians - qd cosh(qd x) { + inline qd cosh(qd x) { return qd(std::cosh(double(x))); } // hyperbolic tangent of an angle of x radians - qd tanh(qd x) { + inline qd tanh(qd x) { return qd(std::tanh(double(x))); } // hyperbolic cotangent of an angle of x radians - qd atanh(qd x) { + inline qd atanh(qd x) { return qd(std::atanh(double(x))); } // hyperbolic cosecant of an angle of x radians - qd acosh(qd x) { + inline qd acosh(qd x) { return qd(std::acosh(double(x))); } // hyperbolic secant of an angle of x radians - qd asinh(qd x) { + inline qd asinh(qd x) { return qd(std::asinh(double(x))); } diff --git a/include/universal/number/qd/math/hypot.hpp b/include/universal/number/qd/math/hypot.hpp index 41969f3d1..8e6d72eb3 100644 --- a/include/universal/number/qd/math/hypot.hpp +++ b/include/universal/number/qd/math/hypot.hpp @@ -8,7 +8,7 @@ namespace sw { namespace universal { - qd hypot(const qd& x, const qd& y) { + inline qd hypot(const qd& x, const qd& y) { return qd(std::hypot(double(x), double(y))); } diff --git a/include/universal/number/qd/math/logarithm.hpp b/include/universal/number/qd/math/logarithm.hpp index 4150281d0..53f2de81c 100644 --- a/include/universal/number/qd/math/logarithm.hpp +++ b/include/universal/number/qd/math/logarithm.hpp @@ -9,7 +9,7 @@ namespace sw { namespace universal { - qd log(const qd& a) { + inline qd log(const qd& a) { if (a.isnan() || a.isinf()) return a; if (a.iszero()) return qd(SpecificValue::infneg); @@ -49,11 +49,11 @@ namespace sw { namespace universal { } /// -/// binary logarithm (base = 2) -/// -/// input -/// binary logarithm of a - qd log2(const qd& a) { + /// binary logarithm (base = 2) + /// + /// input + /// binary logarithm of a + inline qd log2(const qd& a) { if (a.isnan() || a.isinf()) return a; if (a.iszero()) return qd(SpecificValue::infneg); @@ -74,7 +74,7 @@ namespace sw { namespace universal { /// /// input /// binary logarithm of a - qd log10(const qd& a) { + inline qd log10(const qd& a) { if (a.isnan() || a.isinf()) return a; if (a.iszero()) return qd(SpecificValue::infneg); @@ -96,7 +96,7 @@ namespace sw { namespace universal { /// /// input /// binary logarithm of a - qd log1p(const qd& a) { + inline qd log1p(const qd& a) { if (a.isnan() || a.isinf()) return a; if (a.iszero()) return qd(0.0); diff --git a/include/universal/number/qd/math/minmax.hpp b/include/universal/number/qd/math/minmax.hpp index 997ce74ba..be67a2ca6 100644 --- a/include/universal/number/qd/math/minmax.hpp +++ b/include/universal/number/qd/math/minmax.hpp @@ -8,11 +8,11 @@ namespace sw { namespace universal { - qd min(const qd& x, const qd& y) { + inline qd min(const qd& x, const qd& y) { return qd(std::min(double(x), double(y))); } - qd max(const qd& x, const qd& y) { + inline qd max(const qd& x, const qd& y) { return qd(std::max(double(x), double(y))); } diff --git a/include/universal/number/qd/math/next.hpp b/include/universal/number/qd/math/next.hpp index eb4cc0d6c..68b415ddc 100644 --- a/include/universal/number/qd/math/next.hpp +++ b/include/universal/number/qd/math/next.hpp @@ -8,63 +8,31 @@ namespace sw { namespace universal { -/* -Parameters - x Base value. - t Value toward which the return value is approximated. -If both parameters compare equal, the function returns t. + /* + Parameters + x Base value. + t Value toward which the return value is approximated. + If both parameters compare equal, the function returns t. -Return Value - The next representable value after x in the direction of t. + Return Value + The next representable value after x in the direction of t. - If x is the largest finite value representable in the type, - and the result is infinite or not representable, an overflow range error occurs. + If x is the largest finite value representable in the type, + and the result is infinite or not representable, an overflow range error occurs. - If an overflow range error occurs: - - And math_errhandling has MATH_ERRNO set: the global variable errno is set to ERANGE. - - And math_errhandling has MATH_ERREXCEPT set: FE_OVERFLOW is raised. - */ -qd nextafter(qd x, qd target) { - if (x == target) return target; - if (target.isnan()) { - if (x.isneg()) { - --x; - } - else { - ++x; - } - } - else { - if (x > target) { - --x; - } - else { - ++x; - } + If an overflow range error occurs: + - And math_errhandling has MATH_ERRNO set: the global variable errno is set to ERANGE. + - And math_errhandling has MATH_ERREXCEPT set: FE_OVERFLOW is raised. + */ + inline qd nextafter(qd x, qd target) { + + return x; } - return x; -} -qd nexttoward(qd x, qd target) { - if (x == target) return x; - if (target.isnan()) { - if (x.isneg()) { - --x; - } - else { - ++x; - } - } - else { - if (x > target) { - --x; - } - else { - ++x; - } + inline qd nexttoward(qd x, qd target) { + + return x; } - return x; -} }} // namespace sw::universal diff --git a/include/universal/number/qd/math/numerics.hpp b/include/universal/number/qd/math/numerics.hpp index 54317c9a7..e9c17587b 100644 --- a/include/universal/number/qd/math/numerics.hpp +++ b/include/universal/number/qd/math/numerics.hpp @@ -13,7 +13,7 @@ namespace sw { namespace universal { // copysign returns a value with the magnitude of a, and the sign of b - qd copysign(const qd& a, const qd& b) { + inline qd copysign(const qd& a, const qd& b) { auto signA = std::copysign(1.0, a[0]); auto signB = std::copysign(1.0, b[0]); @@ -21,7 +21,7 @@ namespace sw { namespace universal { } // decompose quad-double into a fraction and an exponent - qd frexp(const qd& a, int* pexp) { + inline qd frexp(const qd& a, int* pexp) { double a0 = std::frexp(a[0], pexp); double a1 = std::ldexp(a[1], -*pexp); double a2 = std::ldexp(a[2], -*pexp); @@ -30,7 +30,7 @@ namespace sw { namespace universal { } // recompose quad-double from a fraction and an exponent - qd ldexp(const qd& a, int exponent) { + inline qd ldexp(const qd& a, int exponent) { static_assert(std::numeric_limits< qd >::radix == 2, "CONFIGURATION: qd radix must be 2!"); static_assert(std::numeric_limits< double >::radix == 2, "CONFIGURATION: double radix must be 2!"); diff --git a/include/universal/number/qd/math/pow.hpp b/include/universal/number/qd/math/pow.hpp index d52ae6325..742b030cb 100644 --- a/include/universal/number/qd/math/pow.hpp +++ b/include/universal/number/qd/math/pow.hpp @@ -12,21 +12,21 @@ namespace sw { namespace universal { // fwd reference - qd exp(const qd&); + inline qd exp(const qd&); // power function - qd pow(const qd& a, const qd& b) { + inline qd pow(const qd& a, const qd& b) { return exp(b * log(a)); } // power function of a qd to double - qd pow(const qd& x, double y) { + inline qd pow(const qd& x, double y) { return pow(x, qd(y)); } // Computes the n-th power of a quad-double number. // NOTE: 0^0 causes an error. - qd npwr(const qd& a, int n) { + inline qd npwr(const qd& a, int n) { if (n == 0) { #if QUADDOUBLE_THROW_ARITHMETIC_EXCEPTION if (a.iszero()) throw qd_invalid_argument(); @@ -61,7 +61,7 @@ namespace sw { namespace universal { return s; } - qd pow(const qd& a, int n) { + inline qd pow(const qd& a, int n) { return npwr(a, n); } diff --git a/include/universal/number/qd/math/sqrt.hpp b/include/universal/number/qd/math/sqrt.hpp index 2d6ab09cf..fd93c9cd2 100644 --- a/include/universal/number/qd/math/sqrt.hpp +++ b/include/universal/number/qd/math/sqrt.hpp @@ -66,7 +66,7 @@ namespace sw { namespace universal { /* Computes the n-th root of the quad-double number a. NOTE: n must be a positive integer. NOTE: If n is even, then a must not be negative. */ - qd nroot(const qd& a, int n) { + inline qd nroot(const qd& a, int n) { /* Strategy: Use Newton iteration for the function f(x) = x^(-n) - a diff --git a/include/universal/number/qd/math/trigonometry.hpp b/include/universal/number/qd/math/trigonometry.hpp index caccfa31c..c723f4444 100644 --- a/include/universal/number/qd/math/trigonometry.hpp +++ b/include/universal/number/qd/math/trigonometry.hpp @@ -11,65 +11,56 @@ namespace sw { namespace universal { // value representing an angle expressed in radians - // One radian is equivalent to 180/PI degrees + // One radian is equivalent to 180/PI degrees - // sine of an angle of x radians + // sine of an angle of x radians - qd sin(const qd& x) { + inline qd sin(const qd& x) { return qd(std::sin(double(x))); } - // cosine of an angle of x radians - - qd cos(const qd& x) { + // cosine of an angle of x radians + inline qd cos(const qd& x) { return qd(std::cos(double(x))); } // tangent of an angle of x radians - - qd tan(const qd& x) { + inline qd tan(const qd& x) { return qd(std::tan(double(x))); } // cotangent of an angle of x radians - - qd atan(const qd& x) { + inline qd atan(const qd& x) { return qd(std::atan(double(x))); } // Arc tangent with two parameters - - qd atan2(qd y, const qd& x) { + inline qd atan2(qd y, const qd& x) { return qd(std::atan2(double(y), double(x))); } // cosecant of an angle of x radians - - qd acos(const qd& x) { + inline qd acos(const qd& x) { return qd(std::acos(double(x))); } // secant of an angle of x radians - - qd asin(const qd& x) { + inline qd asin(const qd& x) { return qd(std::asin(double(x))); } // cotangent an angle of x radians - - qd cot(const qd& x) { + inline qd cot(const qd& x) { return qd(std::tan(std::numbers::pi*2 - double(x))); } // secant of an angle of x radians - - qd sec(const qd& x) { + inline qd sec(const qd& x) { return qd(1.0 / std::cos(double(x))); } // cosecant of an angle of x radians - - qd csc(const qd& x) { + inline qd csc(const qd& x) { return qd(1.0 / std::sin(double(x))); } diff --git a/include/universal/number/qd/math/truncate.hpp b/include/universal/number/qd/math/truncate.hpp index bcbcb3236..293494312 100644 --- a/include/universal/number/qd/math/truncate.hpp +++ b/include/universal/number/qd/math/truncate.hpp @@ -9,12 +9,12 @@ namespace sw { namespace universal { // Truncate value by rounding toward zero, returning the nearest integral value that is not larger in magnitude than x - qd trunc(qd x) { + inline qd trunc(qd x) { return qd(std::trunc(double(x))); } // Round to nearest: returns the integral value that is nearest to x, with halfway cases rounded away from zero - qd round(qd x) { + inline qd round(qd x) { return qd(std::round(double(x))); } diff --git a/include/universal/number/qd/qd_impl.hpp b/include/universal/number/qd/qd_impl.hpp index 27bdb4ea4..a77859a41 100644 --- a/include/universal/number/qd/qd_impl.hpp +++ b/include/universal/number/qd/qd_impl.hpp @@ -1372,7 +1372,7 @@ inline qd floor(const qd& a) { } // Round to Nearest integer -qd nint(const qd& a) { +inline qd nint(const qd& a) { double x0{ 0.0 }, x1{ 0.0 }, x2{ 0.0 }, x3{ 0.0 }; x0 = nint(a[0]); @@ -1468,7 +1468,7 @@ inline qd sqr(const qd& a) { } // Computes pow(qd, n), where n is an integer -qd pown(const qd& a, int n) { +inline qd pown(const qd& a, int n) { if (n == 0) return 1.0; @@ -1525,7 +1525,7 @@ inline std::istream& operator>>(std::istream& istr, qd& v) { ////////////////// string operators // parse a decimal ASCII floating-point format and make a quad-double (qd) out of it -bool parse(const std::string& number, qd& value) { +inline bool parse(const std::string& number, qd& value) { char const* p = number.c_str(); // Skip any leading spaces diff --git a/static/appenv/multifile/areals.cpp b/static/appenv/multifile/areals.cpp index c8f66f4a1..1dafe4804 100644 --- a/static/appenv/multifile/areals.cpp +++ b/static/appenv/multifile/areals.cpp @@ -1,6 +1,7 @@ // areals.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/cfloats.cpp b/static/appenv/multifile/cfloats.cpp index b8635705c..7e3c290ef 100644 --- a/static/appenv/multifile/cfloats.cpp +++ b/static/appenv/multifile/cfloats.cpp @@ -1,6 +1,7 @@ // cfloats.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/dbns.cpp b/static/appenv/multifile/dbns.cpp index fbd0d43c5..b95545253 100644 --- a/static/appenv/multifile/dbns.cpp +++ b/static/appenv/multifile/dbns.cpp @@ -1,6 +1,7 @@ // dbns.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/dd.cpp b/static/appenv/multifile/dd.cpp new file mode 100644 index 000000000..fbaea2aab --- /dev/null +++ b/static/appenv/multifile/dd.cpp @@ -0,0 +1,25 @@ +// dds.cpp: compilation test to check arithmetic type usage in application environments +// +// 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 +#include +#include + +using DoubleDouble = sw::universal::dd; + +DoubleDouble ddPolynomial(const std::vector& coef, const DoubleDouble& x) { + using namespace sw::universal; + if (coef.size() < 2) { + std::cerr << "Coefficient set is too small to represent a polynomial\n"; + return DoubleDouble(0); + } + + DoubleDouble v = coef[0]; + for (size_t i = 1; i < coef.size(); ++i) { + v += coef[i] * pow(x, DoubleDouble(i)); + } + return v; +} diff --git a/static/appenv/multifile/fixpnts.cpp b/static/appenv/multifile/fixpnts.cpp index eda8d41e9..a359616fa 100644 --- a/static/appenv/multifile/fixpnts.cpp +++ b/static/appenv/multifile/fixpnts.cpp @@ -1,6 +1,7 @@ // fixpnts.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/integers.cpp b/static/appenv/multifile/integers.cpp index 7173a7a83..31db289d0 100644 --- a/static/appenv/multifile/integers.cpp +++ b/static/appenv/multifile/integers.cpp @@ -1,7 +1,8 @@ // integers.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 @@ -10,7 +11,7 @@ using Integer = sw::universal::integer<8, uint8_t, sw::universal::IntegerNumberType::IntegerNumber>; -Integer integerPolynomial(const std::vector& coef, const Integer& x) { +Integer integerPolynomial(const std::vector& coef, const Integer& x) { using namespace sw::universal; if (coef.size() < 2) { std::cerr << "Coefficient set is too small to represent a polynomial\n"; diff --git a/static/appenv/multifile/logs.cpp b/static/appenv/multifile/logs.cpp index 95ba9755c..4b3eb527a 100644 --- a/static/appenv/multifile/logs.cpp +++ b/static/appenv/multifile/logs.cpp @@ -1,6 +1,7 @@ // logs.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/main.cpp b/static/appenv/multifile/main.cpp index 4005cfa60..480c8cf2a 100644 --- a/static/appenv/multifile/main.cpp +++ b/static/appenv/multifile/main.cpp @@ -1,6 +1,7 @@ // multifile.cpp: compilation test to check arithmetic type usage in application environments // // 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 @@ -18,6 +19,8 @@ #include #include #include +#include +#include #include @@ -38,29 +41,40 @@ #endif // forward references -using Integer = sw::universal::integer<8, uint8_t, sw::universal::IntegerNumberType::IntegerNumber>; -using Fixpnt = sw::universal::fixpnt<8, 4, sw::universal::Saturate, uint8_t>; -using Cfloat = sw::universal::half; -using Posit = sw::universal::posit<8,2>; -using Lns = sw::universal::lns<8, 2>; -using Lns2b = sw::universal::dbns<8, 6>; - -Integer integerPolynomial(const std::vector& coef, const Integer& x); -Fixpnt fixpntPolynomial(const std::vector& coef, const Fixpnt& x); -Cfloat cfloatPolynomial(const std::vector& coef, const Cfloat& x); -Posit positPolynomial(const std::vector& coef, const Posit& x); -Lns lnsPolynomial(const std::vector& coef, const Lns& x); -Lns2b dbnsPolynomial(const std::vector& coef, const Lns2b& x); +using Integer = sw::universal::integer<8, uint8_t, sw::universal::IntegerNumberType::IntegerNumber>; +using Fixpnt = sw::universal::fixpnt<8, 4, sw::universal::Saturate, uint8_t>; +using Cfloat = sw::universal::half; +using Posit = sw::universal::posit<8,2>; +using Lns = sw::universal::lns<8, 2>; +using Lns2b = sw::universal::dbns<8, 6>; +using DoubleDouble = sw::universal::dd; +using QuadDouble = sw::universal::qd; + +Integer integerPolynomial(const std::vector& coef, const Integer& x); +Fixpnt fixpntPolynomial(const std::vector& coef, const Fixpnt& x); +Cfloat cfloatPolynomial(const std::vector& coef, const Cfloat& x); +Posit positPolynomial(const std::vector& coef, const Posit& x); +Lns lnsPolynomial(const std::vector& coef, const Lns& x); +Lns2b dbnsPolynomial(const std::vector& coef, const Lns2b& x); +DoubleDouble ddPolynomial(const std::vector& coef, const DoubleDouble& x); +QuadDouble qdPolynomial(const std::vector& coef, const QuadDouble& x); template void EvaluatePolynomial(const std::vector& coefficients, const NumberType& x) { + std::vector intCoefficients; + for (auto e : coefficients) intCoefficients.push_back(static_cast(e)); + std::vector doubleCoefficients; + for (auto e : coefficients) doubleCoefficients.push_back(double(e)); + std::cout << "x : " << x << '\n'; -// std::cout << "integer : " << integerPolynomial(coefficients, Integer(x)) << '\n'; + std::cout << "integer : " << integerPolynomial(intCoefficients, Integer(x)) << '\n'; std::cout << "fixpnt : " << fixpntPolynomial(coefficients, Fixpnt(x)) << '\n'; std::cout << "cfloat : " << cfloatPolynomial(coefficients, Cfloat(x)) << '\n'; std::cout << "posit : " << positPolynomial(coefficients, Posit(x)) << '\n'; std::cout << "lns : " << lnsPolynomial(coefficients, Lns(x)) << '\n'; - std::cout << "dbns : " << dbnsPolynomial(coefficients, Lns2b(x)) << '\n'; + std::cout << "dbns : " << dbnsPolynomial(coefficients, Lns2b(x)) << '\n'; + std::cout << "double-double: " << ddPolynomial(doubleCoefficients, DoubleDouble(x)) << '\n'; + std::cout << "quad-double : " << qdPolynomial(doubleCoefficients, QuadDouble(x)) << '\n'; } int main() diff --git a/static/appenv/multifile/posits.cpp b/static/appenv/multifile/posits.cpp index 2458771a1..fa5802c0f 100644 --- a/static/appenv/multifile/posits.cpp +++ b/static/appenv/multifile/posits.cpp @@ -1,6 +1,7 @@ // posits.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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 diff --git a/static/appenv/multifile/qd.cpp b/static/appenv/multifile/qd.cpp new file mode 100644 index 000000000..470d29119 --- /dev/null +++ b/static/appenv/multifile/qd.cpp @@ -0,0 +1,25 @@ +// dds.cpp: compilation test to check arithmetic type usage in application environments +// +// 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 +#include +#include + +using QuadDouble = sw::universal::qd; + +QuadDouble qdPolynomial(const std::vector& coef, const QuadDouble& x) { + using namespace sw::universal; + if (coef.size() < 2) { + std::cerr << "Coefficient set is too small to represent a polynomial\n"; + return QuadDouble(0); + } + + QuadDouble v = coef[0]; + for (size_t i = 1; i < coef.size(); ++i) { + v += coef[i] * pow(x, QuadDouble(i)); + } + return v; +} diff --git a/static/appenv/multifile/unums.cpp b/static/appenv/multifile/unums.cpp index d97291c72..f56726f94 100644 --- a/static/appenv/multifile/unums.cpp +++ b/static/appenv/multifile/unums.cpp @@ -1,6 +1,7 @@ // unums.cpp: compilation test to check arithmetic type usage in application environments // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// 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