diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/README.md b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/README.md index 0adc87bb4937..37381a417f68 100644 --- a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/README.md +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/README.md @@ -85,6 +85,109 @@ for ( i = 0; i < 100; i++ ) { + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/kernel_betaincinv.h" +``` + +#### stdlib_base_kernel_betaincinv( a, b, p, q, &out1, &out2 ) + +Computes the inverse of the lower incomplete beta function. + +```c +double out1; +double out2; + +stdlib_base_kernel_betaincinv( 3.0, 3.0, 0.2, 0.8, &out1, &out2 ); +``` + +The function accepts the following arguments: + +- **a**: `[in] double` first function parameter (a positive number). +- **b**: `[in] double` second function parameter (a positive number). +- **p**: `[in] double` probability. +- **q**: `[in] double` probability equal to `1-p`. +- **out1**: `[out] double*` destination pointer to store the function value `y`. +- **out2**: `[out] double*` destination pointer to store `1-y`. + +```c +void stdlib_base_kernel_betaincinv( const double a, const double b, const double p, const double q, double *out1, double *out2 ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/kernel_betaincinv.h" +#include "stdlib/random/base/randu.h" +#include +#include + +int main( void ) { + struct BasePRNGObject *obj = stdlib_base_random_randu_allocate( 0 ); + double out1; + double out2; + int32_t i; + double p; + double a; + double b; + + for ( i = 0; i < 100; i++ ) { + p = stdlib_base_random_randu( obj ); + a = stdlib_base_random_randu( obj ) * 10.0; + b = stdlib_base_random_randu( obj ) * 10.0; + stdlib_base_kernel_betaincinv( a, b, p, 1.0-p, &out1, &out2 ); + printf( "p: %lf, a: %lf, b: %lf, f(p,a,b): %lf\n", p, a, b, out1 ); + } + + stdlib_base_random_randu_free( obj ); +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/benchmark.native.js new file mode 100644 index 000000000000..31909ca7ac8f --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/benchmark.native.js @@ -0,0 +1,68 @@ +/** +* @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 EPS = require( '@stdlib/constants/float64/eps' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var kernelBetaincinv = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( kernelBetaincinv instanceof Error ) +}; + + +// MAIN // + +bench( format( '%s::native', pkg ), opts, function benchmark( assert ) { + var out; + var p; + var a; + var b; + var i; + + p = uniform( 100, 0.0, 1.0 ); + a = uniform( 100, EPS, 1000.0 ); + b = uniform( 100, EPS, 1000.0 ); + + out = [ 0.0, 0.0 ]; + assert.tic(); + for ( i = 0; i < assert.iterations; i++ ) { + out = kernelBetaincinv( a[ i%a.length ], b[ i%b.length ], p[ i%p.length ], 1.0-p[ i%p.length ] ); + if ( isnan( out[ 0 ] ) || isnan( out[ 1 ] ) ) { + assert.fail( 'should not return NaN' ); + } + } + assert.toc(); + if ( isnan( out[ 0 ] ) || isnan( out[ 1 ] ) ) { + assert.fail( 'should not return NaN' ); + } + assert.pass( 'benchmark finished' ); + assert.end(); +}); diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/c/native/Makefile new file mode 100644 index 000000000000..979768abbcec --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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/kernel-betaincinv/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/c/native/benchmark.c new file mode 100644 index 000000000000..078b243fed96 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/benchmark/c/native/benchmark.c @@ -0,0 +1,144 @@ +/** +* @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/kernel_betaincinv.h" +#include "stdlib/constants/float64/eps.h" +#include +#include +#include +#include +#include + +#define NAME "kernel-betaincinv" +#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 p[ 100 ]; + double a[ 100 ]; + double b[ 100 ]; + double elapsed; + double out1; + double out2; + 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, 100.0 ); + b[ i ] = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 100.0 ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + stdlib_base_kernel_betaincinv( a[ i % 100 ], b[ i % 100 ], p[ i % 100 ], 1.0 - p[ i % 100 ], &out1, &out2 ); + if ( out1 != out1 || out2 != out2 ) { + printf( "unexpected results\n" ); + break; + } + } + elapsed = tic() - t; + if ( out1 != out1 || out2 != out2 ) { + printf( "unexpected results\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/kernel-betaincinv/binding.gyp b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/binding.gyp new file mode 100644 index 000000000000..0d6508a12e99 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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/kernel-betaincinv/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/examples/c/Makefile new file mode 100644 index 000000000000..c8f8e9a1517b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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/kernel-betaincinv/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/examples/c/example.c new file mode 100644 index 000000000000..2245db338ad1 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/examples/c/example.c @@ -0,0 +1,42 @@ +/** +* @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/kernel_betaincinv.h" +#include "stdlib/random/base/randu.h" +#include +#include + +int main( void ) { + struct BasePRNGObject *obj = stdlib_base_random_randu_allocate( 0 ); + double out1; + double out2; + int32_t i; + double p; + double a; + double b; + + for ( i = 0; i < 100; i++ ) { + p = stdlib_base_random_randu( obj ); + a = stdlib_base_random_randu( obj ) * 10.0; + b = stdlib_base_random_randu( obj ) * 10.0; + stdlib_base_kernel_betaincinv( a, b, p, 1.0-p, &out1, &out2 ); + printf( "p: %lf, a: %lf, b: %lf, f(p,a,b): %lf\n", p, a, b, out1 ); + } + + stdlib_base_random_randu_free( obj ); +} diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/include.gypi b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/include.gypi new file mode 100644 index 000000000000..bee8d41a2caf --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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': [ + '[ ~0.327, ~0.673 ] +* +* @example +* var out = kernelBetaincinv( 3.0, 3.0, 0.4, 0.6 ); +* // returns [ ~0.446, ~0.554 ] +* +* @example +* var out = kernelBetaincinv( 1.0, 6.0, 0.4, 0.6 ); +* // returns [ ~0.082, ~0.918 ] +* +* @example +* var out = kernelBetaincinv( 1.0, 6.0, 0.8, 0.2 ); +* // returns [ ~0.235, ~0.765 ] +*/ +function kernelBetaincinv( a, b, p, q ) { + var out = new Float64Array( 2 ); + addon( a, b, p, q, out ); + return out; +} + + +// EXPORTS // + +module.exports = kernelBetaincinv; diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/manifest.json b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/manifest.json new file mode 100644 index 000000000000..43cbbc37fba2 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/manifest.json @@ -0,0 +1,160 @@ +{ + "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-float64array", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/acos", + "@stdlib/math/base/special/asin", + "@stdlib/math/base/special/beta", + "@stdlib/math/base/special/cos", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/expm1", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gamma-delta-ratio", + "@stdlib/math/base/special/gammaincinv", + "@stdlib/math/base/special/kernel-betainc", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/log1p", + "@stdlib/math/base/special/max", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/round", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sin", + "@stdlib/math/base/special/sqrt", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/half-pi", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/smallest-normal", + "@stdlib/constants/float64/smallest-subnormal", + "@stdlib/constants/float64/sqrt-two" + ] + }, + { + "task": "benchmark", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/acos", + "@stdlib/math/base/special/asin", + "@stdlib/math/base/special/beta", + "@stdlib/math/base/special/cos", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/expm1", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gamma-delta-ratio", + "@stdlib/math/base/special/gammaincinv", + "@stdlib/math/base/special/kernel-betainc", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/log1p", + "@stdlib/math/base/special/max", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/round", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sin", + "@stdlib/math/base/special/sqrt", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/half-pi", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/smallest-normal", + "@stdlib/constants/float64/smallest-subnormal", + "@stdlib/constants/float64/sqrt-two" + ] + }, + { + "task": "examples", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/random/base/randu", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/acos", + "@stdlib/math/base/special/asin", + "@stdlib/math/base/special/beta", + "@stdlib/math/base/special/cos", + "@stdlib/math/base/special/erfcinv", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/expm1", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gamma-delta-ratio", + "@stdlib/math/base/special/gammaincinv", + "@stdlib/math/base/special/kernel-betainc", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/log1p", + "@stdlib/math/base/special/max", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/round", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sin", + "@stdlib/math/base/special/sqrt", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/half-pi", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/smallest-normal", + "@stdlib/constants/float64/smallest-subnormal", + "@stdlib/constants/float64/sqrt-two" + ] + } + ] +} diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/package.json b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/package.json index 2405a0594d4b..0c0dde56a92a 100644 --- a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/package.json +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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/kernel-betaincinv/scripts/evalpoly.js b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/scripts/evalpoly.js index 108a622d67c6..4320adb596bf 100644 --- a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/scripts/evalpoly.js +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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 // @@ -208,6 +213,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 // /** @@ -218,7 +250,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; function main() { var fpath; var coefs; + var copts; var opts; + var file; var str; var i; @@ -255,6 +289,22 @@ function main() { str = header + compile( coefs[i][0] ); writeFileSync( fpath, str, opts ); } + + copts = { + 'dtype': 'double', + 'name': '' + }; + + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + for ( i = 0; i < coefs.length; i++ ) { + copts.name = 'polyval_'+coefs[i][1]; + str = compileC( coefs[i][0], copts ); + file = insert( file, copts.name, str ); + } + + writeFileSync( fpath, file, opts ); } main(); diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/Makefile b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/Makefile new file mode 100644 index 000000000000..2caf905cedbe --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/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/kernel-betaincinv/src/addon.c b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/addon.c new file mode 100644 index 000000000000..f551699e2f6f --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/addon.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/kernel_betaincinv.h" +#include "stdlib/napi/argv.h" +#include "stdlib/napi/argv_double.h" +#include "stdlib/napi/argv_float64array.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, 5 ); + STDLIB_NAPI_ARGV_DOUBLE( env, a, argv, 0 ); + STDLIB_NAPI_ARGV_DOUBLE( env, b, argv, 1 ); + STDLIB_NAPI_ARGV_DOUBLE( env, p, argv, 2 ); + STDLIB_NAPI_ARGV_DOUBLE( env, q, argv, 3 ); + STDLIB_NAPI_ARGV_FLOAT64ARRAY( env, out, len, argv, 4 ); + stdlib_base_kernel_betaincinv( a, b, p, q, &out[ 0 ], &out[ 1 ] ); + return NULL; +} + +STDLIB_NAPI_MODULE_EXPORT_FCN( addon ) diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/main.c b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/main.c new file mode 100644 index 000000000000..b23188cfc29a --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/src/main.c @@ -0,0 +1,1978 @@ +/** +* @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 original C++ code and copyright notice are from the Boost library ([ibeta_inverse.hpp]{@link http://www.boost.org/doc/libs/1_62_0/boost/math/special_functions/detail/ibeta_inverse.hpp}, [t_distribution_inv.hpp]{@link http://www.boost.org/doc/libs/1_62_0/boost/math/special_functions/detail/t_distribution_inv.hpp}, and [roots.hpp]{@link http://www.boost.org/doc/libs/1_62_0/boost/math/tools/roots.hpp}). The implementation has been modified according to stdlib conventions. +* +* ```text +* Copyright John Maddock 2006. +* Copyright Paul A. Bristow 2007. +* +* Use, modification and distribution are subject to the +* Boost Software License, Version 1.0. (See accompanying file +* LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) +* ``` +*/ + +#include "stdlib/math/base/special/kernel_betaincinv.h" +#include "stdlib/constants/float64/eps.h" +#include "stdlib/constants/float64/half_pi.h" +#include "stdlib/constants/float64/max.h" +#include "stdlib/constants/float64/pi.h" +#include "stdlib/constants/float64/smallest_normal.h" +#include "stdlib/constants/float64/smallest_subnormal.h" +#include "stdlib/constants/float64/sqrt_two.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/math/base/special/acos.h" +#include "stdlib/math/base/special/asin.h" +#include "stdlib/math/base/special/beta.h" +#include "stdlib/math/base/special/cos.h" +#include "stdlib/math/base/special/erfcinv.h" +#include "stdlib/math/base/special/exp.h" +#include "stdlib/math/base/special/expm1.h" +#include "stdlib/math/base/special/floor.h" +#include "stdlib/math/base/special/gamma_delta_ratio.h" +#include "stdlib/math/base/special/gammaincinv.h" +#include "stdlib/math/base/special/kernel_betainc.h" +#include "stdlib/math/base/special/ldexp.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/special/log1p.h" +#include "stdlib/math/base/special/max.h" +#include "stdlib/math/base/special/min.h" +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/round.h" +#include "stdlib/math/base/special/signum.h" +#include "stdlib/math/base/special/sin.h" +#include "stdlib/math/base/special/sqrt.h" +#include +#include + +static const double EPSILON = STDLIB_CONSTANT_FLOAT64_EPS; +static const double HALF_PI = STDLIB_CONSTANT_FLOAT64_HALF_PI; +static const double FLOAT64_MAX = STDLIB_CONSTANT_FLOAT64_MAX; +static const double PI = STDLIB_CONSTANT_FLOAT64_PI; +static const double FLOAT64_MIN_NORM = STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL; +static const double SMALLEST_SUBNORMAL = STDLIB_CONSTANT_FLOAT64_SMALLEST_SUBNORMAL; +static const double SQRT2 = STDLIB_CONSTANT_FLOAT64_SQRT2; + +// Maximum tolerable search bound (see `root_finder.js`): +static const double BIG = STDLIB_CONSTANT_FLOAT64_MAX / 4.0; + +// Degrees-of-freedom threshold (see `inverse_students_t.js`): +static const double DF_THRESHOLD = 0x10000000; // 2^28 + +static const double ONE_THIRD = 1.0 / 3.0; + +// Constant used when finding the inverse of the Student's t-distribution via Hill's approach (JS: `C` in `inverse_students_t.js`): +static const double CONST_C = 0.85498797333834849467655443627193; + +// Number of binary digits which determines the tolerance for Halley iteration (see `main.js`): +static const int32_t DIGITS = 32; + +// Maximum number of iterations when performing Halley iteration (see `main.js`): +static const int32_t MAX_ITERATIONS = 1000; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_co1 + +/** +* 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_co1( const double x ) { + return -1.0 + (x * (-5.0 + (x * 5.0))); +} + +// END: polyval_co1 + +// BEGIN: polyval_co2 + +/** +* 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_co2( const double x ) { + return 1.0 + (x * (21.0 + (x * (-69.0 + (x * 46.0))))); +} + +// END: polyval_co2 + +// BEGIN: polyval_co3 + +/** +* 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_co3( const double x ) { + return 7.0 + (x * (-2.0 + (x * (33.0 + (x * (-62.0 + (x * 31.0))))))); +} + +// END: polyval_co3 + +// BEGIN: polyval_co4 + +/** +* 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_co4( const double x ) { + return 25.0 + (x * (-52.0 + (x * (-17.0 + (x * (88.0 + (x * (-115.0 + (x * 46.0))))))))); +} + +// END: polyval_co4 + +// BEGIN: polyval_co5 + +/** +* 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_co5( const double x ) { + return 7.0 + (x * (12.0 + (x * (-78.0 + (x * 52.0))))); +} + +// END: polyval_co5 + +// BEGIN: polyval_co6 + +/** +* 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_co6( const double x ) { + return -7.0 + (x * (2.0 + (x * (183.0 + (x * (-370.0 + (x * 185.0))))))); +} + +// END: polyval_co6 + +// BEGIN: polyval_co7 + +/** +* 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_co7( const double x ) { + return -533.0 + (x * (776.0 + (x * (-1835.0 + (x * (10240.0 + (x * (-13525.0 + (x * 5410.0))))))))); +} + +// END: polyval_co7 + +// BEGIN: polyval_co8 + +/** +* 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_co8( const double x ) { + return -1579.0 + (x * (3747.0 + (x * (-3372.0 + (x * (-15821.0 + (x * (45588.0 + (x * (-45213.0 + (x * 15071.0))))))))))); +} + +// END: polyval_co8 + +// BEGIN: polyval_co9 + +/** +* 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_co9( const double x ) { + return 449.0 + (x * (-1259.0 + (x * (-769.0 + (x * (6686.0 + (x * (-9260.0 + (x * 3704.0))))))))); +} + +// END: polyval_co9 + +// BEGIN: polyval_co10 + +/** +* 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_co10( const double x ) { + return 63149.0 + (x * (-151557.0 + (x * (140052.0 + (x * (-727469.0 + (x * (2239932.0 + (x * (-2251437.0 + (x * 750479.0))))))))))); +} + +// END: polyval_co10 + +// BEGIN: polyval_co11 + +/** +* 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_co11( const double x ) { + return 29233.0 + (x * (-78755.0 + (x * (105222.0 + (x * (146879.0 + (x * (-1602610.0 + (x * (3195183.0 + (x * (-2554139.0 + (x * 729754.0))))))))))))); +} + +// END: polyval_co11 + +// BEGIN: polyval_co12 + +/** +* 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_co12( const double x ) { + return 1.0 + (x * (-13.0 + (x * 13.0))); +} + +// END: polyval_co12 + +// BEGIN: polyval_co13 + +/** +* 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_co13( const double x ) { + return 1.0 + (x * (21.0 + (x * (-69.0 + (x * 46.0))))); +} + +// END: polyval_co13 + +// BEGIN: polyval_co14 + +/** +* 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_co14( const double x ) { + return 0.16666666666666666 + (x * 0.16666666666666666); +} + +// END: polyval_co14 + +// BEGIN: polyval_co15 + +/** +* 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_co15( const double x ) { + return 0.058333333333333334 + (x * (0.06666666666666667 + (x * 0.008333333333333333))); +} + +// END: polyval_co15 + +// BEGIN: polyval_co16 + +/** +* 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_co16( const double x ) { + return 0.0251984126984127 + (x * (0.026785714285714284 + (x * (0.0017857142857142857 + (x * 0.0001984126984126984))))); +} + +// END: polyval_co16 + +// BEGIN: polyval_co17 + +/** +* 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_co17( const double x ) { + return 0.012039792768959435 + (x * (0.010559964726631394 + (x * (-0.0011078042328042327 + (x * (0.0003747795414462081 + (x * 0.0000027557319223985893))))))); +} + +// END: polyval_co17 + +// BEGIN: polyval_co18 + +/** +* 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_co18( const double x ) { + return 0.003837005972422639 + (x * (0.00610392115600449 + (x * (-0.0016095979637646305 + (x * (0.0005945867404200738 + (x * (-0.00006270542728876062 + (x * 2.505210838544172e-8))))))))); +} + +// END: polyval_co18 + +// BEGIN: polyval_co19 + +/** +* 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_co19( const double x ) { + return 0.0032177478835464946 + (x * (0.0010898206731540065 + (x * (-0.0012579159844784845 + (x * (0.0006908420797309686 + (x * (-0.00016376804137220805 + (x * (0.0000154012654012654 + (x * 1.6059043836821613e-10))))))))))); +} + +// END: polyval_co19 + +// BEGIN: polyval_co20 + +/** +* 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_co20( const double x ) { + return 0.001743826229834001 + (x * (0.00003353097688001788 + (x * (-0.0007624513544032393 + (x * (0.0006451304695145635 + (x * (-0.000249472580470431 + (x * (0.000049255746366361444 + (x * (-0.0000039851014346715405 + (x * 7.647163731819816e-13))))))))))))); +} + +// END: polyval_co20 + +// BEGIN: polyval_co21 + +/** +* 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_co21( const double x ) { + return 0.0009647274732138864 + (x * (-0.0003110108632631878 + (x * (-0.00036307660358786886 + (x * (0.0005140660578834113 + (x * (-0.00029133414466938067 + (x * (0.00009086710793521991 + (x * (-0.000015303004486655377 + (x * (0.0000010914179173496788 + (x * 2.8114572543455206e-15))))))))))))))); +} + +// END: polyval_co21 + +// BEGIN: polyval_co22 + +/** +* 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_co22( const double x ) { + return 0.0005422926281312969 + (x * (-0.0003694266780000966 + (x * (-0.00010230378073700413 + (x * (0.00035764655430568635 + (x * (-0.00028690924218514614 + (x * (0.00012645437628698076 + (x * (-0.000033202652391372056 + (x * (0.000004890304529197534 + (x * (-3.123956959982987e-7 + (x * 8.22063524662433e-18))))))))))))))))); +} + +// END: polyval_co22 + +/* End auto-generated functions. */ + +/** +* Evaluates a polynomial using double-precision floating-point arithmetic. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param c polynomial coefficients sorted in ascending degree +* @param len number of coefficients +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double evalpoly( const double *c, const int32_t len, const double x ) { + int32_t i; + double p; + + i = len; + if ( i < 2 || x == 0.0 ) { + return ( i == 0 ) ? 0.0 : c[ 0 ]; + } + i -= 1; + p = ( c[ i ] * x ) + c[ i-1 ]; + i -= 2; + while ( i >= 0 ) { + p = ( p * x ) + c[ i ]; + i -= 1; + } + return p; +} + +/** +* Helper function used by root finding code to convert `eta` to `x`. +* +* ## Notes +* +* - The function writes the function value and its first derivative to the two-element output array `out`. +* +* @param t first parameter +* @param a second parameter +* @param x function value +* @param out output array +*/ +static void temmeRoots( const double t, const double a, const double x, double *out ) { + double f1; + double f; + double y; + + y = 1.0 - x; + if ( y == 0.0 ) { + out[ 0 ] = -BIG; + out[ 1 ] = -BIG; + return; + } + if ( x == 0.0 ) { + out[ 0 ] = -BIG; + out[ 1 ] = -BIG; + return; + } + f = stdlib_base_ln( x ) + ( a * stdlib_base_ln( y ) ) + t; + f1 = ( 1.0 / x ) - ( a / y ); + out[ 0 ] = f; + out[ 1 ] = f1; + return; +} + +/** +* Calculates the roots of the (upper or lower) incomplete beta function together with its first and second derivatives. +* +* ## Notes +* +* - The function writes the function value and its first and second derivatives to the three-element output array `out`. +* +* @param a function parameter +* @param b function parameter +* @param target probability value +* @param invert boolean indicating whether to find the roots of the upper or lower incomplete beta function +* @param x input value +* @param out output array +*/ +static void ibetaRoots( const double a, const double b, const double target, const bool invert, const double x, double *out ) { + double buf[ 2 ]; + double f1; + double f2; + double xc; + double f; + double y; + + xc = x; + y = 1.0 - xc; + + buf[ 0 ] = 0.0; + buf[ 1 ] = 0.0; + stdlib_base_kernel_betainc( xc, a, b, true, invert, &buf[ 0 ], &buf[ 1 ] ); + f = buf[ 0 ] - target; + f1 = buf[ 1 ]; + if ( invert ) { + f1 = -f1; + } + if ( y == 0.0 ) { + y = FLOAT64_MIN_NORM * 64.0; + } + if ( xc == 0.0 ) { + xc = FLOAT64_MIN_NORM * 64.0; + } + f2 = f1 * ( -( y * a ) + ( ( b - 2.0 ) * xc ) + 1.0 ); + if ( stdlib_base_abs( f2 ) < y * xc * FLOAT64_MAX ) { + f2 /= ( y * xc ); + } + if ( invert ) { + f2 = -f2; + } + // Make sure we don't have a zero derivative: + if ( f1 == 0.0 ) { + f1 = ( ( invert ) ? -1.0 : 1.0 ) * FLOAT64_MIN_NORM * 64.0; + } + out[ 0 ] = f; + out[ 1 ] = f1; + out[ 2 ] = f2; + return; +} + +/** +* Performs root finding via second order Newton-Raphson iteration. +* +* ## Notes +* +* - In the original Boost/JavaScript implementations, the routine accepts a functor returning a two-element array of the function and its first derivative. As the functor is always `temmeRoots`, the functor parameters `t` and `a` are passed directly. +* +* @param t first `temmeRoots` parameter +* @param a second `temmeRoots` parameter +* @param guess initial starting value +* @param minimum minimum possible value for the result, used as initial lower bracket +* @param maximum maximum possible value for the result, used as initial upper bracket +* @param digits desired number of binary digits +* @param maxIter maximum number of iterations +* @return function value +*/ +static double newtonRaphsonIterate( const double t, const double a, const double guess, const double minimum, const double maximum, const int32_t digits, const int32_t maxIter ) { + double res[ 2 ]; + double f0last; + double delta1; + double delta2; + double factor; + double result; + double guessc; + double delta; + int32_t count; + double min; + double max; + double f0; + double f1; + + guessc = guess; + min = minimum; + max = maximum; + + f0 = 0.0; + f0last = 0.0; // cppcheck-suppress unreadVariable + result = guessc; + + factor = stdlib_base_ldexp( 1.0, 1 - digits ); + delta = FLOAT64_MAX; + delta1 = FLOAT64_MAX; + delta2 = FLOAT64_MAX; // cppcheck-suppress unreadVariable + + count = maxIter; + do { + f0last = f0; + delta2 = delta1; + delta1 = delta; + temmeRoots( t, a, result, res ); + f0 = res[ 0 ]; + f1 = res[ 1 ]; + count -= 1; + if ( f0 == 0.0 ) { + break; + } + if ( f1 == 0.0 ) { + // Oops zero derivative!!! + if ( f0last == 0.0 ) { + // Must be the first iteration, pretend that we had a previous one at either min or max: + if ( result == min ) { + guessc = max; + } else { + guessc = min; + } + // NOTE: the original JavaScript implementation assigns the two-element array returned by the functor to `f0last`, so the subsequent `signum` call always evaluates to `1.0`. We faithfully replicate that behavior: + f0last = 1.0; + delta = guessc - result; + } + if ( stdlib_base_signum( f0last ) * stdlib_base_signum( f0 ) < 0.0 ) { + // We've crossed over so move in opposite direction to last step: + if ( delta < 0.0 ) { + delta = ( result - min ) / 2.0; + } else { + delta = ( result - max ) / 2.0; + } + } else if ( delta < 0.0 ) { + delta = ( result - max ) / 2.0; + } else { + delta = ( result - min ) / 2.0; + } + } else { + delta = f0 / f1; + } + if ( stdlib_base_abs( delta * 2.0 ) > stdlib_base_abs( delta2 ) ) { + // Last two steps haven't converged, try bisection: + delta = ( delta > 0.0 ) ? ( result - min ) / 2.0 : ( result - max ) / 2.0; + } + guessc = result; + result -= delta; + if ( result <= min ) { + delta = 0.5 * ( guessc - min ); + result = guessc - delta; + if ( result == min || result == max ) { + break; + } + } else if ( result >= max ) { + delta = 0.5 * ( guessc - max ); + result = guessc - delta; + if ( result == min || result == max ) { + break; + } + } + // Update brackets: + if ( delta > 0.0 ) { + max = guessc; + } else { + min = guessc; + } + } + while ( count && ( stdlib_base_abs( result * factor ) < stdlib_base_abs( delta ) ) ); + + return result; +} + +/** +* Performs root finding via third order Halley iteration. +* +* ## Notes +* +* - In the original Boost/JavaScript implementations, the routine accepts a functor returning a three-element array of the function and its first two derivatives. As the functor is always `ibetaRoots`, the functor parameters `a`, `b`, `target`, and `invert` are passed directly. +* +* @param a function parameter of the incomplete beta function +* @param b function parameter of the incomplete beta function +* @param target probability value +* @param invert boolean indicating whether to find the roots of the upper or lower incomplete beta function +* @param guess initial starting value +* @param minimum minimum possible value for the result, used as initial lower bracket +* @param maximum maximum possible value for the result, used as initial upper bracket +* @param digits desired number of binary digits +* @param maxIter maximum number of iterations +* @return function value +*/ +static double halleyIterate( const double a, const double b, const double target, const bool invert, const double guess, const double minimum, const double maximum, const int32_t digits, const int32_t maxIter ) { + double convergence; + bool outOfBounds; + double res[ 3 ]; + double delta1; + double delta2; + double factor; + double result; + double f0Last; + double guessc; + int32_t count; + double delta; + double denom; + double diff; + double min; + double max; + double num; + double f0; + double f1; + double f2; + + guessc = guess; + min = minimum; + max = maximum; + + f0 = 0.0; + outOfBounds = false; + result = guessc; + factor = stdlib_base_ldexp( 1.0, 1 - digits ); + delta = stdlib_base_max( 10000000.0 * guessc, 10000000.0 ); // Arbitrarily large delta... + f0Last = 0.0; // cppcheck-suppress unreadVariable + delta1 = delta; + delta2 = delta; // cppcheck-suppress unreadVariable + + count = maxIter; + do { + f0Last = f0; + delta2 = delta1; + delta1 = delta; + ibetaRoots( a, b, target, invert, result, res ); + f0 = res[ 0 ]; + f1 = res[ 1 ]; + f2 = res[ 2 ]; + count -= 1; + + if ( f0 == 0.0 ) { + break; + } + if ( f1 == 0.0 ) { + // Oops zero derivative!!! + if ( f0Last == 0.0 ) { + // Must be the first iteration, pretend that we had a previous one at either min or max: + if ( result == min ) { + guessc = max; + } else { + guessc = min; + } + // NOTE: the original JavaScript implementation assigns the three-element array returned by the functor to `f0Last`, so the subsequent `signum` call always evaluates to `1.0`. We faithfully replicate that behavior: + f0Last = 1.0; + delta = guessc - result; + } + if ( stdlib_base_signum( f0Last ) * stdlib_base_signum( f0 ) < 0.0 ) { + // We've crossed over so move in opposite direction to last step: + if ( delta < 0.0 ) { + delta = ( result - min ) / 2.0; + } else { + delta = ( result - max ) / 2.0; + } + // Move in same direction as last step: + } else if ( delta < 0.0 ) { + delta = ( result - max ) / 2.0; + } else { + delta = ( result - min ) / 2.0; + } + } else if ( f2 == 0.0 ) { + delta = f0 / f1; + } else { + denom = 2.0 * f0; + num = ( 2.0 * f1 ) - ( f0 * ( f2 / f1 ) ); + if ( stdlib_base_abs( num ) < 1.0 && ( stdlib_base_abs( denom ) >= stdlib_base_abs( num ) * FLOAT64_MAX ) ) { + // Possible overflow, use Newton step: + delta = f0 / f1; + } else { + delta = denom / num; + } + if ( delta * f1 / f0 < 0.0 ) { + // Probably cancellation error, try a Newton step instead: + delta = f0 / f1; + if ( stdlib_base_abs( delta ) > 2.0 * stdlib_base_abs( guessc ) ) { + delta = ( ( delta < 0.0 ) ? -1.0 : 1.0 ) * 2.0 * stdlib_base_abs( guessc ); + } + } + } + convergence = stdlib_base_abs( delta / delta2 ); + if ( convergence > 0.8 && convergence < 2.0 ) { + // Last two steps haven't converged, try bisection: + delta = ( delta > 0.0 ) ? ( result - min ) / 2.0 : ( result - max ) / 2.0; + if ( stdlib_base_abs( delta ) > result ) { + delta = stdlib_base_signum( delta ) * result; // Protect against huge jumps! + } + // Reset delta2 so that this branch will *not* be taken on the next iteration: + delta2 = delta * 3.0; // cppcheck-suppress unreadVariable + } + guessc = result; + result -= delta; + + // Check for out of bounds step: + if ( result < min ) { + if ( + stdlib_base_abs( min ) < 1.0 && + stdlib_base_abs( result ) > 1.0 && + ( FLOAT64_MAX / stdlib_base_abs( result ) < stdlib_base_abs( min ) ) + ) { + diff = 1000.0; + } else { + diff = result / min; + } + if ( stdlib_base_abs( diff ) < 1.0 ) { + diff = 1.0 / diff; + } + if ( !outOfBounds && diff > 0.0 && diff < 3.0 ) { + // Only a small out of bounds step, let's assume that the result is probably approximately at minimum: + delta = 0.99 * ( guessc - min ); + result = guessc - delta; + outOfBounds = true; // Only take this branch once! + } else { + delta = ( guessc - min ) / 2.0; + result = guessc - delta; + if ( result == min || result == max ) { + break; + } + } + } else if ( result > max ) { + if ( + stdlib_base_abs( max ) < 1.0 && + stdlib_base_abs( result ) > 1.0 && + ( FLOAT64_MAX / stdlib_base_abs( result ) < stdlib_base_abs( max ) ) + ) { + diff = 1000.0; + } else { + diff = result / max; + } + if ( stdlib_base_abs( diff ) < 1.0 ) { + diff = 1.0 / diff; + } + if ( !outOfBounds && diff > 0.0 && diff < 3.0 ) { + // Only a small out of bounds step, let's assume that the result is probably approximately at minimum: + delta = 0.99 * ( guessc - max ); + result = guessc - delta; + outOfBounds = true; // Only take this branch once! + } else { + delta = ( guessc - max ) / 2.0; + result = guessc - delta; + if ( result == min || result == max ) { + break; + } + } + } + // Update brackets: + if ( delta > 0.0 ) { + max = guessc; + } else { + min = guessc; + } + } while ( count && ( stdlib_base_abs( result * factor ) < stdlib_base_abs( delta ) ) ); + + return result; +} + +/** +* Carries out the first method by Temme (described in section 2). +* +* ## References +* +* - Temme, N. M. 1992. "Incomplete Laplace Integrals: Uniform Asymptotic Expansion with Application to the Incomplete Beta Function." _Journal of Computational and Applied Mathematics_ 41 (1–2): 1638–63. doi:[10.1016/0377-0427(92)90244-R](https://doi.org/10.1016/0377-0427(92)90244-R). +* +* @param a function parameter +* @param b function parameter +* @param z function parameter +* @return function value +*/ +static double temme1( const double a, const double b, const double z ) { + double workspace[ 7 ]; + double terms[ 4 ]; + double eta0; + double eta2; + double eta; + double B2; + double B3; + double B; + double c; + + // Get the first approximation for eta from the inverse error function (Eq: 2.9 and 2.10): + eta0 = stdlib_base_erfcinv( 2.0 * z ); + eta0 /= -stdlib_base_sqrt( a / 2.0 ); + + terms[ 0 ] = eta0; + + // Calculate powers: + B = b - a; + B2 = B * B; + B3 = B2 * B; + + // Calculate correction terms: + + // See eq following 2.15: + workspace[ 0 ] = -B * SQRT2 / 2.0; + workspace[ 1 ] = ( 1.0 - ( 2.0 * B ) ) / 8.0; + workspace[ 2 ] = -( B * SQRT2 / 48.0 ); + workspace[ 3 ] = -1.0 / 192.0; + workspace[ 4 ] = -B * SQRT2 / 3840.0; + workspace[ 5 ] = 0.0; + workspace[ 6 ] = 0.0; + terms[ 1 ] = evalpoly( workspace, 7, eta0 ); + + // Eq Following 2.17: + workspace[ 0 ] = B * SQRT2 * ( ( 3.0 * B ) - 2.0 ) / 12.0; + workspace[ 1 ] = ( ( 20.0 * B2 ) - ( 12.0 * B ) + 1.0 ) / 128.0; + workspace[ 2 ] = B * SQRT2 * ( ( 20.0 * B ) - 1.0 ) / 960.0; + workspace[ 3 ] = ( ( 16.0 * B2 ) + ( 30.0 * B ) - 15.0 ) / 4608.0; + workspace[ 4 ] = B * SQRT2 * ( ( 21.0 * B ) + 32.0 ) / 53760.0; + workspace[ 5 ] = ( -( 32.0 * B2 ) + 63.0 ) / 368640.0; + workspace[ 6 ] = -B * SQRT2 * ( ( 120.0 * B ) + 17.0 ) / 25804480.0; + terms[ 2 ] = evalpoly( workspace, 7, eta0 ); + + // Eq Following 2.17: + workspace[ 0 ] = B * SQRT2 * ( ( -75.0 * B2 ) + ( 80.0 * B ) - 16.0 ) / 480.0; + workspace[ 1 ] = ( ( -1080.0 * B3 ) + ( 868.0 * B2 ) - ( 90.0 * B ) - 45.0 ) / 9216.0; + workspace[ 2 ] = B * SQRT2 * ( ( -1190.0 * B2 ) + ( 84.0 * B ) + 373.0 ) / 53760.0; + workspace[ 3 ] = ( ( -2240.0 * B3 ) - ( 2508.0 * B2 ) + ( 2100.0 * B ) - 165.0 ) / 368640.0; + workspace[ 4 ] = 0.0; + workspace[ 5 ] = 0.0; + workspace[ 6 ] = 0.0; + terms[ 3 ] = evalpoly( workspace, 7, eta0 ); + + // Bring them together to get a final estimate for eta: + eta = evalpoly( terms, 4, 1.0 / a ); + + // Now we need to convert eta to the return value `x`, by solving the appropriate quadratic equation: + eta2 = eta * eta; + c = -stdlib_base_exp( -eta2 / 2.0 ); + if ( eta2 == 0.0 ) { + return 0.5; + } + return ( 1.0 + ( eta * stdlib_base_sqrt( ( 1.0 + c ) / eta2 ) ) ) / 2.0; +} + +/** +* Carries out the second method by Temme (described in section 3). +* +* ## References +* +* - Temme, N. M. 1992. "Incomplete Laplace Integrals: Uniform Asymptotic Expansion with Application to the Incomplete Beta Function." _Journal of Computational and Applied Mathematics_ 41 (1–2): 1638–63. doi:[10.1016/0377-0427(92)90244-R](https://doi.org/10.1016/0377-0427(92)90244-R). +* +* @param z function parameter +* @param r function parameter +* @param theta function parameter +* @return function value +*/ +static double temme2( const double z, const double r, const double theta ) { + double workspace[ 6 ]; + double terms[ 4 ]; + double upper; + double lower; + double alpha; + double eta0; + double eta; + double sc7; + double sc6; + double sc5; + double sc4; + double sc3; + double sc2; + double sc; + double lu; + double s2; + double c2; + double c; + double s; + double u; + double x; + + // Get first estimate for eta, see Eq 3.9 and 3.10, but note there is a typo in Eq 3.10: + eta0 = stdlib_base_erfcinv( 2.0 * z ) / ( -stdlib_base_sqrt( r / 2.0 ) ); + + s = stdlib_base_sin( theta ); + c = stdlib_base_cos( theta ); + + // Now we need to perturb eta0 to get eta, which we do by evaluating the polynomial in 1/r at the bottom of page 151, to do this we first need the error terms e1, e2 e3 which we'll fill into the array "terms". Since these terms are themselves polynomials, we'll need another array "workspace" to calculate those... + terms[ 0 ] = eta0; + + // Some powers of sin(theta) cos(theta) that we'll need later: + s2 = s * s; + c2 = c * c; + sc = s * c; + sc2 = sc * sc; + sc3 = sc2 * sc; + sc4 = sc2 * sc2; + sc5 = sc2 * sc3; + sc6 = sc3 * sc3; + sc7 = sc4 * sc3; + + // Calculate e1 and put it in terms[1], see the middle of page 151: + workspace[ 0 ] = ( ( 2.0 * s2 ) - 1.0 ) / ( 3.0 * sc ); + workspace[ 1 ] = -polyval_co1( s2 ) / ( 36.0 * sc2 ); + workspace[ 2 ] = polyval_co2( s2 ) / ( 1620.0 * sc3 ); + workspace[ 3 ] = polyval_co3( s2 ) / ( 6480.0 * sc4 ); + workspace[ 4 ] = polyval_co4( s2 ) / ( 90720.0 * sc5 ); + workspace[ 5 ] = 0.0; + terms[ 1 ] = evalpoly( workspace, 6, eta0 ); + + // Now evaluate e2 and put it in terms[2]: + workspace[ 0 ] = -polyval_co5( s2 ) / ( 405.0 * sc3 ); + workspace[ 1 ] = polyval_co6( s2 ) / ( 2592.0 * sc4 ); + workspace[ 2 ] = -polyval_co7( s2 ) / ( 204120.0 * sc5 ); + workspace[ 3 ] = -polyval_co8( s2 ) / ( 2099520.0 * sc6 ); + workspace[ 4 ] = 0.0; + workspace[ 5 ] = 0.0; + terms[ 2 ] = evalpoly( workspace, 6, eta0 ); + + // And e3, and put it in terms[3]: + workspace[ 0 ] = polyval_co9( s2 ) / ( 102060.0 * sc5 ); + workspace[ 1 ] = -polyval_co10( s2 ) / ( 20995200.0 * sc6 ); + workspace[ 2 ] = polyval_co11( s2 ) / ( 36741600.0 * sc7 ); + workspace[ 3 ] = 0.0; + workspace[ 4 ] = 0.0; + workspace[ 5 ] = 0.0; + terms[ 3 ] = evalpoly( workspace, 6, eta0 ); + + // Bring the correction terms together to evaluate eta; this is the last equation on page 151: + eta = evalpoly( terms, 4, 1.0 / r ); + + // Now that we have eta we need to back solve for x, we seek the value of x that gives eta in Eq 3.2. The two methods used are described in section 5. Begin by defining a few variables we'll need later: + alpha = c / s; + alpha *= alpha; + lu = ( -( eta * eta ) / ( 2.0 * s2 ) ) + stdlib_base_ln( s2 ) + ( c2 * stdlib_base_ln( c2 ) / s2 ); + + // Temme doesn't specify what value to switch on here, but this seems to work pretty well: + if ( stdlib_base_abs( eta ) < 0.7 ) { + // Small eta use the expansion Temme gives in the second equation of section 5, it's a polynomial in eta: + workspace[ 0 ] = s2; + workspace[ 1 ] = sc; + workspace[ 2 ] = ( 1.0 - ( 2.0 * s2 ) ) / 3.0; + workspace[ 3 ] = polyval_co12( s2 ) / ( 36.0 * sc ); + workspace[ 4 ] = polyval_co13( s2 ) / ( 270.0 * sc2 ); + workspace[ 5 ] = 0.0; + x = evalpoly( workspace, 6, eta ); + } else { + // If eta is large we need to solve Eq 3.2 more directly, begin by getting an initial approximation for x from the last equation on page 155, this is a polynomial in u: + u = stdlib_base_exp( lu ); + workspace[ 0 ] = u; + workspace[ 1 ] = alpha; + workspace[ 2 ] = 0.0; + workspace[ 3 ] = 3.0 * alpha * ( ( 3.0 * alpha ) + 1.0 ) / 6.0; + workspace[ 4 ] = 4.0 * alpha * ( ( 4.0 * alpha ) + 1.0 ) * ( ( 4.0 * alpha ) + 2.0 ) / 24.0; + workspace[ 5 ] = 5.0 * alpha * ( ( 5.0 * alpha ) + 1.0 ) * ( ( 5.0 * alpha ) + 2.0 ) * ( ( 5.0 * alpha ) + 3.0 ) / 120.0; + x = evalpoly( workspace, 6, u ); + + // At this point we may or may not have the right answer, Eq-3.2 has two solutions for x for any given eta, however the mapping in 3.2 is 1:1 with the sign of eta and x-sin^2(theta) being the same. So we can check if we have the right root of 3.2, and if not switch x for 1-x. This transformation is motivated by the fact that the distribution is *almost* symmetric so 1-x will be in the right ball park for the solution: + if ( ( x - s2 ) * eta < 0.0 ) { + x = 1.0 - x; + } + } + // The final step is a few Newton-Raphson iterations to clean up our approximation for x, this is pretty cheap in general, and very cheap compared to an incomplete beta evaluation. The limits set on x come from the observation that the sign of eta and x-sin^2(theta) are the same. + if ( eta < 0.0 ) { + lower = 0.0; + upper = s2; + } else { + lower = s2; + upper = 1.0; + } + // If our initial approximation is out of bounds then bisect: + if ( x < lower || x > upper ) { + x = ( lower + upper ) / 2.0; + } + // And iterate: + x = newtonRaphsonIterate( -lu, alpha, x, lower, upper, 32, 100 ); + return x; +} + +/** +* Carries out the third method by Temme (described in section 4). +* +* ## References +* +* - Temme, N. M. 1992. "Incomplete Laplace Integrals: Uniform Asymptotic Expansion with Application to the Incomplete Beta Function." _Journal of Computational and Applied Mathematics_ 41 (1–2): 1638–63. doi:[10.1016/0377-0427(92)90244-R](https://doi.org/10.1016/0377-0427(92)90244-R). +* +* @param a function parameter +* @param b function parameter +* @param p function parameter +* @param q probability equal to `1-p` +* @return function value +*/ +static double temme3( const double a, const double b, const double p, const double q ) { + double cross; + double lower; + double upper; + double eta0; + double eta; + double w10; + double w12; + double w13; + double w14; + double e1; + double e2; + double e3; + double mu; + double d2; + double d3; + double d4; + double w2; + double w3; + double w4; + double w5; + double w6; + double w7; + double w8; + double w9; + double w1; + double d; + double w; + double u; + double x; + + // Begin by getting an initial approximation for the quantity eta from the dominant part of the incomplete beta: + if ( p < q ) { + eta0 = stdlib_base_gammaincinv( p, b, true ); + } else { + eta0 = stdlib_base_gammaincinv( q, b, false ); + } + eta0 /= a; + + // Define the variables and powers we'll need later on: + mu = b / a; + w = stdlib_base_sqrt( 1.0 + mu ); + w2 = w * w; + w3 = w2 * w; + w4 = w2 * w2; + w5 = w3 * w2; + w6 = w3 * w3; + w7 = w4 * w3; + w8 = w4 * w4; + w9 = w5 * w4; + w10 = w5 * w5; + d = eta0 - mu; + d2 = d * d; + d3 = d2 * d; + d4 = d2 * d2; + w1 = w + 1.0; + w12 = w1 * w1; + w13 = w1 * w12; + w14 = w12 * w12; + + // Now we need to compute the perturbation error terms that convert eta0 to eta, these are all polynomials of polynomials. Probably these should be re-written to use tabulated data (see examples above), but it's less of a win in this case as we need to calculate the individual powers for the denominator terms anyway, so we might as well use them for the numerator-polynomials as well. Refer to p154-p155 for the details of these expansions: + e1 = ( w + 2.0 ) * ( w - 1.0 ) / ( 3.0 * w ); + e1 += ( w3 + ( 9.0 * w2 ) + ( 21.0 * w ) + 5.0 ) * d / ( 36.0 * w2 * w1 ); + e1 -= ( w4 - ( 13.0 * w3 ) + ( 69.0 * w2 ) + ( 167.0 * w ) + 46.0 ) * d2 / ( 1620.0 * w12 * w3 ); + e1 -= ( ( 7.0 * w5 ) + ( 21.0 * w4 ) + ( 70.0 * w3 ) + ( 26.0 * w2 ) - ( 93.0 * w ) - 31.0 ) * d3 / ( 6480.0 * w13 * w4 ); + e1 -= ( ( 75.0 * w6 ) + ( 202.0 * w5 ) + ( 188.0 * w4 ) - ( 888.0 * w3 ) - ( 1345.0 * w2 ) + ( 118.0 * w ) + 138.0 ) * d4 / ( 272160.0 * w14 * w5 ); + + e2 = ( ( 28.0 * w4 ) + ( 131.0 * w3 ) + ( 402.0 * w2 ) + ( 581.0 * w ) + 208.0 ) * ( w - 1.0 ) / ( 1620.0 * w1 * w3 ); + e2 -= ( ( 35.0 * w6 ) - ( 154.0 * w5 ) - ( 623.0 * w4 ) - ( 1636.0 * w3 ) - ( 3983.0 * w2 ) - ( 3514.0 * w ) - 925.0 ) * d / ( 12960.0 * w12 * w4 ); + e2 -= ( ( 2132.0 * w7 ) + ( 7915.0 * w6 ) + ( 16821.0 * w5 ) + ( 35066.0 * w4 ) + ( 87490.0 * w3 ) + ( 141183.0 * w2 ) + ( 95993.0 * w ) + 21640.0 ) * d2 / ( 816480.0 * w5 * w13 ); + e2 -= ( ( 11053.0 * w8 ) + ( 53308.0 * w7 ) + ( 117010.0 * w6 ) + ( 163924.0 * w5 ) + ( 116188.0 * w4 ) - ( 258428.0 * w3 ) - ( 677042.0 * w2 ) - ( 481940.0 * w ) - 105497.0 ) * d3 / ( 14696640.0 * w14 * w6 ); + + e3 = -( ( ( 3592.0 * w7 ) + ( 8375.0 * w6 ) - ( 1323.0 * w5 ) - ( 29198.0 * w4 ) - ( 89578.0 * w3 ) - ( 154413.0 * w2 ) - ( 116063.0 * w ) - 29632.0 ) * ( w - 1.0 ) ) / ( 816480.0 * w5 * w12 ); + e3 -= ( ( 442043.0 * w9 ) + ( 2054169.0 * w8 ) + ( 3803094.0 * w7 ) + ( 3470754.0 * w6 ) + ( 2141568.0 * w5 ) - ( 2393568.0 * w4 ) - ( 19904934.0 * w3 ) - ( 34714674.0 * w2 ) - ( 23128299.0 * w ) - 5253353.0 ) * d / ( 146966400.0 * w6 * w13 ); + e3 -= ( ( 116932.0 * w10 ) + ( 819281.0 * w9 ) + ( 2378172.0 * w8 ) + ( 4341330.0 * w7 ) + ( 6806004.0 * w6 ) + ( 10622748.0 * w5 ) + ( 18739500.0 * w4 ) + ( 30651894.0 * w3 ) + ( 30869976.0 * w2 ) + ( 15431867.0 * w ) + 2919016.0 ) * d2 / ( 146966400.0 * w14 * w7 ); + + // Combine eta0 and the error terms to compute eta (Second equation p155): + eta = eta0 + ( e1 / a ) + ( e2 / ( a * a ) ) + ( e3 / ( a * a * a ) ); + + /* + * Now we need to solve Eq 4.2 to obtain x. For any given value of eta there are two solutions to this equation, and since the distribution may be very skewed, these are not related by x ~ 1-x we used when implementing section 3 above. However we know that: + * + * cross < x <= 1 ; iff eta < mu + * x == cross ; iff eta == mu + * 0 <= x < cross ; iff eta > mu + * + * Where cross == 1 / (1 + mu). Many thanks to Prof Temme for clarifying this point. Therefore we'll just jump straight into Newton iterations to solve Eq 4.2 using these bounds, and simple bisection as the first guess, in practice this converges pretty quickly and we only need a few digits correct anyway: + */ + if ( eta <= 0.0 ) { + eta = SMALLEST_SUBNORMAL; + } + u = eta - ( mu * stdlib_base_ln( eta ) ) + ( ( 1.0 + mu ) * stdlib_base_ln( 1.0 + mu ) ) - mu; + cross = 1.0 / ( 1.0 + mu ); + lower = ( eta < mu ) ? cross : 0.0; + upper = ( eta < mu ) ? 1.0 : cross; + x = ( lower + upper ) / 2.0; + return newtonRaphsonIterate( u, mu, x, lower, upper, 32, 100 ); +} + +/** +* Evaluates Student's t quantiles via a method due to Hill. +* +* ## References +* +* - Hill, G. W. 1970. "Algorithm 396: Student's T-Quantiles." _Communications of the ACM_ 13 (10). New York, NY, USA: ACM: 619–20. doi:[10.1145/355598.355600](https://doi.org/10.1145/355598.355600). +* +* @param ndf degrees of freedom +* @param u input probability +* @return function value +*/ +static double inverseStudentsTHill( const double ndf, const double u ) { + double a; + double b; + double c; + double d; + double q; + double x; + double y; + + if ( ndf > 1e20 ) { + return -stdlib_base_erfcinv( 2.0 * u ) * SQRT2; + } + a = 1.0 / ( ndf - 0.5 ); + b = 48.0 / ( a * a ); + c = ( ( ( ( ( 20700.0 * a / b ) - 98.0 ) * a ) - 16.0 ) * a ) + 96.36; + d = ( ( ( ( 94.5 / ( b + c ) ) - 3.0 ) / b ) + 1.0 ) * stdlib_base_sqrt( a * HALF_PI ) * ndf; + y = stdlib_base_pow( d * 2.0 * u, 2.0 / ndf ); + + if ( y > ( 0.05 + a ) ) { + // Asymptotic inverse expansion about normal: + x = -stdlib_base_erfcinv( 2.0 * u ) * SQRT2; + y = x * x; + + if ( ndf < 5.0 ) { + c += 0.3 * ( ndf - 4.5 ) * ( x + 0.6 ); + } + c += ( ( ( ( ( ( ( 0.05 * d * x ) - 5.0 ) * x ) - 7.0 ) * x ) - 2.0 ) * x ) + b; + y = ( ( ( ( ( ( ( ( ( ( 0.4 * y ) + 6.3 ) * y ) + 36.0 ) * y ) + 94.5 ) / c ) - y - 3.0 ) / b ) + 1.0 ) * x; + y = stdlib_base_expm1( a * y * y ); + } else { + y = ( ( 1.0 / ( ( ( ( ndf + 6.0 ) / ( ndf * y ) ) - ( 0.089 * d ) - 0.822 ) * ( ndf + 2.0 ) * 3.0 ) + ( 0.5 / ( ndf + 4.0 ) ) ) * y - 1.0 ) * ( ndf + 1.0 ) / ( ndf + 2.0 ) + ( 1.0 / y ); + } + q = stdlib_base_sqrt( ndf * y ); + return -q; +} + +/** +* Evaluates Student's t quantiles via a body series expansion. Tail and body series are due to Shaw. +* +* ## References +* +* - Shaw, William T. 2006. "Sampling Student's T distribution – use of the inverse cumulative distribution function." _Journal of Computational Finance_ 9 (4): 37–73. [www.mth.kcl.ac.uk/~shaww/web\_page/papers/Tdistribution06.pdf](www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf). +* +* @param df degrees of freedom +* @param u input probability +* @return function value +*/ +static double inverseStudentsTBodySeries( const double df, const double u ) { + double c[ 10 ]; + double idf; + double c0; + double v; + + // Body series for small N, start with Eq 56 of Shaw: + v = stdlib_base_gamma_delta_ratio( df / 2.0, 0.5 ) * stdlib_base_sqrt( df * PI ) * ( u - 0.5 ); + + c0 = 0.0; + + // Figure out what the coefficients are. They depend only on the degrees of freedom (Eq 57 of Shaw): + idf = 1.0 / df; + c[ 0 ] = 1.0; + c[ 1 ] = polyval_co14( idf ); + c[ 2 ] = polyval_co15( idf ); + c[ 3 ] = polyval_co16( idf ); + c[ 4 ] = polyval_co17( idf ); + c[ 5 ] = polyval_co18( idf ); + c[ 6 ] = polyval_co19( idf ); + c[ 7 ] = polyval_co20( idf ); + c[ 8 ] = polyval_co21( idf ); + c[ 9 ] = polyval_co22( idf ); + + // Result is then an odd polynomial in v (see Eq 56 of Shaw)... + return c0 + ( v * evalpoly( c, 10, v * v ) ); +} + +/** +* Evaluates Student's t quantiles via a tail series expansion. Tail and body series are due to Shaw. +* +* ## References +* +* - Shaw, William T. 2006. "Sampling Student's T distribution – use of the inverse cumulative distribution function." _Journal of Computational Finance_ 9 (4): 37–73. [www.mth.kcl.ac.uk/~shaww/web\_page/papers/Tdistribution06.pdf](www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf). +* +* @param df degrees of freedom +* @param v function value +* @return tail value +*/ +static double inverseStudentsTTailSeries( const double df, const double v ) { + double result; + double power; + double d[ 7 ]; + double div; + double np2; + double np4; + double np6; + double rn; + double w; + + // Tail series expansion, see section 6 of Shaw's paper. `w` is calculated using Eq 60: + w = stdlib_base_gamma_delta_ratio( df / 2.0, 0.5 ) * stdlib_base_sqrt( df * PI ) * v; + + // Define some variables: + np2 = df + 2.0; + np4 = df + 4.0; + np6 = df + 6.0; + + d[ 0 ] = 1.0; + d[ 1 ] = -( df + 1.0 ) / ( 2.0 * np2 ); + np2 *= ( df + 2.0 ); + d[ 2 ] = -df * ( df + 1.0 ) * ( df + 3.0 ) / ( 8.0 * np2 * np4 ); + np2 *= df + 2.0; + d[ 3 ] = -df * ( df + 1.0 ) * ( df + 5.0 ) * ( ( ( ( 3.0 * df ) + 7.0 ) * df ) - 2.0 ) / ( 48.0 * np2 * np4 * np6 ); + np2 *= ( df + 2.0 ); + np4 *= ( df + 4.0 ); + d[ 4 ] = -df * ( df + 1.0 ) * ( df + 7.0 ) * ( ( ( ( ( ( ( ( ( 15.0 * df ) + 154.0 ) * df ) + 465.0 ) * df + 286.0 ) * df ) - 336.0 ) * df ) + 64.0 ) / ( 384.0 * np2 * np4 * np6 * ( df + 8.0 ) ); + np2 *= ( df + 2.0 ); + d[ 5 ] = -df * ( df + 1.0 ) * ( df + 3.0 ) * ( df + 9.0 ) * ( ( ( ( ( ( ( ( ( ( ( 35.0 * df ) + 452.0 ) * df ) + 1573.0 ) * df ) + 600.0 ) * df ) - 2020.0 ) * df ) + 928.0 ) * df - 128.0 ) / ( 1280.0 * np2 * np4 * np6 * ( df + 8.0 ) * ( df + 10.0 ) ); + np2 *= ( df + 2.0 ); + np4 *= ( df + 4.0 ); + np6 *= ( df + 6.0 ); + d[ 6 ] = -df * ( df + 1.0 ) * ( df + 11.0 ) * ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( 945.0 * df ) + 31506.0 ) * df ) + 425858.0 ) * df ) + 2980236.0 ) * df ) + 11266745.0 ) * df ) + 20675018.0 ) * df ) + 7747124.0 ) * df ) - 22574632.0 ) * df ) - 8565600.0 ) * df ) + 18108416.0 ) * df ) - 7099392.0 ) * df ) + 884736.0 ) ) / ( 46080.0 * np2 * np4 * np6 * ( df + 8.0 ) * ( df + 10.0 ) * ( df + 12.0 ) ); + + // Now bring everything together to provide the result this is Eq 62 of Shaw: + rn = stdlib_base_sqrt( df ); + div = stdlib_base_pow( rn * w, 1.0 / df ); + power = div * div; + result = evalpoly( d, 7, power ); + result *= rn; + result /= div; + return -result; +} + +/** +* Evaluates Student's t quantiles. +* +* @param df degrees of freedom +* @param u input probability +* @param v probability equal to `1-u` +* @return function value +*/ +static double inverseStudentsT( const double df, const double u, const double v ) { + double crossover; + double tolerance; + double rootAlpha; + bool invert; + double result; + double alpha; + double tmp; + double p0; + double p2; + double p4; + double p5; + double uc; + double vc; + double p; + double r; + double x; + double a; + double b; + + uc = u; + vc = v; + + result = 0.0; + if ( uc > vc ) { + // Function is symmetric, so invert it: + tmp = vc; + vc = uc; + uc = tmp; + invert = true; + } else { + invert = false; + } + if ( stdlib_base_floor( df ) == df && df < 20.0 ) { + // We have integer degrees of freedom, try for the special cases first... + + // NOTE: the original JavaScript implementation computes `ldexp( 1.0, (2.0*53.0)/3.0 )`, in which the fractional exponent is truncated toward zero (i.e., `2^35`). We faithfully replicate that behavior: + tolerance = stdlib_base_ldexp( 1.0, 35 ); + + switch ( (int32_t)stdlib_base_floor( df ) ) { + case 1: + // `df = 1` is the same as the Cauchy distribution, see Shaw Eq 35: + if ( uc == 0.5 ) { + result = 0.0; + } else { + result = -stdlib_base_cos( PI * uc ) / stdlib_base_sin( PI * uc ); + } + break; + case 2: + // `df = 2` has an exact result, see Shaw Eq 36: + result = ( ( 2.0 * uc ) - 1.0 ) / stdlib_base_sqrt( 2.0 * uc * vc ); + break; + case 4: + // `df = 4` has an exact result, see Shaw Eq 38 & 39: + alpha = 4.0 * uc * vc; + rootAlpha = stdlib_base_sqrt( alpha ); + r = 4.0 * stdlib_base_cos( stdlib_base_acos( rootAlpha ) / 3.0 ) / rootAlpha; + x = stdlib_base_sqrt( r - 4.0 ); + result = ( uc - 0.5 < 0.0 ) ? -x : x; + break; + case 6: + // We get numeric overflow in this area: + if ( uc < 1.0e-150 ) { + return ( ( invert ) ? -1.0 : 1.0 ) * inverseStudentsTHill( df, uc ); + } + // Newton-Raphson iteration of a polynomial case, choice of seed value is taken from Shaw's online supplement: + a = 4.0 * ( uc - ( uc * uc ) ); // 1 - 4 * (u - 0.5f) * (u - 0.5f); + b = stdlib_base_pow( a, ONE_THIRD ); + p = 6.0 * ( 1.0 + ( CONST_C * ( ( 1.0 / b ) - 1.0 ) ) ); + do { + p2 = p * p; + p4 = p2 * p2; + p5 = p * p4; + p0 = p; + + // Next term is given by Eq 41: + p = 2.0 * ( ( 8.0 * a * p5 ) - ( 270.0 * p2 ) + 2187.0 ) / ( 5.0 * ( ( 4.0 * a * p4 ) - ( 216.0 * p ) - 243.0 ) ); + } while ( stdlib_base_abs( ( p - p0 ) / p ) > tolerance ); + + // Use Eq 45 to extract the result: + p = stdlib_base_sqrt( p - df ); + result = ( uc - 0.5 < 0.0 ) ? -p : p; + break; + default: + if ( df > DF_THRESHOLD ) { // 2^28 + result = stdlib_base_erfcinv( 2.0 * uc ) * SQRT2; + } else if ( df < 3.0 ) { + // Use a roughly linear scheme to choose between Shaw's tail series and body series: + crossover = 0.2742 - ( df * 0.0242143 ); + if ( uc > crossover ) { + result = inverseStudentsTBodySeries( df, uc ); + } else { + result = inverseStudentsTTailSeries( df, uc ); + } + } else { + // Use Hill's method except in the extreme tails where we use Shaw's tail series. The crossover point is roughly exponential in -df: + crossover = stdlib_base_ldexp( 1.0, (int32_t)stdlib_base_round( df / -0.654 ) ); + if ( uc > crossover ) { + result = inverseStudentsTHill( df, uc ); + } else { + result = inverseStudentsTTailSeries( df, uc ); + } + } + } + } else if ( df > DF_THRESHOLD ) { + result = -stdlib_base_erfcinv( 2.0 * uc ) * SQRT2; + } else if ( df < 3.0 ) { + // Use a roughly linear scheme to choose between Shaw's tail series and body series: + crossover = 0.2742 - ( df * 0.0242143 ); + if ( uc > crossover ) { + result = inverseStudentsTBodySeries( df, uc ); + } else { + result = inverseStudentsTTailSeries( df, uc ); + } + } else { + // Use Hill's method except in the extreme tails where we use Shaw's tail series. The crossover point is roughly exponential in -df: + crossover = stdlib_base_ldexp( 1.0, (int32_t)stdlib_base_round( df / -0.654 ) ); + if ( uc > crossover ) { + result = inverseStudentsTHill( df, uc ); + } else { + result = inverseStudentsTTailSeries( df, uc ); + } + } + return ( invert ) ? -result : result; +} + +/** +* Returns the inverse of the incomplete beta function via the Student t distribution. +* +* ## Notes +* +* - The function writes one minus the returned value to `py`. +* +* @param a function parameter +* @param p probability value +* @param py destination for one minus the returned value +* @return function value +*/ +static double findIBetaInvFromTDist( const double a, const double p, double *py ) { + double df; + double u; + double v; + double t; + + u = p / 2.0; + v = 1.0 - u; + df = a * 2.0; + t = inverseStudentsT( df, u, v ); + *py = t * t / ( df + ( t * t ) ); + return df / ( df + ( t * t ) ); +} + +/** +* Computes the inverse of the lower incomplete beta function. +* +* ## Notes +* +* - The function writes the function value `y` to `out1` and `1-y` to `out2`. +* +* @param a function parameter +* @param b function parameter +* @param p function parameter +* @param q probability equal to `1-p` +* @param out1 destination for the function value +* @param out2 destination for one minus the function value +* +* @example +* double out1; +* double out2; +* stdlib_base_kernel_betaincinv( 3.0, 3.0, 0.2, 0.8, &out1, &out2 ); +*/ +void stdlib_base_kernel_betaincinv( const double a, const double b, const double p, const double q, double *out1, double *out2 ) { + double terms[ 5 ]; + int32_t digits; + bool invert; + double lambda; + double deriv; + double lower; + double theta; + double upper; + double maxv; + double minv; + double bet; + double ppa; + double tmp; + double xs2; + double ap1; + double bm1; + double ac; + double bc; + double pc; + double qc; + double fs; + double fv; + double lx; + double ps; + double xg; + double xs; + double yp; + double a2; + double a3; + double b2; + double r; + double l; + double u; + double x; + double y; + + ac = a; + bc = b; + pc = p; + qc = q; + + // The flag invert is set to true if we swap a for b and p for q, in which case the result has to be subtracted from 1: + invert = false; + + // Handle trivial cases first... + if ( qc == 0.0 ) { + *out1 = 1.0; + *out2 = 0.0; + return; + } + if ( pc == 0.0 ) { + *out1 = 0.0; + *out2 = 1.0; + return; + } + if ( ac == 1.0 ) { + if ( bc == 1.0 ) { + *out1 = pc; + *out2 = 1.0 - pc; + return; + } + // Change things around so we can handle as b == 1 special case below: + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + invert = true; + } + // Depending upon which approximation method we use, we may end up calculating either x or y initially (where y = 1-x): + x = 0.0; + y = 0.0; + + // For some of the methods we can put tighter bounds on the result than simply [0,1]: + lower = 0.0; + upper = 1.0; + + // Student's T with b = 0.5 gets handled as a special case, swap around if the arguments are in the "wrong" order: + if ( ac == 0.5 ) { + if ( bc == 0.5 ) { + x = stdlib_base_sin( pc * HALF_PI ); + x *= x; + y = stdlib_base_sin( qc * HALF_PI ); + y *= y; + *out1 = x; + *out2 = y; + return; + } + if ( bc > 0.5 ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + invert = !invert; + } + } + // Select calculation method for the initial estimate: + if ( bc == 0.5 && ac >= 0.5 && pc != 1.0 ) { + // We have a Student's T distribution: + yp = 0.0; + x = findIBetaInvFromTDist( ac, pc, &yp ); + y = yp; + } + else if ( bc == 1.0 ) { + if ( pc < qc ) { + if ( ac > 1.0 ) { + x = stdlib_base_pow( pc, 1.0 / ac ); + y = -stdlib_base_expm1( stdlib_base_ln( pc ) / ac ); + } else { + x = stdlib_base_pow( pc, 1.0 / ac ); + y = 1.0 - x; + } + } else { + x = stdlib_base_exp( stdlib_base_log1p( -qc ) / ac ); + y = -stdlib_base_expm1( stdlib_base_log1p( -qc ) / ac ); + } + if ( invert ) { + tmp = y; + y = x; + x = tmp; + } + *out1 = x; + *out2 = y; + return; + } + else if ( ac + bc > 5.0 ) { + // When a+b is large then we can use one of Prof Temme's asymptotic expansions, begin by swapping things around so that p < 0.5, we do this to avoid cancellations errors when p is large. + if ( pc > 0.5 ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + invert = !invert; + } + minv = stdlib_base_min( ac, bc ); + maxv = stdlib_base_max( ac, bc ); + if ( ( stdlib_base_sqrt( minv ) > ( maxv - minv ) ) && minv > 5.0 ) { + // When a and b differ by a small amount the curve is quite symmetrical and we can use an error function to approximate the inverse. This is the cheapest of the three Temme expansions, and the calculated value for x will never be much larger than p, so we don't have to worry about cancellation as long as p is small. + x = temme1( ac, bc, pc ); + y = 1.0 - x; + } else { + r = ac + bc; + theta = stdlib_base_asin( stdlib_base_sqrt( ac / r ) ); + lambda = minv / r; + if ( + lambda >= 0.2 && + lambda <= 0.8 && + r >= 10.0 + ) { + // The second error function case is the next cheapest to use, it breaks down when the result is likely to be very small, if `a+b` is also small, but we can use a cheaper expansion there in any case. As before `x` won't be much larger than `p`, so as long as `p` is small we should be free of cancellation error. + ppa = stdlib_base_pow( pc, 1.0 / ac ); + if ( ppa < 0.0025 && ( ac + bc ) < 200.0 ) { + x = ppa * stdlib_base_pow( ac * stdlib_base_beta( ac, bc ), 1.0 / ac ); + } else { + x = temme2( pc, r, theta ); + } + y = 1.0 - x; + } else { + // If we get here then a and b are very different in magnitude and we need to use the third of Temme's methods which involves inverting the incomplete gamma. This is much more expensive than the other methods. We also can only use this method when a > b, which can lead to cancellation errors if we really want y (as we will when x is close to 1), so a different expansion is used in that case. + if ( ac < bc ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + invert = !invert; + } + // Try to compute the easy way first: + bet = 0.0; + if ( bc < 2.0 ) { + bet = stdlib_base_beta( ac, bc ); + } + if ( bet == 0.0 ) { + y = 1.0; + } else { + y = stdlib_base_pow( bc * qc * bet, 1.0 / bc ); + x = 1.0 - y; + } + } + if ( y > 1.0e-5 ) { + x = temme3( ac, bc, pc, qc ); + y = 1.0 - x; + } + } + } + else if ( ac < 1.0 && bc < 1.0 ) { + // Both a and b less than 1, there is a point of inflection at xs: + xs = ( 1.0 - ac ) / ( 2.0 - ac - bc ); + + // Now we need to ensure that we start our iteration from the right side of the inflection point: + stdlib_base_kernel_betainc( xs, ac, bc, true, false, &fv, &deriv ); + fs = fv - pc; + if ( stdlib_base_abs( fs ) / pc < EPSILON * 3.0 ) { + // The result is at the point of inflection, best just return it: + if ( invert ) { + *out1 = 1.0 - xs; + *out2 = xs; + return; + } + *out1 = xs; + *out2 = 1.0 - xs; + return; + } + if ( fs < 0.0 ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + invert = !invert; + xs = 1.0 - xs; + } + xg = stdlib_base_pow( ac * pc * stdlib_base_beta( ac, bc ), 1.0 / ac ); + x = xg / ( 1.0 + xg ); + y = 1.0 / ( 1.0 + xg ); + + // And finally we know that our result is below the inflection point, so set an upper limit on our search: + if ( x > xs ) { + x = xs; + } + upper = xs; + } + else if ( ac > 1.0 && bc > 1.0 ) { + // Small a and b, both greater than 1, there is a point of inflection at xs, and its complement is xs2, we must always start our iteration from the right side of the point of inflection. + xs = ( ac - 1.0 ) / ( ac + bc - 2.0 ); + xs2 = ( bc - 1.0 ) / ( ac + bc - 2.0 ); + stdlib_base_kernel_betainc( xs, ac, bc, true, false, &fv, &deriv ); + ps = fv - pc; + + if ( ps < 0.0 ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + tmp = xs2; + xs2 = xs; // cppcheck-suppress unreadVariable + xs = tmp; + + invert = !invert; + } + // Estimate x and y, using expm1 to get a good estimate for y when it's very small: + lx = stdlib_base_ln( pc * ac * stdlib_base_beta( ac, bc ) ) / ac; + x = stdlib_base_exp( lx ); + y = ( x < 0.9 ) ? 1.0 - x : -stdlib_base_expm1( lx ); + + if ( bc < ac && x < 0.2 ) { + // Under a limited range of circumstances we can improve our estimate for x... + ap1 = ac - 1.0; + bm1 = bc - 1.0; + a2 = ac * ac; + a3 = ac * a2; + b2 = bc * bc; + terms[ 0 ] = 0.0; + terms[ 1 ] = 1.0; + terms[ 2 ] = bm1 / ap1; + ap1 *= ap1; + terms[ 3 ] = bm1 * ( ( 3.0 * ac * bc ) + ( 5.0 * bc ) + a2 - ac - 4.0 ) / ( 2.0 * ( ac + 2.0 ) * ap1 ); + ap1 *= ( ac + 1.0 ); + terms[ 4 ] = bm1 * ( ( 33.0 * ac * b2 ) + ( 31.0 * b2 ) + ( 8.0 * a2 * b2 ) - ( 30.0 * ac * bc ) - ( 47.0 * bc ) + ( 11.0 * a2 * bc ) + ( 6.0 * a3 * bc ) + 18.0 + ( 4.0 * ac ) - a3 + ( a2 * a2 ) - ( 10.0 * a2 ) ); + terms[ 4 ] /= ( 3.0 * ( ac + 3.0 ) * ( ac + 2.0 ) * ap1 ); + x = evalpoly( terms, 5, x ); + } + // Know that result is below the inflection point, so set an upper limit on search... + if ( x > xs ) { + x = xs; + } + upper = xs; + } else { + // Case: ( a <= 1 ) != ( b <= 1 ). If all else fails we get here, only one of a and b is above 1, and a+b is small. Start by swapping things around so that we have a concave curve with b > a and no points of inflection in [0,1]. As long as we expect x to be small then we can use the simple (and cheap) power term to estimate x, but when we expect x to be large then this greatly underestimates x and leaves us trying to iterate "round the corner" which may take almost forever. We could use Temme's inverse gamma function case in that case, this works really rather well (albeit expensively) even though strictly speaking we're outside it's defined range. However it's expensive to compute, and an alternative approach which models the curve as a distorted quarter circle is much cheaper to compute, and still keeps the number of iterations required down to a reasonable level. With thanks to Prof. Temme for this suggestion. + if ( bc < ac ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + invert = !invert; + } + if ( stdlib_base_pow( pc, 1.0 / ac ) < 0.5 ) { + x = stdlib_base_pow( pc * ac * stdlib_base_beta( ac, bc ), 1.0 / ac ); + if ( x == 0.0 ) { + x = FLOAT64_MIN_NORM; + } + y = 1.0 - x; + } + // Case: pow(q, 1/b) < 0.1 + else { + // Model a distorted quarter circle: + y = stdlib_base_pow( 1.0 - stdlib_base_pow( pc, bc * stdlib_base_beta( ac, bc ) ), 1.0 / bc ); + if ( y == 0.0 ) { + y = FLOAT64_MIN_NORM; + } + x = 1.0 - y; + } + } + // Now we have a guess for x (and for y) we can set things up for iteration. If x > 0.5 it pays to swap things round: + if ( x > 0.5 ) { + tmp = bc; + bc = ac; + ac = tmp; + + tmp = qc; + qc = pc; + pc = tmp; + + tmp = y; + y = x; // cppcheck-suppress unreadVariable + x = tmp; + + invert = !invert; + l = 1.0 - upper; + u = 1.0 - lower; + lower = l; + upper = u; + } + // Lower bound for our search: We're not interested in denormalized answers as these tend to take up lots of iterations, given that we can't get accurate derivatives in this area (they tend to be infinite). + if ( lower == 0.0 ) { + if ( invert ) { + // We're not interested in answers smaller than machine epsilon: + lower = EPSILON; + if ( x < lower ) { + x = lower; + } + } else { + lower = FLOAT64_MIN_NORM; + } + if ( x < lower ) { + x = lower; + } + } + // Figure out how many digits to iterate towards: + digits = DIGITS; + if ( x < 1.0e-50 && ( ac < 1.0 || bc < 1.0 ) ) { + // If we're in a region where the first derivative is very large, then we have to take care that the root-finder doesn't terminate prematurely. We'll bump the precision up to avoid this, but we have to take care not to set the precision too high or the last few iterations will just thrash around and convergence may be slow in this case. Try 3/4 of machine epsilon: + digits *= 3; + digits /= 2; + } + // Now iterate, we can use either p or q as the target here depending on which is smaller: + x = halleyIterate( ac, bc, ( ( pc < qc ) ? pc : qc ), pc >= qc, x, lower, upper, digits, MAX_ITERATIONS ); + + // Tidy up, if we "lower" was too high then zero is the best answer we have: + if ( x == lower ) { + x = 0.0; + } + if ( invert ) { + *out1 = 1.0 - x; + *out2 = x; + return; + } + *out1 = x; + *out2 = 1.0 - x; + return; +} diff --git a/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/test/test.native.js new file mode 100644 index 000000000000..e329c9348326 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/kernel-betaincinv/test/test.native.js @@ -0,0 +1,68 @@ +/** +* @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 isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' ); +var tryRequire = require( '@stdlib/utils/try-require' ); + + +// VARIABLES // + +var kernelBetaincinv = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( kernelBetaincinv instanceof Error ) +}; + + +// FIXTURES // + +var data = require( './fixtures/cpp/output.json' ); + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof kernelBetaincinv, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'the function evaluates the inverse of the lower regularized incomplete beta function', opts, function test( t ) { + var expected; + var ULP = 20; + var y; + var i; + var a; + var b; + var p; + + p = data.p; + a = data.a; + b = data.b; + expected = data.y; + for ( i = 0; i < p.length; i++ ) { + y = kernelBetaincinv( a[i], b[i], p[i], 1.0-p[i] )[ 0 ]; + t.strictEqual( isAlmostSameValue( y, expected[ i ], ULP ), true, 'returns expected value' ); + } + t.end(); +});