From d24f74f769ea80f04cf30f9c518495ed8061eec2 Mon Sep 17 00:00:00 2001 From: Karan Anand Date: Sun, 11 May 2025 12:45:21 -0700 Subject: [PATCH 01/14] feat: temp commit --- .../math/base/special/gammainc/src/Makefile | 70 ++ .../math/base/special/gammainc/src/addon.c | 22 + .../math/base/special/gammainc/src/main.c | 679 ++++++++++++++++++ 3 files changed, 771 insertions(+) create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/src/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/src/addon.c create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/src/Makefile b/lib/node_modules/@stdlib/math/base/special/gammainc/src/Makefile new file mode 100644 index 000000000000..7733b6180cb4 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/src/Makefile @@ -0,0 +1,70 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2025 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. +#/ + +# VARIABLES # + +ifndef VERBOSE + QUIET := @ +else + QUIET := +endif + +# Determine the OS ([1][1], [2][2]). +# +# [1]: https://en.wikipedia.org/wiki/Uname#Examples +# [2]: http://stackoverflow.com/a/27776822/2225624 +OS ?= $(shell uname) +ifneq (, $(findstring MINGW,$(OS))) + OS := WINNT +else +ifneq (, $(findstring MSYS,$(OS))) + OS := WINNT +else +ifneq (, $(findstring CYGWIN,$(OS))) + OS := WINNT +else +ifneq (, $(findstring Windows_NT,$(OS))) + OS := WINNT +endif +endif +endif +endif + + +# RULES # + +#/ +# Removes generated files for building an add-on. +# +# @example +# make clean-addon +#/ +clean-addon: + $(QUIET) -rm -f *.o *.node + +.PHONY: clean-addon + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: clean-addon + +.PHONY: clean diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/src/addon.c b/lib/node_modules/@stdlib/math/base/special/gammainc/src/addon.c new file mode 100644 index 000000000000..82be2899cca0 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/src/addon.c @@ -0,0 +1,22 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2025 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/hyp2f1.h" +#include "stdlib/math/base/napi/quaternary.h" + +STDLIB_MATH_BASE_NAPI_MODULE_DDDD_D( stdlib_base_hyp2f1 ) diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c new file mode 100644 index 000000000000..b5553ca9bc82 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c @@ -0,0 +1,679 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2018 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. +* +* +* ## Notice +* +* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified according to stdlib conventions. +* +* ```text +* (C) Copyright John Maddock 2006-7, 2013-14. +* (C) Copyright Paul A. Bristow 2007, 2013-14. +* (C) Copyright Nikhar Agrawal 2013-14. +* (C) Christopher Kormanyos 2013-14. +* +* Use, modification and distribution are subject to the +* Boost Software License, Version 1.0. (See accompanying file +* LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) +* ``` +*/ + +#include "stdlib/math/base/special/gammainc.h" +#include "stdlib/math/base/special/gammaln.h" +#include "stdlib/math/base/special/floor.h" +#include "stdlib/math/base/special/gamma.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/math/base/special/exp.h" +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/constants/float64/sqrt_eps.h" +#include "stdlib/constants/float64/max.h" +#include "stdlib/constants/float64/sqrt_two_pi.h" +#include "stdlib/constants/float64/max_ln.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float64/max_nth_factorial.h" +#include + +static const double MACHEP = 1.11022302462515654042e-16; +static const double EPS = 1.0e-13; +static const double ETHRESH = 1.0e-12; +static const double MAX_ITERATIONS = 10000.0; + +static double hys2f1( double a, double b, const double c, const double x, double *loss ); + +/** +* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value +*/ +static bool upperGammaFraction( const double a, const double z ) { + double f = upperIncompleteGammaFract( a, z ); + return 1.0 / ( z - a + 1.0 + continuedFractionA( f ) ); +} + +/** +* Creates a function to evaluate a series expansion of the upper incomplete gamma fraction. +* +* @param a1 function parameter +* @param z1 function parameter +* @return series function +*/ +static bool upperIncompleteGammaFract( const double a1, const double z1 ) { + var z = z1 - a1 + 1.0; + var a = a1; + var k = 0; + return next; + + /** + * Calculate the next term of the series. + * + * @private + * @returns {Array} series expansion terms + */ + function next() { + k += 1; + z += 2.0; + return [ + k * (a - k), + z + ]; + } +} + +/** +* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value +*/ +static bool upperGammaFraction( const double a, const double z ) { + double f = upperIncompleteGammaFract( a, z ); + return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +} + +/** +* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value +*/ +static bool upperGammaFraction( const double a, const double z ) { + double f = upperIncompleteGammaFract( a, z ); + return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +} + +/** +* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value +*/ +static bool upperGammaFraction( const double a, const double z ) { + double f = upperIncompleteGammaFract( a, z ); + return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +} + +/** +* Evaluates the Gaussian hypergeometric function by two-term recurrence in `a`. +* +* ## Notes +* +* - This function helps prevent losing accuracy in the highly alternating hypergeometric series and allows a and b to be reduced to smaller values. +* +* ## References +* +* - AMS55 #15.2.10 +* +* @param a input value +* @param b input value +* @param c input value +* @param x input value +* @param loss starting loss of significance +* @return function value +*/ +static double hyp2f1ra( const double a, const double b, const double c, const double x, double *loss ) { + double err; + double da; + double f2; + double f1; + double f0; + double t; + double n; + + *loss = 0.0; + err = 0.0; + + if ( ( c < 0.0 && a <= c ) || ( c >= 0.0 && a >= c ) ) { + da = stdlib_base_round( a-c ); + } else { + da = stdlib_base_round( a ); + } + + t = a - da; + if ( stdlib_base_abs( da ) > MAX_ITERATIONS ) { + // Too expensive to compute, so return `NaN`: + *loss = 1.0; + return STDLIB_CONSTANT_FLOAT64_NAN; + } + + if ( da < 0.0 ) { + // Backward recurrence... + f1 = hys2f1( t, b, c, x, &err ); + *loss += err; + f0 = hys2f1( t-1.0, b, c, x, &err ); + *loss += err; + t -= 1.0; + for ( n = 1; n < -da; n++ ) { + f2 = f1; + f1 = f0; + + // eslint-disable-next-line max-len + f0 = -( ( ( ( (2.0*t)-c )-( t*x )+( b*x ) ) * f1 ) + ( ( t*(x-1.0) ) * f2 ) ) / ( c-t ); + t -= 1.0; + } + } else { + // Forward recurrence... + f1 = hys2f1( t, b, c, x, &err ); + *loss += err; + f0 = hys2f1( t+1.0, b, c, x, &err ); + *loss += err; + t += 1.0; + for ( n = 1; n < da; n++ ) { + f2 = f1; + f1 = f0; + + // eslint-disable-next-line max-len + f0 = -( ( ( ( (2.0*t)-c )-( t*x )+( b*x ) ) * f1 ) + ( (c-t)*f2 ) ) / ( t*(x-1.0) ); + t += 1.0; + } + } + return f0; +} + +/** +* Evaluates the power series expansion of Gaussian hypergeometric function. +* +* @param a input value +* @param b input value +* @param c input value +* @param x input value +* @param loss starting loss of significance +* @return function value +*/ +static double hys2f1( double a, double b, const double c, const double x, double *loss ) { + double intFlag; + double umax; + double f; + double g; + double h; + double k; + double m; + double s; + double u; + double i; + + intFlag = 0.0; + + if ( stdlib_base_abs( b ) > stdlib_base_abs( a ) ) { + f = b; + b = a; + a = f; + } + + if ( isNonPositiveInteger( b ) && ( stdlib_base_abs( b ) < stdlib_base_abs( a ) ) ) { + f = b; + b = a; + a = f; + intFlag = 1.0; + } + + // eslint-disable-next-line max-len + if ( ( ( stdlib_base_abs( a ) > stdlib_base_abs( c ) + 1.0 ) || intFlag ) && ( stdlib_base_abs( c-a ) > 2.0 ) && ( stdlib_base_abs( a ) > 2.0 ) ) { + return hyp2f1ra( a, b, c, x, loss ); + } + + i = 0.0; + umax = 0.0; + f = a; + g = b; + h = c; + s = 1.0; + u = 1.0; + k = 0.0; + + do { + if ( stdlib_base_abs( h ) < EPS ) { + *loss = 1.0; + return STDLIB_CONSTANT_FLOAT64_PINF; + } + m = k + 1.0; + u *= ( ( f+k ) * ( g+k ) * x / ( ( h+k ) * m ) ); + s += u; + k = stdlib_base_abs( u ); + if ( k > umax ) { + umax = k; + } + k = m; + i += 1.0; + if ( i > MAX_ITERATIONS ) { + *loss = 1.0; + return s; + } + } while ( s == 0.0 || stdlib_base_abs( u/s ) > MACHEP ); + + // Estimate the relative error due to truncation by the series: + *loss = ( ( MACHEP*umax ) / stdlib_base_abs( s ) ) + ( MACHEP*i ); + return s; +} + +/** +* Applies transformations for `|x|` near unity before performing a power series expansion. +* +* @param a input value +* @param b input value +* @param c input value +* @param x input value +* @param loss starting loss of significance +* @return function value +*/ +static double hyt2f1( const double a, const double b, const double c, const double x, double *loss ) { + double negIntA; + double negIntB; + double sign; + double err1; + double err; + double aid; + double ax; + double id; + double d1; + double d2; + double y1; + double i; + double p; + double q; + double r; + double t; + double y; + double w; + double d; + double e; + double s; + + negIntA = isNonPositiveInteger( a ); + negIntB = isNonPositiveInteger( b ); + err = 0.0; + err1 = 0.0; + s = 1.0 - x; + + if ( x < -0.5 && !( negIntA || negIntB ) ) { + if ( b > a ) { + // Transformation based on AMS55 #15.3.4: + y = stdlib_base_pow( s, -a ) * hys2f1( a, c-b, c, -x/s, &err ); + } else { + // Transformation based on AMS55 #15.3.5: + y = stdlib_base_pow( s, -b ) * hys2f1( c-a, b, c, -x/s, &err ); + } + *loss = err; + return y; + } + + d = c - a - b; + id = stdlib_base_round( d ); + + if ( x > 0.9 && !negIntA && !negIntB ) { + if ( isInteger( d ) == false ) { + // Try the power series first: + y = hys2f1( a, b, c, x, &err ); + if ( err < ETHRESH ) { + *loss = err; + return y; + } + // If the power series fails, then apply AMS55 #15.3.6... + q = hys2f1( a, b, 1.0-d, s, &err ); + sign = 1.0; + w = stdlib_base_gammaln( d ); + sign *= stdlib_base_gammasgn( d ); + w -= stdlib_base_gammaln( c-a ); + sign *= stdlib_base_gammasgn( c-a ); + w -= stdlib_base_gammaln( c-b ); + sign *= stdlib_base_gammasgn( c-b ); + q *= sign * stdlib_base_exp( w ); + + r = stdlib_base_pow( s, d ) * hys2f1( c-a, c-b, d+1.0, s, &err1 ); + sign = 1.0; + w = stdlib_base_gammaln( -d ); + sign *= stdlib_base_gammasgn( -d ); + w -= stdlib_base_gammaln( a ); + sign *= stdlib_base_gammasgn( a ); + w -= stdlib_base_gammaln( b ); + sign *= stdlib_base_gammasgn( b ); + r *= sign * stdlib_base_exp( w ); + + y = q + r; + q = stdlib_base_abs( q ); + r = stdlib_base_abs( r ); + if ( q > r ) { + r = q; + } + err += err1 + ( ( MACHEP*r ) / y ); + y *= stdlib_base_gamma( c ); + } else { + // Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12... + if ( id >= 0.0 ) { + e = d; + d1 = d; + d2 = 0.0; + aid = id; + } else { + e = -d; + d1 = 0.0; + d2 = d; + aid = -id; + } + ax = stdlib_base_ln( s ); + + // eslint-disable-next-line max-len + y = stdlib_base_digamma( 1.0 ) + stdlib_base_digamma( 1.0+e ) - stdlib_base_digamma( a+d1 ) - stdlib_base_digamma( b+d1 ) - ax; + y /= stdlib_base_gamma( e+1.0 ); + p = ( a+d1 ) * ( b+d1 ) * s / stdlib_base_gamma( e+2.0 ); + t = 1.0; + do { + // eslint-disable-next-line max-len + r = stdlib_base_digamma( 1.0+t ) + stdlib_base_digamma( 1.0+t+e ) - stdlib_base_digamma( a+t+d1 ) - stdlib_base_digamma( b+t+d1 ) - ax; + q = p * r; + y += q; + p *= s * ( a+t+d1 ) / ( t+1.0 ); + p *= ( b+t+d1 ) / ( t+1.0+e ); + t += 1.0; + if ( t > MAX_ITERATIONS ) { + *loss = 1.0; + return STDLIB_CONSTANT_FLOAT64_NAN; + } + } while ( y == 0.0 || stdlib_base_abs( q/y ) > EPS ); + + if ( id == 0.0 ) { + y *= stdlib_base_gamma( c ) / ( stdlib_base_gamma( a ) * stdlib_base_gamma( b ) ); + *loss = err; + return y; + } + + y1 = 1.0; + if ( aid != 1.0 ) { + t = 0.0; + p = 1.0; + for ( i = 1; i < aid; i++ ) { + r = 1.0 - e + t; + p *= s * ( a+t+d2 ) * ( b+t+d2 ) / r; + t += 1.0; + p /= t; + y1 += p; + } + } + + p = stdlib_base_gamma( c ); + y1 *= stdlib_base_gamma( e ) * p / ( stdlib_base_gamma( a+d1 ) * stdlib_base_gamma( b+d1 ) ); + y *= p / ( stdlib_base_gamma( a+d2 ) * stdlib_base_gamma( b+d2 ) ); + if ( ( (int)aid & 1 ) != 0 ) { + y = -y; + } + q = stdlib_base_pow( s, id ); + y = ( id > 0.0 ) ? y*q : y1*q; + y += y1; + } + *loss = err; + return y; + } + // Perform power series if no special cases: + y = hys2f1( a, b, c, x, &err ); + *loss = err; + return y; +} + +/** +* Computes the regularized incomplete gamma function. The upper tail is calculated via the modified Lentz's method for computing continued fractions, the lower tail using a power expansion. +* +* ## Notes +* +* - When `a >= FLOAT64_MAX_NTH_FACTORIAL` and computing the non-regularized incomplete gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result. +* +* @param x input value +* @param a input value +* @param regularized input value +* @param upper input value +* @return function value +* +* @example +* double v = stdlib_base_gamma( 1.0, 1.0, 1.0, 0.0 ); +* // returns 1.0 +*/ +double stdlib_base_gammainc( const double x, const double a, const bool regularized, const bool upper ) { + double optimisedInvert; + double evalMethod; + double initValue; + double isHalfInt; + double useTemme; + double isSmallA; + double invert; + double result; + double isInt; + double sigma; + double gam; + double res; + double fa; + double g; + + if ( x < 0.0 || a <= 0.0 ) { + return 0.0 / 0.0; // NaN + } + invert = upper; + result = 0.0; + if ( a >= STDLIB_CONSTANT_FLOAT64_MAX_NTH_FACTORIAL && !regularized ) { + if ( invert && ( a * 4.0 < x ) ) { + // This is method 4 below, done in logs: + result = ( a * stdlib_base_ln(x) ) - x; + result += ln( upperGammaFraction( a, x ) ); + } else if ( !invert && ( a > 4.0 * x ) ) { + // This is method 2 below, done in logs: + result = ( a * ln(x) ) - x; + initValue = 0; + result += ln( lowerGammaSeries( a, x, initValue ) / a ); + } else { + result = igammaFinal( a, x, true, invert ); + if ( result == 0.0 ) { + if ( invert ) { + // Try http://functions.wolfram.com/06.06.06.0039.01: + result = 1.0 + ( 1.0 / (12.0*a) ) + ( 1.0 / (288.0*a*a) ); + result = ln( result ) - a + ( ( a-0.5 ) * ln(a) ); + result += ln( STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); + } else { + // This is method 2 below, done in logs, we're really outside the range of this method, but since the result is almost certainly infinite, we should probably be OK: + result = ( a * ln( x ) ) - x; + initValue = 0.0; + result += ln( lowerGammaSeries( a, x, initValue ) / a ); + } + } else { + result = ln( result ) + gammaln( a ); + } + } + if ( result > STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + return exp( result ); + } + // If no special handling is required then we proceeds as normal: + return igammaFinal( a, x, regularized, invert ); + + // isSmallA = ( a < 30 ) && ( a <= x + 1.0 ) && ( x < MAX_LN ); + // if ( isSmallA ) { + // fa = floor( a ); + // isInt = ( fa == a ); + // isHalfInt = ( isInt ) ? false : ( abs( fa - a ) == 0.5 ); + // } else { + // isInt = isHalfInt = false; + // } + // if ( isInt && x > 0.6 ) { + // // Calculate Q via finite sum: + // invert = !invert; + // evalMethod = 0; + // } else if ( isHalfInt && x > 0.2 ) { + // // Calculate Q via finite sum for half integer a: + // invert = !invert; + // evalMethod = 1; + // } else if ( x < SQRT_EPSILON && a > 1.0 ) { + // evalMethod = 6; + // } else if ( x < 0.5 ) { + // // Changeover criterion chosen to give a changeover at Q ~ 0.33: + // if ( -0.4 / ln( x ) < a ) { + // evalMethod = 2; + // } else { + // evalMethod = 3; + // } + // } else if ( x < 1.1 ) { + // // Changeover here occurs when P ~ 0.75 or Q ~ 0.25: + // if ( x * 0.75 < a ) { + // evalMethod = 2; + // } else { + // evalMethod = 3; + // } + // } else { + // // Begin by testing whether we're in the "bad" zone where the result will be near 0.5 and the usual series and continued fractions are slow to converge: + // useTemme = false; + // if ( regularized && a > 20 ) { + // sigma = abs( (x-a)/a ); + // if ( a > 200 ) { + // // Limit chosen so that we use Temme's expansion only if the result would be larger than about 10^-6. Below that the regular series and continued fractions converge OK, and if we use Temme's method we get increasing errors from the dominant erfc term as it's (inexact) argument increases in magnitude. + // if ( 20 / a > sigma * sigma ) { + // useTemme = true; + // } + // } else if ( sigma < 0.4 ) { + // useTemme = true; + // } + // } + // if ( useTemme ) { + // evalMethod = 5; + // } + // // Regular case where the result will not be too close to 0.5: Changeover occurs at P ~ Q ~ 0.5. Note that series computation of P is about x2 faster than continued fraction calculation of Q, so try and use the CF only when really necessary, especially for small x. + // else if ( x - ( 1.0 / (3.0 * x) ) < a ) { + // evalMethod = 2; + // } else { + // evalMethod = 4; + // invert = !invert; + // } + // } + + // /* eslint-disable default-case */ + // switch ( evalMethod ) { + // case 0: + // result = finiteGammaQ( a, x ); + // if (regularized == false ) { + // result *= gamma( a ); + // } + // break; + // case 1: + // result = finiteHalfGammaQ( a, x ); + // if ( regularized == false ) { + // result *= gamma( a ); + // } + // break; + // case 2: + // // Compute P: + // result = ( regularized ) ? + // regularisedGammaPrefix( a, x ) : + // fullIGammaPrefix( a, x ); + // if ( result !== 0.0 ) { + // initValue = 0.0; + // optimisedInvert = false; + // if ( invert ) { + // initValue = ( regularized ) ? 1.0 : gamma( a ); + // if ( + // regularized || + // result >= 1.0 || + // FLOAT64_MAX * result > initValue + // ) { + // initValue /= result; + // if ( + // regularized || + // a < 1.0 || + // ( FLOAT64_MAX / a > initValue ) + // ) { + // initValue *= -a; + // optimisedInvert = true; + // } + // else { + // initValue = 0.0; + // } + // } + // else { + // initValue = 0.0; + // } + // } + // } + // result *= lowerGammaSeries( a, x, initValue ) / a; + // if ( optimisedInvert ) { + // invert = false; + // result = -result; + // } + // break; + // case 3: + // // Compute Q: + // invert = !invert; + // res = tgammaSmallUpperPart( a, x, invert ); + // result = res[ 0 ]; + // g = res[ 1 ]; + // invert = false; + // if ( regularized ) { + // result /= g; + // } + // break; + // case 4: + // // Compute Q: + // result = ( regularized ) ? + // regularisedGammaPrefix( a, x ) : + // fullIGammaPrefix( a, x ); + // if ( result !== 0 ) { + // result *= upperGammaFraction( a, x ); + // } + // break; + // case 5: + // result = igammaTemmeLarge( a, x ); + // if ( x >= a ) { + // invert = !invert; + // } + // break; + // case 6: + // // Since x is so small that P is necessarily very small too, use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/ + // result = ( regularized ) ? + // pow(x, a) / gamma( a + 1.0 ) : + // pow( x, a ) / a; + // result *= 1.0 - ( a * x / ( a + 1.0 ) ); + // break; + // } + // if ( regularized && result > 1.0 ) { + // result = 1.0; + // } + // if ( invert ) { + // gam = ( regularized ) ? 1.0 : gamma( a ); + // result = gam - result; + // } + // return result; +} From a225e46a4c3a11bef4c15f46ead4344264bd6427 Mon Sep 17 00:00:00 2001 From: Karan Anand Date: Wed, 18 Jun 2025 15:57:22 -0700 Subject: [PATCH 02/14] feat: temp commit --- .../base/special/gammainc/scripts/evalpoly.js | 79 + .../math/base/special/gammainc/src/main.c | 1389 +++++++++++------ 2 files changed, 964 insertions(+), 504 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/scripts/evalpoly.js b/lib/node_modules/@stdlib/math/base/special/gammainc/scripts/evalpoly.js index c262caea673b..41613da2042b 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammainc/scripts/evalpoly.js +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/scripts/evalpoly.js @@ -24,10 +24,15 @@ // MODULES // var resolve = require( 'path' ).resolve; +var readFileSync = require( '@stdlib/fs/read-file' ).sync; var writeFileSync = require( '@stdlib/fs/write-file' ).sync; var currentYear = require( '@stdlib/time/current-year' ); var licenseHeader = require( '@stdlib/_tools/licenses/header' ); var compile = require( '@stdlib/math/base/tools/evalpoly-compile' ); +var compileC = require( '@stdlib/math/base/tools/evalpoly-compile-c' ); +var substringBefore = require( '@stdlib/string/substring-before' ); +var substringAfter = require( '@stdlib/string/substring-after' ); +var format = require( '@stdlib/string/format' ); // VARIABLES // @@ -139,6 +144,33 @@ var header = licenseHeader( 'Apache-2.0', 'js', { header += '\n/* This is a generated file. Do not edit directly. */\n'; +// FUNCTIONS // + +/** +* Inserts a compiled function into file content. +* +* @private +* @param {string} text - source content +* @param {string} id - function identifier +* @param {string} str - function string +* @returns {string} updated content +*/ +function insert( text, id, str ) { + var before; + var after; + var begin; + var end; + + begin = '// BEGIN: '+id; + end = '// END: '+id; + + before = substringBefore( text, begin ); + after = substringAfter( text, end ); + + return format( '%s// BEGIN: %s\n\n%s\n%s%s', before, id, str, end, after ); +} + + // MAIN // /** @@ -148,7 +180,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; */ function main() { var fpath; + var copts; var opts; + var file; var str; opts = { @@ -190,6 +224,51 @@ function main() { fpath = resolve( __dirname, '..', 'lib', 'polyval_c8.js' ); str = header + compile( C8 ); writeFileSync( fpath, str, opts ); + + copts = { + 'dtype': 'double', + 'name': '' + }; + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + copts.name = 'polyval_c0'; + str = compileC( C0, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c1'; + str = compileC( C1, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c2'; + str = compileC( C2, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c3'; + str = compileC( C3, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c4'; + str = compileC( C4, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c5'; + str = compileC( C5, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c6'; + str = compileC( C6, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c7'; + str = compileC( C7, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_c8'; + str = compileC( C8, copts ); + file = insert( file, copts.name, str ); + + writeFileSync( fpath, file, opts ); } main(); diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c index b5553ca9bc82..c54547b0abe2 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c @@ -18,7 +18,7 @@ * * ## Notice * -* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified according to stdlib conventions. +* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/stdlib_base_gamma.hpp}. The implementation has been modified according to stdlib conventions. * * ```text * (C) Copyright John Maddock 2006-7, 2013-14. @@ -33,428 +33,974 @@ */ #include "stdlib/math/base/special/gammainc.h" -#include "stdlib/math/base/special/gammaln.h" -#include "stdlib/math/base/special/floor.h" -#include "stdlib/math/base/special/gamma.h" -#include "stdlib/math/base/special/abs.h" -#include "stdlib/math/base/special/exp.h" -#include "stdlib/math/base/special/pow.h" -#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/special/stdlib_base_gammaln.h" +#include "stdlib/math/base/special/stdlib_base_floor.h" +#include "stdlib/math/base/special/stdlib_base_gamma.h" +#include "stdlib/math/base/special/stdlib_base_abs.h" +#include "stdlib/math/base/special/stdlib_base_exp.h" +#include "stdlib/math/base/special/stdlib_base_pow.h" +#include "stdlib/math/base/special/stdlib_base_ln.h" #include "stdlib/math/base/assert/is_nan.h" #include "stdlib/constants/float64/sqrt_eps.h" #include "stdlib/constants/float64/max.h" #include "stdlib/constants/float64/sqrt_two_pi.h" #include "stdlib/constants/float64/max_ln.h" #include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float64/eps.h" +#include "stdlib/constants/float32/smallest_normal.h" #include "stdlib/constants/float64/max_nth_factorial.h" +#include #include -static const double MACHEP = 1.11022302462515654042e-16; -static const double EPS = 1.0e-13; -static const double ETHRESH = 1.0e-12; -static const double MAX_ITERATIONS = 10000.0; +static const int32_t MAX_ITERATIONS = 1000000; static double hys2f1( double a, double b, const double c, const double x, double *loss ); +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_c0 + /** -* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* Evaluates a polynomial. * -* @param a function parameter -* @param z function parameter -* @return function value +* ## 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_c0( const double x ) { + return -0.3333333333333333 + (x * (0.08333333333333333 + (x * (-0.014814814814814815 + (x * (0.0011574074074074073 + (x * (0.0003527336860670194 + (x * (-0.0001787551440329218 + (x * (0.00003919263178522438 + (x * (-0.0000021854485106799924 + (x * (-0.00000185406221071516 + (x * (8.296711340953087e-7 + (x * (-1.7665952736826078e-7 + (x * (6.707853543401498e-9 + (x * (1.0261809784240309e-8 + (x * (-4.382036018453353e-9 + (x * 9.14769958223679e-10))))))))))))))))))))))))))); +} + +// END: polyval_c0// BEGIN: polyval_c1 + +/** +* 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_c1( const double x ) { + return -0.001851851851851852 + (x * (-0.003472222222222222 + (x * (0.0026455026455026454 + (x * (-0.0009902263374485596 + (x * (0.00020576131687242798 + (x * (-4.018775720164609e-7 + (x * (-0.000018098550334489977 + (x * (0.00000764916091608111 + (x * (-0.0000016120900894563446 + (x * (4.647127802807434e-9 + (x * (1.378633446915721e-7 + (x * (-5.752545603517705e-8 + (x * 1.1951628599778148e-8))))))))))))))))))))))); +} + +// END: polyval_c1// BEGIN: polyval_c2 + +/** +* 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_c2( const double x ) { + return 0.004133597883597883 + (x * (-0.0026813271604938273 + (x * (0.0007716049382716049 + (x * (0.0000020093878600823047 + (x * (-0.00010736653226365161 + (x * (0.000052923448829120125 + (x * (-0.000012760635188618728 + (x * (3.423578734096138e-8 + (x * (0.0000013721957309062932 + (x * (-6.298992138380055e-7 + (x * 1.4280614206064242e-7))))))))))))))))))); +} + +// END: polyval_c2// BEGIN: polyval_c3 + +/** +* 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_c3( const double x ) { + return 0.0006494341563786008 + (x * (0.00022947209362139917 + (x * (-0.0004691894943952557 + (x * (0.00026772063206283885 + (x * (-0.00007561801671883977 + (x * (-2.396505113867297e-7 + (x * (0.000011082654115347302 + (x * (-0.0000056749528269915965 + (x * 0.0000014230900732435883))))))))))))))); +} + +// END: polyval_c3// BEGIN: polyval_c4 + +/** +* 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_c4( const double x ) { + return -0.0008618882909167117 + (x * (0.0007840392217200666 + (x * (-0.0002990724803031902 + (x * (-0.0000014638452578843418 + (x * (0.00006641498215465122 + (x * (-0.00003968365047179435 + (x * 0.000011375726970678419))))))))))); +} + +// END: polyval_c4// BEGIN: polyval_c5 + +/** +* 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_c5( const double x ) { + return -0.00033679855336635813 + (x * (-0.00006972813758365858 + (x * (0.0002772753244959392 + (x * (-0.00019932570516188847 + (x * (0.00006797780477937208 + (x * (1.419062920643967e-7 + (x * (-0.000013594048189768693 + (x * (0.000008018470256334202 + (x * -0.000002291481176508095))))))))))))))); +} + +// END: polyval_c5// BEGIN: polyval_c6 + +/** +* 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_c6( const double x ) { + return 0.0005313079364639922 + (x * (-0.0005921664373536939 + (x * (0.0002708782096718045 + (x * (7.902353232660328e-7 + (x * (-0.00008153969367561969 + (x * (0.0000561168275310625 + (x * -0.000018329116582843375))))))))))); +} + +// END: polyval_c6// BEGIN: polyval_c7 + +/** +* 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_c7( const double x ) { + return 0.00034436760689237765 + (x * (0.00005171790908260592 + (x * (-0.00033493161081142234 + (x * (0.0002812695154763237 + (x * -0.00010976582244684731))))))); +} + +// END: polyval_c7// BEGIN: polyval_c8 + +/** +* 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 bool upperGammaFraction( const double a, const double z ) { - double f = upperIncompleteGammaFract( a, z ); - return 1.0 / ( z - a + 1.0 + continuedFractionA( f ) ); +static double polyval_c8( const double x ) { + return -0.0006526239185953094 + (x * (0.0008394987206720873 + (x * -0.000438297098541721))); } +// END: polyval_c8 + +/* End auto-generated functions. */ + /** -* Creates a function to evaluate a series expansion of the upper incomplete gamma fraction. +* Creates a function to evaluate a series expansion of the upper incomplete stdlib_base_gamma fraction. * * @param a1 function parameter * @param z1 function parameter * @return series function */ -static bool upperIncompleteGammaFract( const double a1, const double z1 ) { - var z = z1 - a1 + 1.0; - var a = a1; - var k = 0; - return next; - - /** - * Calculate the next term of the series. - * - * @private - * @returns {Array} series expansion terms - */ - function next() { - k += 1; - z += 2.0; - return [ - k * (a - k), - z - ]; +static void upperIncompleteGammaFract( const double a, double* z, double* k, double *out ) { + *z += 2.0; + *k += 1.0; + out[ 0 ] = ( *k ) * ( a - (*k) ); + out[ 1 ] = *z; + return; +} + +/** +* Evaluates a continued fraction expansion. +* +* ```text +* a1 +* --------------- +* b1 + a2 +* ---------- +* b2 + a3 +* ----- +* b3 + ... +* ``` +* +* @private +* @param {Function} gen - function giving terms of continued fraction expansion +* @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term +* @param {PositiveInteger} maxIter - maximum number of iterations +* @returns {number} evaluated expansion +*/ +static double continuedFractionA( const double a, const double z, const double factor, const int32_t maxIter ) { + int32_t maxIterC; + double out[ 2 ]; + double delta; + double zc; + double a0; + double C; + double D; + double f; + double k; + + zc = z; + k = 0.0; + upperIncompleteGammaFract( a, &zc, &k, out ); + + f = out[ 1 ]; + a0 = out[ 0 ]; + if ( f == 0.0 ) { + f = STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL; } + C = f; + D = 0.0; + maxIterC = maxIter; + + do { + upperIncompleteGammaFract( a, &zc, &k, out ); + D = out[ 1 ] + ( out[ 0 ] * D ); + if ( D == 0.0 ) { + D = STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL; + } + C = out[ 1 ] + ( out[ 0 ] / C ); + if ( C == 0.0 ) { + C = STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL; + } + D = 1.0 / D; + delta = C * D; + f *= delta; + } while ( out && ( stdlib_base_abs( delta - 1.0 ) > factor ) && --maxIterC ); + + return a0 / f; } /** -* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* Evaluate the lower incomplete stdlib_base_gamma integral via a series expansion and divide by `stdlib_base_gamma(z)` to normalize. * * @param a function parameter * @param z function parameter * @return function value */ -static bool upperGammaFraction( const double a, const double z ) { - double f = upperIncompleteGammaFract( a, z ); - return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +static double upperGammaFraction( const double a, const double z ) { + double f; + + f = z - a + 1.0; + return 1.0 / ( f + continuedFractionA( a, f, STDLIB_CONSTANT_FLOAT64_EPS, MAX_ITERATIONS ) ); } /** -* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* Creates a function to evaluate a series expansion of the upper incomplete stdlib_base_gamma fraction. * -* @param a function parameter -* @param z function parameter -* @return function value +* @param a1 function parameter +* @param z1 function parameter +* @return series function +*/ +static double lowerIncompleteGammaSeries( double* a, const double z, double* result ) { + double r; + + r = *result; + *a += 1.0; + *result *= z / ( *a ); + return r; +} + +/** +* Evaluates a continued fraction expansion. +* +* ```text +* a1 +* --------------- +* b1 + a2 +* ---------- +* b2 + a3 +* ----- +* b3 + ... +* ``` +* +* @private +* @param {Function} gen - function giving terms of continued fraction expansion +* @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term +* @param {PositiveInteger} maxIter - maximum number of iterations +* @returns {number} evaluated expansion */ -static bool upperGammaFraction( const double a, const double z ) { - double f = upperIncompleteGammaFract( a, z ); - return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +static double sumSeries( const double a, const double z, const double initialValue ) { + int32_t counter; + double nextTerm; + double result; + double sum; + double ac; + + ac = a; + result = 1.0; + sum = initialValue; + counter = MAX_ITERATIONS; + + // Repeatedly call function... + do { + nextTerm = lowerIncompleteGammaSeries( &ac, z, &result ); + sum += nextTerm; + } + while ( ( stdlib_base_abs(STDLIB_CONSTANT_FLOAT64_EPS * sum) < stdlib_base_abs(nextTerm) ) && --counter ); + + return sum; } /** -* Evaluate the lower incomplete gamma integral via a series expansion and divide by `gamma(z)` to normalize. +* Evaluate the lower incomplete stdlib_base_gamma integral via a series expansion and divide by `stdlib_base_gamma(z)` to normalize. * * @param a function parameter * @param z function parameter * @return function value */ -static bool upperGammaFraction( const double a, const double z ) { - double f = upperIncompleteGammaFract( a, z ); - return 1.0 / ( z - a + 1.0 + continuedFraction( f ) ); +static double lowerGammaSeries( const double a, const double z, const double initValue ) { + return sumSeries( a, z, initValue ); } /** -* Evaluates the Gaussian hypergeometric function by two-term recurrence in `a`. -* -* ## Notes +* Calculates normalized Q when a is an integer. * -* - This function helps prevent losing accuracy in the highly alternating hypergeometric series and allows a and b to be reduced to smaller values. -* -* ## References -* -* - AMS55 #15.2.10 -* -* @param a input value -* @param b input value -* @param c input value -* @param x input value -* @param loss starting loss of significance -* @return function value +* @private +* @param {integer} a - function parameter +* @param {number} x - function parameter +* @returns {number} upper stdlib_base_gamma fraction */ -static double hyp2f1ra( const double a, const double b, const double c, const double x, double *loss ) { - double err; - double da; - double f2; - double f1; - double f0; - double t; +static double finiteGammaQ( const double a, const double x ) { + double term; + double sum; + double e; double n; - *loss = 0.0; - err = 0.0; + e = stdlib_base_exp( -x ); + sum = e; + if ( sum != 0.0 ) { + term = sum; + for ( n = 1; n < a; ++n ) { + term /= n; + term *= x; + sum += term; + } + } + return sum; +} + +/** +* Calculates normalized Q when a is a half-integer. +* +* @private +* @param {number} a - function parameter +* @param {number} x - function parameter +* @returns {number} upper stdlib_base_gamma fraction +*/ +static double finiteHalfGammaQ( const double a, const double x ) { + double half; + double term; + double sum; + double e; + double n; - if ( ( c < 0.0 && a <= c ) || ( c >= 0.0 && a >= c ) ) { - da = stdlib_base_round( a-c ); - } else { - da = stdlib_base_round( a ); + e = stdlib_base_erfc( stdlib_base_sqrt(x) ); + if ( e != 0 && a > 1.0 ) { + term = stdlib_base_exp( -x ) / stdlib_base_sqrt( PI * x ); + term *= x; + half = 0.5; + term /= half; + sum = term; + for ( n = 2; n < a; ++n ) { + term /= n - half; + term *= x; + sum += term; + } + e += sum; } + return e; +} + +/** +* Computes `(z^a)*(e^-z) / stdlib_base_gamma(a)`. +* +* @private +* @param {number} a - input value +* @param {number} z - input value +* @returns {number} function value +*/ +static double regularisedGammaPrefix( const double a, const double z ) { + // TODO: NEEDS TO BE CHANGED (to be continued here) + double prefix; + double amza; + double agh; + double alz; + double amz; + double sq; + double d; - t = a - da; - if ( stdlib_base_abs( da ) > MAX_ITERATIONS ) { - // Too expensive to compute, so return `NaN`: - *loss = 1.0; - return STDLIB_CONSTANT_FLOAT64_NAN; + agh = a + G - 0.5; + d = ( (z - a) - G + 0.5 ) / agh; + if ( a < 1.0 ) { + // Treat a < 1 as a special case because our Lanczos approximations are optimized against the factorials with a > 1, and for high precision types very small values of `a` can give rather erroneous results for stdlib_base_gamma: + if ( z <= STDLIB_CONSTANT_FLOAT64_MIN_LN ) { + // Use logs, so should be free of cancellation errors: + return stdlib_base_exp( ( a * stdlib_base_ln(z) ) - z - stdlib_base_gammaln( a ) ); + } + // No danger of overflow as stdlib_base_gamma(a) < 1/a for small a, so direct calculation: + return stdlib_base_pow( z, a ) * stdlib_base_exp( -z ) / stdlib_base_gamma( a ); + } + if ( stdlib_base_abs(d*d*a) <= 100.0 && a > 150.0 ) { + // Special case for large a and a ~ z: + prefix = ( a * ( log1p( d ) - d ) ) + ( z * ( 0.5-G ) / agh ); + prefix = stdlib_base_exp( prefix ); } + else { + // General case. Direct computation is most accurate, but use various fallbacks for different parts of the problem domain: + alz = a * stdlib_base_ln(z / agh); + amz = a - z; + if ( + min(alz, amz) <= STDLIB_CONSTANT_FLOAT64_MIN_LN || + max(alz, amz) >= STDLIB_CONSTANT_FLOAT64_MAX_LN + ) { + amza = amz / a; + if ( + min(alz, amz)/2.0 > STDLIB_CONSTANT_FLOAT64_MIN_LN && + max(alz, amz)/2.0 < STDLIB_CONSTANT_FLOAT64_MAX_LN + ) { + // Compute square root of the result and then square it: + sq = stdlib_base_pow( z / agh, a / 2.0 ) * stdlib_base_exp( amz / 2.0 ); + prefix = sq * sq; + } + else if ( + min(alz, amz)/4.0 > STDLIB_CONSTANT_FLOAT64_MIN_LN && + max(alz, amz)/4.0 < STDLIB_CONSTANT_FLOAT64_MAX_LN && + z > a + ) { + // Compute the 4th root of the result then square it twice: + sq = stdlib_base_pow( z / agh, a / 4.0 ) * stdlib_base_exp( amz / 4.0 ); + prefix = sq * sq; + prefix *= prefix; + } + else if ( + amza > STDLIB_CONSTANT_FLOAT64_MIN_LN && + amza < STDLIB_CONSTANT_FLOAT64_MAX_LN + ) { + prefix = stdlib_base_pow( (z * stdlib_base_exp(amza)) / agh, a ); + } + else { + prefix = stdlib_base_exp( alz + amz ); + } + } + else + { + prefix = stdlib_base_pow( z / agh, a ) * stdlib_base_exp( amz ); + } + } + prefix *= stdlib_base_sqrt( agh / E ) / lanczosSumExpGScaled( a ); + return prefix; +} + +/** +* Calculates the power term prefix `(z^a)(e^-z)` used in the non-normalized incomplete gammas. +* +* @private +* @param {number} a - function parameter +* @param {number} z - function parameter +* @returns {number} power term prefix +*/ +static double fullIGammaPrefix( a, z ) { + double prefix; + double alz; - if ( da < 0.0 ) { - // Backward recurrence... - f1 = hys2f1( t, b, c, x, &err ); - *loss += err; - f0 = hys2f1( t-1.0, b, c, x, &err ); - *loss += err; - t -= 1.0; - for ( n = 1; n < -da; n++ ) { - f2 = f1; - f1 = f0; - - // eslint-disable-next-line max-len - f0 = -( ( ( ( (2.0*t)-c )-( t*x )+( b*x ) ) * f1 ) + ( ( t*(x-1.0) ) * f2 ) ) / ( c-t ); - t -= 1.0; + alz = a * stdlib_base_ln( z ); + if ( z >= 1.0 ) { + if ( ( alz < STDLIB_CONSTANT_FLOAT64_MAX_LN ) && ( -z > STDLIB_CONSTANT_FLOAT64_MIN_LN) ) { + prefix = stdlib_base_pow( z, a ) * stdlib_base_exp( -z ); + } else if ( a >= 1.0 ) { + prefix = stdlib_base_pow( z / stdlib_base_exp(z/a), a ); } - } else { - // Forward recurrence... - f1 = hys2f1( t, b, c, x, &err ); - *loss += err; - f0 = hys2f1( t+1.0, b, c, x, &err ); - *loss += err; - t += 1.0; - for ( n = 1; n < da; n++ ) { - f2 = f1; - f1 = f0; - - // eslint-disable-next-line max-len - f0 = -( ( ( ( (2.0*t)-c )-( t*x )+( b*x ) ) * f1 ) + ( (c-t)*f2 ) ) / ( t*(x-1.0) ); - t += 1.0; + else { + prefix = stdlib_base_exp( alz - z ); } } - return f0; + else { + if ( alz > STDLIB_CONSTANT_FLOAT64_MIN_LN ) { + prefix = stdlib_base_pow( z, a ) * stdlib_base_exp( -z ); + } else if ( z/a < STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + prefix = stdlib_base_pow( z / stdlib_base_exp(z/a), a ); + } else { + prefix = stdlib_base_exp( alz - z ); + } + } + return prefix; } /** -* Evaluates the power series expansion of Gaussian hypergeometric function. -* -* @param a input value -* @param b input value -* @param c input value -* @param x input value -* @param loss starting loss of significance -* @return function value +* Creates a function to evaluate a series expansion of the upper incomplete stdlib_base_gamma fraction. +* +* @param a1 function parameter +* @param z1 function parameter +* @return series function */ -static double hys2f1( double a, double b, const double c, const double x, double *loss ) { - double intFlag; - double umax; - double f; - double g; - double h; - double k; - double m; - double s; - double u; - double i; +static double smallGamma2Series( double* a, double* x, int32_t* n, double* result ) { + double r; - intFlag = 0.0; + r = ( *result ) / ( *a ); + *result *= ( *x ); + *n += 1; + *result /= ( *n ); + *a += 1.0; + return r; +} - if ( stdlib_base_abs( b ) > stdlib_base_abs( a ) ) { - f = b; - b = a; - a = f; - } +/** +* Evaluates a continued fraction expansion. +* +* ```text +* a1 +* --------------- +* b1 + a2 +* ---------- +* b2 + a3 +* ----- +* b3 + ... +* ``` +* +* @private +* @param {Function} gen - function giving terms of continued fraction expansion +* @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term +* @param {PositiveInteger} maxIter - maximum number of iterations +* @returns {number} evaluated expansion +*/ +static double sumSeries2( const double a, const double x, const double initialValue ) { + int32_t counter; + double nextTerm; + double result; + double sum; + int32_t n; + double ac; + double xc; + + xc = -x; + ac = a + 1.0; + result = -x; + n = 1; + sum = initialValue; + counter = MAX_ITERATIONS; - if ( isNonPositiveInteger( b ) && ( stdlib_base_abs( b ) < stdlib_base_abs( a ) ) ) { - f = b; - b = a; - a = f; - intFlag = 1.0; + // Repeatedly call function... + do { + nextTerm = smallGamma2Series( &ac, &xc, &n, &result ); + sum += nextTerm; } + while ( ( stdlib_base_abs(STDLIB_CONSTANT_FLOAT64_EPS * sum) < stdlib_base_abs(nextTerm) ) && --counter ); + + return sum; +} + +/** +* Evaluate the lower incomplete stdlib_base_gamma integral via a series expansion and divide by `stdlib_base_gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value +*/ +static void tgammaSmallUpperPart( const double a, const double x, const bool invert, double* out ) { + double initialValue; + double result; + double pgam; + double p; - // eslint-disable-next-line max-len - if ( ( ( stdlib_base_abs( a ) > stdlib_base_abs( c ) + 1.0 ) || intFlag ) && ( stdlib_base_abs( c-a ) > 2.0 ) && ( stdlib_base_abs( a ) > 2.0 ) ) { - return hyp2f1ra( a, b, c, x, loss ); + result = stdlib_base_gamma1pm1( a ); + pgam = ( result + 1.0 ) / a; + p = stdlib_base_powm1( x, a ); + result -= p; + result /= a; + p += 1.0; + initialValue = ( invert ) ? pgam : 0.0; + result = -p * sumSeries2( a, x, ( initialValue-result ) / p ); + if ( invert ) { + result = -result; } + out[ 0 ] = result; + out[ 1 ] = pgam; +} - i = 0.0; - umax = 0.0; - f = a; - g = b; - h = c; - s = 1.0; - u = 1.0; - k = 0.0; +/** +* Creates a function to evaluate a series expansion of the upper incomplete stdlib_base_gamma fraction. +* +* @param a1 function parameter +* @param z1 function parameter +* @return series function +*/ +static double tgammaILargeXSeries( double* a, const double z, double* result ) { + double r; + + r = *result; + *result *= ( *a ) / z; + *a -= 1.0; + return r; +} + +/** +* Evaluates a continued fraction expansion. +* +* ```text +* a1 +* --------------- +* b1 + a2 +* ---------- +* b2 + a3 +* ----- +* b3 + ... +* ``` +* +* @private +* @param {Function} gen - function giving terms of continued fraction expansion +* @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term +* @param {PositiveInteger} maxIter - maximum number of iterations +* @returns {number} evaluated expansion +*/ +static double sumSeries3( const double a, const double z ) { + int32_t counter; + double nextTerm; + double result; + double sum; + double ac; + ac = a; + sum = 0.0; + result = 1.0; + counter = MAX_ITERATIONS; + + // Repeatedly call function... do { - if ( stdlib_base_abs( h ) < EPS ) { - *loss = 1.0; - return STDLIB_CONSTANT_FLOAT64_PINF; - } - m = k + 1.0; - u *= ( ( f+k ) * ( g+k ) * x / ( ( h+k ) * m ) ); - s += u; - k = stdlib_base_abs( u ); - if ( k > umax ) { - umax = k; - } - k = m; - i += 1.0; - if ( i > MAX_ITERATIONS ) { - *loss = 1.0; - return s; - } - } while ( s == 0.0 || stdlib_base_abs( u/s ) > MACHEP ); + nextTerm = tgammaILargeXSeries( &ac, z, &result ); + sum += nextTerm; + } + while ( ( stdlib_base_abs(STDLIB_CONSTANT_FLOAT64_EPS * sum) < stdlib_base_abs(nextTerm) ) && --counter ); - // Estimate the relative error due to truncation by the series: - *loss = ( ( MACHEP*umax ) / stdlib_base_abs( s ) ) + ( MACHEP*i ); - return s; + return sum; } /** -* Applies transformations for `|x|` near unity before performing a power series expansion. -* -* @param a input value -* @param b input value -* @param c input value -* @param x input value -* @param loss starting loss of significance -* @return function value +* Evaluate the lower incomplete stdlib_base_gamma integral via a series expansion and divide by `stdlib_base_gamma(z)` to normalize. +* +* @param a function parameter +* @param z function parameter +* @return function value */ -static double hyt2f1( const double a, const double b, const double c, const double x, double *loss ) { - double negIntA; - double negIntB; - double sign; - double err1; - double err; - double aid; - double ax; - double id; - double d1; - double d2; - double y1; - double i; +static double tgammaILargeX( const double a, const double x ) { + return sumSeries3( a, x ); +} + +/** +* Evaluates a polynomial using double-precision floating-point arithmetic. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param c - polynomial coefficients sorted in ascending degree +* @param x - value at which to evaluate the polynomial +* @returns evaluated polynomial +*/ +static double evalpoly( double *c, const double x ){ + int32_t i; double p; - double q; - double r; - double t; + + i = 10; + if ( x == 0.0 ) { + return c[ 0 ]; + } + i -= 1; + p = ( c[ i ] * x ) + c[ i-1 ]; + i -= 2; + while ( i >= 0 ) { + p = ( p * x ) + c[ i ]; + i -= 1; + } + return p; +} + +/** +* Asymptotic expansions of the incomplete stdlib_base_gamma functions when `a` is large and `x ~ a` (IEEE double precision or 10^-17). +* +* @private +* @param {number} a - function parameter +* @param {number} x - function parameter +* @returns {number} value of asymptotic expansion +*/ +static double igammaTemmeLarge( const double a, const double x ) { + double workspace[ 10 ]; + double result; + double sigma; + double phi; double y; - double w; - double d; - double e; - double s; - - negIntA = isNonPositiveInteger( a ); - negIntB = isNonPositiveInteger( b ); - err = 0.0; - err1 = 0.0; - s = 1.0 - x; - - if ( x < -0.5 && !( negIntA || negIntB ) ) { - if ( b > a ) { - // Transformation based on AMS55 #15.3.4: - y = stdlib_base_pow( s, -a ) * hys2f1( a, c-b, c, -x/s, &err ); - } else { - // Transformation based on AMS55 #15.3.5: - y = stdlib_base_pow( s, -b ) * hys2f1( c-a, b, c, -x/s, &err ); - } - *loss = err; - return y; + double z; + + sigma = ( x-a ) / a; + phi = -stdlib_base_log1pmx( sigma ); + y = a * phi; + z = stdlib_base_sqrt( 2.0 * phi ); + if ( x < a ) { + z = -z; } + workspace[ 0 ] = polyvalC0( z ); + workspace[ 1 ] = polyvalC1( z ); + workspace[ 2 ] = polyvalC2( z ); + workspace[ 3 ] = polyvalC3( z ); + workspace[ 4 ] = polyvalC4( z ); + workspace[ 5 ] = polyvalC5( z ); + workspace[ 6 ] = polyvalC6( z ); + workspace[ 7 ] = polyvalC7( z ); + workspace[ 8 ] = polyvalC8( z ); + workspace[ 9 ] = -0.00059676129019274625; + result = evalpoly( workspace, 1.0/a ); + result *= stdlib_base_exp( -y ) / stdlib_base_sqrt( STDLIB_CONSTANT_FLOAT64_TWO_PI * a ); + if ( x < a ) { + result = -result; + } + result += stdlib_base_erfc( stdlib_base_sqrt(y) ) / 2.0; + return result; +} - d = c - a - b; - id = stdlib_base_round( d ); +/** +* Computes the regularized incomplete stdlib_base_gamma function. The upper tail is calculated via the modified Lentz's method for computing continued fractions, the lower tail using a power expansion. +* +* ## Notes +* +* - When `a >= FLOAT64_MAX_NTH_FACTORIAL` and computing the non-regularized incomplete stdlib_base_gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized stdlib_base_gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result. +* +* @param x input value +* @param a input value +* @param regularized input value +* @param upper input value +* @return function value +* +* @example +* double v = stdlib_base_gammainc( 1.0, 1.0, 1.0, 0.0 ); +* // returns 1.0 +*/ +double igammaFinal( const double a, const double x, const bool regularized, const bool upper ) { + bool optimisedInvert; + int32_t evalMethod; + double initValue; + double out[ 2 ]; + bool isHalfInt; + bool useTemme; + bool isSmallA; + double result; + double sigma; + bool invert; + bool isInt; + double gam; + double fa; + double g; - if ( x > 0.9 && !negIntA && !negIntB ) { - if ( isInteger( d ) == false ) { - // Try the power series first: - y = hys2f1( a, b, c, x, &err ); - if ( err < ETHRESH ) { - *loss = err; - return y; - } - // If the power series fails, then apply AMS55 #15.3.6... - q = hys2f1( a, b, 1.0-d, s, &err ); - sign = 1.0; - w = stdlib_base_gammaln( d ); - sign *= stdlib_base_gammasgn( d ); - w -= stdlib_base_gammaln( c-a ); - sign *= stdlib_base_gammasgn( c-a ); - w -= stdlib_base_gammaln( c-b ); - sign *= stdlib_base_gammasgn( c-b ); - q *= sign * stdlib_base_exp( w ); - - r = stdlib_base_pow( s, d ) * hys2f1( c-a, c-b, d+1.0, s, &err1 ); - sign = 1.0; - w = stdlib_base_gammaln( -d ); - sign *= stdlib_base_gammasgn( -d ); - w -= stdlib_base_gammaln( a ); - sign *= stdlib_base_gammasgn( a ); - w -= stdlib_base_gammaln( b ); - sign *= stdlib_base_gammasgn( b ); - r *= sign * stdlib_base_exp( w ); - - y = q + r; - q = stdlib_base_abs( q ); - r = stdlib_base_abs( r ); - if ( q > r ) { - r = q; - } - err += err1 + ( ( MACHEP*r ) / y ); - y *= stdlib_base_gamma( c ); + result = 0.0; + invert = upper; + isSmallA = ( a < 30 ) && ( a <= x + 1.0 ) && ( x < STDLIB_CONSTANT_FLOAT64_MAX_LN ); + if ( isSmallA ) { + fa = stdlib_base_floor( a ); + isInt = ( fa == a ); + isHalfInt = ( isInt ) ? false : ( stdlib_base_abs( fa - a ) == 0.5 ); + } else { + isInt = false; + isHalfInt = false; + } + if ( isInt && x > 0.6 ) { + // Calculate Q via finite sum: + invert = !invert; + evalMethod = 0; + } else if ( isHalfInt && x > 0.2 ) { + // Calculate Q via finite sum for half integer a: + invert = !invert; + evalMethod = 1; + } else if ( x < STDLIB_CONSTANT_FLOAT64_SQRT_EPS && a > 1.0 ) { + evalMethod = 6; + } else if ( ( x > 1000.0 ) && ( ( a < x ) || ( stdlib_base_abs( a - 50.0 ) / x < 1.0 ) ) ) { + // calculate Q via asymptotic approximation: + invert = !invert; + evalMethod = 7; + } else if ( x < 0.5 ) { + // Changeover criterion chosen to give a changeover at Q ~ 0.33: + if ( -0.4 / stdlib_base_ln( x ) < a ) { + evalMethod = 2; } else { - // Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12... - if ( id >= 0.0 ) { - e = d; - d1 = d; - d2 = 0.0; - aid = id; - } else { - e = -d; - d1 = 0.0; - d2 = d; - aid = -id; - } - ax = stdlib_base_ln( s ); - - // eslint-disable-next-line max-len - y = stdlib_base_digamma( 1.0 ) + stdlib_base_digamma( 1.0+e ) - stdlib_base_digamma( a+d1 ) - stdlib_base_digamma( b+d1 ) - ax; - y /= stdlib_base_gamma( e+1.0 ); - p = ( a+d1 ) * ( b+d1 ) * s / stdlib_base_gamma( e+2.0 ); - t = 1.0; - do { - // eslint-disable-next-line max-len - r = stdlib_base_digamma( 1.0+t ) + stdlib_base_digamma( 1.0+t+e ) - stdlib_base_digamma( a+t+d1 ) - stdlib_base_digamma( b+t+d1 ) - ax; - q = p * r; - y += q; - p *= s * ( a+t+d1 ) / ( t+1.0 ); - p *= ( b+t+d1 ) / ( t+1.0+e ); - t += 1.0; - if ( t > MAX_ITERATIONS ) { - *loss = 1.0; - return STDLIB_CONSTANT_FLOAT64_NAN; + evalMethod = 3; + } + } else if ( x < 1.1 ) { + // Changeover here occurs when P ~ 0.75 or Q ~ 0.25: + if ( x * 0.75 < a ) { + evalMethod = 2; + } else { + evalMethod = 3; + } + } else { + // Begin by testing whether we're in the "bad" zone where the result will be near 0.5 and the usual series and continued fractions are slow to converge: + useTemme = false; + if ( regularized && a > 20.0 ) { + sigma = stdlib_base_abs( (x-a)/a ); + if ( a > 200 ) { + // Limit chosen so that we use Temme's expansion only if the result would be larger than about 10^-6. Below that the regular series and continued fractions converge OK, and if we use Temme's method we get increasing errors from the dominant stdlib_base_erfc term as it's (inexact) argument increases in magnitude. + if ( 20 / a > sigma * sigma ) { + useTemme = true; } - } while ( y == 0.0 || stdlib_base_abs( q/y ) > EPS ); - - if ( id == 0.0 ) { - y *= stdlib_base_gamma( c ) / ( stdlib_base_gamma( a ) * stdlib_base_gamma( b ) ); - *loss = err; - return y; + } else if ( sigma < 0.4 ) { + useTemme = true; } + } + if ( useTemme ) { + evalMethod = 5; + } + // Regular case where the result will not be too close to 0.5: Changeover occurs at P ~ Q ~ 0.5. Note that series computation of P is about x2 faster than continued fraction calculation of Q, so try and use the CF only when really necessary, especially for small x. + else if ( x - ( 1.0 / (3.0 * x) ) < a ) { + evalMethod = 2; + } else { + evalMethod = 4; + invert = !invert; + } + } - y1 = 1.0; - if ( aid != 1.0 ) { - t = 0.0; - p = 1.0; - for ( i = 1; i < aid; i++ ) { - r = 1.0 - e + t; - p *= s * ( a+t+d2 ) * ( b+t+d2 ) / r; - t += 1.0; - p /= t; - y1 += p; + switch ( evalMethod ) { + case 0: + result = finiteGammaQ( a, x ); + if ( regularized == false ) { + result *= stdlib_base_gamma( a ); + } + break; + case 1: + result = finiteHalfGammaQ( a, x ); + if ( regularized == false ) { + result *= stdlib_base_gamma( a ); + } + break; + case 2: + // Compute P: + result = ( regularized ) ? + regularisedGammaPrefix( a, x ) : + fullIGammaPrefix( a, x ); + if ( result != 0.0 ) { + initValue = 0.0; + optimisedInvert = false; + if ( invert ) { + initValue = ( regularized ) ? 1.0 : stdlib_base_gamma( a ); + if ( + regularized || + result >= 1.0 || + STDLIB_CONSTANT_FLOAT64_MAX * result > initValue + ) { + initValue /= result; + if ( + regularized || + a < 1.0 || + ( STDLIB_CONSTANT_FLOAT64_MAX / a > initValue ) + ) { + initValue *= -a; + optimisedInvert = true; + } + else { + initValue = 0.0; + } + } + else { + initValue = 0.0; } } - - p = stdlib_base_gamma( c ); - y1 *= stdlib_base_gamma( e ) * p / ( stdlib_base_gamma( a+d1 ) * stdlib_base_gamma( b+d1 ) ); - y *= p / ( stdlib_base_gamma( a+d2 ) * stdlib_base_gamma( b+d2 ) ); - if ( ( (int)aid & 1 ) != 0 ) { - y = -y; - } - q = stdlib_base_pow( s, id ); - y = ( id > 0.0 ) ? y*q : y1*q; - y += y1; } - *loss = err; - return y; + result *= lowerGammaSeries( a, x, initValue ) / a; + if ( optimisedInvert ) { + invert = false; + result = -result; + } + break; + case 3: + // Compute Q: + invert = !invert; + // it can be changed to return result and option to store g: + tgammaSmallUpperPart( a, x, invert, out ); + result = out[ 0 ]; + g = out[ 1 ]; + invert = false; + if ( regularized ) { + result /= g; + } + break; + case 4: + // Compute Q: + result = ( regularized ) ? + regularisedGammaPrefix( a, x ) : + fullIGammaPrefix( a, x ); + if ( result != 0 ) { + result *= upperGammaFraction( a, x ); + } + break; + case 5: + result = igammaTemmeLarge( a, x ); + if ( x >= a ) { + invert = !invert; + } + break; + case 6: + // Since x is so small that P is necessarily very small too, use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/ + result = ( regularized ) ? + stdlib_base_pow(x, a) / stdlib_base_gamma( a + 1.0 ) : + stdlib_base_pow( x, a ) / a; + result *= 1.0 - ( a * x / ( a + 1.0 ) ); + break; + case 7: + // x is large, so compute Q: + result = ( regularized ) ? + regularisedGammaPrefix( a, x ) : + fullIGammaPrefix( a, x ); + result /= x; + if ( result != 0.0 ) { + result *= tgammaILargeX( a, x ); + } + break; + } + if ( regularized && result > 1.0 ) { + result = 1.0; } - // Perform power series if no special cases: - y = hys2f1( a, b, c, x, &err ); - *loss = err; - return y; + if ( invert ) { + gam = ( regularized ) ? 1.0 : stdlib_base_gamma( a ); + result = gam - result; + } + return result; } + + /** -* Computes the regularized incomplete gamma function. The upper tail is calculated via the modified Lentz's method for computing continued fractions, the lower tail using a power expansion. +* Computes the regularized incomplete stdlib_base_gamma function. The upper tail is calculated via the modified Lentz's method for computing continued fractions, the lower tail using a power expansion. * * ## Notes * -* - When `a >= FLOAT64_MAX_NTH_FACTORIAL` and computing the non-regularized incomplete gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result. +* - When `a >= FLOAT64_MAX_NTH_FACTORIAL` and computing the non-regularized incomplete stdlib_base_gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized stdlib_base_gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result. * * @param x input value * @param a input value @@ -463,24 +1009,13 @@ static double hyt2f1( const double a, const double b, const double c, const doub * @return function value * * @example -* double v = stdlib_base_gamma( 1.0, 1.0, 1.0, 0.0 ); +* double v = stdlib_base_gammainc( 1.0, 1.0, 1.0, 0.0 ); * // returns 1.0 */ double stdlib_base_gammainc( const double x, const double a, const bool regularized, const bool upper ) { - double optimisedInvert; - double evalMethod; double initValue; - double isHalfInt; - double useTemme; - double isSmallA; - double invert; double result; - double isInt; - double sigma; - double gam; - double res; - double fa; - double g; + bool invert; if ( x < 0.0 || a <= 0.0 ) { return 0.0 / 0.0; // NaN @@ -491,189 +1026,35 @@ double stdlib_base_gammainc( const double x, const double a, const bool regulari if ( invert && ( a * 4.0 < x ) ) { // This is method 4 below, done in logs: result = ( a * stdlib_base_ln(x) ) - x; - result += ln( upperGammaFraction( a, x ) ); + result += stdlib_base_ln( upperGammaFraction( a, x ) ); } else if ( !invert && ( a > 4.0 * x ) ) { // This is method 2 below, done in logs: - result = ( a * ln(x) ) - x; - initValue = 0; - result += ln( lowerGammaSeries( a, x, initValue ) / a ); + result = ( a * stdlib_base_ln(x) ) - x; + initValue = 0.0; + result += stdlib_base_ln( lowerGammaSeries( a, x, initValue ) / a ); } else { result = igammaFinal( a, x, true, invert ); if ( result == 0.0 ) { if ( invert ) { // Try http://functions.wolfram.com/06.06.06.0039.01: result = 1.0 + ( 1.0 / (12.0*a) ) + ( 1.0 / (288.0*a*a) ); - result = ln( result ) - a + ( ( a-0.5 ) * ln(a) ); - result += ln( STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); + result = stdlib_base_ln( result ) - a + ( ( a-0.5 ) * stdlib_base_ln(a) ); + result += stdlib_base_ln( STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); } else { // This is method 2 below, done in logs, we're really outside the range of this method, but since the result is almost certainly infinite, we should probably be OK: - result = ( a * ln( x ) ) - x; + result = ( a * stdlib_base_ln( x ) ) - x; initValue = 0.0; - result += ln( lowerGammaSeries( a, x, initValue ) / a ); + result += stdlib_base_ln( lowerGammaSeries( a, x, initValue ) / a ); } } else { - result = ln( result ) + gammaln( a ); + result = stdlib_base_ln( result ) + stdlib_base_gammaln( a ); } } if ( result > STDLIB_CONSTANT_FLOAT64_MAX_LN ) { return STDLIB_CONSTANT_FLOAT64_PINF; } - return exp( result ); + return stdlib_base_exp( result ); } - // If no special handling is required then we proceeds as normal: + // If no special handling is required then we proceed as normal: return igammaFinal( a, x, regularized, invert ); - - // isSmallA = ( a < 30 ) && ( a <= x + 1.0 ) && ( x < MAX_LN ); - // if ( isSmallA ) { - // fa = floor( a ); - // isInt = ( fa == a ); - // isHalfInt = ( isInt ) ? false : ( abs( fa - a ) == 0.5 ); - // } else { - // isInt = isHalfInt = false; - // } - // if ( isInt && x > 0.6 ) { - // // Calculate Q via finite sum: - // invert = !invert; - // evalMethod = 0; - // } else if ( isHalfInt && x > 0.2 ) { - // // Calculate Q via finite sum for half integer a: - // invert = !invert; - // evalMethod = 1; - // } else if ( x < SQRT_EPSILON && a > 1.0 ) { - // evalMethod = 6; - // } else if ( x < 0.5 ) { - // // Changeover criterion chosen to give a changeover at Q ~ 0.33: - // if ( -0.4 / ln( x ) < a ) { - // evalMethod = 2; - // } else { - // evalMethod = 3; - // } - // } else if ( x < 1.1 ) { - // // Changeover here occurs when P ~ 0.75 or Q ~ 0.25: - // if ( x * 0.75 < a ) { - // evalMethod = 2; - // } else { - // evalMethod = 3; - // } - // } else { - // // Begin by testing whether we're in the "bad" zone where the result will be near 0.5 and the usual series and continued fractions are slow to converge: - // useTemme = false; - // if ( regularized && a > 20 ) { - // sigma = abs( (x-a)/a ); - // if ( a > 200 ) { - // // Limit chosen so that we use Temme's expansion only if the result would be larger than about 10^-6. Below that the regular series and continued fractions converge OK, and if we use Temme's method we get increasing errors from the dominant erfc term as it's (inexact) argument increases in magnitude. - // if ( 20 / a > sigma * sigma ) { - // useTemme = true; - // } - // } else if ( sigma < 0.4 ) { - // useTemme = true; - // } - // } - // if ( useTemme ) { - // evalMethod = 5; - // } - // // Regular case where the result will not be too close to 0.5: Changeover occurs at P ~ Q ~ 0.5. Note that series computation of P is about x2 faster than continued fraction calculation of Q, so try and use the CF only when really necessary, especially for small x. - // else if ( x - ( 1.0 / (3.0 * x) ) < a ) { - // evalMethod = 2; - // } else { - // evalMethod = 4; - // invert = !invert; - // } - // } - - // /* eslint-disable default-case */ - // switch ( evalMethod ) { - // case 0: - // result = finiteGammaQ( a, x ); - // if (regularized == false ) { - // result *= gamma( a ); - // } - // break; - // case 1: - // result = finiteHalfGammaQ( a, x ); - // if ( regularized == false ) { - // result *= gamma( a ); - // } - // break; - // case 2: - // // Compute P: - // result = ( regularized ) ? - // regularisedGammaPrefix( a, x ) : - // fullIGammaPrefix( a, x ); - // if ( result !== 0.0 ) { - // initValue = 0.0; - // optimisedInvert = false; - // if ( invert ) { - // initValue = ( regularized ) ? 1.0 : gamma( a ); - // if ( - // regularized || - // result >= 1.0 || - // FLOAT64_MAX * result > initValue - // ) { - // initValue /= result; - // if ( - // regularized || - // a < 1.0 || - // ( FLOAT64_MAX / a > initValue ) - // ) { - // initValue *= -a; - // optimisedInvert = true; - // } - // else { - // initValue = 0.0; - // } - // } - // else { - // initValue = 0.0; - // } - // } - // } - // result *= lowerGammaSeries( a, x, initValue ) / a; - // if ( optimisedInvert ) { - // invert = false; - // result = -result; - // } - // break; - // case 3: - // // Compute Q: - // invert = !invert; - // res = tgammaSmallUpperPart( a, x, invert ); - // result = res[ 0 ]; - // g = res[ 1 ]; - // invert = false; - // if ( regularized ) { - // result /= g; - // } - // break; - // case 4: - // // Compute Q: - // result = ( regularized ) ? - // regularisedGammaPrefix( a, x ) : - // fullIGammaPrefix( a, x ); - // if ( result !== 0 ) { - // result *= upperGammaFraction( a, x ); - // } - // break; - // case 5: - // result = igammaTemmeLarge( a, x ); - // if ( x >= a ) { - // invert = !invert; - // } - // break; - // case 6: - // // Since x is so small that P is necessarily very small too, use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/ - // result = ( regularized ) ? - // pow(x, a) / gamma( a + 1.0 ) : - // pow( x, a ) / a; - // result *= 1.0 - ( a * x / ( a + 1.0 ) ); - // break; - // } - // if ( regularized && result > 1.0 ) { - // result = 1.0; - // } - // if ( invert ) { - // gam = ( regularized ) ? 1.0 : gamma( a ); - // result = gam - result; - // } - // return result; } From 8cc8e26dab3e20820e6f6c581a8a7a183b3de0d3 Mon Sep 17 00:00:00 2001 From: Karan Anand Date: Thu, 19 Jun 2025 00:12:38 -0700 Subject: [PATCH 03/14] feat: add C implementation for `math/base/special/gammainc` --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: passed - task: lint_package_json status: passed - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: passed - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: passed - task: lint_c_examples status: passed - task: lint_c_benchmarks status: passed - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../math/base/special/gammainc/README.md | 100 +++++ .../gammainc/benchmark/benchmark.native.js | 135 ++++++ .../gammainc/benchmark/c/native/Makefile | 146 ++++++ .../gammainc/benchmark/c/native/benchmark.c | 140 ++++++ .../math/base/special/gammainc/binding.gyp | 170 +++++++ .../base/special/gammainc/examples/c/Makefile | 146 ++++++ .../special/gammainc/examples/c/example.c | 40 ++ .../math/base/special/gammainc/include.gypi | 53 +++ .../stdlib/math/base/special/gammainc.h | 40 ++ .../math/base/special/gammainc/lib/native.js | 69 +++ .../math/base/special/gammainc/manifest.json | 160 +++++++ .../math/base/special/gammainc/package.json | 4 +- .../math/base/special/gammainc/src/addon.c | 30 +- .../math/base/special/gammainc/src/main.c | 422 +++++++++--------- .../base/special/gammainc/test/test.native.js | 201 +++++++++ 15 files changed, 1634 insertions(+), 222 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/benchmark/benchmark.native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/benchmark/c/native/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/benchmark/c/native/benchmark.c create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/binding.gyp create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/examples/c/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/examples/c/example.c create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/include.gypi create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/include/stdlib/math/base/special/gammainc.h create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/lib/native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/manifest.json create mode 100644 lib/node_modules/@stdlib/math/base/special/gammainc/test/test.native.js diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/README.md b/lib/node_modules/@stdlib/math/base/special/gammainc/README.md index 6abe4afc36b1..476abde9bce5 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammainc/README.md +++ b/lib/node_modules/@stdlib/math/base/special/gammainc/README.md @@ -173,6 +173,106 @@ logEachMap( 'x: %0.4f, \t s: %0.4f, \t f(x,s): %0.4f', x, s, gammainc ); + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/gammainc.h" +``` + +#### stdlib_base_gammainc( x, a, regularized, upper ) + +Evaluates the [incomplete gamma function][incomplete-gamma-function] for inputs `x` and `a`. The third and fourth parameters of the function can be used to specify whether instead to evaluate the non-regularized and/or upper incomplete gamma functions, respectively. + +```c +double out = stdlib_base_gammainc( 0.0, 1.0, true, false ); +// returns 0.0 + +out = stdlib_base_gammainc( 6.0, 2.0, true, false ); +// returns ~0.9826 +``` + +The function accepts the following arguments: + +- **x**: `[in] double` input value. +- **a**: `[in] double` input value. +- **regularized**: `[in] bool` input value. +- **upper**: `[in] bool` input value. + +```c +double stdlib_base_gammainc( const double x, const double a, const bool regularized, const bool upper ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/gammainc.h" +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + double x; + double a; + double y; + int i; + + for ( i = 0; i < 10; i++ ) { + x = random_uniform( 0.0, 1000.0 ); + a = random_uniform( 1.0, 1000.0 ); + y = stdlib_base_gammainc( x, a, true, false ); + printf( "x: %lf, s: %lf, f(x,s): %lf\n", x, a, y ); + } +} +``` + +
+ + + +
+ + +