From 5100709a769ce441d5135a44a5d4b43197cb0823 Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Thu, 1 Nov 2018 20:31:08 -0400 Subject: [PATCH 1/5] Version 1.7.0 --- README.md | 2 +- docs/source/api/algorithms.rst | 4 +- docs/source/api/basic_functions.rst | 6 +- docs/source/api/special_functions.rst | 16 +- docs/source/index.rst | 2 +- include/gcem_incl/abs.hpp | 1 + include/gcem_incl/acos.hpp | 5 +- include/gcem_incl/acosh.hpp | 4 +- include/gcem_incl/asin.hpp | 5 +- include/gcem_incl/asinh.hpp | 4 +- include/gcem_incl/atan.hpp | 10 +- include/gcem_incl/atanh.hpp | 5 +- include/gcem_incl/beta.hpp | 9 +- include/gcem_incl/binomial_coef.hpp | 53 ++++-- include/gcem_incl/cos.hpp | 5 +- include/gcem_incl/cosh.hpp | 4 +- include/gcem_incl/erf.hpp | 9 +- include/gcem_incl/erf_inv.hpp | 21 ++- include/gcem_incl/exp.hpp | 7 +- include/gcem_incl/factorial.hpp | 4 + include/gcem_incl/find_exponent.hpp | 5 +- include/gcem_incl/find_fraction.hpp | 1 + include/gcem_incl/find_whole.hpp | 1 + include/gcem_incl/floor.hpp | 1 + include/gcem_incl/gcd.hpp | 45 +++-- include/gcem_incl/gcem_options.hpp | 10 +- include/gcem_incl/incomplete_beta.hpp | 33 +++- include/gcem_incl/incomplete_beta_inv.hpp | 44 ++++- include/gcem_incl/incomplete_gamma.hpp | 30 +++- include/gcem_incl/incomplete_gamma_inv.hpp | 38 ++++- include/gcem_incl/is_even.hpp | 1 + include/gcem_incl/is_odd.hpp | 1 + include/gcem_incl/lbeta.hpp | 9 +- include/gcem_incl/lcm.hpp | 32 +++- include/gcem_incl/lgamma.hpp | 8 +- include/gcem_incl/lmgamma.hpp | 26 +-- include/gcem_incl/log.hpp | 10 +- include/gcem_incl/mantissa.hpp | 1 + include/gcem_incl/max.hpp | 11 +- include/gcem_incl/min.hpp | 11 +- include/gcem_incl/pow.hpp | 42 ++++- include/gcem_incl/pow_integral.hpp | 102 ++++++----- include/gcem_incl/sgn.hpp | 1 + include/gcem_incl/sin.hpp | 5 +- include/gcem_incl/sinh.hpp | 4 +- include/gcem_incl/sqrt.hpp | 7 +- include/gcem_incl/tan.hpp | 17 +- include/gcem_incl/tanh.hpp | 6 +- include/gcem_incl/tgamma.hpp | 4 +- tests/gcd.cpp | 52 ++++++ tests/gcem_tests.hpp | 190 +++++++++++---------- tests/incomplete_beta.cpp | 2 +- tests/incomplete_beta_inv.cpp | 2 +- tests/incomplete_gamma.cpp | 14 +- tests/lcm.cpp | 52 ++++++ 55 files changed, 727 insertions(+), 267 deletions(-) create mode 100644 tests/gcd.cpp create mode 100644 tests/lcm.cpp diff --git a/README.md b/README.md index ea1b26e..313a5ae 100644 --- a/README.md +++ b/README.md @@ -116,7 +116,7 @@ GCE-Math functions are written as C++ templates with `constexpr` specifiers, the template constexpr return_t -erf(const T x); +erf(const T x) noexcept; ``` where a set of internal templated ```constexpr``` functions will implement a continued fraction expansion to return a value of type ```return_t```. This output type ('```return_t```') is generally determined by the input type, e.g., ```int```, ```float```, ```double```, ```long double```, etc. When ```T``` is an intergral type, the output will be upgraded to ```return_t = double```, otherwise ```return_t = T```. For types not covered by ```std::is_integral```, recasts should be used. diff --git a/docs/source/api/algorithms.rst b/docs/source/api/algorithms.rst index 6cb3110..d2308db 100644 --- a/docs/source/api/algorithms.rst +++ b/docs/source/api/algorithms.rst @@ -8,9 +8,9 @@ Algorithms =============== .. _gcd-function-reference: -.. doxygenfunction:: gcd(const T, const T) +.. doxygenfunction:: gcd(const T1, const T2) :project: gcem .. _lcm-function-reference: -.. doxygenfunction:: lcm(const T, const T) +.. doxygenfunction:: lcm(const T1, const T2) :project: gcem diff --git a/docs/source/api/basic_functions.rst b/docs/source/api/basic_functions.rst index 10540e9..a0e12af 100644 --- a/docs/source/api/basic_functions.rst +++ b/docs/source/api/basic_functions.rst @@ -24,15 +24,15 @@ Basic functions :project: gcem .. _max-function-reference: -.. doxygenfunction:: max(const T, const T) +.. doxygenfunction:: max(const T1, const T2) :project: gcem .. _min-function-reference: -.. doxygenfunction:: min(const T, const T) +.. doxygenfunction:: min(const T1, const T2) :project: gcem .. _pow-function-reference: -.. doxygenfunction:: pow(const Ta, const Tb) +.. doxygenfunction:: pow(const T1, const T2) :project: gcem .. _sgn-function-reference: diff --git a/docs/source/api/special_functions.rst b/docs/source/api/special_functions.rst index 19b2879..9f0e6bc 100644 --- a/docs/source/api/special_functions.rst +++ b/docs/source/api/special_functions.rst @@ -8,15 +8,15 @@ Special functions =============== .. _binomial-func-ref: -.. doxygenfunction:: binomial_coef(const pT, const pT) +.. doxygenfunction:: binomial_coef(const T1, const T2) :project: gcem .. _beta-function-reference: -.. doxygenfunction:: beta(const T, const T) +.. doxygenfunction:: beta(const T1, const T2) :project: gcem .. _lbeta-func-ref: -.. doxygenfunction:: lbeta(const T, const T) +.. doxygenfunction:: lbeta(const T1, const T2) :project: gcem .. _tgamma-func-ref: @@ -28,7 +28,7 @@ Special functions :project: gcem .. _lmgamma-func-ref: -.. doxygenfunction:: lmgamma(const pT, const eT) +.. doxygenfunction:: lmgamma(const T1, const T2) :project: gcem Incomplete integral functions @@ -39,11 +39,11 @@ Incomplete integral functions :project: gcem .. _ib-func-ref: -.. doxygenfunction:: incomplete_beta(const pT, const pT, const eT) +.. doxygenfunction:: incomplete_beta(const T1, const T2, const T3) :project: gcem .. _ig-func-ref: -.. doxygenfunction:: incomplete_gamma(const pT, const eT) +.. doxygenfunction:: incomplete_gamma(const T1, const T2) :project: gcem Inverse incomplete integral functions @@ -54,9 +54,9 @@ Inverse incomplete integral functions :project: gcem .. _iib-ref: -.. doxygenfunction:: incomplete_beta_inv(const pT, const pT, const eT) +.. doxygenfunction:: incomplete_beta_inv(const T1, const T2, const T3) :project: gcem .. _iig-ref: -.. doxygenfunction:: incomplete_gamma_inv(const pT, const eT) +.. doxygenfunction:: incomplete_gamma_inv(const T1, const T2) :project: gcem diff --git a/docs/source/index.rst b/docs/source/index.rst index 6a6785b..739fb5e 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -47,7 +47,7 @@ GCE-Math functions are written as C++ templates with ``constexpr`` specifiers, t template constexpr return_t - erf(const T x); + erf(const T x) noexcept; where a set of internal templated ``constexpr`` functions will implement a continued fraction expansion to return a value of type ``return_t``. This output type ('``return_t``') is generally determined by the input type, e.g., ``int``, ``float``, ``double``, ``long double``, etc. When ``T`` is an intergral type, the output will be upgraded to ``return_t = double``, otherwise ``return_t = T``. For types not covered by ``std::is_integral``, recasts should be used. diff --git a/include/gcem_incl/abs.hpp b/include/gcem_incl/abs.hpp index bec6bda..9f15f27 100644 --- a/include/gcem_incl/abs.hpp +++ b/include/gcem_incl/abs.hpp @@ -32,6 +32,7 @@ template constexpr T abs(const T x) +noexcept { return( x < T(0) ? - x : x ); } diff --git a/include/gcem_incl/acos.hpp b/include/gcem_incl/acos.hpp index d019b6e..34a0f33 100644 --- a/include/gcem_incl/acos.hpp +++ b/include/gcem_incl/acos.hpp @@ -32,6 +32,7 @@ template constexpr T acos_compute(const T x) +noexcept { return( // only defined on [-1,1] abs(x) > T(1) ? \ @@ -49,6 +50,7 @@ template constexpr T acos_check(const T x) +noexcept { return( x > T(0) ? \ // if @@ -70,8 +72,9 @@ template constexpr return_t acos(const T x) +noexcept { - return internal::acos_check>(x); + return internal::acos_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/acosh.hpp b/include/gcem_incl/acosh.hpp index d8b4dfb..42bea73 100644 --- a/include/gcem_incl/acosh.hpp +++ b/include/gcem_incl/acosh.hpp @@ -32,6 +32,7 @@ template constexpr T acosh_compute(const T x) +noexcept { return( // function defined for x >= 1 x < T(1) ? \ @@ -56,8 +57,9 @@ template constexpr return_t acosh(const T x) +noexcept { - return internal::acosh_compute>(x); + return internal::acosh_compute( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/asin.hpp b/include/gcem_incl/asin.hpp index 0a1bd62..885cf76 100644 --- a/include/gcem_incl/asin.hpp +++ b/include/gcem_incl/asin.hpp @@ -32,6 +32,7 @@ template constexpr T asin_compute(const T x) +noexcept { return( // only defined on [-1,1] x > T(1) ? \ @@ -49,6 +50,7 @@ template constexpr T asin_check(const T x) +noexcept { return( x < T(0) ? - asin_compute(-x) : asin_compute(x) ); } @@ -66,8 +68,9 @@ template constexpr return_t asin(const T x) +noexcept { - return internal::asin_check>(x); + return internal::asin_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/asinh.hpp b/include/gcem_incl/asinh.hpp index 32f1a24..36e09f3 100644 --- a/include/gcem_incl/asinh.hpp +++ b/include/gcem_incl/asinh.hpp @@ -32,6 +32,7 @@ template constexpr T asinh_compute(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -53,8 +54,9 @@ template constexpr return_t asinh(const T x) +noexcept { - return internal::asinh_compute>(x); + return internal::asinh_compute( static_cast>(x) ); } diff --git a/include/gcem_incl/atan.hpp b/include/gcem_incl/atan.hpp index cb37d6f..daad65a 100644 --- a/include/gcem_incl/atan.hpp +++ b/include/gcem_incl/atan.hpp @@ -38,6 +38,7 @@ template constexpr T atan_series_order_calc(const T x, const T x_pow, const uint_t order) +noexcept { return( T(1)/( T((order-1)*4 - 1) * x_pow ) \ - T(1)/( T((order-1)*4 + 1) * x_pow*x) ); @@ -47,6 +48,7 @@ template constexpr T atan_series_order(const T x, const T x_pow, const uint_t order, const uint_t max_order) +noexcept { return( order == 1 ? \ GCEM_HALF_PI - T(1)/x + atan_series_order(x*x,pow(x,3),order+1,max_order) : @@ -62,6 +64,7 @@ template constexpr T atan_series_main(const T x) +noexcept { return( x < T(3) ? atan_series_order(x,x,1U,10U) : // O(1/x^39) x < T(4) ? atan_series_order(x,x,1U,9U) : // O(1/x^35) @@ -80,6 +83,7 @@ template constexpr T atan_cf_recur(const T xx, const uint_t depth, const uint_t max_depth) +noexcept { return( depth < max_depth ? \ // if @@ -92,6 +96,7 @@ template constexpr T atan_cf_main(const T x) +noexcept { return( x < T(0.5) ? x/atan_cf_recur(x*x,1U, 15U ) : x < T(1) ? x/atan_cf_recur(x*x,1U, 25U ) : @@ -106,6 +111,7 @@ template constexpr T atan_begin(const T x) +noexcept { return( x > T(2.5) ? atan_series_main(x) : atan_cf_main(x) ); } @@ -114,6 +120,7 @@ template constexpr T atan_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -137,8 +144,9 @@ template constexpr return_t atan(const T x) +noexcept { - return internal::atan_check>(x); + return internal::atan_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/atanh.hpp b/include/gcem_incl/atanh.hpp index 376340d..272c790 100644 --- a/include/gcem_incl/atanh.hpp +++ b/include/gcem_incl/atanh.hpp @@ -32,6 +32,7 @@ template constexpr T atanh_compute(const T x) +noexcept { return( log( (T(1) + x)/(T(1) - x) ) / T(2) ); } @@ -40,6 +41,7 @@ template constexpr T atanh_check(const T x) +noexcept { return( // function is defined for |x| < 1 T(1) < abs(x) ? \ @@ -66,8 +68,9 @@ template constexpr return_t atanh(const T x) +noexcept { - return internal::atanh_check>(x); + return internal::atanh_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/beta.hpp b/include/gcem_incl/beta.hpp index 97e420a..dd2b39b 100644 --- a/include/gcem_incl/beta.hpp +++ b/include/gcem_incl/beta.hpp @@ -30,12 +30,13 @@ * where \f$ \Gamma \f$ denotes the gamma function. */ -template +template constexpr -return_t -beta(const T a, const T b) +common_return_type_t +beta(const T1 a, const T2 b) +noexcept { - return exp( lbeta>(a,b) ); + return exp( lbeta(a,b) ); } #endif diff --git a/include/gcem_incl/binomial_coef.hpp b/include/gcem_incl/binomial_coef.hpp index b2c033a..34f4ade 100644 --- a/include/gcem_incl/binomial_coef.hpp +++ b/include/gcem_incl/binomial_coef.hpp @@ -24,16 +24,44 @@ namespace internal { -template +template constexpr -eT -binomial_coef_recur(const pT n, const pT k) +T +binomial_coef_recur(const T n, const T k) +noexcept { return( // edge cases - (k == pT(0) || n == k) ? eT(1) : // deals with 0 choose 0 case - n == pT(0) ? eT(0) : + (k == T(0) || n == k) ? T(1) : // deals with 0 choose 0 case + n == T(0) ? T(0) : // else - binomial_coef_recur(n-1,k-1) + binomial_coef_recur(n-1,k) ); + binomial_coef_recur(n-1,k-1) + binomial_coef_recur(n-1,k) ); +} + +template::value>::type* = nullptr> +constexpr +T +binomial_coef_check(const T n, const T k) +noexcept +{ + return binomial_coef_recur(n,k); +} + +template::value>::type* = nullptr> +constexpr +T +binomial_coef_sgn_check(const T n, const T k) +noexcept +{ + return binomial_coef_recur(static_cast(n),static_cast(k)); +} + +template> +constexpr +TC +binomial_coef_type_check(const T1 n, const T2 k) +noexcept +{ + return binomial_coef_sgn_check(static_cast(n),static_cast(k)); } } @@ -48,16 +76,13 @@ binomial_coef_recur(const pT n, const pT k) * also known as '\c n choose \c k '. */ -template +template constexpr -eT -binomial_coef(const pT n, const pT k) +common_type_t +binomial_coef(const T1 n, const T2 k) +noexcept { - return( std::is_integral::value ? \ - // if - internal::binomial_coef_recur(n,k) : - // else - internal::binomial_coef_recur(static_cast(n),static_cast(k)) ); + return internal::binomial_coef_type_check(n,k); } #endif \ No newline at end of file diff --git a/include/gcem_incl/cos.hpp b/include/gcem_incl/cos.hpp index 65575c1..b82fa08 100644 --- a/include/gcem_incl/cos.hpp +++ b/include/gcem_incl/cos.hpp @@ -32,6 +32,7 @@ template constexpr T cos_compute(const T x) +noexcept { return( T(1) - x*x)/(T(1) + x*x ); } @@ -40,6 +41,7 @@ template constexpr T cos_check(const T x) +noexcept { return( // indistinguishable from 0 GCLIM::epsilon() > abs(x) ? @@ -70,8 +72,9 @@ template constexpr return_t cos(const T x) +noexcept { - return internal::cos_check>(x); + return internal::cos_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/cosh.hpp b/include/gcem_incl/cosh.hpp index 57e4de9..0331824 100644 --- a/include/gcem_incl/cosh.hpp +++ b/include/gcem_incl/cosh.hpp @@ -32,6 +32,7 @@ template constexpr T cosh_compute(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -53,8 +54,9 @@ template constexpr return_t cosh(const T x) +noexcept { - return internal::cosh_compute>(x); + return internal::cosh_compute( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/erf.hpp b/include/gcem_incl/erf.hpp index 0837073..476cc60 100644 --- a/include/gcem_incl/erf.hpp +++ b/include/gcem_incl/erf.hpp @@ -35,6 +35,7 @@ template constexpr T erf_cf_large_recur(const T x, const int depth) +noexcept { return( depth < GCEM_ERF_MAX_ITER ? \ // if @@ -47,6 +48,7 @@ template constexpr T erf_cf_large_main(const T x) +noexcept { return( T(1) - T(2) * ( exp(-x*x) / T(GCEM_SQRT_PI) ) \ / erf_cf_large_recur(T(2)*x,1) ); @@ -59,6 +61,7 @@ template constexpr T erf_cf_small_recur(const T xx, const int depth) +noexcept { return( depth < GCEM_ERF_MAX_ITER ? \ // if @@ -72,6 +75,7 @@ template constexpr T erf_cf_small_main(const T x) +noexcept { return( T(2) * x * ( exp(-x*x) / T(GCEM_SQRT_PI) ) \ / erf_cf_small_recur(x*x,1) ); @@ -83,6 +87,7 @@ template constexpr T erf_begin(const T x) +noexcept { return( x > T(2.1) ? \ // if @@ -95,6 +100,7 @@ template constexpr T erf_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -121,8 +127,9 @@ template constexpr return_t erf(const T x) +noexcept { - return internal::erf_check>(x); + return internal::erf_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/erf_inv.hpp b/include/gcem_incl/erf_inv.hpp index ed5d2d3..6981adb 100644 --- a/include/gcem_incl/erf_inv.hpp +++ b/include/gcem_incl/erf_inv.hpp @@ -32,7 +32,7 @@ namespace internal { template -constexpr T erf_inv_decision(const T value, const T p, const T direc, const int iter_count); +constexpr T erf_inv_decision(const T value, const T p, const T direc, const int iter_count) noexcept; // // initial value @@ -43,6 +43,7 @@ template constexpr T erf_inv_initial_val_coef_2(const T a, const T p_term, const int order) +noexcept { return( order == 1 ? T(-0.000200214257L) : order == 2 ? T( 0.000100950558L) + a*p_term : @@ -60,6 +61,7 @@ template constexpr T erf_inv_initial_val_case_2(const T a, const T p_term, const int order) +noexcept { return( order == 9 ? \ // if @@ -72,6 +74,7 @@ template constexpr T erf_inv_initial_val_coef_1(const T a, const T p_term, const int order) +noexcept { return( order == 1 ? T( 2.81022636e-08L) : order == 2 ? T( 3.43273939e-07L) + a*p_term : @@ -89,6 +92,7 @@ template constexpr T erf_inv_initial_val_case_1(const T a, const T p_term, const int order) +noexcept { return( order == 9 ? \ // if @@ -101,6 +105,7 @@ template constexpr T erf_inv_initial_val_int(const T a) +noexcept { return( a < T(5) ? \ // if @@ -113,6 +118,7 @@ template constexpr T erf_inv_initial_val(const T x) +noexcept { return x*erf_inv_initial_val_int( -log( (T(1) - x)*(T(1) + x) ) ); } @@ -124,6 +130,7 @@ template constexpr T erf_inv_err_val(const T value, const T p) +noexcept { // err_val = f(x) return( erf(value) - p ); } @@ -132,6 +139,7 @@ template constexpr T erf_inv_deriv_1(const T value) +noexcept { // derivative of the error function w.r.t. x return( exp( -value*value ) ); } @@ -140,6 +148,7 @@ template constexpr T erf_inv_deriv_2(const T value, const T deriv_1) +noexcept { // second derivative of the error function w.r.t. x return( deriv_1*( -T(2)*value ) ); } @@ -148,6 +157,7 @@ template constexpr T erf_inv_ratio_val_1(const T value, const T p, const T deriv_1) +noexcept { return( erf_inv_err_val(value,p) / deriv_1 ); } @@ -156,6 +166,7 @@ template constexpr T erf_inv_ratio_val_2(const T value, const T deriv_1) +noexcept { return( erf_inv_deriv_2(value,deriv_1) / deriv_1 ); } @@ -164,6 +175,7 @@ template constexpr T erf_inv_halley(const T ratio_val_1, const T ratio_val_2) +noexcept { return( ratio_val_1 / max( T(0.8), min( T(1.2), T(1) - T(0.5)*ratio_val_1*ratio_val_2 ) ) ); } @@ -172,6 +184,7 @@ template constexpr T erf_inv_recur(const T value, const T p, const T deriv_1, const int iter_count) +noexcept { return erf_inv_decision( value, p, erf_inv_halley(erf_inv_ratio_val_1(value,p,deriv_1), @@ -183,6 +196,7 @@ template constexpr T erf_inv_decision(const T value, const T p, const T direc, const int iter_count) +noexcept { return( iter_count < GCEM_ERF_INV_MAX_ITER ? \ // if @@ -195,6 +209,7 @@ template constexpr T erf_inv_recur_begin(const T initial_val, const T p) +noexcept { return erf_inv_recur(initial_val,p,erf_inv_deriv_1(initial_val),1); } @@ -203,6 +218,7 @@ template constexpr T erf_inv_begin(const T p) +noexcept { return erf_inv_recur_begin(erf_inv_initial_val(p),p); } @@ -226,8 +242,9 @@ template constexpr return_t erf_inv(const T p) +noexcept { - return internal::erf_inv_begin>(p); + return internal::erf_inv_begin( static_cast>(p) ); } diff --git a/include/gcem_incl/exp.hpp b/include/gcem_incl/exp.hpp index a593f51..ba8f9ef 100644 --- a/include/gcem_incl/exp.hpp +++ b/include/gcem_incl/exp.hpp @@ -32,6 +32,7 @@ template constexpr T exp_cf_recur(const T x, const int depth) +noexcept { return( depth < GCEM_EXP_MAX_ITER_SMALL ? \ // if @@ -46,6 +47,7 @@ template constexpr T exp_cf(const T x) +noexcept { return( T(1)/exp_cf_recur(x,1) ); } @@ -54,6 +56,7 @@ template constexpr T exp_split(const T x) +noexcept { return( pow_integral(GCEM_E,find_whole(x)) * exp_cf(find_fraction(x)) ); } @@ -62,6 +65,7 @@ template constexpr T exp_check(const T x) +noexcept { return( GCLIM::epsilon() > abs(x) ? \ // if @@ -85,8 +89,9 @@ template constexpr return_t exp(const T x) +noexcept { - return internal::exp_check>(x); + return internal::exp_check( static_cast>(x) ); } // #ifndef GCEM_EXP_TOL diff --git a/include/gcem_incl/factorial.hpp b/include/gcem_incl/factorial.hpp index ce8427a..ab52d9e 100644 --- a/include/gcem_incl/factorial.hpp +++ b/include/gcem_incl/factorial.hpp @@ -34,6 +34,7 @@ template constexpr T factorial_table(const T x) +noexcept { // table for x! when x = {2,...,16} return( x == T(2) ? T(2) : x == T(3) ? T(6) : x == T(4) ? T(24) : x == T(5) ? T(120) : @@ -53,6 +54,7 @@ template::value>::type* constexpr T factorial_recur(const T x) +noexcept { return( x == T(0) ? T(1) : x == T(1) ? x : @@ -68,6 +70,7 @@ template::value>::type* constexpr T factorial_recur(const T x) +noexcept { return tgamma(x + 1); } @@ -87,6 +90,7 @@ template constexpr T factorial(const T x) +noexcept { return internal::factorial_recur(x); } diff --git a/include/gcem_incl/find_exponent.hpp b/include/gcem_incl/find_exponent.hpp index bcc1656..1ac2902 100644 --- a/include/gcem_incl/find_exponent.hpp +++ b/include/gcem_incl/find_exponent.hpp @@ -29,11 +29,12 @@ template constexpr llint_t find_exponent(const T x, const llint_t exponent) +noexcept { return( x < T(1) ? \ - find_exponent(x*T(10),exponent-1) : + find_exponent(x*T(10),exponent - llint_t(1)) : x > T(10) ? \ - find_exponent(x/T(10),exponent+1) : + find_exponent(x/T(10),exponent + llint_t(1)) : // else exponent ); } diff --git a/include/gcem_incl/find_fraction.hpp b/include/gcem_incl/find_fraction.hpp index 086b012..f8313aa 100644 --- a/include/gcem_incl/find_fraction.hpp +++ b/include/gcem_incl/find_fraction.hpp @@ -29,6 +29,7 @@ template constexpr T find_fraction(const T x) +noexcept { return( abs(x - internal::floor(x)) > T(0.5) ? \ // if diff --git a/include/gcem_incl/find_whole.hpp b/include/gcem_incl/find_whole.hpp index 4e48615..76d4e70 100644 --- a/include/gcem_incl/find_whole.hpp +++ b/include/gcem_incl/find_whole.hpp @@ -29,6 +29,7 @@ template constexpr llint_t find_whole(const T x) +noexcept { return( abs(x - internal::floor(x)) > T(0.5) ? \ // if diff --git a/include/gcem_incl/floor.hpp b/include/gcem_incl/floor.hpp index 12f24d8..42ee8b8 100644 --- a/include/gcem_incl/floor.hpp +++ b/include/gcem_incl/floor.hpp @@ -28,6 +28,7 @@ template constexpr T floor(const T x) +noexcept { return T(static_cast(x)); } diff --git a/include/gcem_incl/gcd.hpp b/include/gcem_incl/gcd.hpp index fd60867..7dfc43e 100644 --- a/include/gcem_incl/gcd.hpp +++ b/include/gcem_incl/gcd.hpp @@ -28,8 +28,36 @@ template constexpr T gcd_recur(const T a, const T b) +noexcept { - return b == T(0) ? a : gcd_recur(b, a % b); + return( b == T(0) ? a : gcd_recur(b, a % b) ); +} + +template::value>::type* = nullptr> +constexpr +T +gcd_int_check(const T a, const T b) +noexcept +{ + return gcd_recur(a,b); +} + +template::value>::type* = nullptr> +constexpr +T +gcd_int_check(const T a, const T b) +noexcept +{ + return gcd_recur( static_cast(a), static_cast(b) ); +} + +template> +constexpr +TC +gcd_type_check(const T1 a, const T2 b) +noexcept +{ + return gcd_int_check( static_cast(abs(a)), static_cast(abs(b)) ); } } @@ -42,18 +70,13 @@ gcd_recur(const T a, const T b) * @return the greatest common divisor between integers \c a and \c b using a Euclidean algorithm. */ -template +template constexpr -T -gcd(const T a, const T b) +common_type_t +gcd(const T1 a, const T2 b) +noexcept { - return( (a < T(0) || b < T(0)) ? \ - // sanity check - gcd(abs(a),abs(b)) : - // - std::is_integral::value ? \ - internal::gcd_recur(a,b) : - internal::gcd_recur(a,b) ); + return internal::gcd_type_check(a,b); } #endif diff --git a/include/gcem_incl/gcem_options.hpp b/include/gcem_incl/gcem_options.hpp index 833d29a..27caf32 100644 --- a/include/gcem_incl/gcem_options.hpp +++ b/include/gcem_incl/gcem_options.hpp @@ -30,7 +30,7 @@ #endif #ifndef GCEM_VERSION_MINOR - #define GCEM_VERSION_MINOR 6 + #define GCEM_VERSION_MINOR 7 #endif #ifndef GCEM_VERSION_PATCH @@ -52,6 +52,14 @@ namespace gcem template using return_t = typename std::conditional::value,double,T>::type; + +#if(__cplusplus < 201402L) + template + using common_type_t = typename std::common_type::type; +#endif + + template + using common_return_type_t = return_t>; } // diff --git a/include/gcem_incl/incomplete_beta.hpp b/include/gcem_incl/incomplete_beta.hpp index 0cbee26..924fcfe 100644 --- a/include/gcem_incl/incomplete_beta.hpp +++ b/include/gcem_incl/incomplete_beta.hpp @@ -31,7 +31,7 @@ namespace internal { template -constexpr T incomplete_beta_cf(const T a, const T b, const T z, const T c_j, const T d_j, const T f_j, const int depth); +constexpr T incomplete_beta_cf(const T a, const T b, const T z, const T c_j, const T d_j, const T f_j, const int depth) noexcept; // // coefficients; see eq. 18.5.17b @@ -40,6 +40,7 @@ template constexpr T incomplete_beta_coef_even(const T a, const T b, const T z, const int k) +noexcept { return( -z*(a + k)*(a + b + k)/( (a + 2*k)*(a + 2*k + T(1)) ) ); } @@ -48,6 +49,7 @@ template constexpr T incomplete_beta_coef_odd(const T a, const T b, const T z, const int k) +noexcept { return( z*k*(b - k)/((a + 2*k - T(1))*(a + 2*k)) ); } @@ -56,9 +58,10 @@ template constexpr T incomplete_beta_coef(const T a, const T b, const T z, const int depth) +noexcept { return( !is_odd(depth) ? incomplete_beta_coef_even(a,b,z,depth/2) : - incomplete_beta_coef_odd(a,b,z,(depth+1)/2) ); + incomplete_beta_coef_odd(a,b,z,(depth+1)/2) ); } // @@ -68,6 +71,7 @@ template constexpr T incomplete_beta_c_update(const T a, const T b, const T z, const T c_j, const int depth) +noexcept { return( T(1) + incomplete_beta_coef(a,b,z,depth)/c_j ); } @@ -76,6 +80,7 @@ template constexpr T incomplete_beta_d_update(const T a, const T b, const T z, const T d_j, const int depth) +noexcept { return( T(1) / (T(1) + incomplete_beta_coef(a,b,z,depth)*d_j) ); } @@ -87,6 +92,7 @@ template constexpr T incomplete_beta_decision(const T a, const T b, const T z, const T c_j, const T d_j, const T f_j, const int depth) +noexcept { return( // tolerance check abs(c_j*d_j - T(1)) < GCEM_INCML_BETA_TOL ? f_j*c_j*d_j : @@ -102,6 +108,7 @@ template constexpr T incomplete_beta_cf(const T a, const T b, const T z, const T c_j, const T d_j, const T f_j, const int depth) +noexcept { return incomplete_beta_decision(a,b,z, incomplete_beta_c_update(a,b,z,c_j,depth), @@ -116,6 +123,7 @@ template constexpr T incomplete_beta_begin(const T a, const T b, const T z) +noexcept { return ( (exp(a*log(z) + b*log(T(1)-z) - lbeta(a,b)) / a) * \ incomplete_beta_cf(a,b,z,T(1), @@ -128,6 +136,7 @@ template constexpr T incomplete_beta_check(const T a, const T b, const T z) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > z ? \ @@ -138,6 +147,17 @@ incomplete_beta_check(const T a, const T b, const T z) T(1) - incomplete_beta_begin(b,a,T(1) - z) ); } +template> +constexpr +TC +incomplete_beta_type_check(const T1 a, const T2 b, const T3 p) +noexcept +{ + return incomplete_beta_check(static_cast(a), + static_cast(b), + static_cast(p)); +} + } /** @@ -159,12 +179,13 @@ incomplete_beta_check(const T a, const T b, const T z) * \f[ f_j = c_j d_j f_{j-1} \f] */ -template +template constexpr -eT -incomplete_beta(const pT a, const pT b, const eT z) +common_return_type_t +incomplete_beta(const T1 a, const T2 b, const T3 z) +noexcept { - return internal::incomplete_beta_check(a,b,z); + return internal::incomplete_beta_type_check(a,b,z); } #endif diff --git a/include/gcem_incl/incomplete_beta_inv.hpp b/include/gcem_incl/incomplete_beta_inv.hpp index 3fe1364..6d69cee 100644 --- a/include/gcem_incl/incomplete_beta_inv.hpp +++ b/include/gcem_incl/incomplete_beta_inv.hpp @@ -30,7 +30,7 @@ namespace internal template constexpr T incomplete_beta_inv_decision(const T value, const T alpha_par, const T beta_par, const T p, - const T direc, const T lb_val, const int iter_count); + const T direc, const T lb_val, const int iter_count) noexcept; // // initial value for Halley @@ -42,6 +42,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_tval(const T p) +noexcept { // a > 1.0 return( p > T(0.5) ? \ // if @@ -54,6 +55,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_begin(const T t_val) +noexcept { // internal for a > 1.0 return( t_val - ( T(2.515517) + T(0.802853)*t_val + T(0.010328)*t_val*t_val ) \ / ( T(1) + T(1.432788)*t_val + T(0.189269)*t_val*t_val + T(0.001308)*t_val*t_val*t_val ) ); @@ -63,6 +65,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_ab1(const T alpha_par, const T beta_par) +noexcept { return( T(1)/(2*alpha_par - T(1)) + T(1)/(2*beta_par - T(1)) ); } @@ -71,6 +74,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_ab2(const T alpha_par, const T beta_par) +noexcept { return( T(1)/(2*beta_par - T(1)) - T(1)/(2*alpha_par - T(1)) ); } @@ -79,6 +83,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_h(const T ab_term_1) +noexcept { return( T(2) / ab_term_1 ); } @@ -87,6 +92,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_w(const T value, const T ab_term_2, const T h_term) +noexcept { // return( value * sqrt(h_term + lambda)/h_term - ab_term_2*(lambda + 5.0/6.0 -2.0/(3.0*h_term)) ); return( value * sqrt(h_term + (value*value - T(3))/T(6))/h_term \ @@ -97,6 +103,7 @@ template constexpr T incomplete_beta_inv_initial_val_1_int_end(const T alpha_par, const T beta_par, const T w_term) +noexcept { return( alpha_par / (alpha_par + beta_par*exp(2*w_term)) ); } @@ -105,6 +112,7 @@ template constexpr T incomplete_beta_inv_initial_val_1(const T alpha_par, const T beta_par, const T t_val, const T sgn_term) +noexcept { // a > 1.0 return incomplete_beta_inv_initial_val_1_int_end( alpha_par, beta_par, incomplete_beta_inv_initial_val_1_int_w( @@ -124,6 +132,7 @@ template constexpr T incomplete_beta_inv_initial_val_2_s1(const T alpha_par, const T beta_par) +noexcept { return( pow(alpha_par/(alpha_par+beta_par),alpha_par) / alpha_par ); } @@ -132,6 +141,7 @@ template constexpr T incomplete_beta_inv_initial_val_2_s2(const T alpha_par, const T beta_par) +noexcept { return( pow(beta_par/(alpha_par+beta_par),beta_par) / beta_par ); } @@ -140,6 +150,7 @@ template constexpr T incomplete_beta_inv_initial_val_2(const T alpha_par, const T beta_par, const T p, const T s_1, const T s_2) +noexcept { return( p <= s_1/(s_1 + s_2) ? pow(p*(s_1+s_2)*alpha_par,T(1)/alpha_par) : T(1) - pow(p*(s_1+s_2)*beta_par,T(1)/beta_par) ); @@ -151,6 +162,7 @@ template constexpr T incomplete_beta_inv_initial_val(const T alpha_par, const T beta_par, const T p) +noexcept { return( (alpha_par > T(1) && beta_par > T(1)) ? // if @@ -177,6 +189,7 @@ template constexpr T incomplete_beta_inv_err_val(const T value, const T alpha_par, const T beta_par, const T p) +noexcept { // err_val = f(x) return( incomplete_beta(alpha_par,beta_par,value) - p ); } @@ -185,6 +198,7 @@ template constexpr T incomplete_beta_inv_deriv_1(const T value, const T alpha_par, const T beta_par, const T lb_val) +noexcept { // derivative of the incomplete beta function w.r.t. x return( // indistinguishable from zero or one GCLIM::epsilon() > abs(value) ? \ @@ -199,6 +213,7 @@ template constexpr T incomplete_beta_inv_deriv_2(const T value, const T alpha_par, const T beta_par, const T deriv_1) +noexcept { // second derivative of the incomplete beta function w.r.t. x return( deriv_1*((alpha_par - T(1))/value - (beta_par - T(1))/(T(1) - value)) ); } @@ -207,6 +222,7 @@ template constexpr T incomplete_beta_inv_ratio_val_1(const T value, const T alpha_par, const T beta_par, const T p, const T deriv_1) +noexcept { return( incomplete_beta_inv_err_val(value,alpha_par,beta_par,p) / deriv_1 ); } @@ -215,6 +231,7 @@ template constexpr T incomplete_beta_inv_ratio_val_2(const T value, const T alpha_par, const T beta_par, const T deriv_1) +noexcept { return( incomplete_beta_inv_deriv_2(value,alpha_par,beta_par,deriv_1) / deriv_1 ); } @@ -223,6 +240,7 @@ template constexpr T incomplete_beta_inv_halley(const T ratio_val_1, const T ratio_val_2) +noexcept { return( ratio_val_1 / max( T(0.8), min( T(1.2), T(1) - T(0.5)*ratio_val_1*ratio_val_2 ) ) ); } @@ -232,6 +250,7 @@ constexpr T incomplete_beta_inv_recur(const T value, const T alpha_par, const T beta_par, const T p, const T deriv_1, const T lb_val, const int iter_count) +noexcept { return( // derivative = 0 GCLIM::epsilon() > abs(deriv_1) ? \ @@ -250,6 +269,7 @@ constexpr T incomplete_beta_inv_decision(const T value, const T alpha_par, const T beta_par, const T p, const T direc, const T lb_val, const int iter_count) +noexcept { return( iter_count <= GCEM_INCML_BETA_INV_MAX_ITER ? // if @@ -264,6 +284,7 @@ template constexpr T incomplete_beta_inv_begin(const T initial_val, const T alpha_par, const T beta_par, const T p, const T lb_val) +noexcept { return incomplete_beta_inv_recur(initial_val,alpha_par,beta_par,p, incomplete_beta_inv_deriv_1(initial_val,alpha_par,beta_par,lb_val), @@ -274,6 +295,7 @@ template constexpr T incomplete_beta_inv_check(const T alpha_par, const T beta_par, const T p) +noexcept { return( // indistinguishable from zero or one GCLIM::epsilon() > p ? \ @@ -285,6 +307,17 @@ incomplete_beta_inv_check(const T alpha_par, const T beta_par, const T p) alpha_par,beta_par,p,lbeta(alpha_par,beta_par)) ); } +template> +constexpr +TC +incomplete_beta_inv_type_check(const T1 a, const T2 b, const T3 p) +noexcept +{ + return incomplete_beta_inv_check(static_cast(a), + static_cast(b), + static_cast(p)); +} + } /** @@ -304,12 +337,13 @@ incomplete_beta_inv_check(const T alpha_par, const T beta_par, const T p) * \f[ \frac{\partial^2}{\partial x^2} \left(\frac{\text{B}(x;\alpha,\beta)}{\text{B}(\alpha,\beta)}\right) = \frac{1}{\text{B}(\alpha,\beta)} x^{\alpha-1} (1-x)^{\beta-1} \left( \frac{\alpha-1}{x} - \frac{\beta-1}{1 - x} \right) \f] */ -template +template constexpr -eT -incomplete_beta_inv(const pT a, const pT b, const eT p) +common_type_t +incomplete_beta_inv(const T1 a, const T2 b, const T3 p) +noexcept { - return internal::incomplete_beta_inv_check(a,b,p); + return internal::incomplete_beta_inv_type_check(a,b,p); } #endif diff --git a/include/gcem_incl/incomplete_gamma.hpp b/include/gcem_incl/incomplete_gamma.hpp index 67110a1..7125a04 100644 --- a/include/gcem_incl/incomplete_gamma.hpp +++ b/include/gcem_incl/incomplete_gamma.hpp @@ -34,6 +34,7 @@ template constexpr T incomplete_gamma_quad_inp_vals(const T lb, const T ub, const int counter) +noexcept { return (ub-lb) * gauss_legendre_50_points[counter] / T(2) + (ub + lb) / T(2); } @@ -42,6 +43,7 @@ template constexpr T incomplete_gamma_quad_weight_vals(const T lb, const T ub, const int counter) +noexcept { return (ub-lb) * gauss_legendre_50_weights[counter] / T(2); } @@ -50,6 +52,7 @@ template constexpr T incomplete_gamma_quad_fn(const T x, const T a, const T lg_term) +noexcept { return exp( -x + (a-T(1))*log(x) - lg_term ); } @@ -58,6 +61,7 @@ template constexpr T incomplete_gamma_quad_recur(const T lb, const T ub, const T a, const T lg_term, const int counter) +noexcept { return( counter < 49 ? \ // if @@ -73,6 +77,7 @@ template constexpr T incomplete_gamma_quad_lb(const T a, const T z) +noexcept { return( a > T(1000) ? max(T(0),min(z,a) - 11*sqrt(a)) : // break integration into ranges a > T(800) ? max(T(0),min(z,a) - 11*sqrt(a)) : @@ -92,6 +97,7 @@ template constexpr T incomplete_gamma_quad_ub(const T a, const T z) +noexcept { return( a > T(1000) ? min(z, a + 10*sqrt(a)) : a > T(800) ? min(z, a + 10*sqrt(a)) : @@ -110,6 +116,7 @@ template constexpr T incomplete_gamma_quad(const T a, const T z) +noexcept { return incomplete_gamma_quad_recur(incomplete_gamma_quad_lb(a,z), incomplete_gamma_quad_ub(a,z), a,lgamma(a),0); } @@ -121,6 +128,7 @@ template constexpr T incomplete_gamma_cf_coef(const T a, const T z, const int depth) +noexcept { return( is_odd(depth) ? - (a - 1 + T(depth+1)/T(2)) * z : T(depth)/T(2) * z ); } @@ -129,6 +137,7 @@ template constexpr T incomplete_gamma_cf_recur(const T a, const T z, const int depth) +noexcept { return( depth < GCEM_INCML_GAMMA_MAX_ITER ? \ // if @@ -141,6 +150,7 @@ template constexpr T incomplete_gamma_cf(const T a, const T z) +noexcept { // lower (regularized) incomplete gamma function return( exp(a*log(z) - z) / tgamma(a) / incomplete_gamma_cf_recur(a,z,1) ); } @@ -151,6 +161,7 @@ template constexpr T incomplete_gamma_check(const T a, const T z) +noexcept { return( GCLIM::epsilon() > z ? \ T(0) : @@ -163,6 +174,16 @@ incomplete_gamma_check(const T a, const T z) incomplete_gamma_quad(a,z) ); } +template> +constexpr +TC +incomplete_gamma_type_check(const T1 a, const T2 p) +noexcept +{ + return incomplete_gamma_check(static_cast(a), + static_cast(p)); +} + } /** @@ -180,12 +201,13 @@ incomplete_gamma_check(const T a, const T z) * When \f$ a > 10 \f$, a 50-point Gauss-Legendre quadrature scheme is employed. */ -template +template constexpr -return_t -incomplete_gamma(const pT a, const eT x) +common_return_type_t +incomplete_gamma(const T1 a, const T2 x) +noexcept { - return internal::incomplete_gamma_check>(a,x); + return internal::incomplete_gamma_type_check(a,x); } #endif diff --git a/include/gcem_incl/incomplete_gamma_inv.hpp b/include/gcem_incl/incomplete_gamma_inv.hpp index 2febd97..51070f5 100644 --- a/include/gcem_incl/incomplete_gamma_inv.hpp +++ b/include/gcem_incl/incomplete_gamma_inv.hpp @@ -29,7 +29,7 @@ namespace internal { template -constexpr T incomplete_gamma_inv_decision(const T value, const T a, const T p, const T direc, const T lg_val, const int iter_count); +constexpr T incomplete_gamma_inv_decision(const T value, const T a, const T p, const T direc, const T lg_val, const int iter_count) noexcept; // // initial value for Halley @@ -38,6 +38,7 @@ template constexpr T incomplete_gamma_inv_t_val_1(const T p) +noexcept { // a > 1.0 return( p > T(0.5) ? sqrt(-2*log(T(1) - p)) : sqrt(-2*log(p)) ); } @@ -46,6 +47,7 @@ template constexpr T incomplete_gamma_inv_t_val_2(const T a) +noexcept { // a <= 1.0 return( T(1) - T(0.253) * a - T(0.12) * a*a ); } @@ -56,6 +58,7 @@ template constexpr T incomplete_gamma_inv_initial_val_1_int_begin(const T t_val) +noexcept { // internal for a > 1.0 return( t_val - (T(2.515517L) + T(0.802853L)*t_val + T(0.010328L)*t_val*t_val) \ / (T(1) + T(1.432788L)*t_val + T(0.189269L)*t_val*t_val + T(0.001308L)*t_val*t_val*t_val) ); @@ -65,6 +68,7 @@ template constexpr T incomplete_gamma_inv_initial_val_1_int_end(const T value_inp, const T a) +noexcept { // internal for a > 1.0 return max( T(1E-04), a*pow(T(1) - T(1)/(9*a) - value_inp/(3*sqrt(a)), 3) ); } @@ -73,6 +77,7 @@ template constexpr T incomplete_gamma_inv_initial_val_1(const T a, const T t_val, const T sgn_term) +noexcept { // a > 1.0 return incomplete_gamma_inv_initial_val_1_int_end(sgn_term*incomplete_gamma_inv_initial_val_1_int_begin(t_val), a); } @@ -81,6 +86,7 @@ template constexpr T incomplete_gamma_inv_initial_val_2(const T a, const T p, const T t_val) +noexcept { // a <= 1.0 return( p < t_val ? \ // if @@ -95,6 +101,7 @@ template constexpr T incomplete_gamma_inv_initial_val(const T a, const T p) +noexcept { return( a > T(1) ? \ // if @@ -113,6 +120,7 @@ template constexpr T incomplete_gamma_inv_err_val(const T value, const T a, const T p) +noexcept { // err_val = f(x) return( incomplete_gamma(a,value) - p ); } @@ -121,6 +129,7 @@ template constexpr T incomplete_gamma_inv_deriv_1(const T value, const T a, const T lg_val) +noexcept { // derivative of the incomplete gamma function w.r.t. x return( exp( - value + (a - T(1))*log(value) - lg_val ) ); } @@ -129,6 +138,7 @@ template constexpr T incomplete_gamma_inv_deriv_2(const T value, const T a, const T deriv_1) +noexcept { // second derivative of the incomplete gamma function w.r.t. x return( deriv_1*((a - T(1))/value - T(1)) ); } @@ -137,6 +147,7 @@ template constexpr T incomplete_gamma_inv_ratio_val_1(const T value, const T a, const T p, const T deriv_1) +noexcept { return( incomplete_gamma_inv_err_val(value,a,p) / deriv_1 ); } @@ -145,6 +156,7 @@ template constexpr T incomplete_gamma_inv_ratio_val_2(const T value, const T a, const T deriv_1) +noexcept { return( incomplete_gamma_inv_deriv_2(value,a,deriv_1) / deriv_1 ); } @@ -153,6 +165,7 @@ template constexpr T incomplete_gamma_inv_halley(const T ratio_val_1, const T ratio_val_2) +noexcept { return( ratio_val_1 / max( T(0.8), min( T(1.2), T(1) - T(0.5)*ratio_val_1*ratio_val_2 ) ) ); } @@ -161,6 +174,7 @@ template constexpr T incomplete_gamma_inv_recur(const T value, const T a, const T p, const T deriv_1, const T lg_val, const int iter_count) +noexcept { return incomplete_gamma_inv_decision(value, a, p, incomplete_gamma_inv_halley(incomplete_gamma_inv_ratio_val_1(value,a,p,deriv_1), @@ -172,6 +186,7 @@ template constexpr T incomplete_gamma_inv_decision(const T value, const T a, const T p, const T direc, const T lg_val, const int iter_count) +noexcept { // return( abs(direc) > GCEM_INCML_GAMMA_INV_TOL ? incomplete_gamma_inv_recur(value - direc, a, p, incomplete_gamma_inv_deriv_1(value,a,lg_val), lg_val) : value - direc ); return( iter_count <= GCEM_INCML_GAMMA_INV_MAX_ITER ? \ @@ -187,6 +202,7 @@ template constexpr T incomplete_gamma_inv_begin(const T initial_val, const T a, const T p, const T lg_val) +noexcept { return incomplete_gamma_inv_recur(initial_val,a,p, incomplete_gamma_inv_deriv_1(initial_val,a,lg_val), lg_val,1); @@ -196,6 +212,7 @@ template constexpr T incomplete_gamma_inv_check(const T a, const T p) +noexcept { return( GCLIM::epsilon() > p ? \ T(0) : @@ -210,6 +227,16 @@ incomplete_gamma_inv_check(const T a, const T p) incomplete_gamma_inv_begin(incomplete_gamma_inv_initial_val(a,p),a,p,lgamma(a)) ); } +template> +constexpr +TC +incomplete_gamma_inv_type_check(const T1 a, const T2 p) +noexcept +{ + return incomplete_gamma_inv_check(static_cast(a), + static_cast(p)); +} + } /** @@ -228,12 +255,13 @@ incomplete_gamma_inv_check(const T a, const T p) * \f[ \frac{\partial^2}{\partial x^2} \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1} \exp(-x) \left( \frac{a-1}{x} - 1 \right) \f] */ -template +template constexpr -eT -incomplete_gamma_inv(const pT a, const eT p) +common_return_type_t +incomplete_gamma_inv(const T1 a, const T2 p) +noexcept { - return internal::incomplete_gamma_inv_check(a,p); + return internal::incomplete_gamma_inv_type_check(a,p); } #endif diff --git a/include/gcem_incl/is_even.hpp b/include/gcem_incl/is_even.hpp index ca90689..3ebab5e 100644 --- a/include/gcem_incl/is_even.hpp +++ b/include/gcem_incl/is_even.hpp @@ -28,6 +28,7 @@ constexpr bool is_even(const llint_t x) +noexcept { return !is_odd(x); } diff --git a/include/gcem_incl/is_odd.hpp b/include/gcem_incl/is_odd.hpp index e03c683..11b036c 100644 --- a/include/gcem_incl/is_odd.hpp +++ b/include/gcem_incl/is_odd.hpp @@ -28,6 +28,7 @@ constexpr bool is_odd(const llint_t x) +noexcept { return( x % llint_t(2) == llint_t(0) ? false : true ); } diff --git a/include/gcem_incl/lbeta.hpp b/include/gcem_incl/lbeta.hpp index b886e7b..d2eb777 100644 --- a/include/gcem_incl/lbeta.hpp +++ b/include/gcem_incl/lbeta.hpp @@ -30,12 +30,13 @@ * where \f$ \Gamma \f$ denotes the gamma function. */ -template +template constexpr -T -lbeta(const T a, const T b) +common_return_type_t +lbeta(const T1 a, const T2 b) +noexcept { - return( lgamma(a) + lgamma(b) - lgamma(a+b) ); + return( (lgamma(a) + lgamma(b)) - lgamma(a+b) ); } #endif diff --git a/include/gcem_incl/lcm.hpp b/include/gcem_incl/lcm.hpp index eaa23dc..1861e10 100644 --- a/include/gcem_incl/lcm.hpp +++ b/include/gcem_incl/lcm.hpp @@ -21,6 +21,29 @@ #ifndef _gcem_lcm_HPP #define _gcem_lcm_HPP +namespace internal +{ + +template +constexpr +T +lcm_compute(const T a, const T b) +noexcept +{ + return abs(a * (b / gcd(a,b))); +} + +template> +constexpr +TC +lcm_type_check(const T1 a, const T2 b) +noexcept +{ + return lcm_compute(static_cast(a),static_cast(b)); +} + +} + /** * Compile-time least common multiple (LCM) function * @@ -30,12 +53,13 @@ * where \f$ \text{gcd}(a,b) \f$ denotes the greatest common divisor between \f$ a \f$ and \f$ b \f$. */ -template +template constexpr -T -lcm(const T a, const T b) +common_type_t +lcm(const T1 a, const T2 b) +noexcept { - return abs(a * (b / gcd(a,b))); + return internal::lcm_type_check(a,b); } #endif diff --git a/include/gcem_incl/lgamma.hpp b/include/gcem_incl/lgamma.hpp index 0bda161..84b1416 100644 --- a/include/gcem_incl/lgamma.hpp +++ b/include/gcem_incl/lgamma.hpp @@ -52,6 +52,7 @@ namespace internal constexpr long double lgamma_coef_term(const long double x) +noexcept { return( 0.99999999999999709182L + 57.156235665862923517L / (x+1) \ - 59.597960355475491248L / (x+2) + 14.136097974741747174L / (x+3) \ @@ -67,6 +68,7 @@ template constexpr T lgamma_term_2(const T x) +noexcept { // return( T(GCEM_LOG_SQRT_2PI) + log(T(lgamma_coef_term(x))) ); } @@ -75,6 +77,7 @@ template constexpr T lgamma_term_1(const T x) +noexcept { // note: 607/128 + 0.5 = 5.2421875 return( (x + T(0.5))*log(x + T(5.2421875L)) - (x + T(5.2421875L)) ); } @@ -83,6 +86,7 @@ template constexpr T lgamma_begin(const T x) +noexcept { // returns lngamma(x+1) return( lgamma_term_1(x) + lgamma_term_2(x) ); } @@ -91,6 +95,7 @@ template constexpr T lgamma_check(const T x) +noexcept { return( // indistinguishable from one or <= zero GCLIM::epsilon() > abs(x - T(1)) ? \ @@ -119,8 +124,9 @@ template constexpr return_t lgamma(const T x) +noexcept { - return internal::lgamma_check>(x); + return internal::lgamma_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/lmgamma.hpp b/include/gcem_incl/lmgamma.hpp index 958e778..fc55b1b 100644 --- a/include/gcem_incl/lmgamma.hpp +++ b/include/gcem_incl/lmgamma.hpp @@ -30,18 +30,19 @@ namespace internal // see https://en.wikipedia.org/wiki/Multivariate_gamma_function -template +template constexpr -Ta -lmgamma_recur(const Ta a, const Tb p) +T1 +lmgamma_recur(const T1 a, const T2 p) +noexcept { - return( p == Tb(1) ? \ + return( p == T2(1) ? \ lgamma(a) : - p < Tb(1) ? \ - GCLIM::quiet_NaN() : + p < T2(1) ? \ + GCLIM::quiet_NaN() : // else - Ta(GCEM_LOG_PI) * (p - Ta(1))/Ta(2) \ - + lgamma(a) + lmgamma_recur(a - Ta(0.5),p-1) ); + T1(GCEM_LOG_PI) * (p - T1(1))/T1(2) \ + + lgamma(a) + lmgamma_recur(a - T1(0.5),p-T2(1)) ); } } @@ -56,12 +57,13 @@ lmgamma_recur(const Ta a, const Tb p) * where \f$ \Gamma_1(a) = \Gamma(a) \f$. */ -template +template constexpr -return_t -lmgamma(const pT a, const eT p) +return_t +lmgamma(const T1 a, const T2 p) +noexcept { - return internal::lmgamma_recur(return_t(a),p); + return internal::lmgamma_recur(static_cast>(a),p); } #endif diff --git a/include/gcem_incl/log.hpp b/include/gcem_incl/log.hpp index 8c059f7..a611f6c 100644 --- a/include/gcem_incl/log.hpp +++ b/include/gcem_incl/log.hpp @@ -35,6 +35,7 @@ template constexpr T log_cf_main(const T xx, const int depth) +noexcept { return( depth < GCEM_LOG_MAX_ITER_SMALL ? \ // if @@ -47,6 +48,7 @@ template constexpr T log_cf_begin(const T x) +noexcept { return( T(2)*x/log_cf_main(x*x,1) ); } @@ -55,6 +57,7 @@ template constexpr T log_main(const T x) +noexcept { return( log_cf_begin((x - T(1))/(x + T(1))) ); } @@ -62,6 +65,7 @@ log_main(const T x) constexpr long double log_mantissa_integer(const int x) +noexcept { return( x == 2 ? 0.6931471805599453094172321214581765680755L : x == 3 ? 1.0986122886681096913952452369225257046475L : @@ -79,6 +83,7 @@ template constexpr T log_mantissa(const T x) +noexcept { // divide by the integer part of x, which will be in [1,10], then adjust using tables return( log_main(x/T(static_cast(x))) + T(log_mantissa_integer(static_cast(x))) ); } @@ -87,6 +92,7 @@ template constexpr T log_breakup(const T x) +noexcept { // x = a*b, where b = 10^c return( log_mantissa(mantissa(x)) + T(GCEM_LOG_10)*T(find_exponent(x,0)) ); } @@ -95,6 +101,7 @@ template constexpr T log_check(const T x) +noexcept { return( // x < 0 x < T(0) ? \ @@ -127,8 +134,9 @@ template constexpr return_t log(const T x) +noexcept { - return internal::log_check>(x); + return internal::log_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/mantissa.hpp b/include/gcem_incl/mantissa.hpp index c89efa1..349e3cd 100644 --- a/include/gcem_incl/mantissa.hpp +++ b/include/gcem_incl/mantissa.hpp @@ -29,6 +29,7 @@ template constexpr T mantissa(const T x) +noexcept { return( x < T(1) ? \ mantissa(x*T(10)) : diff --git a/include/gcem_incl/max.hpp b/include/gcem_incl/max.hpp index eaf3f5d..50c9ec7 100644 --- a/include/gcem_incl/max.hpp +++ b/include/gcem_incl/max.hpp @@ -22,19 +22,20 @@ #define _gcem_max_HPP /** - * Compile-time maximum function + * Compile-time pairwise maximum function * * @param x a real-valued input. * @param y a real-valued input. * @return Computes the maximum between \c x and \c y, where \c x and \c y have the same type (e.g., \c int, \c double, etc.) */ -template +template constexpr -T -max(const T x, const T y) +common_type_t +max(const T1 x, const T2 y) +noexcept { - return( x > y ? x : y ); + return( y < x ? x : y ); } #endif diff --git a/include/gcem_incl/min.hpp b/include/gcem_incl/min.hpp index e0efe47..0053faa 100644 --- a/include/gcem_incl/min.hpp +++ b/include/gcem_incl/min.hpp @@ -22,19 +22,20 @@ #define _gcem_min_HPP /** - * Compile-time minimum function + * Compile-time pairwise minimum function * * @param x a real-valued input. * @param y a real-valued input. * @return Computes the minimum between \c x and \c y, where \c x and \c y have the same type (e.g., \c int, \c double, etc.) */ -template +template constexpr -T -min(const T x, const T y) +common_type_t +min(const T1 x, const T2 y) +noexcept { - return( x < y ? x : y ); + return( y > x ? x : y ); } #endif diff --git a/include/gcem_incl/pow.hpp b/include/gcem_incl/pow.hpp index 9fcea2c..c5da8b3 100644 --- a/include/gcem_incl/pow.hpp +++ b/include/gcem_incl/pow.hpp @@ -25,14 +25,43 @@ #ifndef _gcem_pow_HPP #define _gcem_pow_HPP +namespace internal +{ + template constexpr T pow_dbl(const T base, const T exp_term) +noexcept { return exp(exp_term*log(base)); } +template, + typename std::enable_if::value>::type* = nullptr> +constexpr +TC +pow_check(const T1 base, const T2 exp_term) +noexcept +{ + return( base < T1(0) ? \ + GCLIM::quiet_NaN() : + // + pow_dbl(static_cast(base),static_cast(exp_term)) ); +} + +template, + typename std::enable_if::value>::type* = nullptr> +constexpr +TC +pow_check(const T1 base, const T2 exp_term) +noexcept +{ + return pow_integral(base,exp_term); +} + +} + /** * Compile-time power function * @@ -41,16 +70,13 @@ pow_dbl(const T base, const T exp_term) * @return Computes \c base raised to the power \c exp_term. In the case where \c exp_term is integral-valued, recursion by squaring is used, otherwise \f$ \text{base}^{\text{exp\_term}} = e^{\text{exp\_term} \log(\text{base})} \f$ */ -template +template constexpr -Ta -pow(const Ta base, const Tb exp_term) +common_type_t +pow(const T1 base, const T2 exp_term) +noexcept { - return( std::is_integral::value ? \ - // if - pow_integral(base,exp_term) : - // else - pow_dbl(base,Ta(exp_term)) ); + return internal::pow_check(base,exp_term); } #endif diff --git a/include/gcem_incl/pow_integral.hpp b/include/gcem_incl/pow_integral.hpp index 543558e..10a8c90 100644 --- a/include/gcem_incl/pow_integral.hpp +++ b/include/gcem_incl/pow_integral.hpp @@ -28,87 +28,101 @@ namespace internal { -template -constexpr -Ta pow_integral_compute(const Ta base, const Tb exp_term); +template +constexpr T1 pow_integral_compute(const T1 base, const T2 exp_term) noexcept; -// integral-valued powers -// using method described in https://en.wikipedia.org/wiki/Exponentiation_by_squaring +// integral-valued powers using method described in +// https://en.wikipedia.org/wiki/Exponentiation_by_squaring -template +template constexpr -Ta -pow_integral_compute_recur(const Ta base, const Ta val, const Tb exp_term) +T1 +pow_integral_compute_recur(const T1 base, const T1 val, const T2 exp_term) +noexcept { - return( exp_term > Tb(1) ? \ + return( exp_term > T2(1) ? \ (is_odd(exp_term) ? \ pow_integral_compute_recur(base*base,val*base,exp_term/2) : pow_integral_compute_recur(base*base,val,exp_term/2)) : - (exp_term == Tb(1) ? val*base : val) ); + (exp_term == T2(1) ? val*base : val) ); } -template::value>::type* = nullptr> +template::value>::type* = nullptr> constexpr -Ta -pow_integral_sgn_check(const Ta base, const Tb exp_term) +T1 +pow_integral_sgn_check(const T1 base, const T2 exp_term) +noexcept { - return( exp_term < Tb(0) ? \ + return( exp_term < T2(0) ? \ // - Ta(1) / pow_integral_compute(base, - exp_term) : + T1(1) / pow_integral_compute(base, - exp_term) : // - pow_integral_compute_recur(base,Ta(1),exp_term) ); + pow_integral_compute_recur(base,T1(1),exp_term) ); } -template::value>::type* = nullptr> +template::value>::type* = nullptr> constexpr -Ta -pow_integral_sgn_check(const Ta base, const Tb exp_term) +T1 +pow_integral_sgn_check(const T1 base, const T2 exp_term) +noexcept { - return( pow_integral_compute_recur(base,Ta(1),exp_term) ); + return( pow_integral_compute_recur(base,T1(1),exp_term) ); } -template +template constexpr -Ta -pow_integral_compute(const Ta base, const Tb exp_term) +T1 +pow_integral_compute(const T1 base, const T2 exp_term) +noexcept { - return( exp_term == Tb(3) ? \ + return( exp_term == T2(3) ? \ base*base*base : - exp_term == Tb(2) ? \ + exp_term == T2(2) ? \ base*base : - exp_term == Tb(1) ? \ + exp_term == T2(1) ? \ base : - exp_term == Tb(0) ? \ - Ta(1) : + exp_term == T2(0) ? \ + T1(1) : // check for overflow - exp_term == GCLIM::min() ? \ - Ta(0) : - exp_term == GCLIM::max() ? \ - GCLIM::infinity() : + exp_term == GCLIM::min() ? \ + T1(0) : + exp_term == GCLIM::max() ? \ + GCLIM::infinity() : // else pow_integral_sgn_check(base,exp_term) ); } +template::value>::type* = nullptr> +constexpr +T1 +pow_integral_type_check(const T1 base, const T2 exp_term) +noexcept +{ + return pow_integral_compute(base,exp_term); } -// -// main function - -template::value>::type* = nullptr> +template::value>::type* = nullptr> constexpr -return_t -pow_integral(const Ta base, const Tb exp_term) +T1 +pow_integral_type_check(const T1 base, const T2 exp_term) +noexcept { - return internal::pow_integral_compute(return_t(base),exp_term); + // return GCLIM>::quiet_NaN(); + return pow_integral_compute(base,static_cast(exp_term)); } -template::value>::type* = nullptr> +} + +// +// main function + +template constexpr -return_t -pow_integral(const Ta base, const Tb exp_term) +T1 +pow_integral(const T1 base, const T2 exp_term) +noexcept { - // return GCLIM>::quiet_NaN(); - return internal::pow_integral_compute(return_t(base),static_cast(exp_term)); + return internal::pow_integral_type_check(base,exp_term); } #endif diff --git a/include/gcem_incl/sgn.hpp b/include/gcem_incl/sgn.hpp index 614ebad..7c42598 100644 --- a/include/gcem_incl/sgn.hpp +++ b/include/gcem_incl/sgn.hpp @@ -32,6 +32,7 @@ template constexpr int sgn(const T x) +noexcept { return( // positive x > T(0) ? 1 : diff --git a/include/gcem_incl/sin.hpp b/include/gcem_incl/sin.hpp index 4cd5fa2..0335f1f 100644 --- a/include/gcem_incl/sin.hpp +++ b/include/gcem_incl/sin.hpp @@ -34,6 +34,7 @@ template constexpr T sin_compute(const T x) +noexcept { return T(2)*x/(T(1) + x*x); } @@ -42,6 +43,7 @@ template constexpr T sin_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -72,8 +74,9 @@ template constexpr return_t sin(const T x) +noexcept { - return internal::sin_check>(x); + return internal::sin_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/sinh.hpp b/include/gcem_incl/sinh.hpp index 8f7bc00..c6a34f2 100644 --- a/include/gcem_incl/sinh.hpp +++ b/include/gcem_incl/sinh.hpp @@ -32,6 +32,7 @@ template constexpr T sinh_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -53,8 +54,9 @@ template constexpr return_t sinh(const T x) +noexcept { - return internal::sinh_check>(x); + return internal::sinh_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/sqrt.hpp b/include/gcem_incl/sqrt.hpp index 04ba567..5c359d1 100644 --- a/include/gcem_incl/sqrt.hpp +++ b/include/gcem_incl/sqrt.hpp @@ -31,7 +31,8 @@ namespace internal template constexpr T -sqrt_recur(const T x, const T xn, const std::size_t count) +sqrt_recur(const T x, const T xn, const int count) +noexcept { return( abs(xn - x/xn) / (T(1) + xn) < GCLIM::epsilon() ? \ // if @@ -46,6 +47,7 @@ template constexpr T sqrt_check(const T x, const T m_val) +noexcept { return( // negative values x < T(0) ? \ @@ -75,8 +77,9 @@ template constexpr return_t sqrt(const T x) +noexcept { - return internal::sqrt_check>(x,T(1)); + return internal::sqrt_check( static_cast>(x), return_t(1) ); } #endif diff --git a/include/gcem_incl/tan.hpp b/include/gcem_incl/tan.hpp index 64ce2ef..e7476ea 100644 --- a/include/gcem_incl/tan.hpp +++ b/include/gcem_incl/tan.hpp @@ -32,14 +32,16 @@ template constexpr T tan_series_exp_long(const T z) +noexcept { // this is based on a fourth-order expansion of tan(z) using Bernoulli numbers - return( - 1/z + z/3 + pow_integral(z,3)/45 + 2*pow_integral(z,5)/945 + pow_integral(z,7)/4725 ); + return( - 1/z + (z/3 + (pow_integral(z,3)/45 + (2*pow_integral(z,5)/945 + pow_integral(z,7)/4725))) ); } template constexpr T tan_series_exp(const T x) +noexcept { return( GCLIM::epsilon() > abs(x - T(GCEM_HALF_PI)) ? \ // the value tan(pi/2) is somewhat of a convention; @@ -55,6 +57,7 @@ template constexpr T tan_cf_recur(const T xx, const int depth, const int max_depth) +noexcept { return( depth < max_depth ? \ // if @@ -67,6 +70,7 @@ template constexpr T tan_cf_main(const T x) +noexcept { return( (x > T(1.55) && x < T(1.60)) ? \ tan_series_exp(x) : // deals with a singularity at tan(pi/2) @@ -83,6 +87,7 @@ template constexpr T tan_begin(const T x, const int count = 0) +noexcept { // tan(x) = tan(x + pi) return( x > T(GCEM_PI) ? \ // if @@ -96,6 +101,7 @@ template constexpr T tan_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -112,15 +118,20 @@ tan_check(const T x) * Compile-time tangent function * * @param x a real-valued input. - * @return the tangent function using \f[ \tan(x) = \dfrac{x}{1 - \dfrac{x^2}{3 - \dfrac{x^2}{5 - \ddots}}} \f] + * @return the tangent function using + * \f[ \tan(x) = \dfrac{x}{1 - \dfrac{x^2}{3 - \dfrac{x^2}{5 - \ddots}}} \f] + * To deal with a singularity at \f$ \pi / 2 \f$, the following expansion is employed: + * \f[ \tan(x) = - \frac{1}{x-\pi/2} - \sum_{k=1}^\infty \frac{(-1)^k 2^{2k} B_{2k}}{(2k)!} (x - \pi/2)^{2k - 1} \f] + * where \f$ B_n \f$ is the n-th Bernoulli number. */ template constexpr return_t tan(const T x) +noexcept { - return internal::tan_check>(x); + return internal::tan_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/tanh.hpp b/include/gcem_incl/tanh.hpp index 3f13b6c..a6027cf 100644 --- a/include/gcem_incl/tanh.hpp +++ b/include/gcem_incl/tanh.hpp @@ -32,6 +32,7 @@ template constexpr T tanh_cf(const T xx, const int depth) +noexcept { return( depth < GCEM_TANH_MAX_ITER ? \ // if @@ -44,6 +45,7 @@ template constexpr T tanh_begin(const T x) +noexcept { return( x/tanh_cf(x*x,1) ); } @@ -52,6 +54,7 @@ template constexpr T tanh_check(const T x) +noexcept { return( // indistinguishable from zero GCLIM::epsilon() > abs(x) ? \ @@ -75,8 +78,9 @@ template constexpr return_t tanh(const T x) +noexcept { - return internal::tanh_check>(x); + return internal::tanh_check( static_cast>(x) ); } #endif diff --git a/include/gcem_incl/tgamma.hpp b/include/gcem_incl/tgamma.hpp index 9a85628..3e68be2 100644 --- a/include/gcem_incl/tgamma.hpp +++ b/include/gcem_incl/tgamma.hpp @@ -32,6 +32,7 @@ template constexpr T tgamma_check(const T x) +noexcept { return( // indistinguishable from one or zero GCLIM::epsilon() > abs(x - T(1)) ? \ @@ -60,8 +61,9 @@ template constexpr return_t tgamma(const T x) +noexcept { - return internal::tgamma_check>(x); + return internal::tgamma_check( static_cast>(x) ); } #endif diff --git a/tests/gcd.cpp b/tests/gcd.cpp new file mode 100644 index 0000000..7835d99 --- /dev/null +++ b/tests/gcd.cpp @@ -0,0 +1,52 @@ +/*################################################################################ + ## + ## Copyright (C) 2016-2018 Keith O'Hara + ## + ## This file is part of the GCE-Math C++ library. + ## + ## Licensed under the Apache License, Version 2.0 (the "License"); + ## you may not use this file except in compliance with the License. + ## You may obtain a copy of the License at + ## + ## http://www.apache.org/licenses/LICENSE-2.0 + ## + ## Unless required by applicable law or agreed to in writing, software + ## distributed under the License is distributed on an "AS IS" BASIS, + ## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + ## See the License for the specific language governing permissions and + ## limitations under the License. + ## + ################################################################################*/ + +#define TEST_VAL_TYPES int + +#include "gcem_tests.hpp" + +int main() +{ + std::cout << "\n*** begin gcd test ***\n" << std::endl; + + // + + std::function test_fn = gcem::gcd; + std::string test_fn_name = "gcem::gcd"; + + // + + static constexpr int test_vals_1[] = { 12, -4, 42 }; // shape + static constexpr int test_vals_2[] = { 18, 14, 56 }; // x + + static constexpr int expected_vals[] = { 6, 2, 14 }; + + // + + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,0,test_fn,expected_vals[0],true,"",3,18,false,false); + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,1,test_fn,expected_vals[1],true,"",3,18,false,false); + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,2,test_fn,expected_vals[2],false,"",3,18,false,false); + + // + + std::cout << "\n*** end gcd test ***\n" << std::endl; + + return 0; +} diff --git a/tests/gcem_tests.hpp b/tests/gcem_tests.hpp index 411f4ba..50b7cc9 100644 --- a/tests/gcem_tests.hpp +++ b/tests/gcem_tests.hpp @@ -29,10 +29,28 @@ // // tolerance +#ifndef TEST_VAL_TYPES +#define TEST_VAL_TYPES long double +#define TEST_VAL_TYPES_V 1 +#endif + #ifndef TEST_ERR_TOL #define TEST_ERR_TOL 1e-14 #endif +#ifdef TEST_VAL_TYPES_V +#define PRINT_ERR(err_val) \ +{ \ + printf(" error value = %LE.", err_val); \ +} +#else +#define PRINT_ERR(err_val) \ +{ \ + printf(" error value = %d.", err_val); \ +} +#endif + + // // one input @@ -55,16 +73,16 @@ print_test_1(const std::string fn_name, const T val_inp_1, std::function +inline +T +print_test_3(const std::string fn_name, const T val_inp_1, const T val_inp_2, const T val_inp_3, + std::function fn_eval, + const bool new_line=false, const std::string extra_space="", + const int precision_1=2, const int precision_2=18) +{ + T f_val = fn_eval(val_inp_1,val_inp_2,val_inp_3); + std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(precision_1) << fn_name \ + << "(" << val_inp_1 << "," << val_inp_2 << "," << val_inp_3 << ")" << extra_space << " = " \ + << std::setprecision(precision_2) << f_val << std::endl; + if (new_line) std::cout << std::endl; + return f_val; +} + +#define PRINT_TEST_3_EXPECT(fn_name, vals_inp_1, vals_inp_2, vals_inp_3, \ + j, fn_eval, expected_val, \ + new_line, extra_space, \ + precision_1, precision_2, is_inf, is_nan) \ { \ - constexpr long double c_val_1 = vals_inp_1[j]; \ - constexpr long double c_val_2 = vals_inp_2[j]; \ - long double x_val_1 = vals_inp_1[j]; \ - long double x_val_2 = vals_inp_2[j]; \ + constexpr TEST_VAL_TYPES val_1 = vals_inp_1[j]; \ + constexpr TEST_VAL_TYPES val_2 = vals_inp_2[j]; \ + constexpr TEST_VAL_TYPES val_3 = vals_inp_3[j]; \ \ - long double f_val_1 = print_test_2(fn_name_1,c_val_1,c_val_2,fn_eval_1,false,"",precision_1,precision_2); \ - long double f_val_2 = print_test_2(fn_name_2,x_val_1,x_val_2,fn_eval_2,false,extra_space,precision_1,precision_2); \ + TEST_VAL_TYPES f_val = print_test_3(fn_name,val_1,val_2,val_3,fn_eval,false,extra_space,precision_1,precision_2); \ \ bool test_pass = false; \ \ printf(" -"); \ if (!is_nan && !is_inf) { \ - long double err_val = std::abs(f_val_1 - f_val_2); \ - printf(" error value = %LE.", err_val); \ + TEST_VAL_TYPES err_val = std::abs(f_val - expected_val); \ + PRINT_ERR(err_val); \ test_pass = (err_val < TEST_ERR_TOL) ? true : false; \ } else if (is_nan) { \ - if (std::isnan(f_val_1) && std::isnan(f_val_2)) { \ + if (std::isnan(f_val)) { \ test_pass = true; \ } \ } else if (is_inf) { \ - if (std::isinf(f_val_1) && std::isinf(f_val_2)) { \ + if (std::isinf(f_val)) { \ test_pass = true; \ } \ } else { \ @@ -221,61 +297,3 @@ print_test_2(const std::string fn_name, const T val_inp_1, const T val_inp_2, st if (new_line) std::cout << std::endl; \ } \ -// -// three inputs - -template -inline -T -print_test_3(const std::string fn_name, const T val_inp_1, const T val_inp_2, const T val_inp_3, - std::function fn_eval, - const bool new_line=false, const std::string extra_space="", - const int precision_1=2, const int precision_2=18) -{ - T f_val = fn_eval(val_inp_1,val_inp_2,val_inp_3); - std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(precision_1) << fn_name \ - << "(" << val_inp_1 << "," << val_inp_2 << "," << val_inp_3 << ")" << extra_space << " = " \ - << std::setprecision(precision_2) << f_val << std::endl; - if (new_line) std::cout << std::endl; - return f_val; -} - -#define PRINT_TEST_3_EXPECT(fn_name, vals_inp_1, vals_inp_2, vals_inp_3, \ - j, fn_eval, expected_val, \ - new_line, extra_space, \ - precision_1, precision_2, is_inf, is_nan) \ -{ \ - constexpr long double val_1 = vals_inp_1[j]; \ - constexpr long double val_2 = vals_inp_2[j]; \ - constexpr long double val_3 = vals_inp_3[j]; \ - \ - long double f_val = print_test_3(fn_name,val_1,val_2,val_3,fn_eval,false,extra_space,precision_1,precision_2); \ - \ - bool test_pass = false; \ - \ - printf(" -"); \ - if (!is_nan && !is_inf) { \ - long double err_val = std::abs(f_val - expected_val); \ - printf(" error value = %LE.", err_val); \ - test_pass = (err_val < TEST_ERR_TOL) ? true : false; \ - } else if (is_nan) { \ - if (std::isnan(f_val)) { \ - test_pass = true; \ - } \ - } else if (is_inf) { \ - if (std::isinf(f_val)) { \ - test_pass = true; \ - } \ - } else { \ - printf("unknown error!"); \ - } \ - \ - if (test_pass) { \ - std::cout << "\033[32m OK.\033[0m" << std::endl; \ - } else { \ - std::cout << "\033[31m FAIL.\033[0m" << std::endl; \ - } \ - \ - if (new_line) std::cout << std::endl; \ -} \ - diff --git a/tests/incomplete_beta.cpp b/tests/incomplete_beta.cpp index 20e8b8a..2876676 100644 --- a/tests/incomplete_beta.cpp +++ b/tests/incomplete_beta.cpp @@ -26,7 +26,7 @@ int main() // - std::function test_fn = gcem::incomplete_beta; + std::function test_fn = gcem::incomplete_beta; std::string test_fn_name = "gcem::incomplete_beta"; // diff --git a/tests/incomplete_beta_inv.cpp b/tests/incomplete_beta_inv.cpp index cc563fc..749726a 100644 --- a/tests/incomplete_beta_inv.cpp +++ b/tests/incomplete_beta_inv.cpp @@ -26,7 +26,7 @@ int main() // - std::function test_fn = gcem::incomplete_beta_inv; + std::function test_fn = gcem::incomplete_beta_inv; std::string test_fn_name = "gcem::incomplete_beta_inv"; // diff --git a/tests/incomplete_gamma.cpp b/tests/incomplete_gamma.cpp index 0560105..e3efb93 100644 --- a/tests/incomplete_gamma.cpp +++ b/tests/incomplete_gamma.cpp @@ -38,13 +38,13 @@ int main() 1.0L, 1.0L, 3.0L, 5.0L, 9.0L, 9.0L, 0.0L, 0.0L, 11.0L, 18.0L, 18.0L, 19.0L,\ 39.0L, 55.0L, 99.0L, 101.0L, 297.0L, 301.0L, 497.0L, 501.0L, 799.0L, 800.0L, 999.0L, 1001.0L }; // x - static constexpr long double expected_vals[] = { 0.26424111765711527644L, \ - 0.42759329552912134220L, \ - 0.80085172652854419439L, \ - 0.95957231800548714595L, \ - 0.99876590195913317327L, \ - 1.0L, \ - 0.0L,\ + static constexpr long double expected_vals[] = { 0.26424111765711527644L, + 0.42759329552912134220L, + 0.80085172652854419439L, + 0.95957231800548714595L, + 0.99876590195913317327L, + 1.0L, + 0.0L, 0.0L, 0.47974821959920432857L, 0.75410627202774938027L, diff --git a/tests/lcm.cpp b/tests/lcm.cpp new file mode 100644 index 0000000..03ffea4 --- /dev/null +++ b/tests/lcm.cpp @@ -0,0 +1,52 @@ +/*################################################################################ + ## + ## Copyright (C) 2016-2018 Keith O'Hara + ## + ## This file is part of the GCE-Math C++ library. + ## + ## Licensed under the Apache License, Version 2.0 (the "License"); + ## you may not use this file except in compliance with the License. + ## You may obtain a copy of the License at + ## + ## http://www.apache.org/licenses/LICENSE-2.0 + ## + ## Unless required by applicable law or agreed to in writing, software + ## distributed under the License is distributed on an "AS IS" BASIS, + ## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + ## See the License for the specific language governing permissions and + ## limitations under the License. + ## + ################################################################################*/ + +#define TEST_VAL_TYPES int + +#include "gcem_tests.hpp" + +int main() +{ + std::cout << "\n*** begin lcm test ***\n" << std::endl; + + // + + std::function test_fn = gcem::lcm; + std::string test_fn_name = "gcem::lcm"; + + // + + static constexpr int test_vals_1[] = { 3, 12, 17 }; // shape + static constexpr int test_vals_2[] = { 4, 15, 81 }; // x + + static constexpr int expected_vals[] = { 12, 60, 1377 }; + + // + + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,0,test_fn,expected_vals[0],true,"",3,18,false,false); + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,1,test_fn,expected_vals[1],true,"",3,18,false,false); + PRINT_TEST_2_EXPECT(test_fn_name,test_vals_1,test_vals_2,2,test_fn,expected_vals[2],false,"",3,18,false,false); + + // + + std::cout << "\n*** end lcm test ***\n" << std::endl; + + return 0; +} From 83dab5f7464103265bd606230fac7f9dd634d332 Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Thu, 1 Nov 2018 20:35:22 -0400 Subject: [PATCH 2/5] update tests --- .appveyor.yml | 2 ++ tests/CMakeLists.txt | 2 ++ 2 files changed, 4 insertions(+) diff --git a/.appveyor.yml b/.appveyor.yml index bf9ceb7..e536738 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -40,10 +40,12 @@ build_script: - .\erf.test - .\erf_inv.test - .\factorial.test + - .\gcd.test - .\incomplete_beta.test - .\incomplete_beta_inv.test - .\incomplete_gamma.test - .\incomplete_gamma_inv.test + - .\lcm.test - .\lgamma.test - .\log.test - .\other.test diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index cab405b..205396b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -68,10 +68,12 @@ set(GCEM_TESTS erf_inv.cpp exp.cpp factorial.cpp + gcd.cpp incomplete_beta.cpp incomplete_beta_inv.cpp incomplete_gamma.cpp incomplete_gamma_inv.cpp + lcm.cpp lgamma.cpp log.cpp other.cpp From 1921072f7ea7aca8b20279d7298ec29f2530b332 Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Thu, 1 Nov 2018 20:49:10 -0400 Subject: [PATCH 3/5] update tests --- tests/gcem_tests.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/gcem_tests.hpp b/tests/gcem_tests.hpp index 50b7cc9..8349d07 100644 --- a/tests/gcem_tests.hpp +++ b/tests/gcem_tests.hpp @@ -35,7 +35,11 @@ #endif #ifndef TEST_ERR_TOL -#define TEST_ERR_TOL 1e-14 +#ifdef _WIN32 + #define TEST_ERR_TOL 1e-10 +#else + #define TEST_ERR_TOL 1e-14 +#endif #endif #ifdef TEST_VAL_TYPES_V From 65f3ad82fb199b51d9efb595a70ab4a5a0187a63 Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Thu, 1 Nov 2018 21:52:11 -0400 Subject: [PATCH 4/5] satisfy MSVC --- include/gcem_incl/binomial_coef.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/gcem_incl/binomial_coef.hpp b/include/gcem_incl/binomial_coef.hpp index 34f4ade..1b835ad 100644 --- a/include/gcem_incl/binomial_coef.hpp +++ b/include/gcem_incl/binomial_coef.hpp @@ -52,7 +52,7 @@ T binomial_coef_sgn_check(const T n, const T k) noexcept { - return binomial_coef_recur(static_cast(n),static_cast(k)); + return static_cast(binomial_coef_recur(static_cast(n),static_cast(k))); } template> From 32ebef21239fff0575f2676b9b99cc2e337bd80f Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Thu, 1 Nov 2018 21:52:22 -0400 Subject: [PATCH 5/5] update tests --- tests/gcem_tests.hpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/tests/gcem_tests.hpp b/tests/gcem_tests.hpp index 8349d07..181a50f 100644 --- a/tests/gcem_tests.hpp +++ b/tests/gcem_tests.hpp @@ -54,6 +54,13 @@ } #endif +#ifdef TEST_VAL_TYPES_V +#define TEST_IS_INF(val) std::isinf(val) +#define TEST_IS_NAN(val) std::isnan(val) +#else +#define TEST_IS_INF(val) false +#define TEST_IS_NAN(val) false +#endif // // one input @@ -89,11 +96,11 @@ print_test_1(const std::string fn_name, const T val_inp_1, std::function