From 4b88319545e0e48a635d35b7eedf24612cb3c4d8 Mon Sep 17 00:00:00 2001 From: camilo Date: Thu, 25 Jan 2024 12:41:05 -0500 Subject: [PATCH] added ffmath::factorial --- check/qlibs_cpp_test.cpp | 4 +-- src/ffmath.cpp | 67 +++++++++++++++++++++++++++++----------- src/include/ffmath.hpp | 12 +++++++ 3 files changed, 63 insertions(+), 20 deletions(-) diff --git a/check/qlibs_cpp_test.cpp b/check/qlibs_cpp_test.cpp index 33787cd..352c8ba 100644 --- a/check/qlibs_cpp_test.cpp +++ b/check/qlibs_cpp_test.cpp @@ -304,8 +304,8 @@ void test_ffmath(void) cout << ffmath::wrapToPi( -4.5 ) << endl; cout << ffmath::getNan() << endl; cout << ffmath::getInf() << endl; - cout << ffmath::tgamma( 12.5 ) << endl; - cout << ffmath::lgamma( 12.5 ) << endl; + cout << ffmath::tgamma( 2.5 ) << endl; + cout << ffmath::lgamma( 2.5 ) << endl; } void test_mat( void ) diff --git a/src/ffmath.cpp b/src/ffmath.cpp index 20e492e..34f3b0a 100644 --- a/src/ffmath.cpp +++ b/src/ffmath.cpp @@ -5,7 +5,7 @@ using namespace qlibs; template static inline void cast_reinterpret( T1 &f, const T2 &u ) { - static_assert( sizeof(T1)==sizeof(T2), "Types must match sizes" ); + static_assert( sizeof(T1) == sizeof(T2), "Types must match sizes" ); (void)memcpy( &f, &u, sizeof(T2) ); } @@ -36,7 +36,7 @@ bool ffmath::isEqual( const float a, const float b, const float tol ) noexcept { - return ( ffmath::absf( a - b ) <= ffmath::absf( tol ) ); + return ( absolute( a - b ) <= absolute( tol ) ); } /*============================================================================*/ float ffmath::getInf( void ) @@ -98,7 +98,7 @@ float ffmath::sign( float x ) /*============================================================================*/ float ffmath::absf( float x ) { - return ( x < 0.0F ) ? -x : x; + return absolute( x ); } /*============================================================================*/ float ffmath::recip( float x ) @@ -130,7 +130,7 @@ float ffmath::sqrt( float x ) cast_reinterpret( y, x ); y = ( ( y - 0x00800000U ) >> 1U ) + 0x20000000U; cast_reinterpret( z, y ); - retVal = ( ( x/z ) + z ) * 0.5F; + retVal = 0.5F*( ( x/z ) + z ); } return retVal; @@ -565,10 +565,10 @@ float ffmath::nextAfter( float x, float y ) uxi = ( 0U == ay ) ? uyi : ( ( uyi & 0x80000000U ) | 1U ); } else if ( ( ax > ay ) || ( 0U != ( ( uxi^uyi ) & 0x80000000U ) ) ) { - uxi--; + --uxi; } else { - uxi++; + ++uxi; } cast_reinterpret( retVal, uxi ); } @@ -588,12 +588,7 @@ float ffmath::tgamma( float x ) result = getInf(); /* a huge value */ } else if ( classification::FFP_INFINITE == fClass ) { - if ( x > 0.0F ) { - result = getInf(); /* a huge value */ - } - else { - result = getNan(); - } + result = ( x > 0.0F ) ? getInf() : getNan(); } else { bool parity = false; @@ -763,7 +758,7 @@ static float lgamma_positive( float x ) float result; if ( x > 171.624F ) { - result = ffmath::getInf(); /* a huge value */ + result = ffmath::getInf(); } else { float y, corrector, num, den; @@ -853,10 +848,7 @@ float ffmath::lgamma( float x ) if ( classification::FFP_NAN == fClass ) { result = getNan(); } - else if ( classification::FFP_ZERO == fClass ) { - result = getInf(); - } - else if ( classification::FFP_INFINITE == fClass ) { + else if ( ( classification::FFP_ZERO == fClass ) || ( classification::FFP_INFINITE == fClass ) ) { result = getInf(); } else { @@ -877,7 +869,7 @@ float ffmath::lgamma( float x ) float a; a = sin( FFP_PI*isItAnInt ); - result = ffmath::log( FFP_PI/ffmath::absf( a*x ) ) - lgamma_positive( -x ); + result = ffmath::log( FFP_PI/absolute( a*x ) ) - lgamma_positive( -x ); } } } @@ -888,4 +880,43 @@ float ffmath::lgamma( float x ) return result; } +/*============================================================================*/ +float ffmath::factorial( float x ) +{ + constexpr float ft[ 35 ] = { 1.0F, 1.0F, 2.0F, 6.0F, 24.0F, 120.0F, 720.0F, + 5040.0F, 40320.0F, 362880.0F, 3628800.0F, + 39916800.0F, 479001600.0F, 6227020800.0F, + 87178291200.0F, 1307674368000.0F, + 20922789888000.0F, 355687428096000.0F, + 6402373705728001.0F, 121645100408832000.0F, + 2432902008176640000.0F, 51090942171709440000.0F, + 1124000727777607680000.0F, + 25852016738884978212864.0F, + 620448401733239544217600.0F, + 15511210043330988202786816.0F, + 403291461126605719042260992.0F, + 10888869450418351940239884288.0F, + 304888344611713836734530715648.0F, + 8841761993739700772720181510144.0F, + 265252859812191104246398737973248.0F, + 8222838654177921277277005322125312.0F, + 263130836933693591553328612565319680.0F, + 8683317618811885938715673895318323200.0F, + 295232799039604119555149671006000381952.0F, + }; + float y; + + if ( x > 34.0F ) { + y = getInf(); + } + else if ( x >= 0.0F ) { + const auto i = static_cast( x ); + y = ft[ i ]; + } + else { + y = getNan(); + } + + return y; +} /*============================================================================*/ \ No newline at end of file diff --git a/src/include/ffmath.hpp b/src/include/ffmath.hpp index 387d781..c9d27c0 100644 --- a/src/include/ffmath.hpp +++ b/src/include/ffmath.hpp @@ -582,6 +582,17 @@ namespace qlibs { */ float lgamma( float x ); + /** + * @brief Return the factorial of the integer part of @a x. + * @note The argument @a x needs to be positive + * @note This function overflows when @a x > 34 + * @param[in] x The floating point value + * @return Upon successful completion, this function shall return the + * factorial of the integer part of x. If x is non-positive, factorial() shall + * return NaN. If the correct value would cause overflow, lgamma() shall + * return +Inf. + */ + float factorial( float x ); /*cstat -MISRAC++2008-0-1-4_b*/ @@ -615,6 +626,7 @@ namespace qlibs { constexpr float FFP_LN_SQRT_2PI = ( 0.9189385332046727417803297F ); + /*cstat +MISRAC++2008-0-1-4_b*/ /** @}*/