diff --git a/include/qffmath.h b/include/qffmath.h
index 1a006ce..d5e1a70 100644
--- a/include/qffmath.h
+++ b/include/qffmath.h
@@ -91,6 +91,13 @@ extern "C" {
/** @brief The natural logarithm of the square root of 2π given as a single-precision floating-point number */
#define QFFM_LN_SQRT_2PI ( 0.918938533204672669540968854562379419F )
+ /** @brief Constant Euler-Mascheroni given as a single-precision floating-point number*/
+ #define QFFM_GAMMA_E ( 0.577215664901532860606512090082402431F )
+ /** @brief Twice circumference of a circle with diameter 1, ( 2π ) given as a single-precision floating-point number */
+ #define QFFM_2PI ( 6.283185307179586231995926937088370323F )
+ /** @brief The natural logarithm of π ( ln(π) ) given as a single-precision floating-point number */
+ #define QFFM_LN_PI ( 1.144729885849400163877476188645232468F )
+
/** @brief The maximum value of a non-infinite single-precision floating-point number */
#define QFFM_MAXFLOAT ( 3.40282347e+38F )
/** @brief Positive infinity given as a single-precision floating-point number */
@@ -464,6 +471,42 @@ extern "C" {
*/
float qFFMath_ATanh( float x );
+ /**
+ * @brief Wraps angle @a x, in radians, to the interval [−pi, pi] such that
+ * @c pi maps to @c pi and @c -pi maps to @c -pi. In general, odd, positive multiples
+ * of @c pi map to @c pi and odd, negative multiples of @c pi map to @c -pi.
+ * @param x The angle in radians
+ * @return The angle @a x wrapped to the [-pi, pi] interval
+ */
+ float qFFMath_WrapToPi( float x );
+
+ /**
+ * @brief Wraps angle @a x, in radians, to the interval [0, 2*pi] such that
+ * @c 0 maps to @c 0 and @c 2*pi and @c 2*pi maps to @c 2*pi. In general, positive multiples
+ * of @c 2*pi map to @c 2*pi and negative multiples of @c 2*pi map to @c 0.
+ * @param x The angle in radians
+ * @return The angle @a x wrapped to the [0, 2*pi] interval
+ */
+ float qFFMath_WrapTo2Pi( float x );
+
+ /**
+ * @brief Wraps angle @a x, in degrees, to the interval [–180, 180] such
+ * that @c 180 maps to @c 180 and @c -180 maps to @c -180. In general, odd, positive
+ * multiples of @c 180 map to @c 180 and odd, negative multiples of @c 180 map to @c -180.
+ * @param x The angle in degrees
+ * @return The angle @a x wrapped to the [-180, 180] interval
+ */
+ float qFFMath_WrapTo180( float x );
+
+ /**
+ * @brief Wraps angle @a x, in degrees, to the interval [0, 360] such
+ * that @c 0 maps to @c 0 and @c 360 maps to @c 360. In general, positive multiples
+ * of @c 360 map to @c 360 and negative multiples of @c 360 map to zero.
+ * @param x The angle in degrees
+ * @return The angle @a x wrapped to the [0, 360] interval
+ */
+ float qFFMath_WrapTo360( float x );
+
/**
* @brief Computes the error function of @a x.
* @param[in] x The floating point value
@@ -708,6 +751,265 @@ extern "C" {
*/
float qFFMath_Factorial( float x );
+ /**
+ * @brief Computes the associated Laguerre polynomials of the degree @a n,
+ * order @a m, and argument @a x.
+ * @param[in] n The degree of the polynomial, an unsigned integer value
+ * @param[in] m The order of the polynomial, an unsigned integer value
+ * @param[in] x The argument, a floating-point or integer value
+ * @return Upon successful completion, the value of the associated Laguerre
+ * polynomial of @a x shall be returned. If the argument is @c nan, a @c nan is
+ * returned. If @a x is negative, @c nan is returned. If @a n or @a m is greater or
+ * equal to @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Assoc_laguerre( size_t n,
+ size_t m,
+ float x );
+
+ /**
+ * @brief Computes the associated Legendre polynomials of the degree @a n,
+ * order @a m, and argument @a x.
+ * @param[in] n The degree of the polynomial, an unsigned integer value
+ * @param[in] m The order of the polynomial, an unsigned integer value
+ * @param[in] x The argument, a floating-point or integer value
+ * @return Upon successful completion, the value of the associated Legendre
+ * polynomial of x shall be returned. If the argument is @c nan, a @c nan is
+ * returned. If |x| > 1, @c nan is returned due the domain error.
+ * If @a n is greater or equal to @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Assoc_legendre( size_t n,
+ size_t m,
+ float x );
+
+ /**
+ * @brief Computes the Beta function of @a x and @a y.
+ * @param[in] x The argument, a floating-point or integer value
+ * @param[in] y The argument, a floating-point or integer value
+ * @return Upon successful completion, the value of the Beta function of
+ * @a x and @a y. If the argument is @c nan, @c nan is returned. The function
+ * is only required to be defined where both @a x and @a y are greater
+ * than zero, and is allowed to return @c nan otherwise.
+ */
+ float qFFMath_Beta( float x,
+ float y );
+
+ /**
+ * @brief Computes the complete elliptic integral of the first kind of @a k
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @return Upon successful completion, the value of the complete elliptic
+ * integral of the first kind of @a k. If the argument is @c nan, @c nan is
+ * returned. If |k| > 1, NaN is returned due the domain error
+ */
+ float qFFMath_Comp_ellint_1( float k );
+
+ /**
+ * @brief Computes the complete elliptic integral of the second kind of @a k
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @return Upon successful completion, the value of the complete elliptic
+ * integral of the second kind of @a k. If the argument is @c nan, @c nan is
+ * returned. If |k| > 1, @c nan is returned due the domain error
+ */
+ float qFFMath_Comp_ellint_2( float k );
+
+ /**
+ * @brief Computes the complete elliptic integral of the third kind of
+ * @a k and @a nu.
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @param[in] nu Elliptic characteristic as a floating-point value
+ * @return Upon successful completion, the value of the complete elliptic
+ * integral of the third kind of @a k and @a nu. If the argument is @c nan,
+ * @c nan is returned. If |k| > 1, @c nan is returned due the domain error
+ */
+ float qFFMath_Comp_ellint_3( float k,
+ float nu );
+
+ /**
+ * @brief Computes the incomplete elliptic integral of the first kind of
+ * @a k and @a phi
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @param[in] phi Jacobi amplitude as a floating-point value given in radians
+ * @return Upon successful completion, the value of the incomplete elliptic
+ * integral of the first kind of @a k and @a phi. If the argument is @c nan,
+ * @c nan is returned. If |k| > 1, @c nan is returned due the domain error
+ */
+ float qFFMath_Ellint_1( float k,
+ float phi );
+
+ /**
+ * @brief Computes the incomplete elliptic integral of the second kind of
+ * @a k and @a phi
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @param[in] phi Jacobi amplitude as a floating-point value given in radians
+ * @return Upon successful completion, the value of the incomplete elliptic
+ * integral of the second kind of @a k and @a phi. If the argument is @c nan,
+ * @c nan is returned. If |k| > 1, @c nan is returned due the domain error
+ */
+ float qFFMath_Ellint_2( float k,
+ float phi );
+
+ /**
+ * @brief Computes the incomplete elliptic integral of the third kind of
+ * @a k, @a nu and @a phi
+ * @param[in] k Elliptic modulus or eccentricity as a floating-point value
+ * @param[in] nu Elliptic characteristic as a floating-point value
+ * @param[in] phi Jacobi amplitude as a floating-point value given in radians
+ * @return Upon successful completion, the value of the incomplete elliptic
+ * integral of the third kind of @a k, @a nu and @a phi. If the argument is @c nan,
+ * @c nan is returned. If |k| > 1, @c nan is returned due the domain error
+ */
+ float qFFMath_Ellint_3( float k,
+ float nu,
+ float phi );
+
+ /**
+ * @brief Computes the Exponential integral of @a num
+ * @param[in] num A floating-point value
+ * @return Upon successful completion, the value of the exponential integral
+ * of @a num. If the argument is @c nan, @c nan is returned. If the argument
+ * is @c ±0, @c -inf is returned.
+ */
+ float qFFMath_Expint( float num );
+
+ /**
+ * @brief Computes the (physicist's) Hermite polynomials of the degree
+ * @a n and argument @a x
+ * @param[in] n The degree of the polynomial
+ * @param[in] x The argument, a floating-point value
+ * @return Upon successful completion, the value of order-n Hermite polynomial
+ * of @a x. If the argument is @c nan, @c nan is returned. If n>=128,
+ * the behavior is implementation-defined.
+ */
+ float qFFMath_Hermite( size_t n,
+ float x );
+
+ /**
+ * @brief Computes the non-associated Laguerre polynomials of the degree @a n,
+ * and argument @a x.
+ * @param[in] n The degree of the polynomial, an unsigned integer value
+ * @param[in] x The argument, a floating-point or integer value
+ * @return Upon successful completion, the value of the non-associated Laguerre
+ * polynomial of @a x shall be returned. If the argument is @c nan, a @c nan is
+ * returned. If @a x is negative, @c nan is returned. If @a n is greater or
+ * equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Laguerre( size_t n,
+ float x );
+
+ /**
+ * @brief Computes the unassociated Legendre polynomials of the degree @a n,
+ * and argument @a x.
+ * @param[in] n The degree of the polynomial, an unsigned integer value
+ * @param[in] x The argument, a floating-point or integer value
+ * @return Upon successful completion, the value of the unassociated Legendre
+ * polynomial of @a x shall be returned. If the argument is @c nan, a @c nan is
+ * returned. The function is not required to be defined for |x| > 1 .
+ * If @a n is greater or equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Legendre( size_t n,
+ float x );
+
+ /**
+ * @brief Computes the Riemann zeta function of @a s
+ * @param[in] s A floating-point value
+ * @return Upon successful completion, the value of the Riemann zeta function
+ * of @a s. If the argument is @c nan, @c nan is returned.
+ */
+ float qFFMath_Riemann_zeta( float s );
+
+ /**
+ * @brief Computes the spherical Bessel function of the first kind @a n,
+ * and @a x.
+ * @param[in] n The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the spherical Bessel
+ * function of the first kind of @a n and @a x. If the argument is @c nan, a @c nan is
+ * returned. If @a n is greater or equal than @c 128, the behavior is
+ * implementation-defined.
+ */
+ float qFFMath_Sph_bessel( size_t n,
+ float x );
+
+ /**
+ * @brief Computes the spherical Bessel function of the second kind also
+ * known as the spherical Neumann function of @a n and @a x.
+ * @param[in] n The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the spherical Bessel
+ * function of the second kind (spherical Neumann function) of @a n and
+ * @a x. If the argument is @c nan, a @c nan is
+ * returned. If @a n is greater or equal than @c 128, the behavior is
+ * implementation-defined.
+ */
+ float qFFMath_Sph_neumann( size_t n,
+ float x );
+
+ /**
+ * @brief Computes the regular modified cylindrical Bessel function of
+ * @a nu and @a x
+ * @param[in] nu The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the regular modified
+ * cylindrical Bessel function of @a nu and @a x. If the argument is
+ * @c nan, a @c nan is returned. If @a nu is greater or equal than @c 128,
+ * the behavior is implementation-defined.
+ */
+ float qFFMath_Cyl_bessel_i( float nu,
+ float x );
+
+ /**
+ * @brief Computes the cylindrical Bessel function of the first kind of @a nu
+ * and @a x
+ * @param[in] nu The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the irregular modified cylindrical Bessel function
+ * (also known as modified Bessel function of the second kind) of @a nu
+ * and @a x. If the argument is @c nan, a @c nan is returned. If @a nu is
+ * greater or equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Cyl_bessel_j( float nu,
+ float x );
+
+ /**
+ * @brief Computes the irregular modified cylindrical Bessel function
+ * (also known as modified Bessel function of the second kind) of @a nu
+ * and @a x
+ * @param[in] nu The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the irregular modified cylindrical Bessel function
+ * (also known as modified Bessel function of the second kind) of @a nu
+ * and @a x. If the argument is @c nan, a @c nan is returned. If @a nu is
+ * greater or equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Cyl_bessel_k( float nu,
+ float x );
+
+ /**
+ * @brief Computes the cylindrical Neumann function ( also known as Bessel
+ * function of the second kind or Weber function) of @a nu and @a x.
+ * @param[in] nu The order of the function
+ * @param[in] x The argument to the function, a floating-point or integer value
+ * @return Upon successful completion, the value of the cylindrical Neumann
+ * function ( Bessel function of the second kind) of @a nu
+ * and @a x. If the argument is @c nan, a @c nan is returned. If @a nu is
+ * greater or equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Cyl_neumann( float nu,
+ float x );
+
+ /**
+ * @brief 1-3) Computes the spherical associated Legendre function of
+ * degree @a l, order @a m, and polar angle @a theta
+ * @param[in] l The degree
+ * @param[in] m The order
+ * @param[in] theta Polar angle, measured in radians
+ * @return Upon successful completion, the value of the spherical associated
+ * Legendre function (that is, spherical harmonic with sigma = 0) of @a l,
+ * @a m, and @a theta. If the argument @a theta is @c nan, a @c nan is returned.
+ * If @a l is greater or equal than @c 128, the behavior is implementation-defined.
+ */
+ float qFFMath_Sph_legendre( size_t l,
+ size_t m,
+ float theta );
#endif /*#ifdef QLIBS_USE_STD_MATH*/
/** @}*/
diff --git a/qffmath.c b/qffmath.c
index e675f6b..83ef50f 100644
--- a/qffmath.c
+++ b/qffmath.c
@@ -18,6 +18,67 @@
static float qFFMath_CalcCbrt( float x , bool r );
static float lgamma_positive( float x );
+static float poly_laguerre_recursion( size_t n,
+ float alpha,
+ float x );
+static float poly_laguerre_large_n( size_t n,
+ float alpha,
+ float x );
+static float poly_laguerre_hyperg( size_t n,
+ float alpha,
+ float x );
+static float poly_legendre_p( size_t l,
+ float x );
+static float ellint_rf( float x,
+ float y,
+ float z );
+static float ellint_rd( float x,
+ float y,
+ float z );
+static float ellint_rc( float x,
+ float y );
+static float ellint_rj( float x,
+ float y,
+ float z,
+ float p );
+static float expint_E1_series( float x );
+static float expint_E1_asymp( float x );
+static float expint_Ei_asymp( float x );
+static float expint_E1( float x );
+static float expint_En_cont_frac( size_t n,
+ float x );
+static float expint_Ei_series( float x );
+static float expint_Ei( float x );
+static float riemann_zeta_glob( float s );
+static float riemann_zeta_product( float s );
+static void gamma_temme( float mu,
+ float* gam1,
+ float* gam2,
+ float* gam_pl,
+ float* gam_mi);
+static void bessel_jn( float nu,
+ float x,
+ float* j_nu,
+ float* n_nu,
+ float* j_pnu,
+ float* n_pnu );
+static void sph_bessel_jn( size_t n,
+ float x,
+ float* j_n,
+ float* n_n,
+ float* jp_n,
+ float* np_n );
+static void bessel_ik( float nu,
+ float x,
+ float* i_nu,
+ float* k_nu,
+ float* i_pnu,
+ float* k_pnu );
+static float cyl_bessel_ij_series( float nu,
+ float x,
+ float sgn,
+ size_t max_iter );
+
/*============================================================================*/
float _qFFMath_GetAbnormal( const int i )
{
@@ -456,6 +517,26 @@ float qFFMath_ATanh( float x )
return qFFMath_Log( ( 1.0F + x )/( 1.0F - x ) )*0.5F;
}
/*============================================================================*/
+float qFFMath_WrapToPi( float x )
+{
+ return qFFMath_Mod( x + QFFM_PI, QFFM_2PI ) - QFFM_PI;
+}
+/*============================================================================*/
+float qFFMath_WrapTo2Pi( float x )
+{
+ return qFFMath_Mod( x, QFFM_2PI );
+}
+/*============================================================================*/
+float qFFMath_WrapTo180( float x )
+{
+ return qFFMath_Mod( x + 180.0F, 360.0F ) - 180.0F;
+}
+/*============================================================================*/
+float qFFMath_WrapTo360( float x )
+{
+ return qFFMath_Mod( x, 360.0F );
+}
+/*============================================================================*/
float qFFMath_Erf( float x )
{
float retVal;
@@ -1068,4 +1149,1677 @@ float qFFMath_Factorial( float x )
return y;
}
+/*============================================================================*/
+static float poly_laguerre_recursion( size_t n,
+ float alpha,
+ float x )
+{
+ const float l0 = 1.0F;
+ float y;
+
+ if ( 0U == n ) {
+ y = l0;
+ }
+ else {
+ const float l1 = -x + 1.0F + alpha;
+ if ( 1U == n ) {
+ y = l1;
+ }
+ else {
+ float ln2 = l0;
+ float ln1 = l1;
+ float ln = 0.0F;
+ for ( size_t i = 2U ; i <= n ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float nn = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ ln = ( ( ( 2.0F*nn ) - 1.0F ) + alpha - x )*( ln1/nn )
+ - ( ( nn - 1.0F ) + alpha )*( ln2/nn );
+ ln2 = ln1;
+ ln1 = ln;
+ }
+ y = ln;
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float poly_laguerre_large_n( size_t n,
+ float alpha,
+ float x )
+{
+ const float PI_2_SQ = 2.467401100272339498076235031476244330406188964F;
+ /*cstat -CERT-FLP36-C*/
+ const float m = (float)n;
+ /*cstat +CERT-FLP36-C*/
+ const float a = -m;
+ const float b = alpha + 1.0F;
+ const float eta = ( 2.0F*b ) - ( 4.0F*a );
+ const float cos2th = x/eta;
+ const float sin2th = 1.0F - cos2th;
+ const float th = qFFMath_ACos( qFFMath_Sqrt( cos2th ) );
+ const float pre_h = PI_2_SQ*eta*eta*cos2th*sin2th;
+ const float lg_b = qFFMath_LGamma( b + m );
+ const float ln_fact = qFFMath_LGamma( m + 1.0F );
+
+ const float preTerm1 = 0.5F*( 1.0F - b )*qFFMath_Log( 0.25F*x*eta );
+ const float preTerm2 = 0.25F*qFFMath_Log( pre_h );
+ const float lnPre = lg_b - ln_fact + ( 0.5F*x ) + preTerm1 - preTerm2;
+ const float serTerm1 = qFFMath_Sin( QFFM_PI*a );
+ const float th2 = 2.0F*th;
+ const float serTerm2 = qFFMath_Sin( ( 0.25F*eta*( ( th2 ) - qFFMath_Sin( th2 ) ) + QFFM_PI_4 ) );
+
+ return qFFMath_Exp( lnPre )*( serTerm1 + serTerm2 );
+}
+/*============================================================================*/
+static float poly_laguerre_hyperg( size_t n,
+ float alpha,
+ float x )
+{
+ const float b = alpha + 1.0F;
+ const float mx = -x;
+ const float tc_sgn = ( x < 0.0F ) ? 1.0F : ( ( 1 == ( n % 2 ) ) ? -1.0F : 1.0F );
+ const float ax = qFFMath_Abs( x );
+ float tc = 1.0F;
+
+ for ( size_t i = 1U ; i <= n ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float k = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ tc *= ax/k;
+ }
+
+ float term = tc*tc_sgn;
+ float sum = term;
+ const int N = (int)n;
+ for ( int i = ( N - 1 ) ; i >= 0; --i ) {
+ /*cstat -CERT-FLP36-C -MISRAC++2008-5-0-7*/
+ const float k = (float)i;
+ term *= ( b + k )/((float)( N - i ))*( k + 1.0F )/mx;
+ /*cstat +CERT-FLP36-C +MISRAC++2008-5-0-7*/
+ sum += term;
+ }
+
+ return sum;
+}
+/*============================================================================*/
+float qFFMath_Assoc_laguerre( size_t n,
+ size_t m,
+ float x )
+{
+ float y;
+ /*cstat -CERT-FLP36-C*/
+ const float alpha = (float)m;
+ const float N = (float)n;
+ /*cstat +CERT-FLP36-C*/
+ if ( ( x < 0.0F ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( 0U == n ) {
+ y = 1.0F;
+ }
+ else if ( 1U == n ) {
+ y = 1.0F + alpha - x;
+ }
+ else if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ float prod = alpha + 1.0F;
+ for ( size_t i = 2U ; i < n ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float k = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ prod *= ( alpha + k )/k;
+ }
+ y = prod;
+ }
+ else if ( ( alpha > -1.0F ) && ( x < ( ( 2.0F*( alpha + 1.0F ) ) + ( 4.0F*N ) ) ) ) {
+ y = poly_laguerre_large_n( n, alpha, x );
+ }
+ else if ( ( ( x > 0.0F ) && ( alpha < -( N + 1.0F ) ) ) ) {
+ y = poly_laguerre_recursion( n, alpha, x );
+ }
+ else {
+ y = poly_laguerre_hyperg( n, alpha, x );
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float poly_legendre_p( size_t l,
+ float x )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( -1.0F, x ) ) {
+ y = ( 1 == ( l % 2 ) ) ? -1.0F : 1.0F;
+ }
+ else {
+ float p_lm2 = 1.0F;
+
+ if ( 0U == l ) {
+ y = p_lm2;
+ }
+ else {
+ float p_lm1 = x;
+
+ if ( 1U == l ) {
+ y = p_lm1;
+ }
+ else {
+ float p_l = 0.0F;
+
+ for ( size_t i = 2U ; i <= l ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float ll = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ p_l = ( 2.0F*x* p_lm1 ) - p_lm2 - ( ( x*p_lm1 ) - p_lm2)/ll;
+ p_lm2 = p_lm1;
+ p_lm1 = p_l;
+ }
+ y = p_l;
+ }
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Assoc_legendre( size_t n,
+ size_t m,
+ float x )
+{
+ float y;
+ const float phase = 1.0F;
+
+ if ( m > n ) {
+ y = 0.0F;
+ }
+ else if ( qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( 0U == m ) {
+ y = poly_legendre_p( n, x );
+ }
+ else {
+ float p_mm = 1.0F;
+ const float root = qFFMath_Sqrt( 1.0F - x )*qFFMath_Sqrt( 1.0F + x );
+ float fact = 1.0F;
+
+ for ( size_t i = 1U ; i <= m ; ++i ) {
+ p_mm *= phase*fact*root;
+ fact += 2.0F;
+ }
+
+ if ( n == m ) {
+ y = p_mm;
+ }
+ else {
+ /*cstat -CERT-FLP36-C*/
+ const float p_mp1m = ( 2.0F*( (float)m ) + 1.0F )*x*p_mm;
+ /*cstat +CERT-FLP36-C*/
+ if ( n == ( m + 1U ) ) {
+ y = p_mp1m;
+ }
+ else {
+ float p_lm2m = p_mm;
+ float p_lm1m = p_mp1m;
+ float p_lm = 0.0F;
+ /*cstat -CERT-FLP36-C*/
+ const float M = (float)m;
+ for ( size_t i = ( m + 2U ) ; i <= n ; ++i ) {
+ const float j = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ p_lm = ( ( ( 2.0F*j ) - 1.0F )*x*p_lm1m ) - ( ( j + M - 1.0F )*p_lm2m/( j - M ) );
+ p_lm2m = p_lm1m;
+ p_lm1m = p_lm;
+ }
+ y = p_lm;
+ }
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Beta( float x,
+ float y )
+{
+ float result;
+
+ if ( qFFMath_IsNaN( x ) || qFFMath_IsNaN( y ) ) {
+ result = QFFM_NAN;
+ }
+ else {
+ const float bet = qFFMath_LGamma( x ) + qFFMath_LGamma( y ) - qFFMath_LGamma( x + y );
+ result = qFFMath_Exp( bet );
+ }
+
+ return result;
+}
+/*============================================================================*/
+static float ellint_rf( float x,
+ float y,
+ float z )
+{
+ const float loLim = 5.0F*FLT_MIN;
+ float result;
+
+
+ if ( ( x < 0.0F ) || ( y < 0.0F ) || ( z < 0.0F ) ||
+ ( ( x + y ) < loLim ) || ( ( x + z ) < loLim ) || ( ( y + z) < loLim )
+ ) {
+ result = QFFM_NAN;
+ }
+ else {
+ const float c0 = 1.0F/4.0F;
+ const float c1 = 1.0F/24.0F;
+ const float c2 = 1.0F/10.0F;
+ const float c3 = 3.0F/44.0F;
+ const float c4 = 1.0F/14.0F;
+ const float errTol = 0.0024607833005759250678823324F;
+ const float c13 = 1.0F/3.0F;
+
+ float xn = x;
+ float yn = y;
+ float zn = z;
+ float mu, xnDev, ynDev, znDev;
+ const size_t maxIter = 100U;
+ float e2, e3, s;
+
+ for ( size_t iter = 0U ; iter < maxIter ; ++iter ) {
+ float abs_xnDev, abs_ynDev, abs_znDev, lambda, epsilon;
+ float xRoot, yRoot, zRoot;
+
+ mu = ( xn + yn + zn )*c13;
+ xnDev = 2.0F - ( mu + xn )/mu;
+ ynDev = 2.0F - ( mu + yn )/mu;
+ znDev = 2.0F - ( mu + zn )/mu;
+ abs_xnDev = qFFMath_Abs( xnDev );
+ abs_ynDev = qFFMath_Abs( ynDev );
+ abs_znDev = qFFMath_Abs( znDev );
+ epsilon = ( abs_xnDev > abs_ynDev ) ? abs_xnDev : abs_ynDev;
+ epsilon = ( abs_znDev > epsilon ) ? abs_znDev : epsilon;
+ if ( epsilon < errTol ) {
+ break;
+ }
+ xRoot = qFFMath_Sqrt( xn );
+ yRoot = qFFMath_Sqrt( yn );
+ zRoot = qFFMath_Sqrt( zn );
+ lambda = xRoot*( yRoot + zRoot ) + ( yRoot*zRoot );
+ xn = c0*( xn + lambda );
+ yn = c0*( yn + lambda );
+ zn = c0*( zn + lambda );
+ }
+ e2 = xnDev*ynDev;
+ e3 = ( e2*znDev );
+ e2 = e2 - ( znDev*znDev );
+ s = 1.0F + ( ( c1*e2 ) - c2 - ( c3*e3 ) )*e2 + ( c4*e3 );
+ result = s/qFFMath_Sqrt( mu );
+ }
+
+ return result;
+}
+/*============================================================================*/
+static float ellint_rd( float x,
+ float y,
+ float z )
+{
+ const float errTol = 0.0017400365588678507952624663346341549F;
+ const float loLim = 4.103335708781587555782386855921935e-26F;
+ float result;
+
+ if ( ( x < 0.0F ) || ( y < 0.0F ) || ( ( x + y ) < loLim ) || ( z < loLim ) ) {
+ result = QFFM_NAN;
+ }
+ else {
+ const float c0 = 1.0F/4.0F;
+ const float c1 = 3.0F/14.0F;
+ const float c2 = 1.0F/6.0F;
+ const float c3 = 9.0F/22.0F;
+ const float c4 = 3.0F/26.0F;
+ float xn = x;
+ float yn = y;
+ float zn = z;
+ float sigma = 0.0F;
+ float power4 = 1.0F;
+ float mu, xnDev, ynDev, znDev;
+ float ea, eb, ec, ed, ef, s1, s2;
+ const size_t maxIter = 100U;
+
+ for ( size_t iter = 0U ; iter < maxIter ; ++iter ) {
+ float abs_xnDev, abs_ynDev, abs_znDev, lambda, epsilon;
+ float xRoot, yRoot, zRoot;
+
+ mu = ( xn + yn + ( 3.0F*zn ) )*0.2F;
+ xnDev = ( mu - xn )/mu;
+ ynDev = ( mu - yn )/mu;
+ znDev = ( mu - zn )/mu;
+ abs_xnDev = qFFMath_Abs( xnDev );
+ abs_ynDev = qFFMath_Abs( ynDev );
+ abs_znDev = qFFMath_Abs( znDev );
+ epsilon = ( abs_xnDev > abs_ynDev ) ? abs_xnDev : abs_ynDev;
+ epsilon = ( abs_znDev > epsilon ) ? abs_znDev : epsilon;
+ if ( epsilon < errTol ) {
+ break;
+ }
+ xRoot = qFFMath_Sqrt( xn );
+ yRoot = qFFMath_Sqrt( yn );
+ zRoot = qFFMath_Sqrt( zn );
+ lambda = xRoot*( yRoot + zRoot ) + ( yRoot*zRoot );
+ sigma += power4/( zRoot*( zn + lambda ) );
+ power4 *= c0;
+ xn = c0*( xn + lambda );
+ yn = c0*( yn + lambda );
+ zn = c0*( zn + lambda );
+ }
+ ea = xnDev*ynDev;
+ eb = znDev*znDev;
+ ec = ea - eb;
+ ed = ea - ( 6.0F*eb );
+ ef = ed + ec + ec;
+ s1 = ed*( -c1 + ( c3*ed/3.0F ) - ( 1.5F*c4*znDev*ef ) );
+ s2 = znDev*( ( c2*ef ) + znDev*( -( c3*ec ) - ( znDev*c4 ) - ea ) );
+ result = ( 3.0F*sigma ) + power4*qFFMath_RSqrt( mu )*( 1.0F + s1 + s2)/mu;
+ }
+
+ return result;
+}
+/*============================================================================*/
+static float ellint_rc( float x,
+ float y )
+{
+ float result;
+ const float loLim = 5.8774717550000002558112628881984982848919e-38F;
+ const float errTol = 0.049606282877419791144113503378321F;
+
+ if ( ( x < 0.0F ) || ( y < 0.0F ) || ( y < loLim ) ) {
+ result = QFFM_NAN;
+ }
+ else {
+ const float c0 = 1.0F/4.0F;
+ const float c1 = 1.0F/7.0F;
+ const float c2 = 9.0F/22.0F;
+ const float c3 = 3.0F/10.0F;
+ const float c4 = 3.0F/8.0F;
+ const float c13 = 1.0F/3.0F;
+ float xn = x;
+ float yn = y;
+ const size_t maxIter = 100;
+ float mu, s, sn;
+
+ for ( size_t iter = 0U ; iter < maxIter ; ++iter ) {
+ float lambda;
+
+ mu = ( xn + 2.0F*yn )*c13;
+ sn = ( yn + mu )/mu - 2.0F;
+ if ( qFFMath_Abs( sn ) < errTol ){
+ break;
+ }
+ lambda = ( 2.0F*qFFMath_Sqrt( xn )*qFFMath_Sqrt( yn ) ) + yn;
+ xn = c0*( xn + lambda );
+ yn = c0*( yn + lambda );
+ }
+ s = sn*sn*( c3 + sn*( c1 + sn*( c4 + ( sn*c2 ) ) ) );
+ result = ( 1.0F + s )*qFFMath_RSqrt( mu );
+ }
+
+ return result;
+}
+/*============================================================================*/
+static float ellint_rj( float x,
+ float y,
+ float z,
+ float p )
+{
+ float result;
+ const float loLim = 4.103335708781587555782386855921935e-26F;
+ const float errTol = 0.049606282877419791144113503378321F;
+
+ if ( ( x < 0.0F ) || ( y < 0.0F ) || ( z < 0.0F ) || ( ( x + y ) < loLim ) ||
+ ( ( x + z ) < loLim ) || ( ( y + z) < loLim )|| ( p < loLim ) ) {
+ result = QFFM_NAN;
+ }
+ else {
+ const float c0 = 1.0F/4.0F;
+ const float c1 = 3.0F/14.0F;
+ const float c2 = 1.0F/3.0F;
+ const float c3 = 3.0F/22.0F;
+ const float c4 = 3.0F/26.0F;
+ float xn = x;
+ float yn = y;
+ float zn = z;
+ float pn = p;
+ float sigma = 0.0F;
+ float power4 = 1.0F;
+ float mu, xnDev, ynDev, znDev, pnDev;
+ float ea, eb, ec, e2, e3, s1, s2, s3;
+ const size_t maxIter = 100;
+
+ for ( size_t iter = 0U ; iter < maxIter ; ++iter ) {
+ float abs_xnDev, abs_ynDev, abs_znDev, abs_pnDev, lambda, alpha1;
+ float alpha2, beta, epsilon;
+ float xRoot, yRoot, zRoot;
+ mu = 0.2F*( xn + yn + zn + ( 2.0F*pn ) );
+ xnDev = ( mu - xn )/mu;
+ ynDev = ( mu - yn )/mu;
+ znDev = ( mu - zn )/mu;
+ pnDev = ( mu - pn )/mu;
+ abs_xnDev = qFFMath_Abs( xnDev );
+ abs_ynDev = qFFMath_Abs( ynDev );
+ abs_znDev = qFFMath_Abs( znDev );
+ abs_pnDev = qFFMath_Abs( pnDev );
+ epsilon = ( abs_xnDev > abs_ynDev ) ? abs_xnDev : abs_ynDev;
+ epsilon = ( abs_znDev > epsilon ) ? abs_znDev : epsilon;
+ epsilon = ( abs_pnDev > epsilon ) ? abs_pnDev : epsilon;
+ if ( epsilon < errTol ) {
+ break;
+ }
+ xRoot = qFFMath_Sqrt( xn );
+ yRoot = qFFMath_Sqrt( yn );
+ zRoot = qFFMath_Sqrt( zn );
+ lambda = xRoot*( yRoot + zRoot ) + ( yRoot*zRoot );
+ alpha1 = pn*( xRoot + yRoot + zRoot ) + xRoot*yRoot*zRoot;
+ alpha2 = alpha1*alpha1;
+ beta = pn*( pn + lambda )*( pn + lambda );
+ sigma += power4*ellint_rc( alpha2, beta );
+ power4 *= c0;
+ xn = c0*( xn + lambda );
+ yn = c0*( yn + lambda );
+ zn = c0*( zn + lambda );
+ pn = c0*( pn + lambda );
+ }
+ ea = xnDev*( ynDev + znDev ) + ynDev*znDev;
+ eb = xnDev*ynDev*znDev;
+ ec = pnDev*pnDev;
+ e2 = ea - 3.0F*ec;
+ e3 = eb + 2.0F*pnDev*( ea - ec );
+ s1 = 1.0F + e2*( -c1 + 0.75F*c3*e2 - 1.5F*c4*e3 );
+ s2 = eb*( 0.5F*c2 + pnDev*( -c3 - c3 + pnDev*c4 ) );
+ s3 = pnDev*ea*( c2 - ( pnDev*c3 ) ) - ( c2*pnDev*ec );
+ result = 3.0F*sigma + power4*( s1 + s2 + s3)/( mu * qFFMath_Sqrt( mu ) );
+ }
+
+ return result;
+}
+/*============================================================================*/
+float qFFMath_Comp_ellint_1( float k )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( k ) || ( qFFMath_Abs( k ) >= 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ y = ellint_rf( 0.0F, 1.0F - ( k*k ), 1.0F );
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Comp_ellint_2( float k )
+{
+ float y;
+ const float abs_k = qFFMath_Abs( k );
+
+ if ( qFFMath_IsNaN( k ) || ( abs_k > 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, abs_k ) ) {
+ y = 1.0F;
+ }
+ else {
+ const float kk = k*k;
+ const float one_m_kk = 1.0F - kk;
+ const float c13 = 1.0F/3.0F;
+
+ y = ellint_rf( 0.0F, one_m_kk, 1.0F ) - c13*kk*ellint_rd( 0.0F, one_m_kk, 1.0F );
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Comp_ellint_3( float k,
+ float nu )
+{
+ float y;
+ const float abs_k = qFFMath_Abs( k );
+
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( nu ) || ( abs_k > 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, nu ) ) {
+ y = QFFM_INFINITY;
+ }
+ else {
+ const float kk = k*k;
+ const float one_m_kk = 1.0F - kk;
+ const float c13 = 1.0F/3.0F;
+
+ y = ellint_rf( 0.0F, one_m_kk, 1.0F ) + c13*nu*ellint_rj( 0.0F, one_m_kk, 1.0F, 1.0F - nu );
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_1( float k,
+ float phi )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( phi ) || ( qFFMath_Abs(k) > 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ const float n = qFFMath_Floor( phi/QFFM_PI + 0.5F );
+ const float phi_red = phi - n*QFFM_PI;
+ const float s = qFFMath_Sin( phi_red );
+ const float c = qFFMath_Cos( phi_red );
+ const float f = s*ellint_rf( c*c, 1.0F - k*k*s*s, 1.0F );
+ if ( QFFM_FP_ZERO == qFFMath_FPClassify( n ) ) {
+ y = f;
+ }
+ else {
+ y = f + ( 2.0F*n*qFFMath_Comp_ellint_1( k ) );
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_2( float k,
+ float phi )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( phi ) || ( qFFMath_Abs( k ) > 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ const float c13 = 1.0F/3.0F;
+ const float n = qFFMath_Floor( phi/QFFM_PI + 0.5F );
+ const float phi_red = phi - n*QFFM_PI;
+ const float kk = k*k;
+ const float s = qFFMath_Sin( phi_red );
+ const float ss = s*s;
+ const float sss = ss*s;
+ const float c = qFFMath_Cos( phi_red );
+ const float cc = c*c;
+ const float tmp = 1.0F - kk*ss;
+ const float e = s*ellint_rf( cc, tmp, 1.0F )
+ - c13*kk*sss*ellint_rd( cc, tmp, 1.0F );
+
+ if ( QFFM_FP_ZERO == qFFMath_FPClassify( n ) ) {
+ y = e;
+ }
+ else {
+ y = e + ( 2.0F*n*qFFMath_Comp_ellint_2( k ) );
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_3( float k,
+ float nu,
+ float phi )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( phi ) || ( qFFMath_Abs( k ) > 1.0F ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ const float n = qFFMath_Floor( phi/QFFM_PI + 0.5F );
+ const float phi_red = phi - n*QFFM_PI;
+ const float kk = k*k;
+ const float s = qFFMath_Sin( phi_red );
+ const float ss = s*s;
+ const float sss = ss*s;
+ const float c = qFFMath_Cos( phi_red );
+ const float cc = c*c;
+ const float tmp = 1.0F - kk*ss;
+ const float c13 = 1.0F/3.0F;
+ const float pi = s*ellint_rf( cc, tmp, 1.0F ) + c13*nu*sss*ellint_rj( cc, tmp, 1.0F, 1.0F - nu*ss );
+
+ if ( QFFM_FP_ZERO == qFFMath_FPClassify( n ) ) {
+ y = pi;
+ }
+ else {
+ y = pi + ( 2.0F*n*qFFMath_Comp_ellint_3( k, nu ) );
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float expint_E1_series( float x )
+{
+ float term = 1.0F;
+ float eSum = 0.0F;
+ float oSum = 0.0F;
+ const size_t maxIter = 1000U;
+
+ for ( size_t i = 1U; i < maxIter; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float j = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ term *= -x/j;
+ if ( qFFMath_Abs( term ) < FLT_EPSILON ) {
+ break;
+ }
+ if ( term >= 0.0F ) {
+ eSum += term/j;
+ }
+ else {
+ oSum += term/j;
+ }
+ }
+ return - eSum - oSum - QFFM_GAMMA_E - qFFMath_Log( x );
+}
+/*============================================================================*/
+static float expint_E1_asymp( float x )
+{
+ float term = 1.0F;
+ float eSum = 1.0F;
+ float oSum = 0.0F;
+ const size_t maxIter = 1000U;
+
+ for ( size_t i = 1U; i < maxIter; ++i ) {
+ const float prev = term;
+ /*cstat -CERT-FLP36-C*/
+ term *= -( (float)i )/x;
+ /*cstat +CERT-FLP36-C*/
+ if ( qFFMath_Abs( term ) > qFFMath_Abs( prev ) ) {
+ break;
+ }
+ if ( term >= 0.0F ) {
+ eSum += term;
+ }
+ else {
+ oSum += term;
+ }
+ }
+ return qFFMath_Exp( -x )*( eSum + oSum )/x;
+}
+/*============================================================================*/
+static float expint_Ei_asymp( float x )
+{
+ float term = 1.0F;
+ float sum = 1.0F;
+ const size_t maxIter = 1000U;
+
+ for ( size_t i = 1U; i < maxIter; ++i ) {
+ const float prev = term;
+ /*cstat -CERT-FLP36-C*/
+ term *= -( (float)i )/x;
+ /*cstat +CERT-FLP36-C*/
+ if ( ( term < FLT_EPSILON ) || ( term >= prev ) ) {
+ break;
+ }
+ sum += term;
+ }
+
+ return qFFMath_Exp( x )*sum/x;
+}
+/*============================================================================*/
+static float expint_E1( float x )
+{
+ float y;
+
+ if ( x < 0.0F ) {
+ y = -expint_Ei( -x );
+ }
+ else if ( x < 1.0F ) {
+ y = expint_E1_series( x );
+ }
+ else if ( x < 100.F ) {
+ y = expint_En_cont_frac( 1, x );
+ }
+ else {
+ y = expint_E1_asymp( x );
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float expint_En_cont_frac( size_t n,
+ float x )
+{
+ float y = QFFM_NAN;
+ const int maxIter = 1000;
+ const int nm1 = (int)n - 1;
+ /*cstat -CERT-FLP36-C*/
+ float b = x + (float)n;
+ /*cstat +CERT-FLP36-C*/
+ float c = 1.0F/FLT_MIN;
+ float d = 1.0F/b;
+ float h = d;
+
+ for ( int i = 1; i <= maxIter; ++i ) {
+ /*cstat -MISRAC++2008-5-0-7 -CERT-FLP36-C*/
+ const float a = -( (float)( i*( nm1 + i ) ) );
+ /*cstat MISRAC++2008-5-0-7 +CERT-FLP36-C*/
+ b += 2.0F;
+ d = 1.0F/( ( a*d ) + b );
+ c = b + ( a/c );
+ const float del = c*d;
+ h *= del;
+ if ( qFFMath_Abs( del - 1.0F ) < FLT_EPSILON ) {
+ y = h*qFFMath_Exp( -x );
+ break;
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float expint_Ei_series( float x )
+{
+ float term = 1.0F;
+ float sum = 0.0F;
+ const size_t maxIter = 1000U;
+
+ for ( size_t i = 1U; i < maxIter; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float j = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ term *= x/j;
+ sum += term/j;
+ if ( term < ( FLT_EPSILON*sum ) ) {
+ break;
+ }
+ }
+ return QFFM_GAMMA_E + sum + qFFMath_Log( x );
+}
+/*============================================================================*/
+static float expint_Ei( float x )
+{
+ const float logEps = 36.044F;
+ float y;
+
+ if ( x < 0.0F ) {
+ y = -expint_E1( -x );
+ }
+ else if ( x < logEps ) {
+ y = expint_Ei_series( x );
+ }
+ else {
+ y = expint_Ei_asymp( x );
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Expint( float num )
+{
+ return ( qFFMath_IsNaN( num ) ) ? QFFM_NAN : expint_Ei( num );
+}
+/*============================================================================*/
+float qFFMath_Hermite( size_t n,
+ float x )
+{
+ float y = 0.0F;
+
+ if ( qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ const float H_0 = 1.0F;
+
+ if ( 0U == n ) {
+ y= H_0;
+ }
+ else {
+ const float H_1 = 2.0F*x;
+ if ( 1U == n ) {
+ y = H_1;
+ }
+ else {
+ float H_nm1, H_nm2;
+
+ H_nm2 = H_0;
+ H_nm1 = H_1;
+ for ( size_t i = 2U; i <= n; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float j = (float)( i - 1U );
+ /*cstat +CERT-FLP36-C*/
+ y = 2.0F*( ( x*H_nm1 ) - ( j*H_nm2 ) );
+ H_nm2 = H_nm1;
+ H_nm1 = y;
+ }
+ }
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Laguerre( size_t n,
+ float x )
+{
+ return qFFMath_Assoc_laguerre( n, 0, x );
+}
+/*============================================================================*/
+float qFFMath_Legendre( size_t n,
+ float x )
+{
+ float y = 0.0F;
+
+ if ( qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, x ) ) {
+ y = 1.0F;
+ }
+ else if ( qFFMath_IsEqual( -1.0F, x ) ) {
+ y = ( ( n % 2 ) == 1 ) ? -1.0F : 1.0F;
+ }
+ else {
+ float p_lm2 = 1.0F;
+
+ if ( 0 == n ) {
+ y = p_lm2;
+ }
+ else {
+ float p_lm1 = x;
+
+ if ( 1 == n ) {
+ y = p_lm1;
+ }
+ else {
+ for ( size_t ll = 2U ; ll <= n; ++ll ) {
+ /*cstat -CERT-FLP36-C*/
+ const float ll_f = (float)ll;
+ /*cstat +CERT-FLP36-C*/
+ y = ( 2.0F*x*p_lm1 ) - p_lm2 - ( x*p_lm1 - p_lm2)/ll_f;
+ p_lm2 = p_lm1;
+ p_lm1 = y;
+ }
+ }
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float riemann_zeta_glob( float s )
+{
+ const float maxBinCoeff = 86.4982335337F;
+ float zeta = 0.0F;
+ const float ss = s;
+ bool neg = false;
+
+ if ( s < 0.0F ) {
+ s = 1.0F - s;
+ neg = true;
+ }
+ float num = 0.5F;
+ const size_t maxIt = 10000U;
+ /*cstat -MISRAC++2008-6-6-4*/
+ for ( size_t i = 0U ; i < maxIt; ++i ) {
+ bool punt = false;
+ float sgn = 1.0F;
+ float term = 0.0F;
+ for ( size_t j = 0U ; j <= i ; ++j ) {
+ /*cstat -CERT-FLP36-C*/
+ const float ii = (float)i;
+ const float jj = (float)j;
+ /*cstat +CERT-FLP36-C*/
+ float bin_coeff = qFFMath_LGamma( 1.0F + ii ) -
+ qFFMath_LGamma( 1.0F + jj ) -
+ qFFMath_LGamma( 1.0F + ii - jj );
+
+ if ( bin_coeff > maxBinCoeff ) {
+ punt = true;
+ break;
+ }
+ bin_coeff = qFFMath_Exp( bin_coeff );
+ term += sgn*bin_coeff*qFFMath_Pow( 1.0F + jj, -s );
+ sgn *= -1.0F;
+ }
+ if ( punt ) {
+ break;
+ }
+ term *= num;
+ zeta += term;
+ if ( qFFMath_Abs( term/zeta ) < FLT_EPSILON ) {
+ break;
+ }
+ num *= 0.5F;
+ }
+ /*cstat +MISRAC++2008-6-6-4*/
+ zeta /= 1.0F - qFFMath_Pow( 2.0F, 1.0F - s );
+
+ if ( neg ) {
+ zeta *= qFFMath_Pow( 2.0F*QFFM_PI, ss )*
+ qFFMath_Sin( QFFM_PI_2*ss )*
+ qFFMath_Exp( qFFMath_LGamma( s ) )/QFFM_PI;
+ }
+ return zeta;
+}
+/*============================================================================*/
+static float riemann_zeta_product( float s )
+{
+ static const uint8_t prime[ 29 ] = {
+ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
+ 71, 73, 79, 83, 89, 97, 101, 103, 107, 109
+ };
+ static const size_t num_primes = sizeof(prime)/sizeof(float);
+ float zeta = 1.0F;
+
+ for ( size_t i = 0U ; i < num_primes; ++i ) {
+ const float f = 1.0F - qFFMath_Pow( (float)prime[ i ], -s );
+
+ zeta *= f;
+ if ( ( 1.0F - f ) < FLT_EPSILON ) {
+ break;
+ }
+ }
+ zeta = 1.0F/zeta;
+
+ return zeta;
+}
+/*============================================================================*/
+float qFFMath_Riemann_zeta( float s )
+{
+ float z;
+
+ if ( qFFMath_IsNaN( s ) ) {
+ z = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, s ) ) {
+ z = QFFM_INFINITY;
+ }
+ else if ( s < -19.0F ) {
+ z = riemann_zeta_product( 1.0F - s );
+ z *= qFFMath_Pow( 2.0F*QFFM_PI, s )*
+ qFFMath_Sin( QFFM_PI_2*s )*
+ qFFMath_Exp( qFFMath_LGamma( 1.0F - s ) )/
+ QFFM_PI;
+ }
+ else if ( s < 20.0F ) {
+ z = riemann_zeta_glob( s );
+ }
+ else {
+ z = riemann_zeta_product( s );
+ }
+
+ return z;
+}
+/*============================================================================*/
+static void gamma_temme( float mu,
+ float* gam1,
+ float* gam2,
+ float* gam_pl,
+ float* gam_mi)
+{
+
+ gam_pl[ 0 ] = 1.0F/qFFMath_TGamma( 1.0F + mu );
+ gam_mi[ 0 ] = 1.0F/qFFMath_TGamma( 1.0F - mu );
+
+ if ( qFFMath_Abs( mu ) < FLT_EPSILON ) {
+ gam1[ 0 ] = -QFFM_GAMMA_E;
+ }
+ else {
+ gam1[ 0 ] = ( gam_mi[ 0 ] - gam_pl[ 0 ] )/( 2.0F*mu );
+ }
+ gam2[ 0 ] = 0.5F*( gam_mi[ 0 ] + gam_pl[ 0 ] );
+}
+/*============================================================================*/
+static void bessel_jn( float nu,
+ float x,
+ float* j_nu,
+ float* n_nu,
+ float* j_pnu,
+ float* n_pnu )
+{
+ if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ if ( qFFMath_IsEqual( 0.0F, nu ) ) {
+ j_nu[ 0 ] = 1.0F;
+ j_pnu[ 0 ] = 0.0F;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, nu ) ) {
+ j_nu[ 0 ] = 0.0F;
+ j_pnu[ 0 ] = 0.5F;
+ }
+ else {
+ j_nu[ 0 ] = 0.0F;
+ j_pnu[ 0 ] = 0.0F;
+ }
+ n_nu[ 0 ] = -QFFM_INFINITY;
+ n_pnu[ 0 ] = QFFM_INFINITY;
+ }
+ else {
+ const float eps = FLT_EPSILON;
+ const float fp_min = 1.08420217256745973463717809E-19F;
+ const int max_iter = 15000;
+ const float x_min = 2.0F;
+ /*cstat -CERT-FLP34-C*/
+ const int tmp_nl = (int)( nu - x + 1.5F );
+ const int nl = ( x < x_min ) ? (int)( nu + 0.5F )
+ : ( ( tmp_nl > 0 )? tmp_nl : 0 );
+ /*cstat +CERT-FLP34-C -CERT-FLP36-C*/
+ const float mu = nu - (float)nl;
+ /*cstat +CERT-FLP36-C*/
+ const float mu2 = mu*mu;
+ const float xi = 1.0F/x;
+ const float xi2 = 2.0F*xi;
+ const float w = xi2/QFFM_PI;
+ float iSign = 1.0F;
+ float h = nu*xi;
+ if ( h < fp_min ) {
+ h = fp_min;
+ }
+ float b = xi2*nu;
+ float d = 0.0F;
+ float c = h;
+ for ( int i = 1; i <= max_iter; ++i ) {
+ b += xi2;
+ d = b - d;
+ if ( qFFMath_Abs( d ) < fp_min ) {
+ d = fp_min;
+ }
+ c = b - ( 1.0F/c );
+ if ( qFFMath_Abs( c ) < fp_min ) {
+ c = fp_min;
+ }
+ d = 1.0F / d;
+
+ const float del = c * d;
+ h *= del;
+ if ( d < 0.0F ) {
+ iSign = -iSign;
+ }
+ if ( qFFMath_Abs( del - 1.0F ) < eps ) {
+ break;
+ }
+ }
+ float j_nul = iSign*fp_min;
+ float j_pnu_l = h*j_nul;
+ const float j_nul1 = j_nul;
+ const float j_pnu1 = j_pnu_l;
+ float fact = nu*xi;
+
+ for ( int l = nl; l >= 1; --l ) {
+ const float j_nu_temp = ( fact*j_nul ) + j_pnu_l;
+
+ fact -= xi;
+ j_pnu_l = ( fact*j_nu_temp ) - j_nul;
+ j_nul = j_nu_temp;
+ }
+ if ( qFFMath_IsEqual( 0.0F, j_nul ) ) {
+ j_nul = eps;
+ }
+ const float f = j_pnu_l/j_nul;
+ float n_mu, n_nu1, n_pmu, Jmu;
+ if ( x < x_min ) {
+ const float x2 = 0.5F*x;
+ const float pi_mu = QFFM_PI*mu;
+ const float fact_l = ( qFFMath_Abs( pi_mu ) < eps ) ? 1.0F : pi_mu/qFFMath_Sin( pi_mu );
+ d = -qFFMath_Log( x2 );
+ float e = mu*d;
+ const float fact2 = ( qFFMath_Abs( e ) < eps ) ? 1.0F : qFFMath_Sinh(e)/e;
+ float gam1, gam2, gam_pl, gam_mi;
+
+ gamma_temme( mu, &gam1, &gam2, &gam_pl, &gam_mi );
+ float ff = ( 2.0F/QFFM_PI )*fact_l*( gam1*qFFMath_Cosh( e ) + ( gam2*fact2*d ) );
+ e = qFFMath_Exp( e );
+ float p = e/( QFFM_PI*gam_pl );
+ float q = 1.0F/( e*QFFM_PI*gam_mi );
+ const float pi_mu2 = pi_mu / 2.0F;
+ const float fact3 = ( qFFMath_Abs( pi_mu2 ) < eps ) ? 1.0F : qFFMath_Sin( pi_mu2 )/pi_mu2;
+ const float r = QFFM_PI*pi_mu2*fact3*fact3;
+ float sum = ff + ( r*q );
+ float sum1 = p;
+ c = 1.0F;
+ d = -x2*x2;
+ for ( int i = 1; i <= max_iter ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float j = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ ff = ( j*ff + p + q )/( j*j - mu2 );
+ c *= d/j;
+ p /= j - mu;
+ q /= j + mu;
+ const float del = c*( ff + ( r*q ) );
+ sum += del;
+ const float del1 = ( c*p ) - ( j*del );
+ sum1 += del1;
+ if ( qFFMath_Abs( del ) < ( eps*( 1.0F + qFFMath_Abs( sum ) ) ) ) {
+ break;
+ }
+ }
+
+ n_mu = -sum;
+ n_nu1 = -sum1*xi2;
+ n_pmu = ( mu*xi*n_mu ) - n_nu1;
+ Jmu = w/( n_pmu - ( f*n_mu ) );
+ }
+ else {
+ float a = 0.25F - mu2;
+ float q = 1.0F;
+ float p = -xi*0.5F;
+ const float br = 2.0F*x;
+ float bi = 2.0F;
+ float fact_g = a*xi/( ( p*p ) + ( q*q ) );
+ float cr = br + ( q*fact_g );
+ float ci = bi + ( p*fact_g );
+ float den = ( br*br ) + ( bi*bi );
+ float dr = br/den;
+ float di = -bi/den;
+ float dlr = ( cr*dr ) - ( ci*di );
+ float dli = ( cr*di ) + ( ci*dr );
+ float temp = ( p*dlr ) - ( q*dli );
+ q = ( p*dli ) + ( q * dlr );
+ p = temp;
+
+ for ( int i = 2; i <= max_iter; ++i ) {
+ a += (float)( 2*( i - 1 ) );
+ bi += 2.0F;
+ dr = ( a*dr ) + br;
+ di = ( a*di ) + bi;
+ if ( ( qFFMath_Abs( dr ) + qFFMath_Abs( di ) ) < fp_min ) {
+ dr = fp_min;
+ }
+ fact_g = a/( ( cr*cr ) + ( ci*ci ) );
+ cr = br + ( cr*fact_g );
+ ci = bi - ( ci*fact_g );
+ if ( ( qFFMath_Abs( cr ) + qFFMath_Abs( ci ) ) < fp_min ) {
+ cr = fp_min;
+ }
+ den = ( dr*dr ) + ( di*di );
+ dr /= den;
+ di /= -den;
+ dlr = ( cr*dr ) - ( ci*di );
+ dli = ( cr*di ) + ( ci*dr );
+ temp = ( p*dlr ) - ( q*dli );
+ q = ( p*dli ) + ( q*dlr );
+ p = temp;
+ if ( qFFMath_Abs( dlr - 1.0F ) + qFFMath_Abs( dli ) < eps ) {
+ break;
+ }
+ }
+ const float gam = ( p - f )/q;
+
+ Jmu = qFFMath_Sqrt( w/( ( p - f )*gam + q ) );
+ if ( ( Jmu*j_nul ) < 0.0F) {
+ Jmu = -Jmu;
+ }
+ n_mu = gam*Jmu;
+ n_pmu = ( p + ( q/gam ) )*n_mu;
+ n_nu1 = ( mu*xi*n_mu ) - n_pmu;
+ }
+ fact = Jmu/j_nul;
+ j_nu[ 0 ] = fact*j_nul1;
+ j_pnu[ 0 ] = fact*j_pnu1;
+ for ( int i = 1; i <= nl; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float n_nu_temp = ( mu + (float)i )*xi2*n_nu1 - n_mu;
+ /*cstat +CERT-FLP36-C*/
+ n_mu = n_nu1;
+ n_nu1 = n_nu_temp;
+ }
+ n_nu[ 0 ] = n_mu;
+ n_pnu[ 0 ] = ( nu*xi*n_mu ) - n_nu1;
+ }
+}
+/*============================================================================*/
+static void sph_bessel_jn( size_t n,
+ float x,
+ float* j_n,
+ float* n_n,
+ float* jp_n,
+ float* np_n )
+{
+ /*cstat -CERT-FLP36-C*/
+ const float nu = (float)n + 0.5F;
+ /*cstat +CERT-FLP36-C*/
+ const float sqrtpi2 = 1.25331413731550012080617761967005208134651184F;
+ float j_nu, n_nu, jp_nu, np_nu;
+ const float factor = sqrtpi2*qFFMath_RSqrt( x );
+ const float inv_2x = 1.0F/( 2.0F*x );
+
+ bessel_jn( nu, x, &j_nu, &n_nu, &jp_nu, &np_nu );
+ j_n[ 0 ] = factor*j_nu;
+ n_n[ 0 ] = factor*n_nu;
+ jp_n[ 0 ] = ( factor*jp_nu ) - ( j_n[ 0 ]*inv_2x );
+ np_n[ 0 ] = ( factor*np_nu ) - ( n_n[ 0 ]*inv_2x );
+}
+/*============================================================================*/
+float qFFMath_Sph_bessel( size_t n,
+ float x )
+{
+ float y;
+
+ if ( ( x < 0.0F ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ y = ( 0U == n ) ? 1.0F : 0.0F;
+ }
+ else {
+ float j_n, n_n, jp_n, np_n;
+
+ sph_bessel_jn( n, x, &j_n, &n_n, &jp_n, &np_n );
+ y = j_n;
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Sph_neumann( size_t n,
+ float x )
+{
+ float y;
+
+ if ( ( x < 0.0F ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ y = -QFFM_INFINITY;
+ }
+ else {
+ float j_n, n_n, jp_n, np_n;
+
+ sph_bessel_jn( n, x, &j_n, &n_n, &jp_n, &np_n );
+ y = n_n;
+ }
+ return y;
+}
+/*============================================================================*/
+static void bessel_ik( float nu,
+ float x,
+ float* i_nu,
+ float* k_nu,
+ float* i_pnu,
+ float* k_pnu )
+{
+ if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ if ( qFFMath_IsEqual( 0.0F, nu ) ) {
+ i_nu[ 0 ] = 1.0F;
+ i_pnu[ 0 ] = 0.0F;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, x ) ) {
+ i_nu[ 0 ] = 0.0F;
+ i_pnu[ 0 ] = 0.5F;
+ }
+ else {
+ i_nu[ 0 ] = 0.0F;
+ i_pnu[ 0 ] = 0.0F;
+ }
+ k_nu[ 0 ] = QFFM_INFINITY;
+ k_pnu[ 0 ] = -QFFM_INFINITY;
+ }
+ else {
+ const float eps = FLT_EPSILON;
+ const float fp_min = 10.0F*FLT_EPSILON;
+ const int max_iter = 15000;
+ const float x_min = 2.0F;
+ /*cstat -CERT-FLP34-C -CERT-FLP36-C*/
+ const int nl = (int)( nu + 0.5F );
+ const float mu = nu - (float)nl;
+ /*cstat +CERT-FLP34-C +CERT-FLP36-C*/
+ const float mu2 = mu*mu;
+ const float xi = 1.0F/x;
+ const float xi2 = 2.0F*xi;
+ float h = nu*xi;
+
+ if ( h < fp_min ) {
+ h = fp_min;
+ }
+ float b = xi2*nu;
+ float d = 0.0F;
+ float c = h;
+
+ for ( int i = 1; i <= max_iter; ++i ) {
+ b += xi2;
+ d = 1.0F/( b + d );
+ c = b + ( 1.0F/c );
+ const float del = c*d;
+ h *= del;
+ if ( qFFMath_Abs( del - 1.0F ) < eps ) {
+ break;
+ }
+ }
+
+ float i_nul = fp_min;
+ float i_pnu_l = h*i_nul;
+ const float i_nul1 = i_nul;
+ const float i_pnu_1 = i_pnu_l;
+ float fact_m = nu*xi;
+
+ for ( int l = nl ; l >= 1 ; --l ) {
+ const float i_nu_temp = ( fact_m*i_nul ) + i_pnu_l;
+
+ fact_m -= xi;
+ i_pnu_l = ( fact_m*i_nu_temp ) + i_nul;
+ i_nul = i_nu_temp;
+ }
+ const float f = i_pnu_l/i_nul;
+ float Kmu, k_nu1;
+
+ if ( x < x_min ) {
+ const float x2 = 0.5F*x;
+ const float pi_mu = QFFM_PI*mu;
+ const float fact = ( qFFMath_Abs( pi_mu ) < eps ) ? 1.0F : pi_mu/qFFMath_Sin( pi_mu );
+ d = -qFFMath_Log( x2 );
+ float e = mu*d;
+ const float fact2 = ( qFFMath_Abs( e ) < eps ) ? 1.0F : qFFMath_Sinh( e )/e ;
+ float gam1, gam2, gam_pl, gam_mi;
+ gamma_temme(mu, &gam1, &gam2, &gam_pl, &gam_mi);
+ float ff = fact*( ( gam1*qFFMath_Cosh( e ) ) + ( gam2*fact2*d ) );
+ float sum = ff;
+ e = qFFMath_Exp( e );
+ float p = e/( 2.0F*gam_pl );
+ float q = 1.0F/( 2.0F*e*gam_mi );
+ float sum1 = p;
+ c = 1.0F;
+ d = x2*x2;
+
+ for ( int i = 1; i <= max_iter; ++i ) {
+ const float j = (float)i;
+ ff = ( j*ff + p + q )/( j*j - mu2 );
+ c *= d/j;
+ p /= j - mu;
+ q /= j + mu;
+ const float del = c*ff;
+ sum += del;
+ sum1 += c*( p - ( j*ff ) );
+ if ( qFFMath_Abs( del ) < ( eps*qFFMath_Abs( sum ) ) ) {
+ break;
+ }
+ }
+
+ Kmu = sum;
+ k_nu1 = sum1*xi2;
+ }
+ else {
+ float del_h = d;
+ float q1 = 0.0F;
+ float q2 = 1.0F;
+ const float a1 = 0.25F - mu2;
+ float q = a1;
+ float a = -a1;
+ float s = 1.0F + ( q*del_h );
+
+ b = 2.0F*( 1.0F + x );
+ d = 1.0F/b;
+ h = d;
+ c = a1;
+ for ( int i = 2 ; i <= max_iter; ++i) {
+ a -= (float)( 2*( i - 1 ) );
+ c = -a*c/( (float)i );
+ const float q_new = ( q1 - ( b*q2 ) )/a;
+ q1 = q2;
+ q2 = q_new;
+ q += c*q_new;
+ b += 2.0F;
+ d = 1.0F/( b + ( a*d ) );
+ del_h = ( ( b*d ) - 1.0F )*del_h;
+ h += del_h;
+ const float del_s = q*del_h;
+ s += del_s;
+ if ( qFFMath_Abs( del_s/s ) < eps ) {
+ break;
+ }
+ }
+ h = a1*h;
+ Kmu = qFFMath_Sqrt( QFFM_PI/( 2.0F*x ) )*qFFMath_Exp(-x)/s;
+ k_nu1 = Kmu*( mu + x + 0.5F - h )*xi;
+ }
+ const float k_pmu = ( mu*xi*Kmu ) - k_nu1;
+ const float i_num_u = xi/( ( f*Kmu ) - k_pmu );
+
+ i_nu[ 0 ] = i_num_u*i_nul1/i_nul;
+ i_pnu[ 0 ] = i_num_u*i_pnu_1/i_nul;
+ for ( int i = 1 ; i <= nl ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float k_nu_temp = ( ( mu + (float)i )*xi2*k_nu1 ) + Kmu;
+ /*cstat +CERT-FLP36-C*/
+ Kmu = k_nu1;
+ k_nu1 = k_nu_temp;
+ }
+ k_nu[ 0 ] = Kmu;
+ k_pnu[ 0 ] = ( nu*xi*Kmu ) - k_nu1;
+ }
+}
+/*============================================================================*/
+static float cyl_bessel_ij_series( float nu,
+ float x,
+ float sgn,
+ size_t max_iter )
+{
+ float y;
+
+ if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ y = ( qFFMath_IsEqual( 0.0F, nu ) ) ? 1.0F : 0.0F;
+ }
+ else {
+ const float x2 = 0.5F*x;
+ float fact = nu*qFFMath_Log( x2 );
+ float Jn = 1.0F;
+ float term = 1.0F;
+ const float xx4 = sgn*x2*x2;
+
+ fact -= qFFMath_LGamma( nu + 1.0F );
+ fact = qFFMath_Exp( fact );
+ for ( size_t i = 1U ; i < max_iter; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float j = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ term *= xx4/( j*( nu + j ) );
+ Jn += term;
+ if ( qFFMath_Abs( term/Jn ) < FLT_EPSILON ) {
+ break;
+ }
+ }
+ y = fact*Jn;
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Cyl_bessel_i( float nu,
+ float x )
+{
+ float y;
+
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( ( x*x ) < ( 10.0F*( nu + 1.0F ) ) ) {
+ y = cyl_bessel_ij_series( nu, x, 1.0F, 200U );
+ }
+ else {
+ float I_nu, K_nu, Ip_nu, Kp_nu;
+
+ bessel_ik( nu, x, &I_nu, &K_nu, &Ip_nu, &Kp_nu );
+ y = I_nu;
+ }
+ return y;
+}
+/*============================================================================*/
+static void cyl_bessel_jn_asymp( float nu,
+ float x,
+ float * Jnu,
+ float * Nnu)
+{
+ const float mu = 4.0F*nu*nu;
+ const float x8 = 8.0F*x;
+ float P = 0.0F;
+ float Q = 0.0F;
+ float term = 1.0F;
+ const float eps = FLT_EPSILON;
+ size_t i = 0U;
+
+ do {
+ bool epsP, epsQ;
+ float k2_1;
+ /*cstat -MISRAC++2008-0-1-2_b*/
+ float k = (float)i;
+ /*cstat +MISRAC++2008-0-1-2_b*/
+ k2_1 = 2.0F*k - 1.0F;
+ term *= ( i == 0U ) ? 1.0F : -( mu - ( k2_1*k2_1 ) )/( k*x8 );
+ epsP = qFFMath_Abs( term ) < ( eps*qFFMath_Abs( P ) );
+ P += term;
+ ++i;
+ k = (float)i;
+ k2_1 = 2.0F*k - 1.0F;
+ term *= ( mu - ( k2_1*k2_1 ) )/( k*x8 );
+ epsQ = qFFMath_Abs( term ) < ( eps*qFFMath_Abs( Q ) );
+ Q += term;
+ if ( epsP && epsQ && ( k > ( 0.5F*nu ) ) ) {
+ break;
+ }
+ ++i;
+ } while ( i < 1000U );
+ const float chi = x - ( nu + 0.5F )*QFFM_PI_2;
+ const float c = qFFMath_Cos( chi );
+ const float s = qFFMath_Sin( chi );
+ const float coeff = qFFMath_Sqrt( 2.0F/( QFFM_PI*x ) );
+ Jnu[ 0 ] = coeff*( ( c*P ) - ( s*Q ) );
+ Nnu[ 0 ] = coeff*( ( s*P ) + ( c*Q ) );;
+}
+/*============================================================================*/
+float qFFMath_Cyl_bessel_j( float nu,
+ float x )
+{
+ float y;
+
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( ( x*x ) < ( 10.0F*( nu + 1.0F ) ) ) {
+ y = cyl_bessel_ij_series( nu, x, -1.0F, 200U );
+ }
+ else if ( x > 1000.0F ) {
+ float j_nu, n_nu;
+
+ cyl_bessel_jn_asymp( nu, x, &j_nu, &n_nu );
+ y = j_nu;
+ }
+ else {
+ float J_nu, N_nu, Jp_nu, Np_nu;
+
+ bessel_jn( nu, x, &J_nu, &N_nu, &Jp_nu, &Np_nu );
+ y = J_nu;
+ }
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Cyl_bessel_k( float nu,
+ float x )
+{
+ float y;
+
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ float I_nu, K_nu, Ip_nu, Kp_nu;
+
+ bessel_ik( nu, x, &I_nu, &K_nu, &Ip_nu, &Kp_nu );
+ y = K_nu;
+ }
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Cyl_neumann( float nu,
+ float x )
+{
+ float y;
+
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( x > 1000.0F ) {
+ float J_nu, N_nu;
+ cyl_bessel_jn_asymp( nu, x, &J_nu, &N_nu );
+ y = N_nu;
+ }
+ else {
+ float J_nu, N_nu, Jp_nu, Np_nu;
+
+ bessel_jn( nu, x, &J_nu, &N_nu, &Jp_nu, &Np_nu );
+ y = N_nu;
+ }
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Sph_legendre( size_t l,
+ size_t m,
+ float theta )
+{
+ float y;
+
+ if ( qFFMath_IsNaN( theta ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ const float x = qFFMath_Cos( theta );
+ const float pi4 = 4.0F*QFFM_PI;
+
+ if ( m > l ) { //skipcq : CXX-W2041
+ y = 0.0F;
+ }
+ else if ( 0U == m ) {
+ float P = qFFMath_Legendre( l, x );
+ /*cstat -CERT-FLP36-C*/
+ const float fact = qFFMath_Sqrt( (float)( 2U*l + 1U )/pi4 );
+ /*cstat +CERT-FLP36-C*/
+ P *= fact;
+ y = P;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, x ) || qFFMath_IsEqual( -1.0F, x ) ) {
+ y = 0.0F;
+ }
+ else {
+ /*cstat -CERT-FLP36-C*/
+ const float mf = (float)m;
+ const float y_mp1m_factor = x*qFFMath_Sqrt( (float)( 2U*m + 3U ) );
+ /*cstat +CERT-FLP36-C*/
+ const float sgn = ( 1U == ( m % 2U ) ) ? -1.0F : 1.0F;
+ const float ln_circ = qFFMath_Log( 1.0F - ( x*x ) );
+ const float ln_poc_h = qFFMath_LGamma( mf + 0.5F ) - qFFMath_LGamma( mf );
+ const float ln_pre_val = ( -0.25F*QFFM_LN_PI ) + 0.5F*( ln_poc_h + ( mf*ln_circ ) );
+ const float sr = qFFMath_Sqrt( ( 2.0F + ( 1.0F/mf ) )/pi4);
+ float y_mm = sgn*sr*qFFMath_Exp( ln_pre_val );
+ float y_mp1m = y_mp1m_factor*y_mm;
+
+ if ( l == m ) {
+ y = y_mm;
+ }
+ else if ( l == ( m + 1U ) ) {
+ y = y_mp1m;
+ }
+ else {
+ float y_lm = 0.0F;
+
+ for ( size_t ll = ( m + 2U ) ; ll <= l; ++ll ) {
+ /*cstat -CERT-FLP36-C*/
+ const float ll_m_m = (float)( ll - m );
+ const float ll_p_m = (float)( ll + m );
+ const float ll2_p_1 = (float)( ( 2U*ll ) + 1U );
+ const float ll2_m_1 = (float)( ( 2U*ll ) - 1U );
+ const float ll_pm_m1 = (float)( ll + m - 1U );
+ const float ll_mm_m1 = (float)( ll - m - 1U );
+ /*cstat +CERT-FLP36-C*/
+ const float rat1 = ll_m_m/ll_p_m;
+ const float fact1 = qFFMath_Sqrt( rat1*ll2_p_1*ll2_m_1 );
+ /*cstat -CERT-FLP36-C*/
+ const float fact2 = qFFMath_Sqrt( rat1*( ll_mm_m1/ll_pm_m1 )*ll2_p_1/(float)( 2U*ll - 3U ) );
+ /*cstat -CERT-FLP36-C*/
+ y_lm = ( x*y_mp1m*fact1 - ll_pm_m1*y_mm*fact2 )/ll_m_m;
+ y_mm = y_mp1m;
+ y_mp1m = y_lm;
+ }
+ y = y_lm;
+ }
+ }
+ }
+ return y;
+}
#endif /*#ifndef QLIBS_USE_STD_MATH*/
\ No newline at end of file