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/include/qfis.h b/include/qfis.h
index 681cd3f..e914a94 100644
--- a/include/qfis.h
+++ b/include/qfis.h
@@ -114,6 +114,8 @@ extern "C" {
} qFIS_DeFuzzState_t;
/*! @endcond */
+ /*cstat -MISRAC2012-Dir-4.8*/
+
/*! @cond */
typedef struct
{
@@ -148,6 +150,8 @@ extern "C" {
/*! @endcond */
} qFIS_Output_t;
+ /*cstat +MISRAC2012-Dir-4.8*/
+
/*! @cond */
typedef float (*qFIS_MF_Fcn_t)( const qFIS_IO_Base_t * const in, const float *p, const size_t n );
/*! @endcond */
diff --git a/include/qinterp1.h b/include/qinterp1.h
new file mode 100644
index 0000000..ed30e5e
--- /dev/null
+++ b/include/qinterp1.h
@@ -0,0 +1,105 @@
+/*!
+ * @file qinterp1.hp
+ * @author J. Camilo Gomez C.
+ * @version 1.01
+ * @note This file is part of the qLibs distribution.
+ * @brief Class for a set of one-dimensional interpolators
+ **/
+
+#ifndef QCRC_H
+#define QCRC_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+ #include
+ #include
+
+ /** @addtogroup qinterp1 1D Interpolation
+ * @brief One-dimensional interpolation class.
+ * @{
+ */
+
+ /**
+ * @brief An enum with all the available interpolation methods.
+ */
+ typedef enum {
+ QINTERP1_NEXT = 0, /*!< Return the next neighbor.*/
+ QINTERP1_PREVIOUS, /*!< Return the previous neighbor.*/
+ QINTERP1_NEAREST, /*!< Return the nearest neighbor.*/
+ QINTERP1_LINEAR, /*!< Linear interpolation from nearest neighbors.*/
+ QINTERP1_SINE, /*!< Sine interpolation.*/
+ QINTERP1_CUBIC, /*!< Cubic interpolation.*/
+ QINTERP1_HERMITE, /*!< Piecewise cubic Hermite interpolation.*/
+ QINTERP1_SPLINE, /*!< Catmull spline interpolation.*/
+ QINTERP1_CONSTRAINED_SPLINE, /*!< A special kind of spline that doesn't overshoot.*/
+ QINTERP1_MAX,
+ } qInterp1Method_t;
+
+
+
+ typedef struct {
+ float (*method)( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+ const float *xData;
+ const float *yData;
+ size_t dataSize;
+ } qInterp1_t;
+
+ /**
+ * @brief Setup and initialize the 1D interpolation instance.
+ * @param[in] i A pointer to the interpolation instance.
+ * @param[in] xTable An array of size @a sizeTable with the x points sorted in ascending order.
+ * @param[in] yTable An array of size @a sizeTable with the y points.
+ * @param[in] sizeTable The number of points in @a xTable @a yTable
+ * @return 1 on success, otherwise return 0.
+ */
+ int qInterp1_Setup( qInterp1_t * const i,
+ const float * const xTable,
+ const float * const yTable,
+ const size_t sizeTable );
+
+ /**
+ * @brief Set the data table for the 1D interpolation instance.
+ * @param[in] i A pointer to the interpolation instance.
+ * @param[in] xTable An array of size @a sizeTable with the x points sorted in ascending order.
+ * @param[in] yTable An array of size @a sizeTable with the y points.
+ * @param[in] sizeTable The number of points in @a xTable @a yTable
+ * @return 1 on success, otherwise return 0.
+ */
+ int qInterp1_SetData( qInterp1_t * const i,
+ const float * const xTable,
+ const float * const yTable,
+ const size_t sizeTable );
+
+ /**
+ * @brief Specify the interpolation method to use.
+ * @param[in] i A pointer to the interpolation instance.
+ * @param[in] m The interpolation method.
+ * @return 1 on success, otherwise return 0.
+ */
+ int qInterp1_SetMethod( qInterp1_t * const i,
+ const qInterp1Method_t m );
+
+
+ /**
+ * @brief Interpolate input point @a x to determine the value of @a y
+ * at the points @a xi using the current method. If value is beyond
+ * the endpoints, extrapolation is performed using the current method.
+ * @param[in] i A pointer to the interpolation instance.
+ * @param[in] x The input point.
+ * @return @c The interpolated-extrapolated @a y value.
+ */
+ float qInterp1_Get( qInterp1_t * const i,
+ const float x );
+
+ /** @}*/
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
\ No newline at end of file
diff --git a/qffmath.c b/qffmath.c
index e675f6b..fa32bf1 100644
--- a/qffmath.c
+++ b/qffmath.c
@@ -18,6 +18,66 @@
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_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 )
{
@@ -202,26 +262,142 @@ float qFFMath_RCbrt( float x )
/*============================================================================*/
float qFFMath_Round( float x )
{
- x += 12582912.0F;
- x -= 12582912.0F;
- return x;
+ int32_t i0 = 0, j0;
+ float ret = x;
+
+ cast_reinterpret( i0, x, int32_t );
+ /*cstat -MISRAC2012-Rule-10.1_R7 -MISRAC2012-Rule-10.3 */
+ /*cstat -MISRAC2012-Rule-10.1_R6 -ATH-shift-neg -CERT-INT34-C_c */
+ /*cstat -MISRAC2012-Rule-1.3_n -MISRAC2012-Rule-10.4_a*/
+ j0 = ( ( i0 >> 23 ) & 0xFF ) - 0x7F;
+ if ( j0 < 23 ) {
+ if ( j0 < 0 ) {
+ i0 &= (int32_t)0x80000000U;
+ if ( -1 == j0 ) {
+ i0 |= 0x3F800000;
+ }
+ cast_reinterpret( ret, i0, float );
+ }
+ else {
+ const int32_t i = 0x007FFFFF >> j0;
+ if ( 0 != ( i0 & i ) ) {
+ i0 += 0x00400000 >> j0;
+ i0 &= ~i;
+ cast_reinterpret( ret, i0, float );
+ }
+ }
+ }
+ /*cstat +MISRAC2012-Rule-10.1_R7 +MISRAC2012-Rule-10.3 */
+ /*cstat +MISRAC2012-Rule-10.1_R6 +ATH-shift-neg +CERT-INT34-C_c */
+ /*cstat +MISRAC2012-Rule-1.3_n +MISRAC2012-Rule-10.4_a */
+ return ret;
}
/*============================================================================*/
float qFFMath_Floor( float x )
{
- return qFFMath_Round( x - 0.5F );
+ int32_t i0 = 0, j0;
+ float ret = x;
+
+ cast_reinterpret( i0, x, int32_t );
+ /*cstat -MISRAC2012-Rule-10.1_R7 -MISRAC2012-Rule-10.3 */
+ /*cstat -MISRAC2012-Rule-10.1_R6 -ATH-shift-neg -CERT-INT34-C_c */
+ /*cstat -MISRAC2012-Rule-1.3_n -MISRAC2012-Rule-10.4_a*/
+ j0 = ( ( i0 >> 23 ) & 0xFF ) - 0x7F;
+ if ( j0 < 23 ) {
+ if ( j0 < 0 ) {
+ /*cstat -ATH-neg-check-nonneg*/
+ if ( i0 >= 0 ) {
+ i0 = 0;
+ }
+ else if ( 0 != ( i0 & 0x7FFFFFFF ) ) {
+ i0 = (int32_t)0xBF800000U;
+ }
+ else {
+ /*nothing to do here*/
+ }
+ /*cstat +ATH-neg-check-nonneg*/
+ cast_reinterpret( ret, i0, float);
+ }
+ else {
+ const int32_t i = ( 0x007FFFFF ) >> j0;
+ if ( 0 != ( i0 & i ) ) {
+ /*cstat -ATH-neg-check-nonneg*/
+ if ( i0 < 0 ) {
+ i0 += 0x00800000 >> j0;
+ }
+ /*cstat +ATH-neg-check-nonneg*/
+ i0 &= (~i);
+ cast_reinterpret( ret, i0, float);
+ }
+ }
+ }
+ /*cstat +MISRAC2012-Rule-10.1_R7 +MISRAC2012-Rule-10.3 */
+ /*cstat +MISRAC2012-Rule-10.1_R6 +ATH-shift-neg +CERT-INT34-C_c */
+ /*cstat +MISRAC2012-Rule-1.3_n +MISRAC2012-Rule-10.4_a */
+ return ret;
}
/*============================================================================*/
float qFFMath_Ceil( float x )
{
- return qFFMath_Round( x + 0.5F );
+ int32_t i0 = 0, j0;
+ float ret = x;
+
+ cast_reinterpret( i0, x, int32_t );
+ /*cstat -MISRAC2012-Rule-10.1_R7 -MISRAC2012-Rule-10.3 */
+ /*cstat -MISRAC2012-Rule-10.1_R6 -ATH-shift-neg -CERT-INT34-C_c */
+ /*cstat -MISRAC2012-Rule-1.3_n -MISRAC2012-Rule-10.4_a*/
+ j0 = ( ( i0 >> 23 ) & 0xFF ) - 0x7F;
+ if ( j0 < 23 ) {
+ if ( j0 < 0 ) {
+ /*cstat -ATH-neg-check-nonneg*/
+ if ( i0 < 0 ) {
+ i0 = (int32_t)0x80000000U;
+ }
+ else if ( 0 != i0 ) {
+ i0 = (int32_t)0x3F800000U;
+ }
+ else {
+ /*nothing to do here*/
+ }
+ /*cstat +ATH-neg-check-nonneg*/
+ cast_reinterpret( ret, i0, float);
+ }
+ else {
+ const int32_t i = ( 0x007FFFFF ) >> j0;
+ if ( 0 != ( i0 & i ) ) {
+ if ( i0 > 0 ) {
+ i0 += 0x00800000 >> j0;
+ }
+ i0 &= (~i);
+ cast_reinterpret( ret, i0, float);
+ }
+ }
+ }
+ /*cstat +MISRAC2012-Rule-10.1_R7 +MISRAC2012-Rule-10.3 */
+ /*cstat +MISRAC2012-Rule-10.1_R6 +ATH-shift-neg +CERT-INT34-C_c */
+ /*cstat +MISRAC2012-Rule-1.3_n +MISRAC2012-Rule-10.4_a */
+ return ret;
}
/*============================================================================*/
float qFFMath_Trunc( float x )
{
- /*cstat -CERT-FLP36-C -CERT-FLP34-C*/
- return (float)( (int32_t)x );
- /*cstat +CERT-FLP36-C +CERT-FLP34-C*/
+ int32_t i0 = 0, sx, j0;
+ float ret = x;
+
+ cast_reinterpret( i0, x, int32_t );
+ /*cstat -MISRAC2012-Rule-10.1_R7 -MISRAC2012-Rule-10.3 */
+ /*cstat -MISRAC2012-Rule-10.1_R6 -ATH-shift-neg -CERT-INT34-C_c */
+ /*cstat -MISRAC2012-Rule-1.3_n -MISRAC2012-Rule-10.4_a*/
+ sx = i0 & (int32_t)0x80000000U;
+ j0 = ( ( i0 >> 23 ) & 0xFF ) - 0x7F;
+ if ( j0 < 23 ) {
+ const int32_t tmp = ( j0 < 0 ) ? sx : ( sx | ( i0 & ~( 0x007FFFFF >> j0 ) ) );
+ cast_reinterpret( ret, tmp, float );
+ }
+ /*cstat +MISRAC2012-Rule-10.1_R7 +MISRAC2012-Rule-10.3 */
+ /*cstat +MISRAC2012-Rule-10.1_R6 +ATH-shift-neg +CERT-INT34-C_c */
+ /*cstat +MISRAC2012-Rule-1.3_n +MISRAC2012-Rule-10.4_a */
+ return ret;
}
/*============================================================================*/
float qFFMath_Frac( float x )
@@ -238,8 +414,8 @@ float qFFMath_Remainder( float x,
float qFFMath_Mod( float x,
float y )
{
- return ( QFFM_FP_ZERO == qFFMath_FPClassify( x )) ? QFFM_NAN
- : ( x - ( y*qFFMath_Trunc( x/y ) ) );
+ return ( QFFM_FP_ZERO == qFFMath_FPClassify( x ) ) ? QFFM_NAN
+ : ( x - ( y*qFFMath_Trunc( x/y ) ) );
}
/*============================================================================*/
float qFFMath_Sin( float x )
@@ -456,6 +632,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 +1264,1711 @@ 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;
+ /*cstat -MISRAC2012-Rule-14.3_b*/
+ 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;
+ }
+ }
+ /*cstat +MISRAC2012-Rule-14.3_b*/
+ 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 sin_th2 = qFFMath_Sin( th2 );
+ const float serTerm2 = qFFMath_Sin( ( ( 0.25F*eta*( th2 - 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;
+ /*cstat -MISRAC2012-Rule-14.3_b*/
+ const float tc_sgn = ( x < 0.0F ) ? 1.0F :
+ ( ( 1U == ( n % 2U ) ) ? -1.0F : 1.0F );
+ /*cstat +MISRAC2012-Rule-14.3_b*/
+ 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;
+ const int tmp = N - i;
+ term *= ( b + k )/( (float)tmp )*( 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 a = (float)m;
+ const float N = (float)n;
+ /*cstat +CERT-FLP36-C -MISRAC2012-Rule-13.5*/
+ if ( ( x < 0.0F ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( 0U == n ) {
+ y = 1.0F;
+ }
+ else if ( 1U == n ) {
+ y = 1.0F + a - x;
+ }
+ else if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ float prod = a + 1.0F;
+ for ( size_t i = 2U ; i < n ; ++i ) {
+ /*cstat -CERT-FLP36-C*/
+ const float k = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ prod *= ( a + k )/k;
+ }
+ y = prod;
+ }
+ else if ( ( a > -1.0F ) && ( x < ( ( 2.0F*( a + 1.0F ) ) + ( 4.0F*N ) ) ) ) {
+ y = poly_laguerre_large_n( n, a, x );
+ }
+ else if ( ( ( x > 0.0F ) && ( a < -( N + 1.0F ) ) ) ) {
+ y = poly_laguerre_recursion( n, a, x );
+ }
+ else {
+ y = poly_laguerre_hyperg( n, a, x );
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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 = ( 1U == ( l % 2U ) ) ? -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*/
+ const float x_p_lm1 = x*p_lm1;
+ 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;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( qFFMath_IsNaN( x ) || qFFMath_IsNaN( y ) ) { //no side effects here
+ result = QFFM_NAN;
+ }
+ else {
+ const float bet = qFFMath_LGamma( x ) + qFFMath_LGamma( y ) -
+ qFFMath_LGamma( x + y );
+ result = qFFMath_Exp( bet );
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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 t = 1.0F - kk;
+ const float c13 = 1.0F/3.0F;
+
+ y = ellint_rf( 0.0F, t, 1.0F ) - ( c13*kk*ellint_rd( 0.0F, t, 1.0F ) );
+ }
+
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Comp_ellint_3( float k,
+ float nu )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( nu ) || ( qFFMath_Abs( k ) > 1.0F ) ) { //no side effects here
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 1.0F, nu ) ) {
+ y = QFFM_INFINITY;
+ }
+ else {
+ const float t = 1.0F - ( k*k );
+ const float c13 = 1.0F/3.0F;
+
+ y = ellint_rf( 0.0F, t, 1.0F ) + ( c13*nu*ellint_rj( 0.0F, t, 1.0F, 1.0F - nu ) );
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_1( float k,
+ float phi )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( phi ) || ( qFFMath_Abs( k ) > 1.0F ) ) { //no side effects here
+ 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 ) );
+ }
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_2( float k,
+ float phi )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( qFFMath_IsNaN( k ) || qFFMath_IsNaN( phi ) || ( qFFMath_Abs( k ) > 1.0F ) ) { //no side effects here
+ 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 ) );
+ }
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Ellint_3( float k,
+ float nu,
+ float phi )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ 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 ) );
+ }
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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_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 int tmp = i*( nm1 + i );
+ const float a = -( (float)tmp );
+ /*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 ) {
+ x = -x;
+ 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 );
+ }
+ }
+ 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 size_t tmp = i - 1U;
+ const float j = (float)( tmp );
+ /*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 = ( 1U == ( n % 2U ) ) ? -1.0F : 1.0F;
+ }
+ else {
+ float p_lm2 = 1.0F;
+
+ if ( 0U == n ) {
+ y = p_lm2;
+ }
+ else {
+ float p_lm1 = x;
+
+ if ( 1U == 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*/
+ const float x_plm1 = x*p_lm1;
+ y = ( 2.0F*x_plm1 ) - p_lm2 - ( ( x_plm1 - 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
+ };
+ float zeta = 1.0F;
+
+ for ( size_t i = 0U ; i < sizeof(prime) ; ++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 float pre_tmp = nu - x + 1.5F;
+ const int tmp_nl = (int)pre_tmp;
+ const float tmp_1 = nu + 0.5F;
+ const int nl = ( x < x_min ) ? (int)tmp_1
+ : ( ( 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 = 0.0F, gam2 = 0.0F, gam_pl = 0.0F, gam_mi = 0.0F;
+
+ 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 ) {
+ const int tmp = 2*( i - 1 );
+ /*cstat -CERT-FLP36-C*/
+ a += (float)tmp;
+ /*cstat +CERT-FLP36-C*/
+ 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 = 0.0F, n_nu = 0.0F, jp_nu = 0.0F, np_nu = 0.0F;
+ 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;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ 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 = 0.0F, n_n = 0.0F, jp_n = 0.0F, np_n = 0.0F;
+
+ sph_bessel_jn( n, x, &j_n, &n_n, &jp_n, &np_n );
+ y = j_n;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Sph_neumann( size_t n,
+ float x )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( ( x < 0.0F ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( qFFMath_IsEqual( 0.0F, x ) ) {
+ y = -QFFM_INFINITY;
+ }
+ else {
+ float j_n = 0.0F, n_n = 0.0F, jp_n = 0.0F, np_n = 0.0F;
+
+ sph_bessel_jn( n, x, &j_n, &n_n, &jp_n, &np_n );
+ y = n_n;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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 float nl_f = nu + 0.5F;
+ const int nl = (int)nl_f;
+ const float mu = nu - nl_f;
+ /*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 = 0.0F, gam2 = 0.0F, gam_pl = 0.0F, gam_mi = 0.0F;
+ 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 ) {
+ /*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;
+ 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) {
+ /*cstat -CERT-FLP36-C*/
+ const int tmp = 2*( i - 1 );
+ a -= (float)tmp;
+ c = -a*c/( (float)i );
+ /*cstat +CERT-FLP36-C*/
+ 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;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ 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 = 0.0F, K_nu = 0.0F, Ip_nu = 0.0F, Kp_nu = 0.0F;
+
+ bessel_ik( nu, x, &I_nu, &K_nu, &Ip_nu, &Kp_nu );
+ y = I_nu;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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 -CERT-FLP36-C*/
+ float k = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ 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;
+ /*cstat -CERT-FLP36-C*/
+ k = (float)i;
+ /*cstat +CERT-FLP36-C*/
+ 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;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ 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 = 0.0F, n_nu = 0.0F;
+
+ cyl_bessel_jn_asymp( nu, x, &j_nu, &n_nu );
+ y = j_nu;
+ }
+ else {
+ float J_nu = 0.0F, N_nu = 0.0F, Jp_nu = 0.0F, Np_nu = 0.0F;
+
+ bessel_jn( nu, x, &J_nu, &N_nu, &Jp_nu, &Np_nu );
+ y = J_nu;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Cyl_bessel_k( float nu,
+ float x )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else {
+ float I_nu = 0.0F, K_nu = 0.0F, Ip_nu = 0.0F, Kp_nu = 0.0F;
+
+ bessel_ik( nu, x, &I_nu, &K_nu, &Ip_nu, &Kp_nu );
+ y = K_nu;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ return y;
+}
+/*============================================================================*/
+float qFFMath_Cyl_neumann( float nu,
+ float x )
+{
+ float y;
+ /*cstat -MISRAC2012-Rule-13.5*/
+ if ( ( nu < 0.0F ) || ( x < 0.0F ) || qFFMath_IsNaN( nu ) || qFFMath_IsNaN( x ) ) {
+ y = QFFM_NAN;
+ }
+ else if ( x > 1000.0F ) {
+ float J_nu = 0.0F, N_nu = 0.0F;
+ cyl_bessel_jn_asymp( nu, x, &J_nu, &N_nu );
+ y = N_nu;
+ }
+ else {
+ float J_nu = 0.0F, N_nu = 0.0F, Jp_nu = 0.0F, Np_nu = 0.0F;
+
+ bessel_jn( nu, x, &J_nu, &N_nu, &Jp_nu, &Np_nu );
+ y = N_nu;
+ }
+ /*cstat +MISRAC2012-Rule-13.5*/
+ 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 size_t tmp = ( 2U*l ) + 1U;
+ const float fact = qFFMath_Sqrt( ( (float)tmp )/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 size_t tmp = ( 2U*m ) + 3U;
+ const float y_mp1m_factor = x*qFFMath_Sqrt( (float)tmp );
+ /*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 size_t u_ll_m_m = ll - m;
+ const size_t u_ll_p_m = ll + m;
+ const size_t u_ll2_p_1 = ( 2U*ll ) + 1U;
+ const size_t u_ll2_m_1 = ( 2U*ll ) - 1U;
+ const size_t u_ll_pm_m1 = ll + m - 1U;
+ const size_t u_ll_mm_m1 = ll - m - 1U;
+ const size_t u_ll2_m_3 = ( 2U*ll ) - 3U;
+
+ const float ll_m_m = (float)u_ll_m_m;
+ const float ll_p_m = (float)u_ll_p_m;
+ const float ll2_p_1 = (float)u_ll2_p_1;
+ const float ll2_m_1 = (float)u_ll2_m_1;
+ const float ll_pm_m1 = (float)u_ll_pm_m1;
+ const float ll_mm_m1 = (float)u_ll_mm_m1;
+ const float ll2_m_3 = (float)u_ll2_m_3;
+ /*cstat +CERT-FLP36-C*/
+ const float rat1 = ll_m_m/ll_p_m;
+ const float rl2p1 = rat1*ll2_p_1;
+ const float fact1 = qFFMath_Sqrt( rl2p1*ll2_m_1 );
+ /*cstat -CERT-FLP36-C*/
+ const float fact2 = qFFMath_Sqrt( rl2p1*( ll_mm_m1/ll_pm_m1 )/ll2_m_3 );
+ /*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
diff --git a/qfis.c b/qfis.c
index 75937c0..48d3196 100644
--- a/qfis.c
+++ b/qfis.c
@@ -439,15 +439,17 @@ int qFIS_Fuzzify( qFIS_t * const f )
static float qFIS_ParseFuzzValue( qFIS_MF_t * const mfIO,
qFIS_Rules_t index )
{
+ /*cstat -CERT-STR34-C*/
uint8_t neg = ( index < 0 ) ? 1U : 0U ;
+ /*cstat +CERT-STR34-C*/
float y;
if ( 0U != neg ) {
index = -index;
}
- /*cstat -CERT-INT32-C_a*/
+ /*cstat -CERT-INT32-C_a -CERT-STR34-C*/
y = qFIS_Bound( mfIO[ index - 1 ].fx, 0.0F, 1.0F );
- /*cstat +CERT-INT32-C_a*/
+ /*cstat +CERT-INT32-C_a +CERT-STR34-C*/
/*cppcheck-suppress misra-c2012-12.1 */
y = ( 0U != neg ) ? ( 1.0F - y ) : y ;
@@ -555,9 +557,10 @@ static size_t qFIS_InferenceConsequent( struct _qFIS_s * const f,
MFOutIndex -= 1;
if ( f->wi[ f->ruleCount ] > 0.0F ) {
+ /*cstat -CERT-STR34-C*/
qFIS_Output_t *o = &f->output[ outIndex ];
qFIS_MF_t *m = &f->outMF[ MFOutIndex ];
-
+ /*cstat +CERT-STR34-C*/
if ( Mamdani == f->type ) {
float v;
/*cstat -MISRAC2012-Rule-11.3 -CERT-EXP39-C_d*/
@@ -1039,7 +1042,7 @@ static float qFIS_SMF( const qFIS_IO_Base_t * const in,
const float *p,
const size_t n )
{
- float a, b, tmp, y;
+ float a, b, tmp, y = 0.0F;
float x = in[ 0 ].value;
(void)n;
@@ -1060,7 +1063,7 @@ static float qFIS_SMF( const qFIS_IO_Base_t * const in,
y = ( 1.0F - ( 2.0F*tmp*tmp ) );
}
else {
- y = 0.0F;
+ /* Nothing to do here */
}
return y;
@@ -1094,7 +1097,7 @@ static float qFIS_ZMF( const qFIS_IO_Base_t * const in,
const float *p,
const size_t n )
{
- float a, b, tmp, y;
+ float a, b, tmp, y = 0.0F;
float x = in[ 0 ].value;
(void)n;
@@ -1115,7 +1118,7 @@ static float qFIS_ZMF( const qFIS_IO_Base_t * const in,
y = 2.0F*tmp*tmp;
}
else {
- y = 0.0F;
+ /* Nothing to do here */
}
return y;
diff --git a/qfp16.c b/qfp16.c
index d86c759..38c96d1 100644
--- a/qfp16.c
+++ b/qfp16.c
@@ -374,8 +374,8 @@ qFP16_t qFP16_Sqrt( qFP16_t x )
if ( x > 0 ) {
uint32_t bit;
uint8_t n;
-
retValue = 0;
+ /*cstat -ATH-shift-bounds -MISRAC2012-Rule-12.2 -CERT-INT34-C_b*/
/*cppcheck-suppress [ cert-INT31-c, misra-c2012-12.2, misra-c2012-12.1 ]*/
bit = ( 0 != ( x & (qFP16_t)4293918720 ) ) ? ( 1U << 30U ) : ( 1U << 18U );
while ( bit > (uint32_t)x ) {
@@ -418,7 +418,7 @@ qFP16_t qFP16_Sqrt( qFP16_t x )
if ( ( 1U == fp->rounding ) && ( x > retValue ) ) {
++retValue;
}
-
+ /*cstat +ATH-shift-bounds +MISRAC2012-Rule-12.2 +CERT-INT34-C_b*/
return (qFP16_t)retValue;
}
/*============================================================================*/
diff --git a/qinterp1.c b/qinterp1.c
new file mode 100644
index 0000000..dc09c96
--- /dev/null
+++ b/qinterp1.c
@@ -0,0 +1,573 @@
+/*!
+ * @file qinterp1.c
+ * @author J. Camilo Gomez C.
+ * @note This file is part of the qLibs distribution.
+ **/
+
+#include "qinterp1.h"
+#include "qffmath.h"
+
+typedef float (*qInterp1Fcn_t)( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+
+static float qInterp1_next( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_previous( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_nearest( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_linear( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_sine( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_cubic( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_hermite( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_spline( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+static float qInterp1_cSpline( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize );
+
+static float slope( const float * const tx,
+ const float * const ty,
+ const size_t i );
+float firstDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i );
+float leftSecondDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i );
+float rightSecondDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i );
+
+
+/*cstat -CERT-INT30-C_a*/
+
+/*============================================================================*/
+int qInterp1_Setup( qInterp1_t * const i,
+ const float * const xTable,
+ const float * const yTable,
+ const size_t sizeTable )
+{
+ int retVal = 0;
+
+ if ( ( NULL != i ) && ( sizeTable > 2U ) && ( NULL != xTable ) && ( NULL != xTable ) ) {
+ i->xData = xTable;
+ i->yData = yTable;
+ i->dataSize = sizeTable;
+ i->method = &qInterp1_linear;
+ retVal = 1;
+ }
+
+ return retVal;
+}
+/*============================================================================*/
+int qInterp1_SetData( qInterp1_t * const i,
+ const float * const xTable,
+ const float * const yTable,
+ const size_t sizeTable )
+{
+ return qInterp1_Setup( i, xTable, yTable, sizeTable );
+}
+/*============================================================================*/
+int qInterp1_SetMethod( qInterp1_t * const i,
+ const qInterp1Method_t m )
+{
+ int retVal = 0;
+ static const qInterp1Fcn_t im[ QINTERP1_MAX ] = {
+ &qInterp1_next,
+ &qInterp1_previous,
+ &qInterp1_nearest,
+ &qInterp1_linear,
+ &qInterp1_sine,
+ &qInterp1_cubic,
+ &qInterp1_hermite,
+ &qInterp1_spline,
+ &qInterp1_cSpline,
+ };
+
+ if ( ( NULL != i ) && ( m < QINTERP1_MAX ) ) {
+ i->method = im[ m ];
+ retVal = 1;
+ }
+
+ return retVal;
+}
+/*============================================================================*/
+float qInterp1_Get( qInterp1_t * const i,
+ const float x )
+{
+ return i->method( x, i->xData, i->yData, i->dataSize );
+}
+/*============================================================================*/
+static float qInterp1_next( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ size_t nearestIndex = tableSize - 1U;
+
+ for ( size_t i = 0U ; i < ( tableSize - 1U ) ; ++i ) {
+ if ( x < tx[ i + 1U ] ) {
+ nearestIndex = i;
+ break;
+ }
+ }
+
+ if ( x >= tx[ tableSize - 1U ] ) {
+ y = ty[ tableSize - 1U ];
+ }
+ else {
+ y = ty[ nearestIndex + 1U ];
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_previous( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ if ( x <= tx[ 0 ] ) {
+ y = ty[ 0 ];
+ }
+ else {
+ size_t nearestIndex = 0U;
+
+ for ( size_t i = 1U ; i < tableSize; ++i ) {
+ if ( x < tx [ i ] ) {
+ break;
+ }
+ nearestIndex = i;
+ }
+ y = ty[ nearestIndex ];
+ }
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_nearest( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ size_t nearestIndex = 0U;
+ float minDistance = qFFMath_Abs( x - tx[ 0 ] );
+
+ for ( size_t i = 1U ; i < tableSize; ++i ) {
+ const float distance = qFFMath_Abs( x - tx[ i ] );
+
+ if ( distance <= minDistance) {
+ minDistance = distance;
+ nearestIndex = i;
+ }
+ }
+ y = ty[ nearestIndex ];
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_linear( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ if ( x < tx[ 0 ] ) {
+ const float x0 = tx[ 0 ];
+ const float x1 = tx[ 1 ];
+ const float y0 = ty[ 0 ];
+ const float y1 = ty[ 1 ];
+ y = y0 + ( ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ) );
+ }
+ else if ( x > tx[ tableSize - 1U ] ) {
+ const float x0 = tx[ tableSize - 2U ];
+ const float x1 = tx[ tableSize - 1U ];
+ const float y0 = ty[ tableSize - 2U ];
+ const float y1 = ty[ tableSize - 1U ];
+ y = y1 + ( ( ( y0 - y1 )/( x0 - x1 ) )*( x - x1 ) );
+ }
+ else {
+ const int maxIndex = (int)tableSize - 1;
+ for ( int i = 0; i < maxIndex; ++i ) {
+ if ( ( x >= tx[ i ] ) && ( x <= tx[ i + 1 ] ) ) {
+ const float x0 = tx[ i ];
+ const float x1 = tx[ i + 1 ];
+ const float y0 = ty[ i ];
+ const float y1 = ty[ i + 1 ];
+ y = y0 + ( ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ) );
+ break;
+ }
+ }
+ }
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_sine( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ if ( x < tx[ 0 ] ) {
+ const float x0 = tx[ 0 ];
+ const float x1 = tx[ 1 ];
+ const float y0 = ty[ 0 ];
+ const float y1 = ty[ 1 ];
+ const float w = 0.5F - ( 0.5F*qFFMath_Cos( QFFM_PI*( x - x0 )/( x1 - x0 ) ) );
+ y = y0 + ( w*( y1 - y0 ) );
+ }
+ else if ( x > tx[ tableSize - 1U ] ) {
+ const float x0 = tx[ tableSize - 2U ];
+ const float x1 = tx[ tableSize - 1U ];
+ const float y0 = ty[ tableSize - 2U ];
+ const float y1 = ty[ tableSize - 1U ];
+ const float w = 0.5F - ( 0.5F*qFFMath_Cos( QFFM_PI*( x - x1 )/( x0 - x1 ) ) );
+ y = y1 + ( w*( y0 - y1 ) );
+ }
+ else {
+ for ( size_t i = 1; i < tableSize; ++i ) {
+ if ( x <= tx[ i ] ) {
+ const float x0 = tx[ i - 1U ];
+ const float x1 = tx[ i ];
+ const float y0 = ty[ i - 1U ];
+ const float y1 = ty[ i];
+ const float w = 0.5F - ( 0.5F*qFFMath_Cos( QFFM_PI*( x - x0 )/( x1 - x0 ) ) );
+ y = y0 + ( w*( y1 - y0 ) );
+ }
+ }
+ }
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_cubic( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 4U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ if ( x < tx[ 0 ] ) {
+ const float x0 = tx[ 0 ];
+ const float x1 = tx[ 1 ];
+ const float y0 = ty[ 0 ];
+ const float y1 = ty[ 1 ];
+ const float h = x1 - x0;
+ const float t = ( x - x0 )/h;
+ const float t2 = t*t;
+ const float t3 = t2*t;
+
+ y = ( ( ( 2.0F*t3 ) - ( 3.0F*t2 ) + 1.0F )*y0 ) +
+ ( ( t3 - ( 2.0F*t2 ) + t )*h*( y0 - y1 ) ) +
+ ( ( -( 2.0F*t3 ) + ( 3.0F*t2 ) )*y1 ) +
+ ( ( t3 - t2 )*h*( y1 - y0 ) );
+ }
+ else if ( x > tx[ tableSize - 1U ] ) {
+ const float x0 = tx[ tableSize - 2U ];
+ const float x1 = tx[ tableSize - 1U ];
+ const float y0 = ty[ tableSize - 2U ];
+ const float y1 = ty[ tableSize - 1U ];
+ const float h = x1 - x0;
+ const float t = ( x - x1 )/h;
+ const float t2 = t*t;
+ const float t3 = t2*t;
+
+ y = ( ( ( 2.0F*t3 ) - ( 3.0F*t2 ) + 1.0F )*y1 ) +
+ ( ( t3 - ( 2.0F*t2 ) + t )*h*( y0 - ty[ tableSize - 3U ] ) ) +
+ ( ( -( 2.0F*t3 ) + ( 3.0F*t2 ) )*y0 ) +
+ ( ( t3 - t2 )*h*( y1 - y0 ) );
+ }
+ else {
+ for ( size_t i = 1U; i < tableSize; ++i ) {
+ if ( x <= tx[ i ] ) {
+ const float x0 = tx[ i - 1U ];
+ const float x1 = tx[ i ];
+ const float y0 = ty[ i - 1U ];
+ const float y1 = ty[ i ];
+ const float h = x1 - x0;
+ const float t = ( x - x0 )/h;
+ const float t2 = t*t;
+ const float t3 = t2*t;
+ y = ( ( ( 2.0F*t3 ) - ( 3.0F*t2 ) + 1.0F )*y0 ) +
+ ( ( t3 - ( 2.0F*t2 ) + t )*h*( y0 - ty[ i - 2U ] ) ) +
+ ( ( -( 2.0F*t3 ) + ( 3.0F*t2 ) )*y1 ) +
+ ( ( t3 - t2 )*h*( y1 - y0 ) );
+ break;
+ }
+ }
+ }
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_hermite( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 2U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ if ( x < tx[ 0 ] ) {
+ const float x0 = tx[ 0 ];
+ const float x1 = tx[ 1 ];
+ const float y0 = ty[ 0 ];
+ const float y1 = ty[ 1 ];
+ y = y0 + ( ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ) );
+ }
+ else if ( x > tx[ tableSize -1U ] ) {
+ const float x0 = tx[ tableSize - 2U ];
+ const float x1 = tx[ tableSize - 1U ];
+ const float y0 = ty[ tableSize - 2U ];
+ const float y1 = ty[ tableSize - 1U ];
+ y = y1 + ( ( ( y0 - y1 )/( x0 - x1 ) )*( x - x1 ) );
+ }
+ else {
+ y = 0.0F;
+ for ( size_t i = 0U ; i < tableSize; ++i ) {
+ float term = ty[ i ];
+
+ for ( size_t j = 0U ; j < tableSize; ++j ) {
+ if ( i != j ) {
+ term *= ( x - tx[ j ] )/( tx[ i ] - tx[ j ] );
+ }
+ }
+
+ y += term;
+ }
+ }
+ }
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_spline( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 4U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ size_t i = 0U;
+ /* Extrapolation for x beyond the range*/
+ if ( x <= tx[ 0 ] ) {
+ i = 0U;
+ }
+ else if ( x >= tx[ tableSize - 1U ] ) {
+ i = tableSize - 2U; /* Use the last interval for extrapolation*/
+ }
+ else {
+ while ( x >= tx[ i + 1U ] ) {
+ i++;
+ }
+ }
+
+ if ( qFFMath_IsEqual( x , tx[ i + 1U ] ) ) {
+ y = ty[ i + 1U ];
+ }
+ else {
+ const float t = ( x - tx[ i ] )/( tx[ i + 1U ] - tx[ i ] );
+ const float t_2 = t*t;
+ const float t_3 = t_2*t;
+ const float tt_3 = 2.0F*t_3;
+ const float tt_2 = 3.0F*t_2;
+ const float h01 = tt_2 - tt_3;
+ const float h00 = 1.0F - h01;
+ const float h10 = t_3 - ( 2.0F*t_2 ) + t;
+ const float h11 = t_3 - t_2;
+ const float x1_x0 = tx[ i + 1U ] - tx[ i ];
+ const float y0 = ty[ i ];
+ const float y1 = ty[ i + 1U ];
+ float m0;
+ float m1;
+
+ if ( 0U == i ) {
+ m0 = ( ty[ 1 ] - ty[ 0 ] )/( tx[ 1 ] - tx[ 0 ] );
+ m1 = ( ty[ 2 ] - ty[ 0 ] )/( tx[ 2 ] - tx[ 0 ] );
+ }
+ else if ( ( tableSize - 2U ) == i ) {
+ m0 = ( ty[ tableSize - 1U ] - ty[ tableSize - 3U ] )/( tx[ tableSize - 1U ] - tx[ tableSize - 3U ] );
+ m1 = ( ty[ tableSize - 1U ] - ty[ tableSize - 2U ] )/( tx[ tableSize - 1U ] - tx[ tableSize - 2U ] );
+ }
+ else {
+ m0 = slope( tx, ty, i );
+ m1 = slope( tx, ty, i + 1U );
+ }
+ y = ( h00*y0 ) + ( h01*y1 ) + ( h10*x1_x0*m0 ) + ( h11*x1_x0*m1 );
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float qInterp1_cSpline( const float x,
+ const float * const tx,
+ const float * const ty,
+ const size_t tableSize )
+{
+ float y = QFFM_NAN;
+
+ if ( ( tableSize >= 4U ) && ( NULL != tx ) && ( NULL != ty ) ) {
+ size_t i = 0U;
+ /* Extrapolation for x beyond the range*/
+ if ( x <= tx[ 0 ] ) {
+ i = 0U;
+ }
+ else if ( x >= tx[ tableSize - 1U ] ) {
+ i = tableSize - 2U; /* Use the last interval for extrapolation*/
+ }
+ else {
+ while ( x >= tx[ i + 1U ] ) {
+ i++;
+ }
+ }
+
+ if ( qFFMath_IsEqual( x , tx[ i + 1U ] ) ) {
+ y = ty[ i + 1U ];
+ }
+ else {
+ const float x0 = tx[ i ];
+ const float x1 = tx[ i + 1U ];
+ const float y0 = ty[ i ];
+ const float y1 = ty[ i + 1U ];
+ const float fd2i_xl1 = leftSecondDerivate( tx, ty, tableSize - 1U, i + 1U );
+ const float fd2i_x = rightSecondDerivate( tx, ty, tableSize - 1U, i + 1U );
+ const float x0_x1 = x0 - x1;
+ const float x1_2 = x1*x1;
+ const float x1_3 = x1_2*x1;
+ const float x0_2 = x0*x0;
+ const float x0_3 = x0_2*x0;
+ const float inv_x0_x1 = 1.0F/x0_x1;
+ const float d = ( fd2i_x - fd2i_xl1 )*( 0.166666667F*inv_x0_x1 );
+ const float c = ( ( x0*fd2i_xl1 ) - ( x1*fd2i_x ) )*( 0.5F*inv_x0_x1 );
+ const float b = ( y0 - y1 - ( c*( x0_2 - x1_2 ) ) - ( d*( x0_3 - x1_3 ) ) )*inv_x0_x1;
+ const float a = y1 - ( b*x1 ) - ( c*x1_2 ) - ( d*x1_3 );
+ y = a + ( x*( b + ( x*( c + ( x*d ) ) ) ) );
+ }
+ }
+
+ return y;
+}
+/*============================================================================*/
+static float slope( const float * const tx,
+ const float * const ty,
+ const size_t i )
+{
+ float m;
+
+ if ( qFFMath_IsEqual( tx[ i + 1U ], tx[ i - 1U ] ) ) {
+ m = 0.0F;
+ }
+ else {
+ m = ( ty[ i + 1U ] - ty[ i - 1U ] )/( tx[ i + 1U ] - tx[ i - 1U ] );
+ }
+
+ return m;
+}
+/*============================================================================*/
+float firstDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i )
+{
+ float fd1_x;
+
+ if ( 0U == i ) {
+ const float dx = tx[ 1 ] - tx[ 0 ];
+ const float dy = ty[ 1 ] - ty[ 0 ];
+ fd1_x = 1.5F*( dy/dx );
+ fd1_x -= 1.0F/( ( ( tx[ 2 ] - tx[ 0 ] )/( ty[ 2 ] - ty[ 0 ] ) ) + ( dx/dy ) );
+ }
+ else if ( n == i ) {
+ const float dx = tx[ n ] - tx[ n - 1U ];
+ const float dy = ty[ n ] - ty[ n - 1U ];
+ fd1_x = 1.5F*( dy/dx );
+ fd1_x -= 1.0F/( ( ( tx[ n ] - tx[ n - 2U ] )/( ty[ n ] - ty[ n - 2U ] ) ) + ( dx/dy ) );
+ }
+ else {
+ const float tmp1 = ( tx[ i + 1U ] - tx[ i ] )/( ty[ i + 1U ] - ty[ i ] );
+ const float tmp2 = ( tx[ i ] - tx[ i - 1U ] )/( ty[ i ] - ty[ i - 1U ] );
+
+ if ( ( tmp1*tmp2 ) < 0.0F ) {
+ fd1_x = 0.0F;
+ }
+ else {
+ fd1_x = 2.0F/( tmp1 + tmp2 );
+ }
+ }
+ return fd1_x;
+}
+/*============================================================================*/
+float leftSecondDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i )
+{
+ const float fdi_x = firstDerivate( tx, ty, n, i );
+ const float fdi_xl1 = firstDerivate( tx, ty, n, i - 1U );
+ const float xi_delta = tx[ i ] - tx[ i - 1U ];
+ float fd2l_x = -2.0F*( fdi_x + ( 2.0F*fdi_xl1 ) )/xi_delta;
+
+ fd2l_x += 6.0F*( ty[ i ] - ty[ i - 1U ] )/( xi_delta*xi_delta );
+ return fd2l_x;
+}
+/*============================================================================*/
+float rightSecondDerivate( const float * const tx,
+ const float * const ty,
+ const size_t n,
+ const size_t i )
+{
+ const float fdi_x = firstDerivate( tx, ty, n, i );
+ const float fdi_xl1 = firstDerivate( tx, ty, n, i - 1U );
+ const float xi_delta = tx[ i ] - tx[ i - 1U ];
+ float fd2r_x = 2.0F*( ( 2.0F*fdi_x ) + fdi_xl1 )/xi_delta;
+
+ fd2r_x -= 6.0F*( ty[ i ] - ty[ i - 1U ] )/( xi_delta*xi_delta );
+ return fd2r_x;
+}
+/*============================================================================*/
+
+/*cstat +CERT-INT30-C_a*/
\ No newline at end of file