diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c new file mode 100644 index 000000000000..3a7b90f3cd9c --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c @@ -0,0 +1,592 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2022 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "stdlib/math/base/special/gammaln.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/special/exp.h" +#include "stdlib/math/base/special/gamma.h" +#include "stdlib/constants/float64/ln_sqrt_two_pi.h" +#include "stdlib/constants/float64/sqrt_two_pi.h" +#include "stdlib/constants/float32/smallest_normal.h" +#include "stdlib/constants/float32/max.h" + +static const double C6 = 0.30865217988013567769; +static const double A[] ={ + 1.996379051590076518221, + -0.17971032528832887213e-2, + 0.131292857963846713e-4, + -0.2340875228178749e-6, + 0.72291210671127e-8, + -0.3280997607821e-9, + 0.198750709010e-10, + -0.15092141830e-11, + 0.1375340084e-12, + -0.145728923e-13, + 0.17532367e-14, + -0.2351465e-15, + 0.346551e-16, + -0.55471e-17, + 0.9548e-18, + -0.1748e-18, + 0.332e-19, + -0.58e-20 +}; + +/** +* 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 +* +* @private +* @param {number} x - value at which to evaluate the polynomial +* @returns {number} evaluated polynomial +*/ +static double polyval_ak1( double x ) { + if ( x == 0.0 ) { + return 0.0; + } + return 0.0 + (x * (1.0 + (x * (1.0 + (x * (1.5 + (x * (2.6666666666666665 + (x * (5.208333333333333 + (x * 10.8))))))))))); // eslint-disable-line max-len +} +/** +* 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 +* +* @private +* @param {number} x - value at which to evaluate the polynomial +* @returns {number} evaluated polynomial +*/ +static double polyval_ak2( double x ) { + if ( x == 0.0 ) { + return 1.0; + } + return 1.0 + (x * (1.0 + (x * (0.3333333333333333 + (x * (0.027777777777777776 + (x * (-0.003703703703703704 + (x * (0.0002314814814814815 + (x * 0.00005878894767783657))))))))))); // eslint-disable-line max-len +} +/** +* 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 +* +* @private +* @param {number} x - value at which to evaluate the polynomial +* @returns {number} evaluated polynomial +*/ +static double polyval_c( double x ) { + if ( x == 0.0 ) { + return 0.025721014990011306; + } + return 0.025721014990011306 + (x * (0.08247596616699963 + (x * (-0.0025328157302663564 + (x * (0.0006099292666946337 + (x * (-0.00033543297638406 + (x * 0.000250505279903))))))))); // eslint-disable-line max-len +} +/** +* 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 +* +* @private +* @param {number} x - value at which to evaluate the polynomial +* @returns {number} evaluated polynomial +*/ +static double polyval_d( double x ) { + if ( x == 0.0 ) { + return 0.08333333333333333; + } + return 0.08333333333333333 + (x * (-0.002777777777777778 + (x * (0.0007936507936507937 + (x * -0.0005952380952380953))))); // eslint-disable-line max-len +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak0bk0( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.3333333333438; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.3333333333438 + (x * (-0.2070740359969 + (x * (-0.05041806657154 + (x * (-0.004923635739372 + (x * -0.00004293658292782))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.7045554412463 + (x * (0.2118190062224 + (x * (0.03048648397436 + (x * 0.001605037988091))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = -0.00004293658292782 + (x * (-0.004923635739372 + (x * (-0.05041806657154 + (x * (-0.2070740359969 + (x * -0.3333333333438))))))); // eslint-disable-line max-len + s2 = 0.001605037988091 + (x * (0.03048648397436 + (x * (0.2118190062224 + (x * (0.7045554412463 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak1bk1( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0172847633523; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0172847633523 + (x * (-0.0159372646475 + (x * (-0.00464910887221 + (x * (-0.00060683488776 + (x * -0.00000614830384279))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.764050615669 + (x * (0.297143406325 + (x * (0.0579490176079 + (x * 0.00574558524851))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = -0.00000614830384279 + (x * (-0.00060683488776 + (x * (-0.00464910887221 + (x * (-0.0159372646475 + (x * -0.0172847633523))))))); // eslint-disable-line max-len + s2 = 0.00574558524851 + (x * (0.0579490176079 + (x * (0.297143406325 + (x * (0.764050615669 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak2bk2( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0172839517431; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0172839517431 + (x * (-0.0146362417966 + (x * (-0.00357406772616 + (x * (-0.000391032032692 + (x * 0.00000249634036069))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.690560400696 + (x * (0.249962384741 + (x * (0.0443843438769 + (x * 0.00424073217211))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = 0.00000249634036069 + (x * (-0.000391032032692 + (x * (-0.00357406772616 + (x * (-0.0146362417966 + (x * -0.0172839517431))))))); // eslint-disable-line max-len + s2 = 0.00424073217211 + (x * (0.0443843438769 + (x * (0.249962384741 + (x * (0.690560400696 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak3bk3( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.99994466948; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.99994466948 + (x * (104.649839762 + (x * (857.204033806 + (x * (731.901559577 + (x * 45.5174411671))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (104.526456943 + (x * (823.313447808 + (x * (3119.93802124 + (x * 3970.03311219))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = 45.5174411671 + (x * (731.901559577 + (x * (857.204033806 + (x * (104.649839762 + (x * 0.99994466948))))))); // eslint-disable-line max-len + s2 = 3970.03311219 + (x * (3119.93802124 + (x * (823.313447808 + (x * (104.526456943 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak4bk4( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.0495346498136; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.0495346498136 + (x * (0.0299521337141 + (x * (0.00688296911516 + (x * (0.000512634846317 + (x * -0.0000201411722031))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.759803615283 + (x * (0.261547111595 + (x * (0.0464854522477 + (x * 0.00403751193496))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = -0.0000201411722031 + (x * (0.000512634846317 + (x * (0.00688296911516 + (x * (0.0299521337141 + (x * 0.0495346498136))))))); // eslint-disable-line max-len + s2 = 0.00403751193496 + (x * (0.0464854522477 + (x * (0.261547111595 + (x * (0.759803615283 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak5bk5( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.00452313583942; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.00452313583942 + (x * (0.00120744920113 + (x * (-0.0000789724156582 + (x * (-0.0000504476066942 + (x * -0.00000535770949796))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.912203410349 + (x * (0.405368773071 + (x * (0.0901638932349 + (x * 0.00948935714996))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = -0.00000535770949796 + (x * (-0.0000504476066942 + (x * (-0.0000789724156582 + (x * (0.00120744920113 + (x * 0.00452313583942))))))); // eslint-disable-line max-len + s2 = 0.00948935714996 + (x * (0.0901638932349 + (x * (0.405368773071 + (x * (0.912203410349 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak6bk6( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.00439937562904; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.00439937562904 + (x * (0.000487225670639 + (x * (-0.000128470657374 + (x * (0.00000529110969589 + (x * 1.5716677175e-7))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (0.794435257415 + (x * (0.333094721709 + (x * (0.0703527806143 + (x * 0.00806110846078))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = 1.5716677175e-7 + (x * (0.00000529110969589 + (x * (-0.000128470657374 + (x * (0.000487225670639 + (x * 0.00439937562904))))))); // eslint-disable-line max-len + s2 = 0.00806110846078 + (x * (0.0703527806143 + (x * (0.333094721709 + (x * (0.794435257415 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak7bk7( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0011481191232; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0011481191232 + (x * (-0.112850923276 + (x * (1.51623048511 + (x * (-0.218472031183 + (x * 0.0730002451555))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (14.2482206905 + (x * (69.7360396285 + (x * (218.938950816 + (x * 277.067027185))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = 0.0730002451555 + (x * (-0.218472031183 + (x * (1.51623048511 + (x * (-0.112850923276 + (x * -0.0011481191232))))))); // eslint-disable-line max-len + s2 = 277.067027185 + (x * (218.938950816 + (x * (69.7360396285 + (x * (14.2482206905 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} +/** +* Evaluates a rational function, i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @private +* @param {number} x - value at which to evaluate the rational function +* @returns {number} evaluated rational function +*/ +static double rational_ak8bk8( double x ) { + double ax; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.000145727889667; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.000145727889667 + (x * (-0.290806748131 + (x * (-13.308504545 + (x * (199.722374056 + (x * -11.4311378756))))))); // eslint-disable-line max-len + s2 = 1.0 + (x * (139.612587808 + (x * (2189.01116348 + (x * (7115.24019009 + (x * 45574.6081453))))))); // eslint-disable-line max-len + } else { + x = 1.0 / x; + s1 = -11.4311378756 + (x * (199.722374056 + (x * (-13.308504545 + (x * (-0.290806748131 + (x * -0.000145727889667))))))); // eslint-disable-line max-len + s2 = 45574.6081453 + (x * (7115.24019009 + (x * (2189.01116348 + (x * (139.612587808 + (x * 1.0))))))); // eslint-disable-line max-len + } + return s1 / s2; +} + +/** +* Computes the sum of a Chebyshev polynomial. +* +* @private +* @param {PositiveInteger} n - degree of polynomial +* @param {number} t - input value +* @returns {number} Chebyshev sum +*/ +static double chepolsum( double n, double t ) { + CBLAS_INT k; + double tt; + double u0; + double u1; + double u2; + + u0 = 0.0; + u1 = 0.0; + tt = t + t; + k = n; + do { + u2 = u1; + u1 = u0; + u0 = ( tt*u1 ) - u2 + A[ k ]; + k -= 1; + } while ( k >= 0 ); + return ( u0-u2 ) / 2.0; +} +/** +* Computes the Stirling series corresponding to asymptotic series for the logarithm of the gamma function. +* +* ```tex +* \frac{1}{12x}-\frac{1}{360x^3}\ldots; x \ge 3 +* ``` +* +* @private +* @param {number} x - input value +* @returns {number} function value +*/ +static double stirling( double x ) { + double z; + if ( x < STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL ) { + return STDLIB_CONSTANT_FLOAT32_MAX; + } + if ( x < 1.0 ) { + return stdlib_base_gammaln( x+1.0 ) - ( (x+0.5) * stdlib_base_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI; + } + if ( x < 2.0 ) { + return stdlib_base_gammaln( x ) - ( (x-0.5) * stdlib_base_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI; + } + if ( x < 3.0 ) { + return stdlib_base_gammaln( x-1.0 ) - ( (x-0.5) * stdlib_base_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI + stdlib_base_ln( x-1.0 ); // eslint-disable-line max-len + } + if ( x < 12.0 ) { + z = ( 18.0/( x*x ) ) - 1.0; + return chepolsum( 17, z ) / ( 12.0*x ); + } + z = 1.0 / ( x * x ); + if ( x < 1000.0 ) { + return polyval_c( z ) / ( C6+z ) / x; + } + return polyval_d( z ) / x; +} +/** +* Evaluates the `eps2` function. +* +* @private +* @param {number} eta - eta value +* @returns {number} function value +*/ +static double eps2( double eta ) { + double lnmeta; + double x; + if ( eta < -5.0 ) { + x = eta * eta; + lnmeta = stdlib_base_ln( -eta ); + return ( 12.0 - x - ( 6.0*( lnmeta*lnmeta ) ) ) / ( 12.0*x*eta ); + } + if ( eta < -2.0 ) { + return rational_ak1bk1( eta ); + } + if ( eta < 2.0 ) { + return rational_ak2bk2( eta ); + } + if ( eta < 1000.0 ) { + x = 1.0 / eta; + return rational_ak3bk3( eta ) / ( -12.0*eta ); + } + return -1.0 / ( 12.0 * eta ); +} +/** +* Evaluates the `eps3` function. +* +* @private +* @param {number} eta - eta value +* @returns {number} function value +*/ +static double eps3( double eta ) { + double x; + double y; + + if ( eta < -8.0 ) { + x = eta * eta; + y = stdlib_base_ln( -eta ) / eta; + return ( -30.0 + ( eta*y*( (6.0*x*y*y)-12.0+x ) ) ) / ( 12.0*eta*x*x ); + } + if ( eta < -4.0 ) { + return rational_ak4bk4( eta ) / ( eta*eta ); + } + if ( eta < -2.0 ) { + return rational_ak5bk5( eta ); + } + if ( eta < 2.0 ) { + return rational_ak6bk6( eta ); + } + if ( eta < 10.0 ) { + x = 1.0 / eta; + return rational_ak7bk7( x ) / ( eta*eta ); + } + if ( eta < 100.0 ) { + x = 1.0 / eta; + return rational_ak8bk8( x ) / ( eta*eta ); + } + return -stdlib_base_ln( eta ) / ( 12.0*eta*eta*eta ); +} +/** +* Computes the regulated gamma function. +* +* @private +* @param {number} x - input value +* @returns {number} function value +*/ +static double gamstar( double x ) { + if ( x >= 3.0 ) { + return stdlib_base_exp( stirling(x) ); + } + if ( x > 0.0 ) { + return stdlib_base_gamma(x) / ( stdlib_base_exp( -x + ( ( x-0.5 ) * stdlib_base_ln(x) ) ) * STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); + } + // Case: x <= 0.0 + return STDLIB_CONSTANT_FLOAT32_MAX; +} +