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..f867574c65ca 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/README.md @@ -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,109 @@ logEachMap( 'p: %0.4f, \t a: %0.4f, \t P^(-1)(p, a): %0.4f', p, a, gammaincinv ) + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### 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 "stdlib/constants/float64/eps.h" +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + double p; + double a; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + p = random_uniform( 0.0, 1.0 ); + a = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 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/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/benchmark.native.js new file mode 100644 index 000000000000..3f103847e1af --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/benchmark.native.js @@ -0,0 +1,64 @@ +/** +* @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 bench = require( '@stdlib/bench' ); +var uniform = require( '@stdlib/random/array/uniform' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var format = require( '@stdlib/string/format' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var gammaincinv = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( gammaincinv instanceof Error ) +}; + + +// MAIN // + +bench( format( '%s::native', pkg ), opts, function benchmark( b ) { + var x; + var y; + var z; + var i; + + x = uniform( 100, 0.0, 1.0 ); + y = uniform( 100, 0.1, 1000.0 ); + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + z = gammaincinv( x[ i%x.length ], y[ i%y.length ] ); + if ( isnan( z ) ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( isnan( z ) ) { + 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/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..966b86ae5d77 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/benchmark/c/native/benchmark.c @@ -0,0 +1,142 @@ +/** +* @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/constants/float64/eps.h" +#include +#include +#include +#include +#include +#include + +#define NAME "gammaincinv" +#define ITERATIONS 1000000 +#define REPEATS 3 + +/** +* Prints the TAP version. +*/ +static void print_version( void ) { + printf( "TAP version 13\n" ); +} + +/** +* Prints the TAP summary. +* +* @param total total number of tests +* @param passing total number of passing tests +*/ +static void print_summary( int total, int passing ) { + printf( "#\n" ); + printf( "1..%d\n", total ); // TAP plan + printf( "# total %d\n", total ); + printf( "# pass %d\n", passing ); + printf( "#\n" ); + printf( "# ok\n" ); +} + +/** +* Prints benchmarks results. +* +* @param elapsed elapsed time in seconds +*/ +static void print_results( double elapsed ) { + double rate = (double)ITERATIONS / elapsed; + printf( " ---\n" ); + printf( " iterations: %d\n", ITERATIONS ); + printf( " elapsed: %0.9f\n", elapsed ); + printf( " rate: %0.9f\n", rate ); + printf( " ...\n" ); +} + +/** +* Returns a clock time. +* +* @return clock time +*/ +static double tic( void ) { + struct timeval now; + gettimeofday( &now, NULL ); + return (double)now.tv_sec + (double)now.tv_usec/1.0e6; +} + +/** +* Generates a random number on the interval [min,max). +* +* @param min minimum value (inclusive) +* @param max maximum value (exclusive) +* @return random number +*/ +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +/** +* Runs a benchmark. +* +* @return elapsed time in seconds +*/ +static double benchmark( void ) { + double elapsed; + double p[ 100 ]; + double a[ 100 ]; + double y; + double t; + int i; + + for ( i = 0; i < 100; i++ ) { + p[ i ] = random_uniform( 0.0, 1.0 ); + a[ i ] = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 20.0 ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_gammaincinv( p[ i%100 ], a[ i%100 ], false ); + if ( y != y ) { + printf( "should not return NaN\n" ); + break; + } + } + elapsed = tic() - t; + if ( y != y ) { + printf( "should not return NaN\n" ); + } + return elapsed; +} + +/** +* Main execution sequence. +*/ +int main( void ) { + double elapsed; + int i; + + // Use the current time to seed the random number generator: + srand( time( NULL ) ); + + print_version(); + for ( i = 0; i < REPEATS; i++ ) { + printf( "# c::native::%s\n", NAME ); + elapsed = benchmark(); + print_results( elapsed ); + printf( "ok %d benchmark finished\n", i+1 ); + } + print_summary( REPEATS, REPEATS ); +} 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..31f312581a61 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/examples/c/example.c @@ -0,0 +1,41 @@ +/** +* @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/constants/float64/eps.h" +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + double p; + double a; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + p = random_uniform( 0.0, 1.0 ); + a = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 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..d4367104c552 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/lib/native.js @@ -0,0 +1,51 @@ +/** +* @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 +* +* @example +* var v = gammaincinv( 0.5, 2.0 ); +* // returns ~1.678 +*/ +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..ac2aae108470 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/manifest.json @@ -0,0 +1,126 @@ +{ + "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/ln-sqrt-two-pi", + "@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/ln-sqrt-two-pi", + "@stdlib/constants/float64/sqrt-two-pi", + "@stdlib/constants/float32/max", + "@stdlib/constants/float64/two-pi", + "@stdlib/constants/float64/eps", + "@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/ln-sqrt-two-pi", + "@stdlib/constants/float64/sqrt-two-pi", + "@stdlib/constants/float32/max", + "@stdlib/constants/float64/two-pi", + "@stdlib/constants/float64/eps", + "@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/scripts/evalpoly.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalpoly.js index 7f27b9817f0e..b80a7af42895 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalpoly.js +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalpoly.js @@ -24,10 +24,15 @@ // MODULES // var resolve = require( 'path' ).resolve; +var readFileSync = require( '@stdlib/fs/read-file' ).sync; var writeFileSync = require( '@stdlib/fs/write-file' ).sync; var currentYear = require( '@stdlib/time/current-year' ); +var substringBefore = require( '@stdlib/string/substring-before' ); +var substringAfter = require( '@stdlib/string/substring-after' ); +var format = require( '@stdlib/string/format' ); var licenseHeader = require( '@stdlib/_tools/licenses/header' ); var compile = require( '@stdlib/math/base/tools/evalpoly-compile' ); +var compileC = require( '@stdlib/math/base/tools/evalpoly-compile-c' ); // VARIABLES // @@ -74,6 +79,33 @@ var header = licenseHeader( 'Apache-2.0', 'js', { header += '\n/* This is a generated file. Do not edit directly. */\n'; +// FUNCTIONS // + +/** +* Inserts a compiled function into file content. +* +* @private +* @param {string} text - source content +* @param {string} id - function identifier +* @param {string} str - function string +* @returns {string} updated content +*/ +function insert( text, id, str ) { + var before; + var after; + var begin; + var end; + + begin = '// BEGIN: '+id; + end = '// END: '+id; + + before = substringBefore( text, begin ); + after = substringAfter( text, end ); + + return format( '%s// BEGIN: %s\n\n%s\n%s%s', before, id, str, end, after ); +} + + // MAIN // /** @@ -83,7 +115,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; */ function main() { var fpath; + var copts; var opts; + var file; var str; opts = { @@ -105,6 +139,32 @@ function main() { fpath = resolve( __dirname, '..', 'lib', 'polyval_ak2.js' ); str = header + compile( AK2 ); writeFileSync( fpath, str, opts ); + + copts = { + 'dtype': 'double', + 'name': '' + }; + + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + copts.name = 'polyval_c'; + str = compileC( C, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_d'; + str = compileC( D, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_ak1'; + str = compileC( AK1, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_ak2'; + str = compileC( AK2, copts ); + file = insert( file, copts.name, str ); + + writeFileSync( fpath, file, opts ); } main(); diff --git a/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalrational.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalrational.js index 670b90cef310..bed682fe5e63 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalrational.js +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/scripts/evalrational.js @@ -24,10 +24,15 @@ // MODULES // var resolve = require( 'path' ).resolve; +var readFileSync = require( '@stdlib/fs/read-file' ).sync; var writeFileSync = require( '@stdlib/fs/write-file' ).sync; var currentYear = require( '@stdlib/time/current-year' ); +var substringBefore = require( '@stdlib/string/substring-before' ); +var substringAfter = require( '@stdlib/string/substring-after' ); +var format = require( '@stdlib/string/format' ); var licenseHeader = require( '@stdlib/_tools/licenses/header' ); var compile = require( '@stdlib/math/base/tools/evalrational-compile' ); +var compileC = require( '@stdlib/math/base/tools/evalrational-compile-c' ); // VARIABLES // @@ -176,6 +181,33 @@ var header = licenseHeader( 'Apache-2.0', 'js', { header += '\n/* This is a generated file. Do not edit directly. */\n'; +// FUNCTIONS // + +/** +* Inserts a compiled function into file content. +* +* @private +* @param {string} text - source content +* @param {string} id - function identifier +* @param {string} str - function string +* @returns {string} updated content +*/ +function insert( text, id, str ) { + var before; + var after; + var begin; + var end; + + begin = '// BEGIN: '+id; + end = '// END: '+id; + + before = substringBefore( text, begin ); + after = substringAfter( text, end ); + + return format( '%s// BEGIN: %s\n\n%s\n%s%s', before, id, str, end, after ); +} + + // MAIN // /** @@ -185,7 +217,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; */ function main() { var fpath; + var copts; var opts; + var file; var str; opts = { @@ -227,6 +261,52 @@ function main() { fpath = resolve( __dirname, '..', 'lib', 'rational_ak8bk8.js' ); str = header + compile( AK8, BK8 ); writeFileSync( fpath, str, opts ); + + copts = { + 'dtype': 'double', + 'name': '' + }; + + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + copts.name = 'rational_ak0bk0'; + str = compileC( AK0, BK0, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak1bk1'; + str = compileC( AK1, BK1, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak2bk2'; + str = compileC( AK2, BK2, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak3bk3'; + str = compileC( AK3, BK3, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak4bk4'; + str = compileC( AK4, BK4, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak5bk5'; + str = compileC( AK5, BK5, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak6bk6'; + str = compileC( AK6, BK6, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak7bk7'; + str = compileC( AK7, BK7, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'rational_ak8bk8'; + str = compileC( AK8, BK8, copts ); + file = insert( file, copts.name, str ); + + writeFileSync( fpath, file, opts ); } main(); 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..2ef048ff8fbe --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/src/main.c @@ -0,0 +1,1111 @@ +/** +* @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. +* +* +* ## Notice +* +* The implementation was translated from the Fortran module by +* +* ```text +* ---------------------------------------------------------------------- +* Authors: +* Amparo Gil (U. Cantabria, Santander, Spain) +* e-mail: amparo.gil@unican.es +* Javier Segura (U. Cantabria, Santander, Spain) +* e-mail: javier.segura@unican.es +* Nico M. Temme (CWI, Amsterdam, The Netherlands) +* e-mail: nico.temme@cwi.nl +* --------------------------------------------------------------------- +* ``` +*/ + +#include "stdlib/math/base/special/gammaincinv.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/special/gammainc.h" +#include "stdlib/math/base/special/erfcinv.h" +#include "stdlib/math/base/special/gammaln.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/float32/smallest_normal.h" +#include "stdlib/constants/float64/ln_sqrt_two_pi.h" +#include "stdlib/constants/float64/sqrt_two_pi.h" +#include "stdlib/constants/float64/two_pi.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float32/max.h" +#include +#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; +static const double LN_SQRT_TWO_PI = STDLIB_CONSTANT_FLOAT64_LN_SQRT_TWO_PI; + +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; +static const double ONEO24 = 0.0416666666666666666666666666667; +static const double ONEO120 = 0.00833333333333333333333333333333; +static const double THRESHOLD = 1.0e-8; +static const double C6 = 0.30865217988013567769; + +// Chebyshev polynomial coefficients... +static const double A[ 18 ] = { + 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 +}; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_ak1 + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_ak1( const double x ) { + return 0.0 + (x * (1.0 + (x * (1.0 + (x * (1.5 + (x * (2.6666666666666665 + (x * (5.208333333333333 + (x * 10.8))))))))))); +} + +// END: polyval_ak1 + +// BEGIN: polyval_ak2 + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_ak2( const double x ) { + return 1.0 + (x * (1.0 + (x * (0.3333333333333333 + (x * (0.027777777777777776 + (x * (-0.003703703703703704 + (x * (0.0002314814814814815 + (x * 0.00005878894767783657))))))))))); +} + +// END: polyval_ak2 + +// BEGIN: polyval_c + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_c( const double x ) { + return 0.025721014990011306 + (x * (0.08247596616699963 + (x * (-0.0025328157302663564 + (x * (0.0006099292666946337 + (x * (-0.00033543297638406 + (x * 0.000250505279903))))))))); +} + +// END: polyval_c + +// BEGIN: polyval_d + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_d( const double x ) { + return 0.08333333333333333 + (x * (-0.002777777777777778 + (x * (0.0007936507936507937 + (x * -0.0005952380952380953))))); +} + +// END: polyval_d + +// BEGIN: rational_ak0bk0 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak0bk0( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.3333333333438; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.3333333333438 + (x * (-0.2070740359969 + (x * (-0.05041806657154 + (x * (-0.004923635739372 + (x * -0.00004293658292782))))))); + s2 = 1.0 + (x * (0.7045554412463 + (x * (0.2118190062224 + (x * (0.03048648397436 + (x * 0.001605037988091))))))); + } else { + ix = 1.0 / x; + s1 = -0.00004293658292782 + (ix * (-0.004923635739372 + (ix * (-0.05041806657154 + (ix * (-0.2070740359969 + (ix * -0.3333333333438))))))); + s2 = 0.001605037988091 + (ix * (0.03048648397436 + (ix * (0.2118190062224 + (ix * (0.7045554412463 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak0bk0 + +// BEGIN: rational_ak1bk1 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak1bk1( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0172847633523; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0172847633523 + (x * (-0.0159372646475 + (x * (-0.00464910887221 + (x * (-0.00060683488776 + (x * -0.00000614830384279))))))); + s2 = 1.0 + (x * (0.764050615669 + (x * (0.297143406325 + (x * (0.0579490176079 + (x * 0.00574558524851))))))); + } else { + ix = 1.0 / x; + s1 = -0.00000614830384279 + (ix * (-0.00060683488776 + (ix * (-0.00464910887221 + (ix * (-0.0159372646475 + (ix * -0.0172847633523))))))); + s2 = 0.00574558524851 + (ix * (0.0579490176079 + (ix * (0.297143406325 + (ix * (0.764050615669 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak1bk1 + +// BEGIN: rational_ak2bk2 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak2bk2( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0172839517431; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0172839517431 + (x * (-0.0146362417966 + (x * (-0.00357406772616 + (x * (-0.000391032032692 + (x * 0.00000249634036069))))))); + s2 = 1.0 + (x * (0.690560400696 + (x * (0.249962384741 + (x * (0.0443843438769 + (x * 0.00424073217211))))))); + } else { + ix = 1.0 / x; + s1 = 0.00000249634036069 + (ix * (-0.000391032032692 + (ix * (-0.00357406772616 + (ix * (-0.0146362417966 + (ix * -0.0172839517431))))))); + s2 = 0.00424073217211 + (ix * (0.0443843438769 + (ix * (0.249962384741 + (ix * (0.690560400696 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak2bk2 + +// BEGIN: rational_ak3bk3 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak3bk3( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.99994466948; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.99994466948 + (x * (104.649839762 + (x * (857.204033806 + (x * (731.901559577 + (x * 45.5174411671))))))); + s2 = 1.0 + (x * (104.526456943 + (x * (823.313447808 + (x * (3119.93802124 + (x * 3970.03311219))))))); + } else { + ix = 1.0 / x; + s1 = 45.5174411671 + (ix * (731.901559577 + (ix * (857.204033806 + (ix * (104.649839762 + (ix * 0.99994466948))))))); + s2 = 3970.03311219 + (ix * (3119.93802124 + (ix * (823.313447808 + (ix * (104.526456943 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak3bk3 + +// BEGIN: rational_ak4bk4 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak4bk4( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.0495346498136; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.0495346498136 + (x * (0.0299521337141 + (x * (0.00688296911516 + (x * (0.000512634846317 + (x * -0.0000201411722031))))))); + s2 = 1.0 + (x * (0.759803615283 + (x * (0.261547111595 + (x * (0.0464854522477 + (x * 0.00403751193496))))))); + } else { + ix = 1.0 / x; + s1 = -0.0000201411722031 + (ix * (0.000512634846317 + (ix * (0.00688296911516 + (ix * (0.0299521337141 + (ix * 0.0495346498136))))))); + s2 = 0.00403751193496 + (ix * (0.0464854522477 + (ix * (0.261547111595 + (ix * (0.759803615283 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak4bk4 + +// BEGIN: rational_ak5bk5 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak5bk5( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.00452313583942; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.00452313583942 + (x * (0.00120744920113 + (x * (-0.0000789724156582 + (x * (-0.0000504476066942 + (x * -0.00000535770949796))))))); + s2 = 1.0 + (x * (0.912203410349 + (x * (0.405368773071 + (x * (0.0901638932349 + (x * 0.00948935714996))))))); + } else { + ix = 1.0 / x; + s1 = -0.00000535770949796 + (ix * (-0.0000504476066942 + (ix * (-0.0000789724156582 + (ix * (0.00120744920113 + (ix * 0.00452313583942))))))); + s2 = 0.00948935714996 + (ix * (0.0901638932349 + (ix * (0.405368773071 + (ix * (0.912203410349 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak5bk5 + +// BEGIN: rational_ak6bk6 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak6bk6( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return 0.00439937562904; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = 0.00439937562904 + (x * (0.000487225670639 + (x * (-0.000128470657374 + (x * (0.00000529110969589 + (x * 1.5716677175e-7))))))); + s2 = 1.0 + (x * (0.794435257415 + (x * (0.333094721709 + (x * (0.0703527806143 + (x * 0.00806110846078))))))); + } else { + ix = 1.0 / x; + s1 = 1.5716677175e-7 + (ix * (0.00000529110969589 + (ix * (-0.000128470657374 + (ix * (0.000487225670639 + (ix * 0.00439937562904))))))); + s2 = 0.00806110846078 + (ix * (0.0703527806143 + (ix * (0.333094721709 + (ix * (0.794435257415 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak6bk6 + +// BEGIN: rational_ak7bk7 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak7bk7( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.0011481191232; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.0011481191232 + (x * (-0.112850923276 + (x * (1.51623048511 + (x * (-0.218472031183 + (x * 0.0730002451555))))))); + s2 = 1.0 + (x * (14.2482206905 + (x * (69.7360396285 + (x * (218.938950816 + (x * 277.067027185))))))); + } else { + ix = 1.0 / x; + s1 = 0.0730002451555 + (ix * (-0.218472031183 + (ix * (1.51623048511 + (ix * (-0.112850923276 + (ix * -0.0011481191232))))))); + s2 = 277.067027185 + (ix * (218.938950816 + (ix * (69.7360396285 + (ix * (14.2482206905 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak7bk7 + +// BEGIN: rational_ak8bk8 + +/** +* Evaluates a rational function (i.e., the ratio of two polynomials described by the coefficients stored in \\(P\\) and \\(Q\\)). +* +* ## Notes +* +* - Coefficients should be sorted in ascending degree. +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the rational function +* @return evaluated rational function +*/ +static double rational_ak8bk8( const double x ) { + double ax; + double ix; + double s1; + double s2; + if ( x == 0.0 ) { + return -0.000145727889667; + } + if ( x < 0.0 ) { + ax = -x; + } else { + ax = x; + } + if ( ax <= 1.0 ) { + s1 = -0.000145727889667 + (x * (-0.290806748131 + (x * (-13.308504545 + (x * (199.722374056 + (x * -11.4311378756))))))); + s2 = 1.0 + (x * (139.612587808 + (x * (2189.01116348 + (x * (7115.24019009 + (x * 45574.6081453))))))); + } else { + ix = 1.0 / x; + s1 = -11.4311378756 + (ix * (199.722374056 + (ix * (-13.308504545 + (ix * (-0.290806748131 + (ix * -0.000145727889667))))))); + s2 = 45574.6081453 + (ix * (7115.24019009 + (ix * (2189.01116348 + (ix * (139.612587808 + (ix * 1.0))))))); + } + return s1 / s2; +} + +// END: rational_ak8bk8 + +/* End auto-generated functions. */ + +/** +* Computes the sum of a Chebyshev polynomial. +* +* @param n degree of polynomial +* @param t input value +* @return Chebyshev sum +*/ +static double chepolsum( const int32_t n, const double t ) { + double tt; + double u0; + double u1; + double u2; + int32_t 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; +} + +/** +* Computes the Stirling series corresponding to asymptotic series for the logarithm of the gamma function. +* +* ```tex +* \frac{1}{12x}-\frac{1}{360x^3}\ldots; x \ge 3 +* ``` +* +* @param x input value +* @return function value +*/ +static double stirling( const double x ) { + double z; + + if ( x < SMALLEST_NORMAL_FLOAT32 ) { + return MAX_FLOAT32; + } + if ( x < 1.0 ) { + return stdlib_base_gammaln( x+1.0 ) - ( ( x+0.5 ) * stdlib_base_ln( x ) ) + x - LN_SQRT_TWO_PI; + } + if ( x < 2.0 ) { + return stdlib_base_gammaln( x ) - ( ( x-0.5 ) * stdlib_base_ln( x ) ) + x - LN_SQRT_TWO_PI; + } + if ( x < 3.0 ) { + return stdlib_base_gammaln( x-1.0 ) - ( ( x-0.5 ) * stdlib_base_ln( x ) ) + x - LN_SQRT_TWO_PI + stdlib_base_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; +} + +/** +* Computes the regulated gamma function. +* +* @param x input value +* @return function value +*/ +static double gamstar( const 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 ); + } + // Case: x <= 0.0 + return MAX_FLOAT32; +} + +/** +* Returns the positive number satisfying \\( \eta^2/2=\lambda-1-\ln(\lambda) \\) with \\( \operatorname{sign}(\lambda-1)=\operatorname{sign}(\eta) \\). +* +* @param eta eta value +* @return value satisfying equation +*/ +static double lambdaeta( const double eta ) { + double AK[ 6 ]; + double L2; + double L3; + double L4; + double L5; + double la; + double L; + double q; + double r; + double s; + + 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 = polyval_ak1( r ); + } else if ( eta < 1.0 ) { + r = eta; + la = polyval_ak2( 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 * ( AK[ 0 ] + ( r * ( AK[ 1 ] + ( r * ( AK[ 2 ] + ( r * ( AK[ 3 ] + ( r * ( AK[ 4 ] + ( r*AK[ 5 ] ) ) ) ) ) ) ) ) ) ) ); + } + 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; +} + +/** +* Evaluates the `eps1` function. +* +* @param eta eta value +* @return function value +*/ +static double eps1( const double eta ) { + double la; + + if ( stdlib_base_abs( eta ) < 1.0 ) { + return rational_ak0bk0( eta ); + } + la = lambdaeta( eta ); + return stdlib_base_ln( eta / ( la-1.0 ) ) / eta; +} + +/** +* Evaluates the `eps2` function. +* +* @param eta eta value +* @return function value +*/ +static double eps2( const double eta ) { + double lnmeta; + double x; + + if ( eta < -5.0 ) { + x = eta * eta; + lnmeta = stdlib_base_ln( -eta ); + return ( 12.0 - x - ( 6.0*( lnmeta*lnmeta ) ) ) / ( 12.0*x*eta ); + } + if ( eta < -2.0 ) { + return rational_ak1bk1( eta ); + } + if ( eta < 2.0 ) { + return rational_ak2bk2( eta ); + } + if ( eta < 1000.0 ) { + return rational_ak3bk3( eta ) / ( -12.0*eta ); + } + return -1.0 / ( 12.0*eta ); +} + +/** +* Evaluates the `eps3` function. +* +* @param eta eta value +* @return function value +*/ +static double eps3( const double eta ) { + double x; + double y; + + if ( eta < -8.0 ) { + x = eta * eta; + y = stdlib_base_ln( -eta ) / eta; + return ( -30.0 + ( eta*y*( (6.0*x*y*y) - 12.0 + x ) ) ) / ( 12.0*eta*x*x ); + } + if ( eta < -4.0 ) { + return rational_ak4bk4( eta ) / ( eta*eta ); + } + if ( eta < -2.0 ) { + return rational_ak5bk5( eta ); + } + if ( eta < 2.0 ) { + return rational_ak6bk6( eta ); + } + if ( eta < 10.0 ) { + x = 1.0 / eta; + return rational_ak7bk7( x ) / ( eta*eta ); + } + if ( eta < 100.0 ) { + x = 1.0 / eta; + return rational_ak8bk8( x ) / ( eta*eta ); + } + return -stdlib_base_ln( eta ) / ( 12.0*eta*eta*eta ); +} + +/** +* Implementation of the high order Newton-like method. +* +* @param x0 initial value +* @param a scale parameter +* @param m indicator +* @param p probability value +* @param q probability value +* @param lgama logarithm of scale parameter +* @param invfp one over `fp` +* @param pcase boolean indicating whether `p < 0.5` +* @return function value of the inverse +*/ +static double higherNewton( double x0, const double a, const int32_t 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; + int32_t 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 ) ) { + // Overflow problems in one or more steps of the computation; the initial approximation to the root is returned... + return xini; + } + r = stdlib_base_exp( dlnr ); + } else { + r = -invfp * x; + } + if ( pcase ) { + // Call: gammainc( x, s[, regularized = true ][, upper = false ] ) + px = stdlib_base_gammainc( x, a, true, false ); + ck0 = -r * ( px - p ); + } else { + // Call: gammainc( x, s[, regularized = true ][, upper = true ] ) + 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 ) ); + if ( x == 0.0 || stdlib_base_is_nan( x ) ) { + xr = 0.0; + } else { + xr = x; + } + return xr; +} + +/** +* 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`. +* +* ## Notes +* +* - The equation is inverted with `min(p,q)`. +* +* @param a scale value of incomplete gamma function +* @param p probability value +* @param q probability value +* @return solution of the equations `P(a,xr) = p` and `Q(a,xr) = q` where `a` is a positive parameter +*/ +static double compute( const double a, const double p, const double q ) { + double CK[ 5 ]; + double ap1inv; + double invfp; + double lgama; + 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 p6; + double p5; + double x0; + double a2; + double L2; + double L3; + double L4; + double b2; + double b3; + double p3; + double a4; + double fp; + double p4; + double p2; + double a3; + double xr; + double ck; + double b; + double L; + int32_t i; + int32_t k; + int32_t m; + double r; + double s; + double t; + double y; + + invfp = 0.0; + lgama = 0.0; + if ( p < HALF ) { + 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*HALF ) + ( p3*ONEO3 ) + ( p4*ONEO4 ) + ( p5*ONEO5 ) + ( p6*ONEO6 ); + } 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 < HALF ) { + 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( ONEO5 * ( 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 ] = HALF * ( ( 3.0*a ) + 5.0 ) / ( ap12*ap2 ); + CK[ 3 ] = ONEO3 * ( 31.0 + (8.0*a2) + (33.0*a) ) / ( ap13*ap2*ap3 ); + CK[ 4 ] = ONEO24 * ( 2888.0 + (1179.0*a3) + (125.0*a4) + (3971.0*a2) + (5661.0*a) ) / ( ap14*ap22*ap3*( a+4.0 ) ); + x0 = r * ( CK[ 0 ] + ( r * ( CK[ 1 ] + ( r * ( CK[ 2 ] + ( r * ( CK[ 3 ] + ( r*CK[ 4 ] ) ) ) ) ) ) ) ); + 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 ) * 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*( CK[ 0 ] + ( r * ( CK[ 1 ] + ( r * ( CK[ 2 ] + ( r * ( CK[ 3 ] + ( r*CK[ 4 ] ) ) ) ) ) ) ) ) ); + } else { + r = 1.0 / x0; + ck = L - 1.0; + t = L - ( b*r*ck ); + if ( t < x0 ) { + x0 -= t; + } + } + lgama = stdlib_base_gammaln( a ); + k = 1; + } + } + if ( ( stdlib_base_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_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*HALF ); + 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( -HALF*a*y*y ) / ( gamstar( a ) ); + invfp = 1.0 / fp; + } else { + // Overflow problems in one or more steps of the computation... + return 0.0 / 0.0; // NaN + } + } + 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`. +* +* ## Method +* +* The present code uses different methods of computation depending on the values of the input values: Taylor, asymptotic expansions and high-order Newton methods. +* +* ## Notes +* +* - The claimed accuracy obtained using this inversion routine is near `1e-12`. +* +* ## References +* +* - A. Gil, J. Segura and N.M. Temme, GammaCHI: a package for the inversion and computation of the gamma and chi-square distribution functions (central and noncentral). Computer Physics Commun +* - A. Gil, J. Segura and N.M. Temme. Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function ratios. SIAM J Sci Comput. (2012) 34(6), A2965-A2981 +* +* @param p probability value +* @param a scale parameter +* @param upper 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` +* @return function value of the inverse +* +* @example +* double y = stdlib_base_gammaincinv( 0.5, 2.0, false ); +* // returns ~1.678 +*/ +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.js b/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.js index 9a09c68a8156..48beadd93c19 100644 --- a/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.js +++ b/lib/node_modules/@stdlib/math/base/special/gammaincinv/test/test.js @@ -25,6 +25,7 @@ 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 abs = require( '@stdlib/math/base/special/abs' ); +var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' ); var PINF = require( '@stdlib/constants/float64/pinf' ); var gammaincinv = require( './../lib' ); @@ -116,54 +117,54 @@ tape( 'the function returns `NaN` if provided `p` outside the interval `[0,1]`', t.end(); }); -tape( 'the function returns `0` if provided `p = 1`', function test( t ) { +tape( 'the function returns `+Infinity` if provided `p = 1`', function test( t ) { var val; var a; var i; - a = incrspace( 0, 100, 10 ); + a = incrspace( 10, 100, 10 ); for ( i = 0; i < a.length; i++ ) { val = gammaincinv( 1, a[ i ] ); - t.notOk( val === 0, 'returns expected value' ); + t.strictEqual( val, PINF, 'returns expected value' ); } t.end(); }); -tape( 'the function returns `+Infinity` if provided `p = 0`', function test( t ) { +tape( 'the function returns `0` if provided `p = 0`', function test( t ) { var val; var a; var i; - a = incrspace( 0, 100, 10 ); + a = incrspace( 10, 100, 10 ); for ( i = 0; i < a.length; i++ ) { val = gammaincinv( 0, a[ i ] ); - t.notOk( val === PINF, 'returns expected value' ); + t.strictEqual( val, 0.0, 'returns expected value' ); } t.end(); }); -tape( 'the function returns `0` if provided `p = 0` when `upper = true`', function test( t ) { +tape( 'the function returns `+Infinity` if provided `p = 0` when `upper = true`', function test( t ) { var val; var a; var i; - a = incrspace( 0, 100, 10 ); + a = incrspace( 10, 100, 10 ); for ( i = 0; i < a.length; i++ ) { val = gammaincinv( 0, a[ i ], true ); - t.notOk( val === 0, 'returns expected value' ); + t.strictEqual( val, PINF, 'returns expected value' ); } t.end(); }); -tape( 'the function returns `+Infinity` if provided `p = 1` when `upper = true`', function test( t ) { +tape( 'the function returns `0` if provided `p = 1` when `upper = true`', function test( t ) { var val; var a; var i; - a = incrspace( 0, 100, 10 ); + a = incrspace( 10, 100, 10 ); for ( i = 0; i < a.length; i++ ) { val = gammaincinv( 1, a[ i ], true ); - t.notOk( val === PINF, 'returns expected value' ); + t.strictEqual( val, 0.0, 'returns expected value' ); } t.end(); }); @@ -195,7 +196,7 @@ tape( 'the function inverts the upper incomplete gamma function when `upper=true var b1; var b2; var i; - for ( i = 0; i < arg1.length; i++ ) { + for ( i = 0; i < upperArg1.length; i++ ) { actual = gammaincinv( upperArg1[ i ], upperArg2[ i ], true ); b1 = isfinite( actual ); @@ -203,10 +204,10 @@ tape( 'the function inverts the upper incomplete gamma function when `upper=true t.strictEqual( b1, b2, 'returned result is ' + ( (b2) ? 'finite' : 'not finite' ) ); b1 = isnan( actual ); - b2 = isnan( expected[ i ] ); + b2 = isnan( upperExpected[ i ] ); t.strictEqual( b1, b2, 'returned result is ' + ( (b1) ? '' : 'not' ) + ' NaN' ); if ( !b1 || !b2 ) { - t.ok( abs( actual - upperExpected[ i ] ) < 5e-14, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + upperExpected[ i ] + '.' ); + t.strictEqual( isAlmostSameValue( actual, upperExpected[ i ], 40.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + upperExpected[ i ] + '.' ); } } t.end(); 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..1549d00a26fd --- /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 `+Infinity` if provided `p = 1`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 10, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 1, a[ i ] ); + t.strictEqual( val, PINF, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `0` if provided `p = 0`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 10, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 0, a[ i ] ); + t.strictEqual( val, 0.0, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `+Infinity` if provided `p = 0` when `upper = true`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 10, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 0, a[ i ], true ); + t.strictEqual( val, PINF, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function returns `0` if provided `p = 1` when `upper = true`', opts, function test( t ) { + var val; + var a; + var i; + + a = incrspace( 10, 100, 10 ); + for ( i = 0; i < a.length; i++ ) { + val = gammaincinv( 1, a[ i ], true ); + t.strictEqual( val, 0.0, '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 ], 25.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 < upperArg1.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( upperExpected[ 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 ], 400.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 ], 50.0 ), true, 'returned result is within tolerance. actual: ' + actual + '; expected: ' + largeExpected[ i ] + '.' ); + } + } + t.end(); +});