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..64c2dc236129 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md @@ -2,7 +2,7 @@ @license Apache-2.0 -Copyright (c) 2018 The Stdlib Authors. +Copyright (c) 2026 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. @@ -70,7 +70,7 @@ Again, for given `p` and `a` the function returns the `x` which satisfies `p = Q var gammaincinv = require( '@stdlib/math/base/special/gammaincinv' ); ``` -#### gammaincinv( p, s\[, upper ] ) +#### gammaincinv( p, a\[, upper ] ) Inverts the regularized incomplete gamma function. Contrary to the more commonly used definition, in this implementation the first argument is the probability `p` and the second argument is the scale factor `a`. By default, the function inverts the _lower_ regularized incomplete gamma function, `P(x,a)`. To invert the _upper_ function instead, i.e. `Q(x,a)`, set the `upper` argument to `true`. @@ -136,6 +136,92 @@ logEachMap( 'p: %0.4f, \t a: %0.4f, \t P^(-1)(p, a): %0.4f', p, a, gammaincinv ) + +
+ +## C API + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/gammaincinv.h" +``` + +#### stdlib_base_gammaincinv( p, a, upper ) + +Inverts the regularized incomplete gamma function. + +```c +double y = stdlib_base_gammaincinv( 0.5, 2.0, false ); +// returns ~1.678 + +y = stdlib_base_gammaincinv( 0.1, 10.0, false ); +// returns ~6.221 + +y = stdlib_base_gammaincinv( 0.75, 3.0, true ); +// returns ~1.727 +``` + +The function accepts the following arguments: + +- `p`: `[in] double` input probability. +- `a`: `[in] double` scale parameter. +- `upper`: `[in] bool` boolean indicating if the function should invert the upper tail of the incomplete gamma function. + +```c +double stdlib_base_gammaincinv( const double p, const double a, const bool upper ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/gammaincinv.h" +#include +#include + +int main( void ) { + double p; + double a; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + p = (double)rand() / (double)RAND_MAX; + a = ( ( (double)rand() / (double)RAND_MAX ) * 20.0 ); + y = stdlib_base_gammaincinv( p, a, false ); + printf( "p: %lf, a: %lf, y: %lf\n", p, a, y ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/Makefile new file mode 100644 index 000000000000..979768abbcec --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2026 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 + +# Define the program used for compiling C source files: +ifdef C_COMPILER + CC := $(C_COMPILER) +else + CC := gcc +endif + +# Define the command-line options when compiling C files: +CFLAGS ?= \ + -std=c99 \ + -O3 \ + -Wall \ + -pedantic + +# Determine whether to generate position independent code ([1][1], [2][2]). +# +# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options +# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option +ifeq ($(OS), WINNT) + fPIC ?= +else + fPIC ?= -fPIC +endif + +# List of includes (e.g., `-I /foo/bar -I /beep/boop/include`): +INCLUDE ?= + +# List of source files: +SOURCE_FILES ?= + +# List of libraries (e.g., `-lopenblas -lpthread`): +LIBRARIES ?= + +# List of library paths (e.g., `-L /foo/bar -L /beep/boop`): +LIBPATH ?= + +# List of C targets: +c_targets := benchmark.out + + +# RULES # + +#/ +# Compiles source files. +# +# @param {string} [C_COMPILER] - C compiler (e.g., `gcc`) +# @param {string} [CFLAGS] - C compiler options +# @param {(string|void)} [fPIC] - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} [INCLUDE] - list of includes (e.g., `-I /foo/bar -I /beep/boop/include`) +# @param {string} [SOURCE_FILES] - list of source files +# @param {string} [LIBPATH] - list of library paths (e.g., `-L /foo/bar -L /beep/boop`) +# @param {string} [LIBRARIES] - list of libraries (e.g., `-lopenblas -lpthread`) +# +# @example +# make +# +# @example +# make all +#/ +all: $(c_targets) + +.PHONY: all + +#/ +# Compiles C source files. +# +# @private +# @param {string} CC - C compiler (e.g., `gcc`) +# @param {string} CFLAGS - C compiler options +# @param {(string|void)} fPIC - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} INCLUDE - list of includes (e.g., `-I /foo/bar`) +# @param {string} SOURCE_FILES - list of source files +# @param {string} LIBPATH - list of library paths (e.g., `-L /foo/bar`) +# @param {string} LIBRARIES - list of libraries (e.g., `-lopenblas`) +#/ +$(c_targets): %.out: %.c + $(QUIET) $(CC) $(CFLAGS) $(fPIC) $(INCLUDE) -o $@ $(SOURCE_FILES) $< $(LIBPATH) -lm $(LIBRARIES) + +#/ +# Runs compiled benchmarks. +# +# @example +# make run +#/ +run: $(c_targets) + $(QUIET) ./$< + +.PHONY: run + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: + $(QUIET) -rm -f *.o *.out + +.PHONY: clean diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/benchmark.c new file mode 100644 index 000000000000..8c1bb8c7efd2 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/benchmark.c @@ -0,0 +1,44 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 +#include + +#define MAX_ITER 1000000 + +int main( void ) { + double p, a, y; + int i; + clock_t start, end; + double cpu_time_used; + + start = clock(); + for ( i = 0; i < MAX_ITER; i++ ) { + p = (double)rand() / (double)RAND_MAX; + a = ( ( (double)rand() / (double)RAND_MAX ) * 20.0 ); + y = stdlib_base_gammaincinv( p, a, false ); + if ( y != y ) { + printf("Error"); + } + } + end = clock(); + cpu_time_used = ( (double) (end - start) ) / CLOCKS_PER_SEC; + printf("Elapsed time: %f seconds\n", cpu_time_used); +} 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..0d6508a12e99 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/binding.gyp @@ -0,0 +1,170 @@ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. + +# A `.gyp` file for building a Node.js native add-on. +# +# [1]: https://gyp.gsrc.io/docs/InputFormatReference.md +# [2]: https://gyp.gsrc.io/docs/UserDocumentation.md +{ + # List of files to include in this file: + 'includes': [ + './include.gypi', + ], + + # Define variables to be used throughout the configuration for all targets: + 'variables': { + # Target name should match the add-on export name: + 'addon_target_name%': 'addon', + + # Set variables based on the host OS: + 'conditions': [ + [ + 'OS=="win"', + { + # Define the object file suffix: + 'obj': 'obj', + }, + { + # Define the object file suffix: + 'obj': 'o', + } + ], # end condition (OS=="win") + ], # end conditions + }, # end variables + + # Define compile targets: + 'targets': [ + + # Target to generate an add-on: + { + # The target name should match the add-on export name: + 'target_name': '<(addon_target_name)', + + # Define dependencies: + 'dependencies': [], + + # Define directories which contain relevant include headers: + 'include_dirs': [ + # Local include directory: + '<@(include_dirs)', + ], + + # List of source files: + 'sources': [ + '<@(src_files)', + ], + + # Settings which should be applied when a target's object files are used as linker input: + 'link_settings': { + # Define libraries: + 'libraries': [ + '<@(libraries)', + ], + + # Define library directories: + 'library_dirs': [ + '<@(library_dirs)', + ], + }, + + # C/C++ compiler flags: + 'cflags': [ + # Enable commonly used warning options: + '-Wall', + + # Aggressive optimization: + '-O3', + ], + + # C specific compiler flags: + 'cflags_c': [ + # Specify the C standard to which a program is expected to conform: + '-std=c99', + ], + + # C++ specific compiler flags: + 'cflags_cpp': [ + # Specify the C++ standard to which a program is expected to conform: + '-std=c++11', + ], + + # Linker flags: + 'ldflags': [], + + # Apply conditions based on the host OS: + 'conditions': [ + [ + 'OS=="mac"', + { + # Linker flags: + 'ldflags': [ + '-undefined dynamic_lookup', + '-Wl,-no-pie', + '-Wl,-search_paths_first', + ], + }, + ], # end condition (OS=="mac") + [ + 'OS!="win"', + { + # C/C++ flags: + 'cflags': [ + # Generate platform-independent code: + '-fPIC', + ], + }, + ], # end condition (OS!="win") + ], # end conditions + }, # end target <(addon_target_name) + + # Target to copy a generated add-on to a standard location: + { + 'target_name': 'copy_addon', + + # Declare that the output of this target is not linked: + 'type': 'none', + + # Define dependencies: + 'dependencies': [ + # Require that the add-on be generated before building this target: + '<(addon_target_name)', + ], + + # Define a list of actions: + 'actions': [ + { + 'action_name': 'copy_addon', + 'message': 'Copying addon...', + + # Explicitly list the inputs in the command-line invocation below: + 'inputs': [], + + # Declare the expected outputs: + 'outputs': [ + '<(addon_output_dir)/<(addon_target_name).node', + ], + + # Define the command-line invocation: + 'action': [ + 'cp', + '<(PRODUCT_DIR)/<(addon_target_name).node', + '<(addon_output_dir)/<(addon_target_name).node', + ], + }, + ], # end actions + }, # end target copy_addon + ], # end targets +} 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..c8f8e9a1517b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2026 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 + +# Define the program used for compiling C source files: +ifdef C_COMPILER + CC := $(C_COMPILER) +else + CC := gcc +endif + +# Define the command-line options when compiling C files: +CFLAGS ?= \ + -std=c99 \ + -O3 \ + -Wall \ + -pedantic + +# Determine whether to generate position independent code ([1][1], [2][2]). +# +# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options +# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option +ifeq ($(OS), WINNT) + fPIC ?= +else + fPIC ?= -fPIC +endif + +# List of includes (e.g., `-I /foo/bar -I /beep/boop/include`): +INCLUDE ?= + +# List of source files: +SOURCE_FILES ?= + +# List of libraries (e.g., `-lopenblas -lpthread`): +LIBRARIES ?= + +# List of library paths (e.g., `-L /foo/bar -L /beep/boop`): +LIBPATH ?= + +# List of C targets: +c_targets := example.out + + +# RULES # + +#/ +# Compiles source files. +# +# @param {string} [C_COMPILER] - C compiler (e.g., `gcc`) +# @param {string} [CFLAGS] - C compiler options +# @param {(string|void)} [fPIC] - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} [INCLUDE] - list of includes (e.g., `-I /foo/bar -I /beep/boop/include`) +# @param {string} [SOURCE_FILES] - list of source files +# @param {string} [LIBPATH] - list of library paths (e.g., `-L /foo/bar -L /beep/boop`) +# @param {string} [LIBRARIES] - list of libraries (e.g., `-lopenblas -lpthread`) +# +# @example +# make +# +# @example +# make all +#/ +all: $(c_targets) + +.PHONY: all + +#/ +# Compiles C source files. +# +# @private +# @param {string} CC - C compiler (e.g., `gcc`) +# @param {string} CFLAGS - C compiler options +# @param {(string|void)} fPIC - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} INCLUDE - list of includes (e.g., `-I /foo/bar`) +# @param {string} SOURCE_FILES - list of source files +# @param {string} LIBPATH - list of library paths (e.g., `-L /foo/bar`) +# @param {string} LIBRARIES - list of libraries (e.g., `-lopenblas`) +#/ +$(c_targets): %.out: %.c + $(QUIET) $(CC) $(CFLAGS) $(fPIC) $(INCLUDE) -o $@ $(SOURCE_FILES) $< $(LIBPATH) -lm $(LIBRARIES) + +#/ +# Runs compiled examples. +# +# @example +# make run +#/ +run: $(c_targets) + $(QUIET) ./$< + +.PHONY: run + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: + $(QUIET) -rm -f *.o *.out + +.PHONY: clean 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..d808b83a57f3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/example.c @@ -0,0 +1,35 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 ) { + double p; + double a; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + p = (double)rand() / (double)RAND_MAX; + a = ( ( (double)rand() / (double)RAND_MAX ) * 20.0 ); + y = stdlib_base_gammaincinv( p, a, false ); + printf( "p: %lf, a: %lf, y: %lf\n", p, a, y ); + } +} diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/include.gypi b/lib/node_modules/@stdlib/math/base/special/gammaincinv/include.gypi new file mode 100644 index 000000000000..bee8d41a2caf --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/include.gypi @@ -0,0 +1,53 @@ +# @license Apache-2.0 +# +# Copyright (c) 2026 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. + +# A GYP include file for building a Node.js native add-on. +# +# Main documentation: +# +# [1]: https://gyp.gsrc.io/docs/InputFormatReference.md +# [2]: https://gyp.gsrc.io/docs/UserDocumentation.md +{ + # Define variables to be used throughout the configuration for all targets: + 'variables': { + # Source directory: + 'src_dir': './src', + + # Include directories: + 'include_dirs': [ + ' + +/* +* If C++, prevent name mangling so that the compiler emits a binary file having undecorated names, thus mirroring the behavior of a C compiler. +*/ +#ifdef __cplusplus +extern "C" { +#endif + +/** +* Computes the inverse of the incomplete gamma function. +*/ +double stdlib_base_gammaincinv( const double p, const double a, const bool upper ); + +#ifdef __cplusplus +} +#endif + +#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..5e6ee0b99593 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/lib/native.js @@ -0,0 +1,47 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 addon = require( './../src/addon.node' ); + + +// MAIN // + +/** +* Inverts the lower gamma function; i.e., computes `xr` such that `P(a,xr) = p`. +* +* @private +* @param {Probability} p - probability value +* @param {number} a - scale parameter +* @param {boolean} [upper=false] - boolean indicating if the function should invert the upper tail of the incomplete gamma function instead; i.e., compute `xr` such that `Q(a,xr) = p`. +* @returns {number} function value of the inverse +*/ +function gammaincinv( p, a, upper ) { + if ( arguments.length < 3 ) { + upper = false; + } + return addon( p, a, upper ); +} + + +// EXPORTS // + +module.exports = gammaincinv; diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/manifest.json b/lib/node_modules/@stdlib/math/base/special/gammaincinv/manifest.json new file mode 100644 index 000000000000..823e44ef5362 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/manifest.json @@ -0,0 +1,121 @@ +{ + "options": { + "task": "build" + }, + "fields": [ + { + "field": "src", + "resolve": true, + "relative": true + }, + { + "field": "include", + "resolve": true, + "relative": true + }, + { + "field": "libraries", + "resolve": false, + "relative": false + }, + { + "field": "libpath", + "resolve": true, + "relative": false + } + ], + "confs": [ + { + "task": "build", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/napi/export", + "@stdlib/napi/argv", + "@stdlib/napi/argv-double", + "@stdlib/napi/argv-bool", + "@stdlib/napi/create-double", + "@stdlib/math/base/assert/is-nan", + "@stdlib/constants/float32/smallest-normal", + "@stdlib/constants/float64/pinf", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/gamma", + "@stdlib/math/base/special/sqrt", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/ln", + "@stdlib/constants/float64/sqrt-two-pi", + "@stdlib/constants/float32/max", + "@stdlib/constants/float64/two-pi", + "@stdlib/math/base/special/gammainc" + ] + }, + { + "task": "benchmark", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/assert/is-nan", + "@stdlib/constants/float32/smallest-normal", + "@stdlib/constants/float64/pinf", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/gamma", + "@stdlib/math/base/special/sqrt", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/ln", + "@stdlib/constants/float64/sqrt-two-pi", + "@stdlib/constants/float32/max", + "@stdlib/constants/float64/two-pi", + "@stdlib/math/base/special/gammainc" + ] + }, + { + "task": "examples", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/assert/is-nan", + "@stdlib/constants/float32/smallest-normal", + "@stdlib/constants/float64/pinf", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/gamma", + "@stdlib/math/base/special/sqrt", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/ln", + "@stdlib/constants/float64/sqrt-two-pi", + "@stdlib/constants/float32/max", + "@stdlib/constants/float64/two-pi", + "@stdlib/math/base/special/gammainc" + ] + } + ] +} 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..5344be32824e 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/package.json @@ -14,12 +14,15 @@ } ], "main": "./lib", + "gypfile": true, "directories": { "benchmark": "./benchmark", "doc": "./docs", "example": "./examples", + "include": "./include", "lib": "./lib", "scripts": "./scripts", + "src": "./src", "test": "./test" }, "types": "./docs/types", diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/Makefile b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/Makefile new file mode 100644 index 000000000000..2caf905cedbe --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/Makefile @@ -0,0 +1,70 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2026 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/gammaincinv/src/addon.c b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/addon.c new file mode 100644 index 000000000000..ab56f1350e8d --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/addon.c @@ -0,0 +1,43 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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/napi/argv.h" +#include "stdlib/napi/argv_double.h" +#include "stdlib/napi/argv_bool.h" +#include "stdlib/napi/create_double.h" +#include "stdlib/napi/export.h" +#include + +/** +* Receives JavaScript callback invocation data. +* +* @param env environment under which the function is invoked +* @param info callback data +* @return Node-API value +*/ +static napi_value addon( napi_env env, napi_callback_info info ) { + STDLIB_NAPI_ARGV( env, info, argv, argc, 3 ); + STDLIB_NAPI_ARGV_DOUBLE( env, p, argv, 0 ); + STDLIB_NAPI_ARGV_DOUBLE( env, a, argv, 1 ); + STDLIB_NAPI_ARGV_BOOL( env, upper, argv, 2 ); + STDLIB_NAPI_CREATE_DOUBLE( env, stdlib_base_gammaincinv( p, a, upper ), out ); + return out; +} + +STDLIB_NAPI_MODULE_EXPORT_FCN( addon ) 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..efc9018262ec --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c @@ -0,0 +1,596 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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/assert/is_nan.h" +#include "stdlib/constants/float32/smallest_normal.h" +#include "stdlib/constants/float64/pinf.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/min.h" +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/constants/float64/sqrt_two_pi.h" +#include "stdlib/constants/float32/max.h" +#include "stdlib/constants/float64/two_pi.h" +#include "stdlib/math/base/special/gammainc.h" +#include + +static const double PINF = STDLIB_CONSTANT_FLOAT64_PINF; +static const float SMALLEST_NORMAL_FLOAT32 = STDLIB_CONSTANT_FLOAT32_SMALLEST_NORMAL; +static const double SQRT_TWO_PI = STDLIB_CONSTANT_FLOAT64_SQRT_TWO_PI; +static const float MAX_FLOAT32 = STDLIB_CONSTANT_FLOAT32_MAX; +static const double TWO_PI = STDLIB_CONSTANT_FLOAT64_TWO_PI; + +// Forward declarations +static double compute( double a, double p, double q ); +static double higherNewton( double x0, double a, int m, double p, double q, double lgama, double invfp, bool pcase ); +static double lambdaeta( double eta ); +static double gamstar( double x ); +static double stirling( double x ); +static double eps1( double eta ); +static double eps2( double eta ); +static double eps3( double eta ); +static double evalpoly( const double c[], int n, double x ); +static double rateval( double eta ); +static double rateval1( double eta ); +static double rateval2( double eta ); +static double rateval3( double eta ); +static double rational1( double eta ); +static double rational2( double eta ); +static double rational3( double eta ); +static double rational4( double x ); +static double rational5( double x ); +static double polyvalAK1( double r ); +static double polyvalAK2( double r ); + +/** +* Implementation of the high order Newton-like method. +*/ +static double higherNewton( double x0, double a, int m, double p, double q, double lgama, double invfp, 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.0; + n = 1; + a2 = a * a; + xini = x0; + do { + x = x0; + x2 = x * x; + if ( m == 0 ) { + dlnr = ( ( 1.0 - a ) * stdlib_base_ln( x ) ) + x + lgama; + if ( dlnr > stdlib_base_ln(MAX_FLOAT32) ) { + return xini; + } + r = stdlib_base_exp( dlnr ); + } else { + r = -invfp * x; + } + if ( pcase ) { + px = stdlib_base_gammainc( x, a, true, false ); + ck0 = -r * ( px - p ); + } else { + qx = stdlib_base_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.0*x2) - (4.0*x*a) + (4.0*x) + (2.0*a2) - (3.0*a) + 1.0 ) / x2; + ck2 /= 6.0; + x0 = x + ( r * ( 1.0 + ( r * ( ck1 + (r*ck2) ) ) ) ); + } else { + x0 = x + r; + } + t = stdlib_base_abs( ( x/x0 ) - 1.0 ); + n += 1; + x = x0; + if ( x < 0.0 ) { + x = xini; + n = 100; + } + } while ( ( ( t > 2e-14 ) && ( n < 35 ) ) ); + xr = x > 0.0 ? x : 0.0; + return xr; +} + +/** +* Returns the positive number satisfying \( \eta^2/2=\lambda-1-\ln(\lambda) \) with \( \operatorname{sign}(\lambda-1)=\operatorname{sign}(\eta) \). +*/ +static double lambdaeta( double eta ) { + double L2, L3, L4, L5; + double la; + double L; + double q, r, s; + static const double THRESHOLD = 1.0e-8; + static const double ONEO12 = 0.08333333333333333; + static const double ONEO120 = 0.008333333333333333; + double AK[6]; + + s = eta * eta * 0.5; + if ( eta == 0.0 ) { + la = 0.0; + } + else if ( eta < -1.0 ) { + r = stdlib_base_exp( -1.0 - s ); + la = polyvalAK1( r ); + } + else if ( eta < 1.0 ) { + r = eta; + la = polyvalAK2( r ); + } + else { + r = 11.0 + s; + L = stdlib_base_ln( r ); + la = r + L; + r = 1.0 / r; + L2 = L * L; + L3 = L2 * L; + L4 = L3 * L; + L5 = L4 * L; + AK[0] = 1.0; + AK[1] = ( 2.0-L ) * 0.5; + AK[2] = ( ( -9.0*L ) + 6.0 + ( 2.0*L2 ) ) / 6.0; + AK[3] = -( (3.0*L3)+ (36.0*L) - (22.0*L2) - 12.0 ) * ONEO12; + AK[4] = ( 60.0 + (350.0*L2) - (300.0*L) - (125.0*L3) + (12.0*L4) ) / 60.0; + AK[5] = -(-120.0 - (274.0*L4) + (900.0*L) - (1700.0*L2) + (1125.0*L3) + (20.0*L5)) * ONEO120; + la += ( L * r * evalpoly(AK, 6, r) ); + } + if ( + ( eta > -3.5 && eta < -0.03 ) || + ( eta > 0.03 && eta < 40.0 ) + ) { + r = 1.0; + q = la; + do { + la = q * ( s+stdlib_base_ln(q) ) / ( q-1.0 ); + r = stdlib_base_abs( ( q/la ) - 1.0 ); + q = la; + } while ( r > THRESHOLD ); + } + return la; +} + +static double polyvalAK1( double r ) { + static const double ak1[] = { + 3.183098861837907e-01, 2.705633337418241e-02, 1.239335416533316e-03, + 2.583338384581978e-05, -3.101829763300332e-07, -1.588693874366961e-08, + -1.378786850299244e-09 + }; + return evalpoly(ak1, 7, r); +} + +static double polyvalAK2( double r ) { + static const double ak2[] = { + 1.0, 1.0, 0.3333333333333333, 0.0, -0.06666666666666667, + -0.03555555555555555, 0.030476190476190474, 0.021547619047619047, + -0.02230842230842231, -0.018617195767195767, 0.020296245962459623 + }; + return evalpoly(ak2, 11, r); +} + +/** +* Computes the regulated gamma function. +*/ +static double gamstar( double x ) { + if ( x >= 3.0 ) { + return stdlib_base_exp( stirling(x) ); + } + if ( x > 0.0 ) { + return stdlib_base_gamma(x) / ( stdlib_base_exp( -x + ( ( x-0.5 ) * stdlib_base_ln(x) ) ) * SQRT_TWO_PI ); + } + return MAX_FLOAT32; +} + +static double stirling( double x ) { + static const double C[] = { + 0.08333333333333333, -0.002777777777777778, 0.0007936507936507936, + -0.0002755731922398589, 0.0001309838779518041, -0.00008129533379969938 + }; + double y = 1.0 / x; + double y2 = y * y; + return evalpoly(C, 6, y2) * y; +} + +/** +* Evaluates the `eps1` function. +*/ +static double eps1( double eta ) { + double la; + if ( stdlib_base_abs( eta ) < 1.0 ) { + return rateval( eta ); + } + la = lambdaeta( eta ); + return stdlib_base_ln( eta / ( la - 1.0 ) ) / eta; +} + +static double rateval( double eta ) { + static const double ak0[] = { + -0.3333333333333333, -0.08333333333333333, 0.014814814814814815, + -0.0011574074074074074, -0.0003527336860670194, 0.0001787551440329218 + }; + static const double bk0[] = { + 1.0, -0.75, 0.3541666666666667, -0.10677083333333333, + 0.021085069444444443, -0.0028741319444444446 + }; + double p = evalpoly(ak0, 6, eta); + double q = evalpoly(bk0, 6, eta); + return p / q; +} + +/** +* Evaluates the `eps2` function. +*/ +static double eps2( double eta ) { + double lnmeta; + double x; + if ( eta < -5.0 ) { + x = eta * eta; + lnmeta = stdlib_base_ln( -eta ); + return ( 12.0 - x - ( 6.0*( lnmeta*lnmeta ) ) ) / ( 12.0*x*eta ); + } + if ( eta < -2.0 ) { + return rateval1( eta ); + } + if ( eta < 2.0 ) { + return rateval2( eta ); + } + if ( eta < 1000.0 ) { + return rateval3( eta ) / ( -12.0*eta ); + } + return -1.0 / ( 12.0 * eta ); +} + +static double rateval1( double eta ) { + static const double ak1[] = { + -0.001851851851851852, 0.003472222222222222, -0.0026455026455026454, + 0.0009902263374485596 + }; + static const double bk1[] = { 1.0, -1.5, 0.875, -0.21875 }; + double p = evalpoly(ak1, 4, eta); + double q = evalpoly(bk1, 4, eta); + return p / q; +} + +static double rateval2( double eta ) { + static const double ak2[] = { + -0.041666666666666664, 0.0, 0.00496031746031746, 0.0, + -0.0006231481481481482 + }; + static const double bk2[] = { 1.0, 0.0, -0.4166666666666667, 0.0, 0.05753968253968254 }; + double p = evalpoly(ak2, 5, eta); + double q = evalpoly(bk2, 5, eta); + return p / q; +} + +static double rateval3( double eta ) { + static const double ak3[] = { + 1.0, -0.15, 0.03, -0.007291666666666667, 0.0020416666666666668 + }; + return evalpoly(ak3, 5, 1.0/eta); +} + +/** +* Evaluates the `eps3` function. +*/ +static double eps3( double eta ) { + double x; + double y; + + if ( eta < -8.0 ) { + x = eta * eta; + y = stdlib_base_ln( -eta ) / eta; + return ( -30.0 + ( eta*y*( (6.0*x*y*y)-12.0+x ) ) ) / ( 12.0*eta*x*x ); + } + if ( eta < -4.0 ) { + return rational1( eta ) / ( eta*eta ); + } + if ( eta < -2.0 ) { + return rational2( eta ); + } + if ( eta < 2.0 ) { + return rational3( eta ); + } + if ( eta < 10.0 ) { + x = 1.0 / eta; + return rational4( x ) / ( eta*eta ); + } + if ( eta < 100.0 ) { + x = 1.0 / eta; + return rational5( x ) / ( eta*eta ); + } + return -stdlib_base_ln( eta ) / ( 12.0*eta*eta*eta ); +} + +static double rational1( double eta ) { + static const double ak4[] = { -0.0006494341563786008, 0.00022947209362139917 }; + static const double bk4[] = { 1.0, -1.0, 0.25 }; + double p = evalpoly(ak4, 2, eta); + double q = evalpoly(bk4, 3, eta); + return p / q; +} + +static double rational2( double eta ) { + static const double ak5[] = { + -0.0008618882909167117, -0.0007840392217200666, -0.0002990724803031902 + }; + static const double bk5[] = { 1.0, -1.5, 0.875, -0.21875 }; + double p = evalpoly(ak5, 3, eta); + double q = evalpoly(bk5, 4, eta); + return p / q; +} + +static double rational3( double eta ) { + static const double ak6[] = { + 0.005952380952380952, 0.0, -0.0011904761904761906, 0.0, + 0.00024195767195767195 + }; + static const double bk6[] = { + 1.0, 0.0, -0.8333333333333334, 0.0, 0.20634920634920634 + }; + double p = evalpoly(ak6, 5, eta); + double q = evalpoly(bk6, 5, eta); + return p / q; +} + +static double rational4( double x ) { + static const double ak7[] = { + -0.0005313079364639922, 0.0005921664373536939, -0.0002708782096718045 + }; + static const double bk7[] = { 1.0, 0.25, 0.025, 0.0015625 }; + double p = evalpoly(ak7, 3, x); + double q = evalpoly(bk7, 4, x); + return p / q; +} + +static double rational5( double x ) { + static const double ak8[] = { + -0.00034436760689237765, 0.00005171790908260592, 0.00033493161081142234 + }; + static const double bk8[] = { 1.0, 0.15, 0.009375, -0.0001953125 }; + double p = evalpoly(ak8, 3, x); + double q = evalpoly(bk8, 4, x); + return p / q; +} + +static double evalpoly( const double c[], int n, double x ) { + double r = c[n-1]; + for (int i = n - 2; i >= 0; --i) { + r = r * x + c[i]; + } + return r; +} + +/** +* Computes `x` in the equations `P(a,xr) = p` and `Q(a,xr) = q`, where `a` is a positive parameter and `p` and `q` satisfy `p+q = 1`. +*/ +static double compute( double a, double p, double q ) { + double ap1inv; + double invfp; + double lgama; + bool pcase; + double porq; + double ainv; + double logr; + double ap22, ap14, ap13, ap12; + double vgam, vmin, xini, ap1, ap2, ap3, eta, p6, p5, x0, a2, L2, L3, L4, b2, b3, p3, a4, fp, p4, p2, a3, xr, b, L; + int i, k, m; + double r, s, t, y; + double CK[ 5 ]; + + invfp = 0.0; + + if ( p < 0.5 ) { + pcase = true; + porq = p; + s = -1.0; + } else { + pcase = false; + porq = q; + s = 1.0; + } + k = 0; + if ( stdlib_base_abs( a-1.0 ) < 1.0e-4 ) { + m = 0; + if ( pcase ) { + if ( p < 1.0e-3 ) { + p2 = p * p; + p3 = p2 * p; + p4 = p3 * p; + p5 = p4 * p; + p6 = p5 * p; + x0 = p + ( p2*0.5 ) + ( p3*(1.0/3.0) ) + ( p4*0.25 ) + ( p5*0.2 ) + ( p6*(1.0/6.0) ); + } else { + x0 = -stdlib_base_ln( 1.0-p ); + } + } else { + x0 = -stdlib_base_ln( q ); + } + if ( a == 1.0 ) { + k = 2; + xr = x0; + } else { + lgama = stdlib_base_gammaln( a ); + k = 1; + } + } + if ( q < 1.0e-30 && a < 0.5 ) { + m = 0; + x0 = -stdlib_base_ln( q*stdlib_base_gamma(a) ) + ( ( a-1.0 ) * stdlib_base_ln( -stdlib_base_ln( q*stdlib_base_gamma(a) ) ) ); + k = 1; + lgama = stdlib_base_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_gammaln( a+1.0 ) + stdlib_base_ln( p ) ) * ainv; + x0 = stdlib_base_exp( x0 ); + xini = x0; + for ( i = 0; i < 10; i++ ) { + x0 = xini * stdlib_base_exp( x0*ainv ) * stdlib_base_pow( 1.0-( x0*ap1inv ), ainv ); + } + k = 1; + lgama = stdlib_base_gammaln( a ); + } + + logr = (1.0/a) * ( stdlib_base_ln(p) + stdlib_base_gammaln( a+1.0 ) ); + if ( ( logr < stdlib_base_ln( 0.2 * ( 1.0+a ) ) ) && ( k == 0 ) ) { + r = stdlib_base_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 ] = 0.5 * ( ( 3.0*a ) + 5.0 ) / ( ap12*ap2 ); + CK[ 3 ] = (1.0/3.0) * ( 31.0 + (8.0*a2) + (33.0*a) ) / ( ap13*ap2*ap3 ); + CK[ 4 ] = (1.0/24.0) * ( 2888.0 + (1179.0*a3) + (125.0*a4) + (3971.0*a2) + (5661.0*a) ) / ( ap14*ap22*ap3*( a+4.0 ) ); + x0 = r * evalpoly( CK, 5, r ); + lgama = stdlib_base_gammaln( a ); + k = 1; + } + if ( ( a < 10.0 ) && ( k == 0 ) ) { + vgam = stdlib_base_sqrt( a ) / ( gamstar(a)*SQRT_TWO_PI ); + vmin = stdlib_base_min( 0.02, vgam ); + if ( q < vmin ) { + m = 0; + b = 1.0 - a; + b2 = b * b; + b3 = b2 * b; + eta = stdlib_base_sqrt( -2.0/a * stdlib_base_ln( q/vgam ) ); + x0 = a * lambdaeta( eta ); + L = stdlib_base_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 ) * 0.5; + 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) ) * (1.0/6.0); + 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) ) ) * (1.0/12.0); + CK[ 4 ] = 0.0; + x0 = x0 - L + ( b*r*evalpoly( CK, 5, r ) ); + } else { + r = 1.0 / x0; + t = L - (b*r*(L-1.0)); + if ( t < x0 ) { + x0 -= t; + } + } + lgama = stdlib_base_gammaln( a ); + k = 1; + } + } + if ( ( stdlib_base_abs( porq-0.5 ) < 1.0e-5 ) && ( k == 0 ) ) { + m = 0; + ainv = 1.0 / a; + x0 = a - (1.0/3.0) + ( ( 0.019753086419753086 + + ( 0.0072114442484812855 * ainv ) ) * ainv ); + lgama = stdlib_base_gammaln( a ); + k = 1; + } + if ( ( a < 1.0 ) && ( k == 0 ) ) { + m = 0; + if ( pcase ) { + x0 = stdlib_base_exp( (1.0/a) * ( stdlib_base_ln(porq) + stdlib_base_gammaln(a+1.0) ) ); + } else { + x0 = stdlib_base_exp( (1.0/a) * ( stdlib_base_ln(1.0-porq) + stdlib_base_gammaln(a+1.0) ) ); + } + lgama = stdlib_base_gammaln( a ); + k = 1; + } + if ( k == 0 ) { + m = 1; + ainv = 1.0 / a; + r = stdlib_base_erfcinv( 2.0 * porq ); + eta = s * r / stdlib_base_sqrt( a*0.5 ); + if ( r < MAX_FLOAT32 ) { + eta += ( eps1(eta) + ( (eps2(eta)+(eps3(eta)*ainv))*ainv ) ) * ainv; + x0 = a * lambdaeta(eta); + y = eta; + fp = -stdlib_base_sqrt( a/TWO_PI ) * stdlib_base_exp( -0.5*a*y*y ) / ( gamstar(a) ); + invfp = 1.0 / fp; + } else { + return 0.0/0.0; + } + } + if ( k < 2 ) { + xr = higherNewton( x0, a, m, p, q, lgama, invfp, pcase ); + } + return xr; +} + +/** +* Inverts the lower gamma function; i.e., computes `xr` such that `P(a,xr) = p`. +*/ +double stdlib_base_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 < SMALLEST_NORMAL_FLOAT32 ) { + 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 ) { + if ( p == 0.0 ) { + return 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 PINF; + } + return compute( a, p, 1.0 - p ); +} 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..d69ad86d982b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.native.js @@ -0,0 +1,267 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 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 resolve = require( 'path' ).resolve; +var tape = require( 'tape' ); +var incrspace = require( '@stdlib/array/base/incrspace' ); +var isfinite = require( '@stdlib/math/base/assert/is-finite' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); + + +// VARIABLES // + +var gammaincinv = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( gammaincinv instanceof Error ) +}; + + +// FIXTURES // + +var arg1 = require( './fixtures/r/arg1.json' ); +var arg2 = require( './fixtures/r/arg2.json' ); +var expected = require( './fixtures/r/expected.json' ); +var smallArg1 = require( './fixtures/r/small_arg1.json' ); +var smallArg2 = require( './fixtures/r/small_arg2.json' ); +var smallExpected = require( './fixtures/r/small_expected.json' ); +var largeArg1 = require( './fixtures/r/large_arg1.json' ); +var largeArg2 = require( './fixtures/r/large_arg2.json' ); +var largeExpected = require( './fixtures/r/large_expected.json' ); +var upperArg1 = require( './fixtures/r/upper_arg1.json' ); +var upperArg2 = require( './fixtures/r/upper_arg2.json' ); +var upperExpected = require( './fixtures/r/upper_expected.json' ); + +var i; +var v; +for ( i = 0; i < expected.length; i++ ) { + v = expected[ i ]; + if ( v === 'Inf' ) { + expected[ i ] = PINF; + } + else if ( v === 'NaN' ) { + expected[ i ] = NaN; + } +} +for ( i = 0; i < smallExpected.length; i++ ) { + v = smallExpected[ i ]; + if ( v === 'Inf' ) { + smallExpected[ i ] = PINF; + } + else if ( v === 'NaN' ) { + smallExpected[ i ] = NaN; + } +} +for ( i = 0; i < largeExpected.length; i++ ) { + v = largeExpected[ i ]; + if ( v === 'Inf' ) { + largeExpected[ i ] = PINF; + } + else if ( v === 'NaN' ) { + largeExpected[ i ] = NaN; + } +} +for ( i = 0; i < upperExpected.length; i++ ) { + v = upperExpected[ i ]; + if ( v === 'Inf' ) { + upperExpected[ i ] = PINF; + } + else if ( v === 'NaN' ) { + upperExpected[ i ] = NaN; + } +} + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof gammaincinv, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'the function returns `NaN` if provided a `NaN`', opts, function test( t ) { + var val; + + val = gammaincinv( NaN, 2 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = gammaincinv( 0.5, NaN ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'the function returns `NaN` if provided `p` outside the interval `[0,1]`', opts, function test( t ) { + var val; + + val = gammaincinv( 1.2, 2 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = gammaincinv( -0.1, 2 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + t.end(); +}); + +tape( 'the function returns `0` if provided `p = 1`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 0, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 1, a[ i ] ); + t.notOk( val === 0, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `+Infinity` if provided `p = 0`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 0, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 0, a[ i ] ); + t.notOk( val === PINF, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `0` if provided `p = 0` when `upper = true`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 0, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 0, a[ i ], true ); + t.notOk( val === 0, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `+Infinity` if provided `p = 1` when `upper = true`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 0, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 1, a[ i ], true ); + t.notOk( val === PINF, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function inverts the lower incomplete gamma function', opts, function test( t ) { + var actual; + var b1; + var b2; + var i; + for ( i = 0; i < arg1.length; i++ ) { + actual = gammaincinv( arg1[ i ], arg2[ i ] ); + + b1 = isfinite( actual ); + b2 = isfinite( expected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b2) ? 'finite' : 'not finite' ) ); + + b1 = isnan( actual ); + b2 = isnan( expected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b1) ? '' : 'not' ) + ' NaN' ); + if ( !b1 || !b2 ) { + t.strictEqual( isAlmostSameValue( actual, expected[ i ], 40.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + expected[ i ] + '.' ); + } + } + t.end(); +}); + +tape( 'the function inverts the upper incomplete gamma function when `upper=true`', opts, function test( t ) { + var actual; + var b1; + var b2; + var i; + for ( i = 0; i < arg1.length; i++ ) { + actual = gammaincinv( upperArg1[ i ], upperArg2[ i ], true ); + + b1 = isfinite( actual ); + b2 = isfinite( upperExpected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b2) ? 'finite' : 'not finite' ) ); + + b1 = isnan( actual ); + b2 = isnan( expected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b1) ? '' : 'not' ) + ' NaN' ); + if ( !b1 || !b2 ) { + t.strictEqual( isAlmostSameValue( actual, upperExpected[ i ], 40.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + upperExpected[ i ] + '.' ); + } + } + t.end(); +}); + +tape( 'the function inverts the lower incomplete gamma function for `p` smaller than `1e-4`', opts, function test( t ) { + var actual; + var b1; + var b2; + var i; + for ( i = 0; i < smallArg2.length; i++ ) { + actual = gammaincinv( smallArg1[ i ], smallArg2[ i ] ); + + b1 = isfinite( actual ); + b2 = isfinite( smallExpected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b2) ? 'finite' : 'not finite' ) ); + + b1 = isnan( actual ); + b2 = isnan( smallExpected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b1) ? '' : 'not' ) + ' NaN' ); + if ( !b1 || !b2 ) { + t.strictEqual( isAlmostSameValue( actual, smallExpected[ i ], 500.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + smallExpected[ i ] + '.' ); + } + } + t.end(); +}); + +tape( 'the function inverts the lower incomplete gamma function for `p` larger than `0.9`', opts, function test( t ) { + var actual; + var b1; + var b2; + var i; + + for ( i = 0; i < largeArg1.length; i++ ) { + actual = gammaincinv( largeArg1[ i ], largeArg2[ i ] ); + + b1 = isfinite( actual ); + b2 = isfinite( largeExpected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b2) ? 'finite' : 'not finite' ) ); + + b1 = isnan( actual ); + b2 = isnan( largeExpected[ i ] ); + t.strictEqual( b1, b2, 'returned result is ' + ( (b1) ? '' : 'not' ) + ' NaN' ); + if ( !b1 || !b2 ) { + t.strictEqual( isAlmostSameValue( actual, largeExpected[ i ], 100.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + largeExpected[ i ] + '.' ); + } + } + t.end(); +});