diff --git a/lib/node_modules/@stdlib/stats/base/dists/truncated-normal/pdf/lib/truncated.c b/lib/node_modules/@stdlib/stats/base/dists/truncated-normal/pdf/lib/truncated.c new file mode 100644 index 000000000000..56f851d971c1 --- /dev/null +++ b/lib/node_modules/@stdlib/stats/base/dists/truncated-normal/pdf/lib/truncated.c @@ -0,0 +1,1455 @@ +#include "stdlib/math/base/assert/is_nan.h" +#include +#include "stdlib/math/base/special/exp.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/special/trunc.h" +#include "stdlib/constants/float64/ninf.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/math/base/special/ldexp.h" +#include + +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float64/high_word_abs_mask.h" +#include "stdlib/number/float64/base/get_high_word.h" +#include "stdlib/math/base/assert/is_odd.h" +#include "stdlib/math/base/special/copysign.h" +#include "stdlib/constants/float64/ninf.h" +#include "stdlib/number/float64/base/set_high_word.h" +#include "stdlib/number/float64/base/set_low_word.h" +#include "stdlib/math/base/special/ldexp.h" +#include "stdlib/constants/float64/ln_two.h" +#include "stdlib/constants/float64/exponent_bias.h" +#include "stdlib/constants/float64/high_word_significand_mask.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/assert/is_infinite.h" +#include "stdlib/math/base/assert/is_integer.h" +#include "stdlib/math/base/special/sqrt.h" +#include "stdlib/number/float64/base/to_words.h" +#include "stdlib/constants/float64/num_high_word_significand_bits.h" + +#include "stdlib/math/base/special/sqrt.h" +#include + +#include "stdlib/math/base/special/erfc.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/special/exp.h" +#include "stdlib/number/float64/base/set_low_word.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float64/ninf.h" +#include + +/** +* Tests if a double-precision floating-point numeric value is NaN. +* +* @param x number +* @return boolean indicating if an input value is NaN +*/ +bool stdlib_base_is_nan( const double x ) { + return ( x != x ); +} + +static const double LN2_HI = 6.93147180369123816490e-01; +static const double LN2_LO = 1.90821492927058770002e-10; +static const double LOG2_E = 1.44269504088896338700e+00; +static const double EXP_OVERFLOW = 7.09782712893383973096e+02; +static const double EXP_UNDERFLOW = -7.45133219101941108420e+02; +static const double NEARZERO = 1.0 / (1 << 28); // 2^-28 +static const double NEG_NEARZERO = -NEARZERO; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_p + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_p( const double x ) { + return 0.16666666666666602 + (x * (-0.0027777777777015593 + (x * (0.00006613756321437934 + (x * (-0.0000016533902205465252 + (x * 4.1381367970572385e-8))))))); +} + +// END: polyval_p + +/* End auto-generated functions. */ + +/** +* Computes \\(e^{r} 2^k\\) where \\(r = \mathrm{hi} - \mathrm{lo}\\) and \\(|r| \leq \ln(2)/2\\). +* +* @param hi upper bound +* @param lo lower bound +* @param k power of 2 +* @return function value +*/ +static double expmulti( const double hi, const double lo, const int32_t k ) { + double r; + double t; + double c; + double y; + + r = hi - lo; + t = r * r; + c = r - (t * polyval_p( t )); + y = 1.0 - ( lo - ( ( r * c ) / ( 2.0 - c ) ) - hi); + + return stdlib_base_ldexp( y, k ); +} + +/** +* Evaluates the natural exponential function. +* +* ## Method +* +* 1. We reduce \\( x \\) to an \\( r \\) so that \\( |r| \leq 0.5 \cdot \ln(2) \approx 0.34658 \\). Given \\( x \\), we find an \\( r \\) and integer \\( k \\) such that +* +* tex +* \begin{align*} +* x &= k \cdot \ln(2) + r \\ +* |r| &\leq 0.5 \cdot \ln(2) +* \end{align*} +* +* +* +* +* \\( r \\) can be represented as \\( r = \mathrm{hi} - \mathrm{lo} \\) for better accuracy. +* +* +* +* 2. We approximate of \\( e^{r} \\) by a special rational function on the interval \\(\[0,0.34658]\\): +* +* tex +* \begin{align*} +* R\left(r^2\right) &= r \cdot \frac{ e^{r}+1 }{ e^{r}-1 } \\ +* &= 2 + \frac{r^2}{6} - \frac{r^4}{360} + \ldots +* \end{align*} +* +* +* We use a special Remes algorithm on \\(\[0,0.34658]\\) to generate a polynomial of degree \\(5\\) to approximate \\(R\\). The maximum error of this polynomial approximation is bounded by \\(2^{-59}\\). In other words, +* +* tex +* R(z) \sim 2 + P_1 z + P_2 z^2 + P_3 z^3 + P_4 z^4 + P_5 z^5 +* +* +* where \\( z = r^2 \\) and +* +* tex +* \left| 2 + P_1 z + \ldots + P_5 z^5 - R(z) \right| \leq 2^{-59} +* +* +* +* +* The values of \\( P_1 \\) to \\( P_5 \\) are listed in the source code. +* +* +* +* The computation of \\( e^{r} \\) thus becomes +* +* tex +* \begin{align*} +* e^{r} &= 1 + \frac{2r}{R-r} \\ +* &= 1 + r + \frac{r \cdot R_1(r)}{2 - R_1(r)}\ \text{for better accuracy} +* \end{align*} +* +* +* where +* +* tex +* R_1(r) = r - P_1\ r^2 + P_2\ r^4 + \ldots + P_5\ r^{10} +* +* +* 3. We scale back to obtain \\( e^{x} \\). From step 1, we have +* +* tex +* e^{x} = 2^k e^{r} +* +* +* +* ## Special Cases +* +* tex +* \begin{align*} +* e^\infty &= \infty \\ +* e^{-\infty} &= 0 \\ +* e^{\mathrm{NaN}} &= \mathrm{NaN} \\ +* e^0 &= 1\ \mathrm{is\ exact\ for\ finite\ argument\ only} +* \end{align*} +* +* +* ## Notes +* +* - According to an error analysis, the error is always less than \\(1\\) ulp (unit in the last place). +* +* - For an IEEE double, +* +* - if \\(x > 7.09782712893383973096\mbox{e+}02\\), then \\(e^{x}\\) overflows +* - if \\(x < -7.45133219101941108420\mbox{e+}02\\), then \\(e^{x}\\) underflows +* +* - The hexadecimal values included in the source code are the intended ones for the used constants. Decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the intended hexadecimal values. +* +* @param x input value +* @return output value +* +* @example +* double out = stdlib_base_exp( 0.0 ); +* // returns 1.0 +*/ +double stdlib_base_exp( const double x ) { + double hi; + double lo; + int32_t k; + + if ( stdlib_base_is_nan( x ) || x == STDLIB_CONSTANT_FLOAT64_PINF ) { + return x; + } + if ( x == STDLIB_CONSTANT_FLOAT64_NINF ) { + return 0.0; + } + if ( x > EXP_OVERFLOW ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( x < EXP_UNDERFLOW ) { + return 0.0; + } + if ( x > NEG_NEARZERO && x < NEARZERO ) { + return 1.0 + x; + } + // Reduce and compute r = hi - lo for extra precision... + if ( x < 0.0 ) { + k = stdlib_base_trunc( ( LOG2_E * x ) - 0.5 ); + } else { + k = stdlib_base_trunc( ( LOG2_E * x ) + 0.5 ); + } + hi = x - ( k * LN2_HI ); + lo = k * LN2_LO; + + return expmulti( hi, lo, k ); +} + +// 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1 +static const int32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; + +static const double HUGE_VALUE = 1.0e300; + +static const double TINY_VALUE = 1.0e-300; + +// 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 +static const int32_t HIGH_MIN_NORMAL_EXP = 0x00100000; + +// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 +static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; + +// High: LN2 +static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 + +// Low: LN2 +static const double LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39 + +// 1/LN2 +static const double INV_LN2 = 1.44269504088896338700e+00; // 0x3FF71547, 0x652B82FE + +// High (24 bits): 1/LN2 +static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60000000 + +// Low: 1/LN2 +static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 + +// 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 +static const int32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; + +// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 +static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; + +// 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 +static const int32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; + +// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 +static const int32_t HIGH_SIGNIFICAND_HALF = 0x00080000; + +static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 + +// 2/(3*LN2) +static const double CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD + +// (float)CP +static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 + +// Low: CP_HI +static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 + +static const double BP [2] = { 1.0, 1.5 }; + +static const double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 + +static const double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 + +// 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 +static const int32_t HIGH_BIASED_EXP_31 = 0x41e00000; + +// 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 +static const int32_t HIGH_BIASED_EXP_64 = 0x43f00000; + +// 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 +static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; + +// 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 +static const int32_t HIGH_1075 = 0x4090cc00; + +// 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 +static const uint32_t HIGH_NEG_1075 = 0xc090cc00; + +static const int32_t HIGH_NUM_NONSIGN_BITS = 31; + +// -(1024-log2(ovfl+.5ulp)) +static const double OVT = 8.0085662595372944372e-17; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_l + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_l( const double x ) { + return 0.5999999999999946 + (x * (0.4285714285785502 + (x * (0.33333332981837743 + (x * (0.272728123808534 + (x * (0.23066074577556175 + (x * 0.20697501780033842))))))))); +} + +// END: polyval_l + +// BEGIN: polyval_w + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_w( const double x ) { + return 0.5 + (x * (-0.3333333333333333 + (x * 0.25))); +} + +// END: polyval_w + +// BEGIN: polyval_p + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_p( const double x ) { + return 0.16666666666666602 + (x * (-0.0027777777777015593 + (x * (0.00006613756321437934 + (x * (-0.0000016533902205465252 + (x * 4.1381367970572385e-8))))))); +} + +// END: polyval_p + +/* End auto-generated functions. */ + +/** +* Evaluates the exponential function when \\( y = \pm \infty\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_infinite( -1.0, 1.0/0.0 ); +* // returns NaN +* +* @example +* double out = y_is_infinite( 1.5, -1.0/0.0 ); +* // returns 0.0 +*/ +static double y_is_infinite( const double x, const double y ) { + if ( x == -1.0 ) { + return 0.0 / 0.0; // NaN + } + if ( x == 1.0 ) { + return 1.0; + } + // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) + if ( ( stdlib_base_abs( x ) < 1.0 ) == ( y == STDLIB_CONSTANT_FLOAT64_PINF ) ) { + return 0.0; + } + // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Evaluates the exponential function when \\(|y| > 2^64\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_huge( 9.0, 3.6893488147419103e19 ); +* // returns Infinity +* +* @example +* double out = y_is_huge( -3.14, -3.6893488147419103e19 ); +* // returns 0.0 +*/ +static double y_is_huge( const double x, const double y ) { + uint32_t ahx; + uint32_t hx; + + stdlib_base_float64_get_high_word( x, &hx ); + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + if ( ahx <= HIGH_MAX_NEAR_UNITY ) { + if ( y < 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; + } + // x has a biased exponent greater than or equal to 0... + + if ( y > 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; +} + +/** +* Evaluates the exponential function when \\(|x| = 0\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = x_is_zero( 0.0, 2 ); +* // returns 0.0 +* +* @example +* double out = x_is_zero( 0.0, -9 ); +* // returns Infinity +*/ +static double x_is_zero( const double x, const double y ) { + if ( y == STDLIB_CONSTANT_FLOAT64_NINF ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( y == STDLIB_CONSTANT_FLOAT64_PINF ) { + return 0.0; + } + if ( y > 0.0 ) { + if ( stdlib_base_is_odd( y ) ) { + return x; // handles +-0 + } + return 0.0; + } + // y < 0.0 + if ( stdlib_base_is_odd( y ) ) { + return stdlib_base_copysign( STDLIB_CONSTANT_FLOAT64_PINF, x ); // handles +-0 + } + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Computes \\(2^{\mathrm{hp} + \mathrm{lp}\\). +* +* @param j high word of hp + lp +* @param hp first power summand +* @param lp second power summand +* @return function value +* +* @example +* double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); +* // returns ~0.79 +*/ +static double pow2( uint32_t j, const double hp, const double lp ) { + uint32_t tmp; + double hpc; + int32_t jc; + int32_t nc; + int32_t i; + int32_t k; + int32_t n; + double t1; + double t; + double r; + double u; + double v; + double w; + double z; + + hpc = hp; + jc = (int32_t)j; + i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + k = ( ( i >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); + n = 0; + nc = (int32_t)n; + + // |z| > 0.5, set n = z+0.5 + if ( i > HIGH_BIASED_EXP_NEG_1 ) { + n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ); + k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n + tmp = ( ( n & ~( HIGH_SIGNIFICAND_MASK >> k ) ) ); + t = 0.0; + stdlib_base_float64_set_high_word( tmp, &t ); + n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS - k ) ); + nc = (int32_t)n; + if ( jc < 0 ) { + nc = -nc; + } + hpc -= t; + } + t = lp + hpc; + stdlib_base_float64_set_low_word( 0, &t ); + u = t * LN2_HI; + v = ( ( lp - ( t - hpc ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); + z = u + v; + w = v - ( z - u ); + t = z * z; + t1 = z - ( t * polyval_p( t ) ); + r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); + z = 1.0 - ( r - z ); + stdlib_base_float64_get_high_word( z, &j ); + j = j + ( nc << STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ); + jc = (int32_t)j; + + // Check for subnormal output... + if ( ( jc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) <= 0 ) { + z = stdlib_base_ldexp( z, nc ); + } else { + stdlib_base_float64_set_high_word( j, &z ); + } + return z; +} + +/** +* Computes \\(\operatorname{log}(x)\\) assuming \\(|1-x|\\) is small and using the approximation \\(x - x^2/2 + x^3/3 - x^4/4\\). +* +* @param ax absolute value of x +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* logx( 9.0, &o1, &o2 ); +*/ +void logx( const double ax, double *o1, double *o2 ) { + double t2; + double t1; + double t; + double w; + double u; + double v; + + t = ax - 1.0; // t has 20 trailing zeros + w = t * t * polyval_w( t ); + u = INV_LN2_HI * t; // INV_LN2_HI has 21 significant bits + v = ( t * INV_LN2_LO ) - ( w * INV_LN2 ); + t1 = u + v; + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = v - (t1 - u); + + *o1 = t1; + *o2 = t2; +} + +/** +* Computes \\(\operatorname{log2}(ax)\\). +* +* @param ax absolute value of x +* @param ahx high word of ax +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* log2ax( 9.0, 1075970048, &o1, &o2 ); +*/ +void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { + uint32_t ahxcc; + uint32_t tmp; + int32_t ahxc; + double axc; + int32_t n; + int32_t j; + int32_t k; + double ss; // hs + ls + double s2; // ss squared + double hs; + double ls; + double ht; + double lt; + double bp; // BP constant + double dp; // DP constant + double hp; + double lp; + double hz; + double lz; + double t1; + double t2; + double t; + double r; + double u; + double v; + + n = 0; + axc = ax; + ahxc = ahx; + ahxcc = (uint32_t)ahx; + // Check if x is subnormal... + if ( ahxc < HIGH_MIN_NORMAL_EXP ) { + axc *= TWO53; + n -= 53; + stdlib_base_float64_get_high_word( axc, &ahxcc ); + ahxc = (int32_t)ahxcc; + } + // Extract the unbiased exponent of x: + n += ( ( ahxc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); + + // Isolate the significand bits of x: + j = ( ahxc & HIGH_SIGNIFICAND_MASK ); + + // Normalize ahx by setting the (biased) exponent to 1023: + ahxc = ( j|HIGH_BIASED_EXP_0 ); + + // Determine the interval of |x| by comparing significand bits... + + // |x| < sqrt(3/2) + if ( j <= 0x3988E ) { // 0 00000000000 00111001100010001110 + k = 0; + } + // |x| < sqrt(3) + else if ( j < 0xBB67A ) { // 0 00000000000 10111011011001111010 + k = 1; + } + // |x| >= sqrt(3) + else { + k = 0; + n += 1; + ahxc -= HIGH_MIN_NORMAL_EXP; + } + ahxcc = (uint32_t)ahxc; + + // Load the normalized high word into |x|: + stdlib_base_float64_set_high_word( ahxcc, &axc ); + + // Compute ss = hs + ls = (x-1)/(x+1) or (x-1.5)/(x+1.5): + bp = BP[ k ]; // BP[0] = 1.0, BP[1] = 1.5 + u = axc - bp; // (x-1) || (x-1.5) + v = 1.0 / (axc + bp); // 1/(x+1) || 1/(x+1.5) + ss = u * v; + hs = ss; + stdlib_base_float64_set_low_word( 0, &hs ); // set all low word (less significant significand) bits to 0s + + // Compute ht = ax + bp (via manipulation, i.e., bit flipping, of the high word): + tmp = ( ( ahxc >> 1 ) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; + tmp += ( k << 18 ); // (k<<18) can be considered the word equivalent of 1.0 or 1.5 + ht = 0.0; + stdlib_base_float64_set_high_word( tmp, &ht ); + lt = axc - ( ht - bp ); + ls = v * ( ( u - ( hs * ht ) ) - ( hs * lt ) ); + + // Compute log(ax)... + + s2 = ss * ss; + r = s2 * s2 * polyval_l( s2 ); + r += ls * (hs + ss); + s2 = hs * hs; + ht = 3.0 + s2 + r; + stdlib_base_float64_set_low_word( 0, &ht ); + lt = r - ( ( ht - 3.0 ) - s2 ); + + // u+v = ss*(1+...): + u = hs * ht; + v = ( ls * ht ) + ( lt * ss ); + + // 2/(3LN2) * (ss+...): + hp = u + v; + stdlib_base_float64_set_low_word( 0, &hp ); + lp = v - ( hp - u ); + hz = CP_HI * hp; // CP_HI+CP_LO = 2/(3*LN2) + lz = ( CP_LO * hp ) + ( lp*CP ) + DP_LO[ k ]; + + // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz + dp = DP_HI[ k ]; + t = n; + t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); + + *o1 = t1; + *o2 = t2; +} + +/** +* Evaluates the exponential function. +* +* @param x base +* @param y exponent +* @return function value +* +* @example +* double out = stdlib_base_pow( 2.0, 3.0 ); +* // returns 8.0 +* +* @example +* double out = stdlib_base_pow( 4.0, 0.5 ); +* // returns 2.0 +*/ +double stdlib_base_pow( const double x, const double y ) { + uint32_t hx; // high word x + uint32_t lx; // low word x + uint32_t hy; // high word y + uint32_t ly; // low word y + int32_t ahx; // absolute value high word x + int32_t ahy; // absolute value high word y + uint32_t i; + uint32_t j; + int32_t ic; + int32_t jc; + int32_t sx; // sign x + int32_t sy; // sign y + double ax; // absolute value x + double y1; + double hp; + double lp; + double t1; + double t2; + double xc; + double z; // y prime + + xc = x; + if ( stdlib_base_is_nan( xc ) || stdlib_base_is_nan( y ) ) { + return 0.0/0.0; // NaN + } + + // Split y into high and low words: + stdlib_base_float64_to_words( y, &hy, &ly ); + + // Special cases y... + if ( ly == 0 ) { + if ( y == 0.0 ) { + return 1.0; + } + if ( y == 1.0 ) { + return xc; + } + if ( y == -1.0 ) { + return 1.0 / xc; + } + if ( y == 0.5 ) { + return stdlib_base_sqrt( xc ); + } + if ( y == -0.5 ) { + return 1.0 / stdlib_base_sqrt( xc ); + } + if ( y == 2.0 ) { + return xc * xc; + } + if ( y == 3.0 ) { + return xc * xc * xc; + } + if ( y == 4.0 ) { + xc *= xc; + return xc * xc; + } + if ( stdlib_base_is_infinite( y ) ) { + return y_is_infinite( xc, y ); + } + } + // Split x into high and low words: + stdlib_base_float64_to_words( xc, &hx, &lx ); + + // Special cases x... + if ( lx == 0 ) { + if ( hx == 0 ) { + return x_is_zero( xc, y ); + } + if ( xc == 1.0 ) { + return 1.0; + } + if ( + xc == -1.0 && + stdlib_base_is_odd( y ) + ) { + return -1.0; + } + if ( stdlib_base_is_infinite( xc ) ) { + if ( xc == STDLIB_CONSTANT_FLOAT64_NINF ) { + // stdlib_base_pow( 1/x, -y ) + return stdlib_base_pow( -0.0, -y ); + } + if ( y < 0.0 ) { + return 0.0; + } + return STDLIB_CONSTANT_FLOAT64_PINF; + } + } + if ( + xc < 0.0 && + stdlib_base_is_integer( y ) == false + ) { + // Signal NaN... + return 0.0/0.0; + } + ax = stdlib_base_abs( x ); + + // Remove the sign bits (i.e., get absolute values): + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + ahy = ( hy & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + // Extract the sign bits: + sx = ( hx >> HIGH_NUM_NONSIGN_BITS ); + sy = ( hy >> HIGH_NUM_NONSIGN_BITS ); + + // Determine the sign of the result... + if ( sx && stdlib_base_is_odd( y ) ) { + sx = -1; + } else { + sx = 1; + } + // Case 1: |y| is huge... + + // |y| > 2^31 + if ( ahy > HIGH_BIASED_EXP_31 ) { + // |y| > 2^64, then must over- or underflow... + if ( ahy > HIGH_BIASED_EXP_64 ) { + return y_is_huge( xc, y ); + } + // Over- or underflow if x is not close to unity... + if ( ahx < HIGH_MAX_NEAR_UNITY ) { + // y < 0 + if ( sy == 1 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( ahx > HIGH_BIASED_EXP_0 ) { + // y > 0 + if ( sy == 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + // At this point, |1-x| is tiny (<= 2^-20). Suffice to compute log(x) by x - x^2/2 + x^3/3 - x^4/4. + logx( ax, &t1, &t2 ); + } + // Case 2: |y| is not huge... + else { + log2ax( ax, ahx, &t1, &t2 ); + } + + // Split y into y1 + y2 and compute (y1+y2) * (t1+t2)... + y1 = y; + stdlib_base_float64_set_low_word( 0, &y1 ); + lp = ( ( y - y1 ) * t1 ) + ( y * t2 ); + hp = y1 * t1; + z = lp + hp; + + // Note: can be more performant to use getHighWord and getLowWord directly, but using toWords looks cleaner. + stdlib_base_float64_to_words( z, &j, &i ); + jc = (int32_t)j; + ic = (int32_t)i; + + // z >= 1024 + if ( jc >= HIGH_BIASED_EXP_10 ) { + // z > 1024 + if ( ( ( jc - HIGH_BIASED_EXP_10 )|ic ) != 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + if ( ( lp + OVT ) > ( z - hp ) ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + } + // z <= -1075 + else if ( ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { + // z < -1075 + if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( lp <= ( z - hp ) ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + } + // Compute 2^(hp+lp)... + z = pow2( j, hp, lp ); + return sx * z; +} + +double stdlib_base_sqrt( const double x ) { + return sqrt( x ); +} + +static const double TINY = 1.0e-300; + +// 2*-56 = 1/(2*56) = 1/72057594037927940 +static const double SMALL = 1.3877787807814457e-17; + +static const double ERX = 8.45062911510467529297e-1; // 0x3FEB0AC1, 0x60000000 + +static const double PPC = 1.28379167095512558561e-1; // 0x3FC06EBA, 0x8214DB68 +static const double QQC = 1.0; + +static const double PAC = -2.36211856075265944077e-3; // 0xBF6359B8, 0xBEF77538 +static const double QAC = 1.0; + +static const double RAC = -9.86494403484714822705e-3; // 0xBF843412, 0x600D6435 +static const double SAC = 1.0; + +static const double RBC = -9.86494292470009928597e-3; // 0xBF843412, 0x39E86F4A +static const double SBC = 1.0; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_pa + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_pa( const double x ) { + return 0.41485611868374833 + (x * (-0.3722078760357013 + (x * (0.31834661990116175 + (x * (-0.11089469428239668 + (x * (0.035478304325618236 + (x * -0.002166375594868791))))))))); +} + +// END: polyval_pa + +// BEGIN: polyval_pp + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_pp( const double x ) { + return -0.3250421072470015 + (x * (-0.02848174957559851 + (x * (-0.005770270296489442 + (x * -0.000023763016656650163))))); +} + +// END: polyval_pp + +// BEGIN: polyval_qa + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_qa( const double x ) { + return 0.10642088040084423 + (x * (0.540397917702171 + (x * (0.07182865441419627 + (x * (0.12617121980876164 + (x * (0.01363708391202905 + (x * 0.011984499846799107))))))))); +} + +// END: polyval_qa + +// BEGIN: polyval_qq + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_qq( const double x ) { + return 0.39791722395915535 + (x * (0.0650222499887673 + (x * (0.005081306281875766 + (x * (0.00013249473800432164 + (x * -0.000003960228278775368))))))); +} + +// END: polyval_qq + +// BEGIN: polyval_ra + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_ra( const double x ) { + return -0.6938585727071818 + (x * (-10.558626225323291 + (x * (-62.375332450326006 + (x * (-162.39666946257347 + (x * (-184.60509290671104 + (x * (-81.2874355063066 + (x * -9.814329344169145))))))))))); +} + +// END: polyval_ra + +// BEGIN: polyval_rb + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_rb( const double x ) { + return -0.799283237680523 + (x * (-17.757954917754752 + (x * (-160.63638485582192 + (x * (-637.5664433683896 + (x * (-1025.0951316110772 + (x * -483.5191916086514))))))))); +} + +// END: polyval_rb + +// BEGIN: polyval_sa + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_sa( const double x ) { + return 19.651271667439257 + (x * (137.65775414351904 + (x * (434.56587747522923 + (x * (645.3872717332679 + (x * (429.00814002756783 + (x * (108.63500554177944 + (x * (6.570249770319282 + (x * -0.0604244152148581))))))))))))); +} + +// END: polyval_sa + +// BEGIN: polyval_sb + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_sb( const double x ) { + return 30.33806074348246 + (x * (325.7925129965739 + (x * (1536.729586084437 + (x * (3199.8582195085955 + (x * (2553.0504064331644 + (x * (474.52854120695537 + (x * -22.44095244658582))))))))))); +} + +// END: polyval_sb + +/* End auto-generated functions. */ + +/** +* Evaluates the complementary error function. +* +* tex +* \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int^{x}_{0} e^{-t^2}\ \mathrm{dt} +* +* +* Note that +* +* tex +* \begin{align*} +* \operatorname{erfc}(x) &= 1 - \operatorname{erf}(x) \\ +* \operatorname{erf}(-x) &= -\operatorname{erf}(x) \\ +* \operatorname{erfc}(-x) &= 2 - \operatorname{erfc}(x) +* \end{align*} +* +* +* ## Method +* +* 1. For \\(|x| \in [0, 0.84375)\\), +* +* tex +* \operatorname{erf}(x) = x + x \cdot \operatorname{R}(x^2) +* +* +* and +* +* tex +* \operatorname{erfc}(x) = \begin{cases} +* 1 - \operatorname{erf}(x) & \textrm{if}\ x \in (-.84375,0.25) \\ +* 0.5 + ((0.5-x)-x \mathrm{R}) & \textrm{if}\ x \in [0.25,0.84375) +* \end{cases} +* +* +* where \\(R = P/Q\\) and where \\(P\\) is an odd polynomial of degree \\(8\\) and \\(Q\\) is an odd polynomial of degree \\(10\\). +* +* tex +* \biggl| \mathrm{R} - \frac{\operatorname{erf}(x)-x}{x} \biggr| \leq 2^{-57.90} +* +* +* +* +* The formula is derived by noting +* +* tex +* \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\biggl(x - \frac{x^3}{3} + \frac{x^5}{10} - \frac{x^7}{42} + \ldots \biggr) +* +* +* and that +* +* tex +* \frac{2}{\sqrt{\pi}} = 1.128379167095512573896158903121545171688 +* +* +* is close to unity. The interval is chosen because the fix point of \\(\operatorname{erf}(x)\\) is near \\(0.6174\\) (i.e., \\(\operatorname{erf(x)} = x\\) when \\(x\\) is near \\(0.6174\\)), and, by some experiment, \\(0.84375\\) is chosen to guarantee the error is less than one ulp for \\(\operatorname{erf}(x)\\). +* +* +* +* 2. For \\(|x| \in [0.84375,1.25)\\), let \\(s = |x|-1\\), and \\(c = 0.84506291151\\) rounded to single (\\(24\\) bits) +* +* tex +* \operatorname{erf}(x) = \operatorname{sign}(x) \cdot \biggl(c + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)}\biggr) +* +* +* and +* +* tex +* \operatorname{erfc}(x) = \begin{cases} +* (1-c) - \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)} & \textrm{if}\ x > 0 \\ +* 1 + \biggl(c + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)}\biggr) & \textrm{if}\ x < 0 +* \end{cases} +* +* +* where +* +* tex +* \biggl|\frac{\mathrm{P1}}{\mathrm{Q1}} - (\operatorname{erf}(|x|)-c)\biggr| \leq 2^{-59.06} +* +* +* +* +* Here, we use the Taylor series expansion at \\(x = 1\\) +* +* tex +* \begin{align*} +* \operatorname{erf}(1+s) &= \operatorname{erf}(1) + s\cdot \operatorname{poly}(s) \\ +* &= 0.845.. + \frac{\operatorname{P1}(s)}{\operatorname{Q1}(s)} +* \end{align*} +* +* +* using a rational approximation to approximate +* +* tex +* \operatorname{erf}(1+s) - (c = (\mathrm{single})0.84506291151) +* +* +* +* +* Note that, for \\(x \in [0.84375,1.25)\\), \\(|\mathrm{P1}/\mathrm{Q1}| < 0.078\\), where +* +* - \\(\operatorname{P1}(s)\\) is a degree \\(6\\) polynomial in \\(s\\) +* - \\(\operatorname{Q1}(s)\\) is a degree \\(6\\) polynomial in \\(s\\) +* +* 3. For \\(x \in [1.25,1/0.35)\\), +* +* tex +* \begin{align*} +* \operatorname{erfc}(x) &= \frac{1}{x}e^{-x^2-0.5625+(\mathrm{R1}/\mathrm{S1})} \\ +* \operatorname{erf}(x) &= 1 - \operatorname{erfc}(x) +* \end{align*} +* +* +* where +* +* - \\(\operatorname{R1}(z)\\) is a degree \\(7\\) polynomial in \\(z\\), where \\(z = 1/x^2\\) +* - \\(\operatorname{S1}(z)\\) is a degree \\(8\\) polynomial in \\(z\\) +* +* 4. For \\(x \in [1/0.35,28)\\), +* +* tex +* \operatorname{erfc}(x) = \begin{cases} +* \frac{1}{x} e^{-x^2-0.5625+(\mathrm{R2}/\mathrm{S2})} & \textrm{if}\ x > 0 \\ +* 2.0 - \frac{1}{x} e^{-x^2-0.5625+(\mathrm{R2}/\mathrm{S2})} & \textrm{if}\ -6 < x < 0 \\ +* 2.0 - \mathrm{tiny} & \textrm{if}\ x \leq -6 +* \end{cases} +* +* +* and +* +* tex +* \operatorname{erf}(x) = \begin{cases} +* \operatorname{sign}(x) \cdot (1.0 - \operatorname{erfc}(x)) & \textrm{if}\ x < 6 \\ +* \operatorname{sign}(x) \cdot (1.0 - \mathrm{tiny}) & \textrm{otherwise} +* \end{cases} +* +* +* where +* +* - \\(\operatorname{R2}(z)\\) is a degree \\(6\\) polynomial in \\(z\\), where \\(z = 1/x^2\\) +* - \\(\operatorname{S2}(z)\\) is a degree \\(7\\) polynomial in \\(z\\) +* +* 5. For \\(x \in [28, \infty)\\), +* +* tex +* \begin{align*} +* \operatorname{erf}(x) &= \operatorname{sign}(x) \cdot (1 - \mathrm{tiny}) & \textrm{(raise inexact)} +* \end{align*} +* +* +* and +* +* tex +* \operatorname{erfc}(x) = \begin{cases} +* \mathrm{tiny} \cdot \mathrm{tiny} & \textrm{if}\ x > 0\ \textrm{(raise underflow)} \\ +* 2 - \mathrm{tiny} & \textrm{if}\ x < 0 +* \end{cases} +* +* +* +* ## Special Cases +* +* tex +* \begin{align*} +* \operatorname{erf}(0) &= 0 \\ +* \operatorname{erf}(-0) &= -0 \\ +* \operatorname{erf}(\infty) &= 1 \\ +* \operatorname{erf}(-\infty) &= -1 \\ +* \operatorname{erfc}(0) &= 1 \\ +* \operatorname{erfc}(\infty) &= 0 \\ +* \operatorname{erfc}(-\infty) &= 2 \\ +* \operatorname{erf}(\mathrm{NaN}) &= \mathrm{NaN} \\ +* \operatorname{erfc}(\mathrm{NaN}) &= \mathrm{NaN} +* \end{align*} +* +* +* +* ## Notes +* +* - To compute \\(\exp(-x^2-0.5625+(\mathrm{R}/\mathrm{S}))\\), let \\(s\\) be a single precision number and \\(s := x\\); then +* +* tex +* -x^2 = -s^2 + (s-x)(s+x) +* +* +* and +* +* tex +* e^{-x^2-0.5626+(\mathrm{R}/\mathrm{S})} = e^{-s^2-0.5625} e^{(s-x)(s+x)+(\mathrm{R}/\mathrm{S})} +* +* +* - #4 and #5 make use of the asymptotic series +* +* tex +* \operatorname{erfc}(x) \approx \frac{e^{-x^2}}{x\sqrt{\pi}} (1 + \operatorname{poly}(1/x^2)) +* +* +* We use a rational approximation to approximate +* +* tex +* g(s) = f(1/x^2) = \ln(\operatorname{erfc}(x) \cdot x) - x^2 + 0.5625 +* +* +* - The error bound for \\(\mathrm{R1}/\mathrm{S1}\\) is +* +* tex +* |\mathrm{R1}/\mathrm{S1} - f(x)| < 2^{-62.57} +* +* +* and for \\(\mathrm{R2}/\mathrm{S2}\\) is +* +* tex +* |\mathrm{R2}/\mathrm{S2} - f(x)| < 2^{-61.52} +* +* +* @param x input value +* @return output value +*S +* @example +* double out = stdlib_base_erfc( 0.0 ); +* // returns 1.0 +*/ +double stdlib_base_erfc( const double x ) { + int32_t sign; + double ax; + double z; + double r; + double s; + double y; + double p; + double q; + + // Special case: NaN + if ( stdlib_base_is_nan( x ) ) { + return 0.0 / 0.0; // NaN + } + // Special case: +infinity + if ( x == STDLIB_CONSTANT_FLOAT64_PINF ) { + return 0.0; + } + // Special case: -infinity + if ( x == STDLIB_CONSTANT_FLOAT64_NINF ) { + return 2.0; + } + // Special case: +-0 + if ( x == 0.0 ) { + return 1.0; + } + if ( x < 0.0 ) { + sign = 1; + ax = -x; + } else { + sign = 0; + ax = x; + } + // |x| < 0.84375 + if ( ax < 0.84375 ) { + if ( ax < SMALL ) { + return 1.0 - x; // raise inexact + } + z = x * x; + r = PPC + ( z * polyval_pp( z ) ); + s = QQC + ( z * polyval_qq( z ) ); + y = r / s; + + // x < 1/4 + if ( x < 0.25 ) { + return 1.0 - ( x + ( x * y ) ); + } + r = x * y; + r += x - 0.5; + return 0.5 - r; + } + // 0.84375 <= |x| < 1.25 + if ( ax < 1.25 ) { + s = ax - 1.0; + p = PAC + ( s * polyval_pa( s ) ); + q = QAC + ( s * polyval_qa( s ) ); + if ( sign == 1 ) { + return 1.0 + ERX + ( p / q ); + } + return 1.0 - ERX - ( p / q ); + } + // |x| < 28 + if ( ax < 28.0 ) { + s = 1.0 / ( ax * ax ); + + // |x| < 1/0.35 ~ 2.857143 + if ( ax < 2.857142857142857 ) { + r = RAC + ( s * polyval_ra( s ) ); + s = SAC + ( s * polyval_sa( s ) ); + } + // |x| >= 1/0.35 ~ 2.857143 + else { + // x < -6 + if ( x < -6.0 ) { + return 2.0 - TINY; // raise inexact + } + r = RBC + ( s * polyval_rb( s ) ); + s = SBC + ( s * polyval_sb( s ) ); + } + z = ax; + stdlib_base_float64_set_low_word( 0, &z ); // pseudo-single (20-bit) precision x + r = stdlib_base_exp( -( z * z ) - 0.5625 ) * stdlib_base_exp( ( ( z - ax ) * ( z + ax ) ) + ( r / s ) ); + if ( sign == 1 ) { + return 2.0 - ( r / ax ); + } + return r / ax; + } + if ( sign == 1 ) { + return 2.0 - TINY; // raise inexact; ~2 + } + return TINY * TINY; // raise inexact; ~0 +} + +#ifndef STDLIB_CONSTANTS_FLOAT64_PI_H +#define STDLIB_CONSTANTS_FLOAT64_PI_H + +/** +* Macro for π. +*/ +#define STDLIB_CONSTANT_FLOAT64_PI 3.141592653589793238 + +#endif // !STDLIB_CONSTANTS_FLOAT64_PI_H + + +double cdf(double x, double mu) { + double denom, xc; + if (stdlib_base_is_nan(x) || stdlib_base_is_nan(mu)) { + return NAN; + } + denom = stdlib_base_sqrt(2.0); + xc = x - mu; + return 0.5 * stdlib_base_erfc(-xc / denom); +} + +var normalCDF = cdf( 0.0, 1.0 ); + +double pdf(double x, double a, double b, double mu, double sigma) { + double s2x2, A, B, C; + + if (stdlib_base_is_nan(x) || stdlib_base_is_nan(a) || stdlib_base_is_nan(b) || sigma <= 0.0 || a >= b) { + return NAN; + } + if (x < a || x > b) { + return 0.0; + } + + s2x2 = 2.0 * stdlib_base_pow(sigma, 2.0); + A = 1.0 / (stdlib_base_sqrt(s2x2 * STDLIB_CONSTANT_FLOAT64_PI)); + B = -1.0 / (s2x2); + C = normalCDF((b - mu) / sigma) - normalCDF((a - mu) / sigma); + + return A * stdlib_base_exp(B * stdlib_base_pow(x - mu, 2.0)) / C; +} diff --git a/lib/node_modules/@stdlib/stats/base/smeanors/benchmark/c/benchmark.length.c b/lib/node_modules/@stdlib/stats/base/smeanors/benchmark/c/benchmark.length.c index 6d0d385a7436..3d852ca11823 100644 --- a/lib/node_modules/@stdlib/stats/base/smeanors/benchmark/c/benchmark.length.c +++ b/lib/node_modules/@stdlib/stats/base/smeanors/benchmark/c/benchmark.length.c @@ -16,7 +16,7 @@ * limitations under the License. */ -#include "stdlib/stats/base/smeanors.h" +#include "stdlib/stats/base/smeanors/include/stdlib/stats/base" #include #include #include diff --git a/lib/node_modules/@stdlib/strided/base/dmskmap/examples/c/example.c b/lib/node_modules/@stdlib/strided/base/dmskmap/examples/c/example.c index b7f829860183..e91615871552 100644 --- a/lib/node_modules/@stdlib/strided/base/dmskmap/examples/c/example.c +++ b/lib/node_modules/@stdlib/strided/base/dmskmap/examples/c/example.c @@ -28,10 +28,10 @@ static double scale( const double x ) { int main( void ) { // Create an input strided array: - double X[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + const double X[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; // Create a mask strided array: - uint8_t M[] = { 0, 0, 1, 0, 0, 1 }; + const uint8_t M[] = { 0, 0, 1, 0, 0, 1 }; // Create an output strided array: double Y[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };