diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md b/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md index 7426de70ddd1..c7c676f38ac8 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md @@ -138,6 +138,64 @@ logEachMap( 'p: %0.4f, \t a: %0.4f, \t P^(-1)(p, a): %0.4f', p, a, gammaincinv ) +
+ +## C API + + +
+
+#include "stdlib/math/base/special/gammaincinv.h"
+#include 
+
+double stdlib_base_special_gammaincinv( const double p, const double a, const bool upper )
+
+
+ +

+Inverts the lower incomplete gamma function. +

+ +
+
+double y = stdlib_base_special_gammaincinv( 0.5, 1.0, false );
+// returns ~0.693
+
+y = stdlib_base_special_gammaincinv( 0.1, 1.0, true );
+// returns ~2.303
+
+
+ +The function accepts the following arguments: + +* p: `[in] double` probability value. +* a: `[in] double` scale parameter. +* upper: `[in] bool` boolean indicating if the function should invert the upper tail of the incomplete gamma function instead. + +
+
+double stdlib_base_special_gammaincinv( const double p, const double a, const bool upper );
+
+
+ +

+ +

+ + +
+
+$ stdlib repl
+> #include "stdlib/math/base/special/gammaincinv.h"
+> #include 
+>
+> double y = stdlib_base_special_gammaincinv( 0.8, 2.0, false );
+> printf( "y = %f\n", y );
+> // => y = 1.213038...
+
+
+ +
* * * diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/benchmark.native.js new file mode 100644 index 000000000000..446d9e073a48 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/benchmark.native.js @@ -0,0 +1,77 @@ +/** +* @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. +*/ + +'use strict'; + +// MODULES // + +var bench = require( '@stdlib/bench' ); +var randu = require( '@stdlib/random/base/randu' ); +var pow = require( '@stdlib/math/base/special/pow' ); +var isFunction = require( '@stdlib/assert/is-function' ); +var pkg = require( './../package.json' ).name; +var gammaincinv = require( './../lib/native.js' ); + + +// MAIN // + +bench( pkg+'::native', function benchmark( b ) { + var x; + var a; + var y; + var i; + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + x = ( randu()*1.0 ); + a = ( randu()*10.0 ); + y = gammaincinv( x, a, false ); + if ( y !== y ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( y !== y ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); + +bench( pkg+'::native:upper', function benchmark( b ) { + var x; + var a; + var y; + var i; + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + x = ( randu()*1.0 ); + a = ( randu()*10.0 ); + y = gammaincinv( x, a, true ); + if ( y !== y ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( y !== y ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/Makefile b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/Makefile new file mode 100644 index 000000000000..4f938a34c498 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/Makefile @@ -0,0 +1,88 @@ +# +# @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. +# + +# Define the Node.js executable: +NODE_EXE ?= node + +# Get the project root directory: +PROJECT_ROOT_DIR := $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammaincinv' ), '..', '..', '..', '..', '..' ) )" ) + +# Define the path to the stdlib make includes: +stdlib_make_path := $(PROJECT_ROOT_DIR)/tools/make/lib + +include $(stdlib_make_path)/Makefile.c.benchmark + +# Define the C source file: +SRC := $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c + +# Define the include directories: +INCLUDE_DIRS := \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaincinv/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/assert/is-nan/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/abs/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/erfcinv/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/exp/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gamma/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammainc/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaln/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/ln/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/pow/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/sqrt/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float32/max/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float32/smallest-normal/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/ln-sqrt-two-pi/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/pinf/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/sqrt-two-pi/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/two-pi/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/randu/include + +# Add commands for including C source files from dependencies: +include $(stdlib_make_path)/c.dependencies + +# Add a command to include the C source file for `is_nan`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/assert/is-nan' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `abs`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/abs' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `erfcinv`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/erfcinv' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `exp`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/exp' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `gamma`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gamma' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `gammainc`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammainc' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `gammaln`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammaln' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `ln`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/ln' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `pow`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/pow' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `sqrt`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/sqrt' ), '..', '..', 'make', 'c.src.mk' ) )" ) + +# Add a command to include the C source file for `randu`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/randu' ), '..', '..', 'make', 'c.src.mk' ) )" ) diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/benchmark.c b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/benchmark.c new file mode 100644 index 000000000000..dff384f15d52 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/benchmark.c @@ -0,0 +1,113 @@ +/** +* @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/gammaincinv.h" +#include "stdlib/math/base/special/randu.h" +#include "stdlib/math/base/special/pow.h" +#include +#include +#include +#include + +// MACROS // + +#define NAME "gammaincinv" +#define ITERATIONS 1000000 +#define REPEATS 3 + +// FUNCTIONS // + +/** +* Returns the current time in nanoseconds. +*/ +static uint64_t tic( void ) { + struct timeval tv; + gettimeofday( &tv, NULL ); + return (uint64_t)tv.tv_sec * 1000000000L + (uint64_t)tv.tv_usec * 1000L; +} + +/** +* Prints the TAP version. +*/ +static void print_version( void ) { + printf( "TAP version 13\n" ); +} + +/** +* Prints the TAP plan. +*/ +static void print_plan( const int n ) { + printf( "1..%d\n", n ); +} + +/** +* Prints the benchmark results. +*/ +static void print_results( const char *name, const uint64_t elapsed, const int iterations ) { + double rate = (double)iterations / ( (double)elapsed / 1000000000.0 ); + printf( " ---\n" ); + printf( " iterations: %d\n", iterations ); + printf( " elapsed: %llu\n", elapsed ); + printf( " rate: %0.2f\n", rate ); + printf( " ...\n" ); +} + +// MAIN // + +int main( void ) { + double x[ ITERATIONS ]; + double a[ ITERATIONS ]; + double y; + int i; + int j; + uint64_t elapsed; + uint64_t t; + + // Use the same values as the JS benchmark: + for ( i = 0; i < ITERATIONS; i++ ) { + x[ i ] = stdlib_base_randu() * 1.0; + a[ i ] = stdlib_base_randu() * 10.0; + } + + print_version(); + print_plan( REPEATS ); + + for ( j = 0; j < REPEATS; j++ ) { + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_special_gammaincinv( x[ i ], a[ i ], false ); + if ( y != y ) { + printf( "should not return NaN\n" ); + } + } + elapsed = tic() - t; + printf( "ok %d benchmark( %s, %s, %s ) - lower\n", j+1, NAME, "number", "number" ); + print_results( NAME, elapsed, ITERATIONS ); + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_special_gammaincinv( x[ i ], a[ i ], true ); + if ( y != y ) { + printf( "should not return NaN\n" ); + } + } + elapsed = tic() - t; + printf( "ok %d benchmark( %s, %s, %s ) - upper\n", j+2, NAME, "number", "number" ); + print_results( NAME, elapsed, ITERATIONS ); + } +} diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/binding.gyp b/lib/node_modules/@stdlib/math/base/special/gammaincinv/binding.gyp new file mode 100644 index 000000000000..747f6e159058 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/binding.gyp @@ -0,0 +1,66 @@ +# +# @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. +# + +{ + "targets": [ + { + "target_name": "math.base.special.gammaincinv", + "include_dirs": [ + "include", # Self + "../../assert/is-nan/include", # Dependency includes... + "../abs/include", + "../erfcinv/include", + "../exp/include", + "../gamma/include", + "../gammainc/include", + "../gammaln/include", + "../ln/include", + "../pow/include", + "../sqrt/include", + "../../../../constants/float32/max/include", + "../../../../constants/float32/smallest-normal/include", + "../../../../constants/float64/ln-sqrt-two-pi/include", + "../../../../constants/float64/pinf/include", + "../../../../constants/float64/sqrt-two-pi/include", + "../../../../constants/float64/two-pi/include" + ], + "sources": [ + "src/main.c", + "src/stdlib_base_special_gammaincinv.h" + ], + "dependencies": [ + "../../assert/is-nan/binding.gyp", # Dependency gyp files... + "../abs/binding.gyp", + "../erfcinv/binding.gyp", + "../exp/binding.gyp", + "../gamma/binding.gyp", + "../gammainc/binding.gyp", + "../gammaln/binding.gyp", + "../ln/binding.gyp", + "../pow/binding.gyp", + "../sqrt/binding.gyp", + "../../../../constants/float32/max/binding.gyp", + "../../../../constants/float32/smallest-normal/binding.gyp", + "../../../../constants/float64/ln-sqrt-two-pi/binding.gyp", + "../../../../constants/float64/pinf/binding.gyp", + "../../../../constants/float64/sqrt-two-pi/binding.gyp", + "../../../../constants/float64/two-pi/binding.gyp" + ] + } + ] +} diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/Makefile new file mode 100644 index 000000000000..2e5f67f3f780 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/Makefile @@ -0,0 +1,84 @@ +# +# @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. +# + +# Define the Node.js executable: +NODE_EXE ?= node + +# Get the project root directory: +PROJECT_ROOT_DIR := $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammaincinv' ), '..', '..', '..', '..', '..' ) )" ) + +# Define the path to the stdlib make includes: +stdlib_make_path := $(PROJECT_ROOT_DIR)/tools/make/lib + +include $(stdlib_make_path)/Makefile.c.example + +# Define the C source file: +SRC := $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c + +# Define the include directories: +INCLUDE_DIRS := \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaincinv/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/assert/is-nan/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/abs/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/erfcinv/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/exp/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gamma/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammainc/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/gammaln/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/ln/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/pow/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/math/base/special/sqrt/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float32/max/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float32/smallest-normal/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/ln-sqrt-two-pi/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/pinf/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/sqrt-two-pi/include \ + $(PROJECT_ROOT_DIR)/lib/node_modules/@stdlib/constants/float64/two-pi/include + +# Add commands for including C source files from dependencies: +include $(stdlib_make_path)/c.dependencies + +# Add a command to include the C source file for `is_nan`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/assert/is-nan' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `abs`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/abs' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `erfcinv`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/erfcinv' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `exp`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/exp' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `gamma`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gamma' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `gammainc`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammainc' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `gammaln`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/gammaln' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `ln`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/ln' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `pow`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/pow' ), '..', '..', 'make' ) )" )/c.src.mk + +# Add a command to include the C source file for `sqrt`: +include $(shell $(NODE_EXE) -e "console.log( require('path').join( require.resolve( '@stdlib/math/base/special/sqrt' ), '..', '..', 'make' ) )" )/c.src.mk diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/example.c new file mode 100644 index 000000000000..dff0273de5a8 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/example.c @@ -0,0 +1,45 @@ +/** +* @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/gammaincinv.h" +#include +#include + +int main( void ) { + const double a[] = { 0.1, 0.5, 1.0, 2.0, 3.0 }; + const double p[] = { 0.1, 0.3, 0.5, 0.7, 0.9 }; + + double y; + int i; + int j; + + for ( i = 0; i < 5; i++ ) { + for ( j = 0; j < 5; j++ ) { + y = stdlib_base_special_gammaincinv( p[ j ], a[ i ], false ); + printf( "gammaincinv( %0.2f, %0.2f, false ) = %0.16f\n", p[ j ], a[ i ], y ); + } + } + + printf( "\n" ); + for ( i = 0; i < 5; i++ ) { + for ( j = 0; j < 5; j++ ) { + y = stdlib_base_special_gammaincinv( p[ j ], a[ i ], true ); + printf( "gammaincinv( %0.2f, %0.2f, true ) = %0.16f\n", p[ j ], a[ i ], y ); + } + } +} diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/include/stdlib/math/base/special/gammaincinv.h b/lib/node_modules/@stdlib/math/base/special/gammaincinv/include/stdlib/math/base/special/gammaincinv.h new file mode 100644 index 000000000000..fedff5399759 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/include/stdlib/math/base/special/gammaincinv.h @@ -0,0 +1,42 @@ +/** +* @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. +*/ + +#ifndef STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H +#define STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H + +// MODULES // + +#include + +// FUNCTIONS // + +/** +* Inverts the lower incomplete gamma function. +* +* @param p - probability value +* @param a - scale parameter +* @param upper - boolean indicating if the function should invert the upper tail +* @return function value +* +* @example +* double y = stdlib_base_special_gammaincinv( 0.5, 1.0, false ); +* // returns ~0.693 +*/ +double stdlib_base_special_gammaincinv( const double p, const double a, const bool upper ); + +#endif // !STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/lib/native.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/lib/native.js new file mode 100644 index 000000000000..a488c4aad254 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/lib/native.js @@ -0,0 +1,59 @@ +/** +* @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. +*/ + +'use strict'; + +// MODULES // + +var join = require( 'path' ).join; +var tryRequire = require( '@stdlib/utils/try-require' ); +var isError = require( '@stdlib/assert/is-error' ); + + +// MAIN // + +// Attempt to load the native addon: +var addon = tryRequire( join( __dirname, '..', 'build', 'Release', 'math.base.special.gammaincinv.node' ) ); // eslint-disable-line stdlib/no-dynamic-require + +// If we failed to load the addon, raise an error: +if ( isError( addon ) ) { + throw new Error( 'failed to load native addon. Error: ' + addon.message ); +} + +/** +* Inverts the lower incomplete gamma function. +* +* @param {Probability} p - probability value +* @param {number} a - scale parameter +* @param {boolean} [upper=false] - boolean indicating if the function should invert the upper tail +* @returns {number} function value +* +* @example +* var y = gammaincinv( 0.5, 1.0, false ); +* // returns ~0.693 +*/ +function gammaincinv( p, a, upper ) { + // Directly call the loaded native addon function + // The native function likely expects (double, double, bool) + return addon( p, a, upper === true ); +} + + +// EXPORTS // + +module.exports = gammaincinv; diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json b/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json index 3bc5f2084385..0e1ebe8f5006 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json @@ -14,10 +14,12 @@ } ], "main": "./lib", + "gypfile": true, "directories": { "benchmark": "./benchmark", "doc": "./docs", "example": "./examples", + "include": "./include", "lib": "./lib", "scripts": "./scripts", "test": "./test" 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..61be12603aff --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c @@ -0,0 +1,831 @@ +/** +* @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/gammaincinv.h" +#include "stdlib/math/base/special/gammainc.h" +#include "stdlib/math/base/special/gammaln.h" +#include "stdlib/math/base/special/erfcinv.h" +#include "stdlib/math/base/special/gamma.h" +#include "stdlib/math/base/special/sqrt.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/pinf.h" +#include "stdlib/constants/float64/sqrt_two_pi.h" +#include "stdlib/constants/float64/two_pi.h" +#include "stdlib/constants/float64/ln_sqrt_two_pi.h" +#include "stdlib/constants/float32/max.h" +#include "stdlib/constants/float32/smallest_normal.h" +#include + +// MODULES // + +// In C, we replace stdlib's evalpoly with our own static helpers based on the JS files. +// We also translate all other helper functions. + +// From: evalpoly/lib/main.js (Horner's method) +static double evalpoly( const double *c, const int len, const double x ) { + double p; + int i; + + if ( len < 2 || x == 0.0 ) { + if ( len == 0 ) { + return 0.0; + } + return c[ 0 ]; + } + i = len - 1; + p = ( c[ i ] * x ) + c[ i - 1 ]; + i -= 2; + while ( i >= 0 ) { + p = ( p * x ) + c[ i ]; + i -= 1; + } + return p; +} + +// From: polyval_ak1.js +static double polyval_ak1( const 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))))))))))); + static const double C[] = { + 0.0, + 1.0, + 1.0, + 1.5, + 2.6666666666666665, + 5.208333333333333, + 10.8 + }; + return evalpoly( C, 7, x ); +} + +// From: polyval_ak2.js +static double polyval_ak2( const 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))))))))))); + static const double C[] = { + 1.0, + 1.0, + 0.3333333333333333, + 0.027777777777777776, + -0.003703703703703704, + 0.0002314814814814815, + 0.00005878894767783657 + }; + return evalpoly( C, 7, x ); +} + +// From: rational_ak0bk0.js +static double rational_ak0bk0( const double x ) { + double s1; + double s2; + static const double P[] = { + -0.3333333333438, + -0.2070740359969, + -0.05041806657154, + -0.004923635739372, + -0.00004293658292782 + }; + static const double Q[] = { + 1.0, + 0.7045554412463, + 0.2118190062224, + 0.03048648397436, + 0.001605037988091 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak1bk1.js +static double rational_ak1bk1( const double x ) { + double s1; + double s2; + static const double P[] = { + -0.0172847633523, + -0.0159372646475, + -0.00464910887221, + -0.00060683488776, + -0.00000614830384279 + }; + static const double Q[] = { + 1.0, + 0.764050615669, + 0.297143406325, + 0.0579490176079, + 0.00574558524851 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak2bk2.js (same as rational_ak3bk3.js) +static double rational_ak2bk2( const double x ) { + double s1; + double s2; + static const double P[] = { + -0.0172839517431, + -0.0146362417966, + -0.00357406772616, + -0.000391032032692, + 0.00000249634036069 + }; + static const double Q[] = { + 1.0, + 0.690560400696, + 0.249962384741, + 0.0443843438769, + 0.00424073217211 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak4bk4.js +static double rational_ak4bk4( const double x ) { + double s1; + double s2; + static const double P[] = { + 0.0495346498136, + 0.0299521337141, + 0.00688296911516, + 0.000512634846317, + -0.0000201411722031 + }; + static const double Q[] = { + 1.0, + 0.759803615283, + 0.261547111595, + 0.0464854522477, + 0.00403751193496 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak5bk5.js +static double rational_ak5bk5( const double x ) { + double s1; + double s2; + static const double P[] = { + 0.00452313583942, + 0.00120744920113, + -0.0000789724156582, + -0.0000504476066942, + -0.00000535770949796 + }; + static const double Q[] = { + 1.0, + 0.912203410349, + 0.405368773071, + 0.0901638932349, + 0.00948935714996 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak6bk6.js +static double rational_ak6bk6( const double x ) { + double s1; + double s2; + static const double P[] = { + 0.00439937562904, + 0.000487225670639, + -0.000128470657374, + 0.00000529110969589, + 1.5716677175e-7 + }; + static const double Q[] = { + 1.0, + 0.794435257415, + 0.333094721709, + 0.0703527806143, + 0.00806110846078 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak7bk7.js +static double rational_ak7bk7( const double x ) { + double s1; + double s2; + static const double P[] = { + -0.0011481191232, + -0.112850923276, + 1.51623048511, + -0.218472031183, + 0.0730002451555 + }; + static const double Q[] = { + 1.0, + 14.2482206905, + 69.7360396285, + 218.938950816, + 277.067027185 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: rational_ak8bk8.js +static double rational_ak8bk8( const double x ) { + double s1; + double s2; + static const double P[] = { + -0.000145727889667, + -0.290806748131, + -13.308504545, + 199.722374056, + -11.4311378756 + }; + static const double Q[] = { + 1.0, + 139.612587808, + 2189.01116348, + 7115.24019009, + 45574.6081453 + }; + s1 = evalpoly( P, 5, x ); + s2 = evalpoly( Q, 5, x ); + return s1 / s2; +} + +// From: chepolsum.js +static double chepolsum( const int n, const double t ) { + 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 + }; + double tt; + double u0; + double u1; + double u2; + int k; + + 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; +} + +// From: polyval_c.js +static double polyval_c( const 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))))))))); + static const double C[] = { + 0.025721014990011306, + 0.08247596616699963, + -0.0025328157302663564, + 0.0006099292666946337, + -0.00033543297638406, + 0.000250505279903 + }; + return evalpoly( C, 6, x ); +} + +// From: polyval_d.js +static double polyval_d( const double x ) { + if ( x == 0.0 ) { + return 0.08333333333333333; + } + // return 0.08333333333333333 + (x * (-0.002777777777777778 + (x * (0.0007936507936507937 + (x * -0.0005952380952380953))))); + static const double C[] = { + 0.08333333333333333, + -0.002777777777777778, + 0.0007936507936507937, + -0.0005952380952380953 + }; + return evalpoly( C, 4, x ); +} + +// From: stirling.js +static double stirling( const double x ) { + double z; + static const double C6 = 0.30865217988013567769; + + if ( x < STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL ) { + return STDLIB_CONSTANT_FLOAT32_MAX; + } + if ( x < 1.0 ) { + return stdlib_base_special_gammaln( x+1.0 ) - ( (x+0.5) * stdlib_base_special_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI; + } + if ( x < 2.0 ) { + return stdlib_base_special_gammaln( x ) - ( (x-0.5) * stdlib_base_special_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI; + } + if ( x < 3.0 ) { + return stdlib_base_special_gammaln( x-1.0 ) - ( (x-0.5) * stdlib_base_special_ln(x) ) + x - STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI + stdlib_base_special_ln( x-1.0 ); + } + 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; +} + +// From: gamstar.js +static double gamstar( const double x ) { + if ( x >= 3.0 ) { + return stdlib_base_special_exp( stirling(x) ); + } + if ( x > 0.0 ) { + return stdlib_base_special_gamma(x) / ( stdlib_base_special_exp( -x + ( ( x-0.5 ) * stdlib_base_special_ln(x) ) ) * STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); + } + // Case: x <= 0.0 + return STDLIB_CONSTANT_FLOAT32_MAX; +} + +// From: lambdaeta.js +static double lambdaeta( const double eta ) { + double L2; + double L3; + double L4; + double L5; + double la; + double L; + double q; + double r; + double s; + + static const double THRESHOLD = 1.0e-8; + static const double ONEO12 = 0.0833333333333333333333333333333; + static const double ONEO120 = 0.00833333333333333333333333333333; + + // Polynomial coefficient workspace: + double AK[] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; + + s = eta * eta * 0.5; + if ( eta == 0.0 ) { + la = 0.0; + } + else if ( eta < -1.0 ) { + r = stdlib_base_special_exp( -1.0 - s ); + la = polyval_ak1( r ); + } + else if ( eta < 1.0 ) { + r = eta; + la = polyval_ak2( r ); + } + else { + r = 11.0 + s; + L = stdlib_base_special_ln( r ); + la = r + L; + r = 1.0 / r; + L2 = L * L; + L3 = L2 * L; + L4 = L3 * L; + L5 = L4 * L; + AK[ 1 ] = ( 2.0-L ) * 0.5; + AK[ 2 ] = ( ( -9.0*L ) + 6.0 + ( 2.0*L2 ) ) / 6.0; + AK[ 3 ] = -( (3*L3)+ (36*L) - (22*L2) - 12 ) * ONEO12; + AK[ 4 ] = ( 60.0 + (350.0*L2) - (300.0*L) - (125.0*L3) + (12.0*L4) ) / 60.0; + AK[ 5 ] = -(-120 - (274*L4) + (900*L) - (1700*L2) + (1125*L3) + (20*L5)) * ONEO120; + la += ( L * r * evalpoly( AK, 6, r ) ); + } + r = 1.0; + if ( + ( eta > -3.5 && eta < -0.03 ) || + ( eta > 0.03 && eta < 40.0 ) + ) { + r = 1.0; + q = la; + do { + la = q * ( s+stdlib_base_special_ln(q) ) / ( q-1.0 ); + r = stdlib_base_special_abs( ( q/la ) - 1.0 ); + q = la; + } while ( r > THRESHOLD ); + } + return la; +} + +// From: eps1.js +static double eps1( const double eta ) { + double la; + if ( stdlib_base_special_abs( eta ) < 1.0 ) { + return rational_ak0bk0( eta ); + } + la = lambdaeta( eta ); + return stdlib_base_special_ln( eta / ( la - 1.0 ) ) / eta; +} + +// From: eps2.js +static double eps2( const double eta ) { + double lnmeta; + double x; + if ( eta < -5.0 ) { + x = eta * eta; + lnmeta = stdlib_base_special_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 ); // rational_ak3bk3 is identical + } + if ( eta < 1000.0 ) { + // x = 1.0 / eta; // This variable 'x' is not used in the JS logic + return rational_ak2bk2( eta ) / ( -12.0*eta ); // Re-using rational_ak2bk2 as it's same as ak3bk3 + } + return -1.0 / ( 12.0 * eta ); +} + +// From: eps3.js +static double eps3( const double eta ) { + double x; + double y; + + if ( eta < -8.0 ) { + x = eta * eta; + y = stdlib_base_special_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 ); // rational_ak6bk6 is provided twice + } + 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_special_ln( eta ) / ( 12.0*eta*eta*eta ); +} + +// From: higher_newton.js +static double higherNewton( double x0, const double a, const int m, const double p, const double q, const double lgama, const double invfp, const bool pcase ) { + double dlnr; + double xini; + double ck0; + double ck1; + double ck2; + double a2; + double x2; + double px; + double qx; + double xr; + double t; + int n; + double r; + double x; + + x = x0; + t = 1; + n = 1; + a2 = a * a; + xini = x0; + do { + x = x0; + x2 = x * x; + if ( m == 0 ) { + dlnr = ( ( 1.0-a ) * stdlib_base_special_ln( x ) ) + x + lgama; + if ( dlnr > stdlib_base_special_ln( STDLIB_CONSTANT_FLOAT32_MAX ) ) { + // Warning: overflow... + return xini; + } + r = stdlib_base_special_exp( dlnr ); + } else { + r = -invfp * x; + } + if ( pcase ) { + // Call: gammainc( x, s, lower ) + px = stdlib_base_special_gammainc( x, a, true, false ); + ck0 = -r * ( px - p ); + } else { + // Call: gammainc( x, s, upper ) + qx = stdlib_base_special_gammainc( x, a, true, true ); + ck0 = r * ( qx - q ); + } + r = ck0; + if ( ( p > 1e-120 ) || ( n > 1 ) ) { + ck1 = 0.5 * ( x - a + 1.0 ) / x; + ck2 = ( (2*x2) - (4*x*a) + (4*x) + (2*a2) - (3*a) + 1 ) / x2; + ck2 /= 6.0; + x0 = x + ( r * ( 1.0 + ( r * ( ck1 + (r*ck2) ) ) ) ); + } else { + x0 = x + r; + } + t = stdlib_base_special_abs( ( x/x0 ) - 1.0 ); + n += 1; + x = x0; + if ( x < 0 ) { + x = xini; + n = 100; + } + } while ( ( ( t > 2e-14 ) && ( n < 35 ) ) ); + + xr = x; + return xr; +} + +// From: compute.js +static double compute( const double a, const double p, const double q ) { + double ap1inv; + double invfp = 0.0; // Initialize to avoid uninitialized variable warning + double lgama = 0.0; // Initialize + bool pcase; + double porq; + double ainv; + double logr; + double ap22; + double ap14; + double ap13; + double ap12; + double vgam; + double vmin; + double xini; + double ap1; + double ap2; + double ap3; + double eta; + double x0; + double a2; + double L2; + double L3; + double L4; + double b2; + double b3; + double a4; + double fp; + double a3; + double xr = 0.0; // Initialize + double b; + double L; + int i; + int k; + int m = 0; // Initialize + double r; + double s; + double t; + double y; + + // Coefficient workspace: + double CK[] = { 0.0, 0.0, 0.0, 0.0, 0.0 }; + + static const double HALF = 0.5; + static const double ONEO3 = 0.333333333333333333333333333333; + static const double ONEO4 = 0.25; + static const double ONEO5 = 0.2; + static const double ONEO6 = 0.166666666666666666666666666667; + static const double ONEO12 = 0.0833333333333333333333333333333; + + if ( p < HALF ) { + pcase = true; + porq = p; + s = -1.0; + } else { + pcase = false; + porq = q; + s = 1.0; + } + k = 0; + if ( stdlib_base_special_abs( a-1.0 ) < 1.0e-4 ) { + m = 0; + if ( pcase ) { + if ( p < 1.0e-3 ) { + // Use polyval_ak1 coefficients for the Taylor expansion + const double P_TAYLOR[] = { 0.0, 1.0, 0.5, ONEO3, ONEO4, ONEO5, ONEO6 }; + x0 = evalpoly( P_TAYLOR, 7, p ); + } else { + x0 = -stdlib_base_special_ln( 1.0-p ); + } + } else { + x0 = -stdlib_base_special_ln( q ); + } + if ( a == 1.0 ) { + k = 2; + xr = x0; + } else { + lgama = stdlib_base_special_gammaln( a ); + k = 1; + } + } + if ( q < 1.0e-30 && a < HALF ) { + m = 0; + double gamma_a = stdlib_base_special_gamma(a); + x0 = -stdlib_base_special_ln( q*gamma_a ) + ( ( a-1.0 ) * stdlib_base_special_ln( -stdlib_base_special_ln( q*gamma_a ) ) ); + k = 1; + lgama = stdlib_base_special_gammaln( a ); + } + if ( a > 1.0 && a < 500.0 && p < 1.0e-80 ) { + m = 0; + ainv = 1.0 / a; + ap1inv = 1.0 / ( a+1.0 ); + x0 = ( stdlib_base_special_gammaln( a+1.0 ) + stdlib_base_special_ln( p ) ) * ainv; + x0 = stdlib_base_special_exp( x0 ); + xini = x0; + for ( i = 0; i < 10; i++ ) { + x0 = xini * stdlib_base_special_exp( x0*ainv ) * stdlib_base_special_pow( 1.0-( x0*ap1inv ), ainv ); + } + k = 1; + lgama = stdlib_base_special_gammaln( a ); + } + + logr = (1.0/a) * ( stdlib_base_special_ln(p) + stdlib_base_special_gammaln( a+1.0 ) ); + if ( ( logr < stdlib_base_special_ln( ONEO5 * ( 1.0+a ) ) ) && ( k == 0 ) ) { + r = stdlib_base_special_exp( logr ); + m = 0; + a2 = a * a; + a3 = a2 * a; + a4 = a3 * a; + ap1 = a + 1.0; + ap12 = ap1 * ap1; + ap13 = ap1 * ap12; + ap14 = ap12 * ap12; + ap2 = a + 2.0; + ap22 = ap2 * ap2; + ap3 = a + 3.0; + CK[ 0 ] = 1.0; + CK[ 1 ] = 1.0 / ap1; + CK[ 2 ] = HALF * ( ( 3.0*a ) + 5.0 ) / ( ap12*ap2 ); + CK[ 3 ] = ONEO3 * ( 31.0 + (8.0*a2) + (33.0*a) ) / ( ap13*ap2*ap3 ); + CK[ 4 ] = ( 2888.0 + (1179.0*a3) + (125.0*a4) + (3971.0*a2) + (5661.0*a) ) / ( 24.0 * ap14*ap22*ap3*( a+4.0 ) ); + x0 = r * evalpoly( CK, 5, r ); + lgama = stdlib_base_special_gammaln( a ); + k = 1; + } + if ( ( a < 10.0 ) && ( k == 0 ) ) { + vgam = stdlib_base_special_sqrt( a ) / ( gamstar(a)*STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI ); + vmin = ( 0.02 < vgam ) ? 0.02 : vgam; // min( 0.02, vgam ) + if ( q < vmin ) { + m = 0; + b = 1.0 - a; + b2 = b * b; + b3 = b2 * b; + eta = stdlib_base_special_sqrt( -2.0/a * stdlib_base_special_ln( q/vgam ) ); + x0 = a * lambdaeta( eta ); + L = stdlib_base_special_ln( x0 ); + if ( x0 > 5.0 ) { + L2 = L * L; + L3 = L2 * L; + L4 = L3 * L; + r = 1.0 / x0; + CK[ 0 ] = L - 1.0; + CK[ 1 ] = ( (3.0*b) - (2.0*b*L) + L2 - ( 2.0*L ) + 2.0 ) * HALF; + CK[ 2 ] =( (24.0*b*L) - (11.0*b2) - (24.0*b) - (6.0*L2) + (12.0*L) - 12.0 - (9.0*b*L2) + (6.0*b2*L) + (2.0*L3) ) * ONEO6; + CK[ 3 ] = ( (-12.0*b3*L) + (8.04*b*L2) - (114.0*b2*L) + (72.0+(36.0*L2)) + (((3.0*L4)-(72.0*L)+162.0) * (b-(168.0*b*L))) - ((12.0*L3)+(25.0*b3)) - ( (22.0*b*L3)+(36.0*b2*L2)+(120.0*b2) ) ) * ONEO12; + CK[ 4 ] = 0.0; + x0 = x0 - L + ( b*r*evalpoly( CK, 5, r ) ); + } else { + r = 1.0 / x0; + double ck = L - 1.0; + t = L - (b*r*ck); + if ( t < x0 ) { + x0 -= t; + } + } + lgama = stdlib_base_special_gammaln( a ); + k = 1; + } + } + if ( ( stdlib_base_special_abs( porq-HALF ) < 1.0e-5 ) && ( k == 0 ) ) { + m = 0; + ainv = 1.0 / a; + x0 = a - ONEO3 + ( ( 0.0197530864197530864197530864198 + + ( 0.00721144424848128551832255535959*ainv ) ) * ainv ); + lgama = stdlib_base_special_gammaln( a ); + k = 1; + } + if ( ( a < 1.0 ) && ( k == 0 ) ) { + m = 0; + if ( pcase ) { + x0 = stdlib_base_special_exp( (1.0/a) * ( stdlib_base_special_ln(porq) + stdlib_base_special_gammaln(a+1.0) ) ); + } else { + x0 = stdlib_base_special_exp( (1.0/a) * ( stdlib_base_special_ln(1.0-porq) + stdlib_base_special_gammaln(a+1.0) ) ); + } + lgama = stdlib_base_special_gammaln( a ); + k = 1; + } + if ( k == 0 ) { + m = 1; + ainv = 1.0 / a; + r = stdlib_base_special_erfcinv( 2.0 * porq ); + eta = s * r / stdlib_base_special_sqrt( a*HALF ); + if ( r < STDLIB_CONSTANT_FLOAT32_MAX ) { + eta += ( eps1(eta) + ( (eps2(eta)+(eps3(eta)*ainv))*ainv ) ) * ainv; + x0 = a * lambdaeta(eta); + y = eta; + fp = -stdlib_base_special_sqrt( a/STDLIB_CONSTANT_FLOAT64_TWO_PI ) * stdlib_base_special_exp( -HALF*a*y*y ) / ( gamstar(a) ); + invfp = 1.0 / fp; + } else { + // Warning: Overflow... + return 0.0/0.0; // NaN + } + } + if ( k < 2 ) { + xr = higherNewton( x0, a, m, p, q, lgama, invfp, pcase ); + } + return xr; +} + + +// MAIN // + +/** +* Inverts the lower incomplete gamma function. +* +* @param p - probability value +* @param a - scale parameter +* @param upper - boolean indicating if the function should invert the upper tail +* @return function value +* +* @example +* double y = stdlib_base_special_gammaincinv( 0.5, 1.0, false ); +* // returns ~0.693 +*/ +double stdlib_base_special_gammaincinv( const double p, const double a, const bool upper ) { + if ( stdlib_base_is_nan( p ) || stdlib_base_is_nan( a ) ) { + return 0.0/0.0; // NaN + } + if ( a < STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL ) { + return 0.0/0.0; // NaN + } + if ( p > 1.0 || p < 0.0 ) { + return 0.0/0.0; // NaN + } + // Case: invert upper gamma function + if ( upper == true ) { + if ( p == 0.0 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( p == 1.0 ) { + return 0.0; + } + return compute( a, 1.0-p, p ); + } + // Default: invert lower gamma function + if ( p == 0.0 ) { + return 0.0; + } + if ( p == 1.0 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + return compute( a, p, 1.0-p ); +} diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/stdlib_base_special_gammaincinv.h b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/stdlib_base_special_gammaincinv.h new file mode 100644 index 000000000000..a83561eea061 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/stdlib_base_special_gammaincinv.h @@ -0,0 +1,26 @@ +/** +* @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. +*/ + +#ifndef STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H +#define STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H + +// MODULES // + +#include "stdlib/math/base/special/gammaincinv.h" + +#endif // !STDLIB_MATH_BASE_SPECIAL_GAMMAINCINV_H diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.native.js new file mode 100644 index 000000000000..22dd3201622a --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.native.js @@ -0,0 +1,92 @@ +/** +* @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. +*/ + +'use strict'; + +// MODULES // + +var tape = require( 'tape' ); +var IS_NODE = require( '@stdlib/assert/is-node' ); +var IS_LITTLE_ENDIAN = require( '@stdlib/assert/is-little-endian' ); +var gammaincinv = require( './../lib/native.js' ); + + +// VARIABLES // + +var FLG = IS_NODE && ( process.platform !== 'win32' || process.arch !== 'ia32' ); + + +// FIXTURES // + +var arg1 = require( './fixtures/r/arg1.json' ); +var arg2 = require( './fixtures/r/arg2.json' ); +var expected_lower = require( './fixtures/r/expected.json' ); +var upperArg1 = require( './fixtures/r/upper_arg1.json' ); +var upperArg2 = require( './fixtures/r/upper_arg2.json' ); +var expected_upper = require( './fixtures/r/upper_expected.json' ); + + +// TESTS // + +tape( 'main export is a function', function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof gammaincinv, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'if node and not Windows 32-bit, the function is the native function', function test( t ) { + if ( FLG && IS_LITTLE_ENDIAN ) { + t.strictEqual( gammaincinv.name, 'gammaincinv', 'function name is `gammaincinv`' ); + } else { + t.ok( true, 'not node or is Windows 32-bit' ); + } + t.end(); +}); + +tape( 'the function computes the inverse incomplete gamma function (lower)', function test( t ) { + var expected; + var actual; + var p; + var a; + var i; + + p = arg1; + a = arg2; + + for ( i = 0; i < p.length; i++ ) { + actual = gammaincinv( p[i], a[i], false ); + t.inDelta( actual, expected_lower[i], 1e-12, 'p: ' + p[i] + ', a: ' + a[i] + '. Actual: ' + actual + '. Expected: ' + expected_lower[i] + '.' ); + } + t.end(); +}); + +tape( 'the function computes the inverse incomplete gamma function (upper)', function test( t ) { + var actual; + var p; + var a; + var i; + + p = upperArg1; + a = upperArg2; + + for ( i = 0; i < p.length; i++ ) { + actual = gammaincinv( p[i], a[i], true ); + t.inDelta( actual, expected_upper[i], 1e-12, 'p: ' + p[i] + ', a: ' + a[i] + '. Actual: ' + actual + '. Expected: ' + expected_upper[i] + '.' ); + } + t.end(); +});