From ac3ec1f58ce40a9108df12de7eebcb56d8a6cc68 Mon Sep 17 00:00:00 2001 From: Philipp Burckhardt Date: Fri, 10 Jul 2026 20:24:11 -0500 Subject: [PATCH] feat: add C implementation for `math/base/special/polygamma` --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown_pkg_readmes status: passed - task: lint_markdown_docs status: na - task: lint_markdown status: na - task: lint_package_json status: passed - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: passed - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: missing_dependencies - task: lint_c_examples status: missing_dependencies - task: lint_c_benchmarks status: missing_dependencies - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../math/base/special/polygamma/README.md | 100 +++ .../polygamma/benchmark/benchmark.native.js | 65 ++ .../polygamma/benchmark/c/native/Makefile | 146 ++++ .../polygamma/benchmark/c/native/benchmark.c | 142 ++++ .../math/base/special/polygamma/binding.gyp | 170 ++++ .../special/polygamma/examples/c/Makefile | 146 ++++ .../special/polygamma/examples/c/example.c | 42 + .../math/base/special/polygamma/include.gypi | 53 ++ .../stdlib/math/base/special/polygamma.h | 40 + .../math/base/special/polygamma/lib/native.js | 67 ++ .../math/base/special/polygamma/manifest.json | 145 ++++ .../math/base/special/polygamma/package.json | 3 + .../special/polygamma/scripts/evalpoly.js | 84 ++ .../math/base/special/polygamma/src/Makefile | 70 ++ .../math/base/special/polygamma/src/addon.c | 22 + .../math/base/special/polygamma/src/main.c | 752 ++++++++++++++++++ .../special/polygamma/test/test.native.js | 267 +++++++ 17 files changed, 2314 insertions(+) create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/benchmark.native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/benchmark.c create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/binding.gyp create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/examples/c/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/examples/c/example.c create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/include.gypi create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/include/stdlib/math/base/special/polygamma.h create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/lib/native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/manifest.json create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/src/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/src/addon.c create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/src/main.c create mode 100644 lib/node_modules/@stdlib/math/base/special/polygamma/test/test.native.js diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/README.md b/lib/node_modules/@stdlib/math/base/special/polygamma/README.md index 8163af4edca0..e23d3b320e2e 100644 --- a/lib/node_modules/@stdlib/math/base/special/polygamma/README.md +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/README.md @@ -123,6 +123,106 @@ logEachMap( 'ψ^(%d)(x = %0.4f): %0.4f', n, x, polygamma ); + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/polygamma.h" +``` + +#### stdlib_base_polygamma( n, x ) + +Evaluates the [polygamma function][polygamma-function] of order `n`; i.e., the `(n+1)`th derivative of the [natural logarithm][@stdlib/math/base/special/ln] of the [gamma function][@stdlib/math/base/special/gamma]. + +```c +double out = stdlib_base_polygamma( 3, 1.2 ); +// returns ~3.245 + +out = stdlib_base_polygamma( 5, 1.2 ); +// returns ~41.39 +``` + +The function accepts the following arguments: + +- **n**: `[in] int32_t` order of derivative (nonnegative integer). +- **x**: `[in] double` input value. + +```c +double stdlib_base_polygamma( const int32_t n, const double x ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/polygamma.h" +#include "stdlib/math/base/special/floor.h" +#include +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + int32_t n; + double x; + double y; + int i; + + for ( i = 0; i < 25; i++ ) { + n = (int32_t)stdlib_base_floor( random_uniform( 0.0, 51.0 ) ); + x = random_uniform( -50.0, 50.0 ); + y = stdlib_base_polygamma( n, x ); + printf( "polygamma(%d, %lf) = %lf\n", n, x, y ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/benchmark.native.js new file mode 100644 index 000000000000..3782e103065b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/benchmark.native.js @@ -0,0 +1,65 @@ +/** +* @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 randu = require( '@stdlib/random/base/randu' ); +var round = require( '@stdlib/math/base/special/round' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var format = require( '@stdlib/string/format' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var polygamma = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( polygamma instanceof Error ) +}; + + +// MAIN // + +bench( format( '%s::native', pkg ), opts, function benchmark( b ) { + var n; + var x; + var y; + var i; + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + n = round( randu()*40.0 ); + x = ( randu()*100.0 ) + EPS; + y = polygamma( n, x ); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/Makefile new file mode 100644 index 000000000000..979768abbcec --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/benchmark.c new file mode 100644 index 000000000000..bf60871e1e5c --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/benchmark/c/native/benchmark.c @@ -0,0 +1,142 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#include "stdlib/math/base/special/polygamma.h" +#include "stdlib/math/base/special/round.h" +#include "stdlib/constants/float64/eps.h" +#include +#include +#include +#include +#include + +#define NAME "polygamma" +#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 ) { + int32_t n[ 100 ]; + double x[ 100 ]; + double elapsed; + double t; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + n[ i ] = (int32_t)stdlib_base_round( random_uniform( 0.0, 40.0 ) ); + x[ i ] = random_uniform( STDLIB_CONSTANT_FLOAT64_EPS, 100.0 + STDLIB_CONSTANT_FLOAT64_EPS ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_polygamma( n[ i % 100 ], x[ i % 100 ] ); + if ( y != y ) { + printf( "should not return NaN\n" ); + break; + } + } + elapsed = tic() - t; + if ( y != y ) { + printf( "should not return NaN\n" ); + } + return elapsed; +} + +/** +* Main execution sequence. +*/ +int main( void ) { + double elapsed; + int i; + + // Use the current time to seed the random number generator: + srand( time( NULL ) ); + + print_version(); + for ( i = 0; i < REPEATS; i++ ) { + printf( "# c::native::%s\n", NAME ); + elapsed = benchmark(); + print_results( elapsed ); + printf( "ok %d benchmark finished\n", i+1 ); + } + print_summary( REPEATS, REPEATS ); +} diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/binding.gyp b/lib/node_modules/@stdlib/math/base/special/polygamma/binding.gyp new file mode 100644 index 000000000000..0d6508a12e99 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/polygamma/examples/c/Makefile new file mode 100644 index 000000000000..c8f8e9a1517b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/polygamma/examples/c/example.c new file mode 100644 index 000000000000..5732e0227ad3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma.h" +#include "stdlib/math/base/special/floor.h" +#include +#include +#include + +static double random_uniform( const double min, const double max ) { + double v = (double)rand() / ( (double)RAND_MAX + 1.0 ); + return min + ( v*(max-min) ); +} + +int main( void ) { + int32_t n; + double x; + double y; + int i; + + for ( i = 0; i < 25; i++ ) { + n = (int32_t)stdlib_base_floor( random_uniform( 0.0, 51.0 ) ); + x = random_uniform( -50.0, 50.0 ); + y = stdlib_base_polygamma( n, x ); + printf( "polygamma(%d, %lf) = %lf\n", n, x, y ); + } +} diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/include.gypi b/lib/node_modules/@stdlib/math/base/special/polygamma/include.gypi new file mode 100644 index 000000000000..bee8d41a2caf --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/include.gypi @@ -0,0 +1,53 @@ +# @license Apache-2.0 +# +# Copyright (c) 2026 The Stdlib Authors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# A GYP include file for building a Node.js native add-on. +# +# Main documentation: +# +# [1]: https://gyp.gsrc.io/docs/InputFormatReference.md +# [2]: https://gyp.gsrc.io/docs/UserDocumentation.md +{ + # Define variables to be used throughout the configuration for all targets: + 'variables': { + # Source directory: + 'src_dir': './src', + + # Include directories: + 'include_dirs': [ + ' + +/* +* If C++, prevent name mangling so that the compiler emits a binary file having undecorated names, thus mirroring the behavior of a C compiler. +*/ +#ifdef __cplusplus +extern "C" { +#endif + +/** +* Evaluates the polygamma function. +*/ +double stdlib_base_polygamma( const int32_t n, const double x ); + +#ifdef __cplusplus +} +#endif + +#endif // !STDLIB_MATH_BASE_SPECIAL_POLYGAMMA_H diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/lib/native.js b/lib/node_modules/@stdlib/math/base/special/polygamma/lib/native.js new file mode 100644 index 000000000000..39c9cd3acbe3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/lib/native.js @@ -0,0 +1,67 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +'use strict'; + +// MODULES // + +var addon = require( './../src/addon.node' ); + + +// MAIN // + +/** +* Evaluates the polygamma function. +* +* @private +* @param {NonNegativeInteger} n - order of derivative +* @param {number} x - input value +* @returns {number} (n+1)'th derivative +* +* @example +* var v = polygamma( 3, 1.2 ); +* // returns ~3.245 +* +* @example +* var v = polygamma( 5, 1.2 ); +* // returns ~41.39 +* +* @example +* var v = polygamma( 3, -4.9 ); +* // returns ~60014.239 +* +* @example +* var v = polygamma( -1, 5.3 ); +* // returns NaN +* +* @example +* var v = polygamma( 2, -2.0 ); +* // returns NaN +* +* @example +* var v = polygamma( 1, NaN ); +* // returns NaN +*/ +function polygamma( n, x ) { + return addon( n, x ); +} + + +// EXPORTS // + +module.exports = polygamma; diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/manifest.json b/lib/node_modules/@stdlib/math/base/special/polygamma/manifest.json new file mode 100644 index 000000000000..b55c79545a76 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/manifest.json @@ -0,0 +1,145 @@ +{ + "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/math/base/napi/binary", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/bernoulli", + "@stdlib/math/base/special/cospi", + "@stdlib/math/base/special/digamma", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/factorial", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/riemann-zeta", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sinpi", + "@stdlib/math/base/special/trigamma", + "@stdlib/math/base/special/trunc", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/ln-pi", + "@stdlib/constants/float64/ln-two", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/max-ln", + "@stdlib/constants/float64/ninf", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/pinf" + ] + }, + { + "task": "benchmark", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/bernoulli", + "@stdlib/math/base/special/cospi", + "@stdlib/math/base/special/digamma", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/factorial", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/riemann-zeta", + "@stdlib/math/base/special/round", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sinpi", + "@stdlib/math/base/special/trigamma", + "@stdlib/math/base/special/trunc", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/ln-pi", + "@stdlib/constants/float64/ln-two", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/max-ln", + "@stdlib/constants/float64/ninf", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/pinf" + ] + }, + { + "task": "examples", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/bernoulli", + "@stdlib/math/base/special/cospi", + "@stdlib/math/base/special/digamma", + "@stdlib/math/base/special/exp", + "@stdlib/math/base/special/factorial", + "@stdlib/math/base/special/floor", + "@stdlib/math/base/special/gammaln", + "@stdlib/math/base/special/ldexp", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/min", + "@stdlib/math/base/special/pow", + "@stdlib/math/base/special/riemann-zeta", + "@stdlib/math/base/special/signum", + "@stdlib/math/base/special/sinpi", + "@stdlib/math/base/special/trigamma", + "@stdlib/math/base/special/trunc", + "@stdlib/constants/float64/eps", + "@stdlib/constants/float64/ln-pi", + "@stdlib/constants/float64/ln-two", + "@stdlib/constants/float64/max", + "@stdlib/constants/float64/max-ln", + "@stdlib/constants/float64/ninf", + "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/pinf" + ] + } + ] +} diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/package.json b/lib/node_modules/@stdlib/math/base/special/polygamma/package.json index bc582aa672e3..c2f88e2d06a7 100644 --- a/lib/node_modules/@stdlib/math/base/special/polygamma/package.json +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma/scripts/evalpoly.js b/lib/node_modules/@stdlib/math/base/special/polygamma/scripts/evalpoly.js index 768215aeb8ff..e37122d63346 100644 --- a/lib/node_modules/@stdlib/math/base/special/polygamma/scripts/evalpoly.js +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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 // @@ -102,6 +107,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 // /** @@ -111,7 +143,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; */ function main() { var fpath; + var copts; var opts; + var file; var str; opts = { @@ -157,6 +191,56 @@ function main() { fpath = resolve( __dirname, '..', 'lib', 'polyval_p12.js' ); str = header + compile( P12 ); writeFileSync( fpath, str, opts ); + + copts = { + 'dtype': 'double', + 'name': '' + }; + + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + copts.name = 'polyval_p3'; + str = compileC( P3, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p4'; + str = compileC( P4, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p5'; + str = compileC( P5, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p6'; + str = compileC( P6, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p7'; + str = compileC( P7, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p8'; + str = compileC( P8, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p9'; + str = compileC( P9, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p10'; + str = compileC( P10, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p11'; + str = compileC( P11, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_p12'; + str = compileC( P12, copts ); + file = insert( file, copts.name, str ); + + writeFileSync( fpath, file, opts ); } main(); diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/src/Makefile b/lib/node_modules/@stdlib/math/base/special/polygamma/src/Makefile new file mode 100644 index 000000000000..2caf905cedbe --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/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/polygamma/src/addon.c b/lib/node_modules/@stdlib/math/base/special/polygamma/src/addon.c new file mode 100644 index 000000000000..89590b7984b5 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/src/addon.c @@ -0,0 +1,22 @@ +/** +* @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/polygamma.h" +#include "stdlib/math/base/napi/binary.h" + +STDLIB_MATH_BASE_NAPI_MODULE_ID_D( stdlib_base_polygamma ) diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/src/main.c b/lib/node_modules/@stdlib/math/base/special/polygamma/src/main.c new file mode 100644 index 000000000000..837c3381da93 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/src/main.c @@ -0,0 +1,752 @@ +/** +* @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]{@link http://www.boost.org/doc/libs/1_65_0/boost/math/special_functions/detail/polygamma.hpp}. The implementation follows the original but has been reformatted and modified for use in stdlib. +* +* ```text +* (C) Copyright Nikhar Agrawal 2013. +* (C) Copyright Christopher Kormanyos 2013. +* (C) Copyright John Maddock 2014. +* (C) Copyright Paul Bristow 2013. +* +* 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/polygamma.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/math/base/special/bernoulli.h" +#include "stdlib/math/base/special/cospi.h" +#include "stdlib/math/base/special/digamma.h" +#include "stdlib/math/base/special/exp.h" +#include "stdlib/math/base/special/factorial.h" +#include "stdlib/math/base/special/floor.h" +#include "stdlib/math/base/special/gammaln.h" +#include "stdlib/math/base/special/ldexp.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/special/min.h" +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/riemann_zeta.h" +#include "stdlib/math/base/special/signum.h" +#include "stdlib/math/base/special/sinpi.h" +#include "stdlib/math/base/special/trigamma.h" +#include "stdlib/math/base/special/trunc.h" +#include "stdlib/constants/float64/eps.h" +#include "stdlib/constants/float64/ln_pi.h" +#include "stdlib/constants/float64/ln_two.h" +#include "stdlib/constants/float64/max.h" +#include "stdlib/constants/float64/max_ln.h" +#include "stdlib/constants/float64/ninf.h" +#include "stdlib/constants/float64/pi.h" +#include "stdlib/constants/float64/pinf.h" +#include +#include + +// NOTE: threshold expressions of the form `0.4 * DIGITS_BASE10` must remain computed products (`7.6000000000000005`), as the computed product differs from the literal `7.6` by one ULP: +static const double DIGITS_BASE10 = 19.0; +static const int32_t MAX_SERIES_ITERATIONS = 1000000; + +// NOTE: the following bound is intentionally `172` (per the Boost implementation), not the `170` overflow bound of the factorial function itself: +static const int32_t MAX_FACTORIAL = 172; + +// π raised to powers two to twelve (obtained from Wolfram Alpha): +static const double PI2 = 9.869604401089358; +static const double PI3 = 31.00627668029982; +static const double PI4 = 97.40909103400244; +static const double PI5 = 306.01968478528147; +static const double PI6 = 961.3891935753045; +static const double PI7 = 3020.2932277767923; +static const double PI8 = 9488.531016070574; +static const double PI9 = 29809.09933344621; +static const double PI10 = 93648.04747608303; +static const double PI11 = 294204.0179738906; +static const double PI12 = 924269.1815233742; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_p3 + +/** +* 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_p3( const double x ) { + return -2.0 + (x * -4.0); +} + +// END: polyval_p3 + +// BEGIN: polyval_p4 + +/** +* 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_p4( const double x ) { + return 16.0 + (x * 8.0); +} + +// END: polyval_p4 + +// BEGIN: polyval_p5 + +/** +* 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_p5( const double x ) { + return -16.0 + (x * (-88.0 + (x * -16.0))); +} + +// END: polyval_p5 + +// BEGIN: polyval_p6 + +/** +* 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_p6( const double x ) { + return 272.0 + (x * (416.0 + (x * 32.0))); +} + +// END: polyval_p6 + +// BEGIN: polyval_p7 + +/** +* 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_p7( const double x ) { + return -272.0 + (x * (-2880.0 + (x * (-1824.0 + (x * -64.0))))); +} + +// END: polyval_p7 + +// BEGIN: polyval_p8 + +/** +* 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_p8( const double x ) { + return 7936.0 + (x * (24576.0 + (x * (7680.0 + (x * 128.0))))); +} + +// END: polyval_p8 + +// BEGIN: polyval_p9 + +/** +* 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_p9( const double x ) { + return -7936.0 + (x * (-137216.0 + (x * (-185856.0 + (x * (-31616.0 + (x * -256.0))))))); +} + +// END: polyval_p9 + +// BEGIN: polyval_p10 + +/** +* 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_p10( const double x ) { + return 353792.0 + (x * (1841152.0 + (x * (1304832.0 + (x * (128512.0 + (x * 512.0))))))); +} + +// END: polyval_p10 + +// BEGIN: polyval_p11 + +/** +* 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_p11( const double x ) { + return -353792.0 + (x * (-9061376.0 + (x * (-21253376.0 + (x * (-8728576.0 + (x * (-518656.0 + (x * -1024.0))))))))); +} + +// END: polyval_p11 + +// BEGIN: polyval_p12 + +/** +* 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_p12( const double x ) { + return 22368256.0 + (x * (175627264.0 + (x * (222398464.0 + (x * (56520704.0 + (x * (2084864.0 + (x * 2048.0))))))))); +} + +// END: polyval_p12 + +/* End auto-generated functions. */ + +/** +* Evaluates a polynomial having an arbitrary number of coefficients. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation and mirrors `@stdlib/math/base/tools/evalpoly`, including the short-circuit when `x` equals zero. The short-circuit must be retained, as, for half-integer reflection arguments, the evaluation point is exactly zero, and plain Horner evaluation could otherwise return NaN when provided non-finite coefficients. +* +* [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; + + if ( len < 2 || x == 0.0 ) { + return ( len == 0 ) ? 0.0 : c[ 0 ]; + } + i = len - 1; + p = ( c[ i ] * x ) + c[ i-1 ]; + i -= 2; + while ( i >= 0 ) { + p = ( p * x ) + c[ i ]; + i -= 1; + } + return p; +} + +/** +* Evaluates the n'th derivative of `cot(πx)` at `x`. +* +* @param n derivative to evaluate +* @param x input value +* @param xc one minus `x` +* @return n'th derivative +*/ +static double polycotpi( const int32_t n, const double x, const double xc ) { + double powTerms; + double *wsp; + double *cur; + double *nxt; + double *tmp; + double sum; + double out; + double s; + double c; + int32_t noffset; // offset for next row + int32_t offset; // 1 if the first cos power is 0; otherwise 0 + int32_t ncols; // how many entries there are in the current row + int32_t mcols; // how many entries there will be in the next row + int32_t len; + int32_t mo; // largest order of the polynomial of cos terms + int32_t so; // order of the sin term + int32_t co; // order of the cosine term in entry "j" + int32_t i; + int32_t j; + int32_t k; + + // Use the argument having the smaller absolute value for improved accuracy (mathematically, `sinpi( x ) == sinpi( xc )`, as `x + xc == 1`): + s = ( stdlib_base_abs( x ) < stdlib_base_abs( xc ) ) ? stdlib_base_sinpi( x ) : stdlib_base_sinpi( xc ); + c = stdlib_base_cospi( x ); + switch ( n ) { + case 1: + return -STDLIB_CONSTANT_FLOAT64_PI / ( s * s ); + case 2: + return 2.0 * PI2 * c / stdlib_base_pow( s, 3.0 ); + case 3: + return PI3 * polyval_p3( c * c ) / stdlib_base_pow( s, 4.0 ); + case 4: + return PI4 * c * polyval_p4( c * c ) / stdlib_base_pow( s, 5.0 ); + case 5: + return PI5 * polyval_p5( c * c ) / stdlib_base_pow( s, 6.0 ); + case 6: + return PI6 * c * polyval_p6( c * c ) / stdlib_base_pow( s, 7.0 ); + case 7: + return PI7 * polyval_p7( c * c ) / stdlib_base_pow( s, 8.0 ); + case 8: + return PI8 * c * polyval_p8( c * c ) / stdlib_base_pow( s, 9.0 ); + case 9: + return PI9 * polyval_p9( c * c ) / stdlib_base_pow( s, 10.0 ); + case 10: + return PI10 * c * polyval_p10( c * c ) / stdlib_base_pow( s, 11.0 ); + case 11: + return PI11 * polyval_p11( c * c ) / stdlib_base_pow( s, 12.0 ); + case 12: + return PI12 * c * polyval_p12( c * c ) / stdlib_base_pow( s, 13.0 ); + default: + break; + } + // We'll have to compute the coefficients up to `n`; however, if the final evaluation would have too many terms, just bail out right away: + if ( ( (double)n / 2.0 ) > (double)MAX_SERIES_ITERATIONS ) { + return 0.0 / 0.0; // NaN + } + /* + * Compute the coefficients of the n'th derivative via the recurrence relation + * + * C[k-1, m+1] -= k * C[k, m] + * C[k+1, m+1] += (k-m-1) * C[k, m] + * + * seeded with `C[0, 1] = -1`, where row `m` stores only the non-zero coefficients (every other term in the polynomials is zero, so row `m` stores the coefficient of `cos^(2j + (m-1)&1)` at index `j`). + * + * The JavaScript implementation lazily grows a module-level memoization table (and Boost, similarly, a vector of vectors). As the recurrence only ever reads row `i` when writing row `i+1`, we instead recompute the final row on each call using two heap-allocated "ping-pong" row buffers, which is bit-identical to evaluating a fully memoized table (the recurrence is deterministic) and differs only in performance. Heap allocation is required, as row lengths grow as `n/2` (given the semantic bound `n <= 2000000`, a row may hold ~10^6 doubles (~8 MB), which cannot be stack-allocated), and a mutable static workspace would not be thread-safe. + */ + len = ( ( n - 1 ) / 2 ) + 1; // length of the final row + wsp = (double *)malloc( (size_t)( 2 * len ) * sizeof( double ) ); + if ( wsp == NULL ) { + return 0.0 / 0.0; // NaN + } + cur = wsp; + nxt = wsp + len; + cur[ 0 ] = -1.0; // row 0 (first derivative): single coefficient `-1.0` + + // Build rows `1..n-1` (note: row `r` has length `floor(r/2)+1 <= len`, and reads of the current row (indices `0..ncols`) touch only entries written when that row was built, so uninitialized memory is never read): + for ( i = 0; i < n - 1; i++ ) { + offset = i & 1; + so = i + 2; + mo = so - 1; + ncols = ( mo - offset ) / 2; + noffset = ( offset ) ? 0 : 1; + mcols = ( mo + 1 - noffset ) / 2; + + // Zero-fill the next row: + for ( j = 0; j <= mcols; j++ ) { + nxt[ j ] = 0.0; + } + for ( j = 0; j <= ncols; j++ ) { + co = ( 2 * j ) + offset; + k = ( co + 1 ) / 2; + nxt[ k ] += ( (double)( co - so ) * cur[ j ] ) / (double)( so - 1 ); + if ( co != 0 ) { + k = ( co - 1 ) / 2; + nxt[ k ] += ( -(double)co * cur[ j ] ) / (double)( so - 1 ); + } + } + // Swap the row buffers: + tmp = cur; + cur = nxt; + nxt = tmp; + } + sum = evalpoly( cur, len, c * c ); + free( wsp ); + if ( ( n - 1 ) & 1 ) { + sum *= c; // the first coefficient is order 1, so the row is really an odd polynomial + } + if ( sum == 0.0 ) { + return sum; // NOTE: may be either positive or negative zero + } + // The remaining terms are computed using logs since the powers and factorials get real large real quick: + powTerms = (double)n * STDLIB_CONSTANT_FLOAT64_LN_PI; + if ( s == 0.0 ) { + return ( sum >= 0.0 ) ? STDLIB_CONSTANT_FLOAT64_PINF : STDLIB_CONSTANT_FLOAT64_NINF; + } + powTerms -= stdlib_base_ln( stdlib_base_abs( s ) ) * ( (double)n + 1.0 ); + powTerms += stdlib_base_gammaln( (double)n ) + stdlib_base_ln( stdlib_base_abs( sum ) ); + + if ( powTerms > STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + return ( sum >= 0.0 ) ? STDLIB_CONSTANT_FLOAT64_PINF : STDLIB_CONSTANT_FLOAT64_NINF; + } + out = stdlib_base_exp( powTerms ) * stdlib_base_signum( sum ); + if ( s < 0.0 && ( ( n + 1 ) & 1 ) ) { + out *= -1.0; + } + return out; +} + +/** +* Evaluates the polygamma function for large values of `x` such as for `x > 400`. +* +* ## References +* +* - [Wolfram Research: Series representations][@wolfram:polygamma2:series] +* +* [@wolfram:polygamma2:series]: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/06/02/0001/ +* +* @param n derivative to evaluate +* @param x input value +* @return (n+1)'th derivative +*/ +static double atinfinityplus( const int32_t n, const double x ) { + double partTerm; // value of current term excluding the Bernoulli number part + double xsquared; + double term; // value of current term to be added to sum + double sum; // current value of accumulated sum + double nlx; + double k2; + int32_t k; + + if ( (double)n + x == x ) { + // If `x` is very large, just concentrate on the first part of the expression and use logs: + if ( n == 1 ) { + // NOTE: while unreachable via the main dispatch function (which delegates `n == 1` to `stdlib_base_trigamma`), the following branch is retained for fidelity with the original (Boost) implementation: + return 1.0 / x; + } + nlx = (double)n * stdlib_base_ln( x ); + if ( nlx < STDLIB_CONSTANT_FLOAT64_MAX_LN && n < MAX_FACTORIAL ) { + return ( ( n & 1 ) ? 1.0 : -1.0 ) * stdlib_base_factorial( (double)n - 1.0 ) * stdlib_base_pow( x, -(double)n ); + } + return ( ( n & 1 ) ? 1.0 : -1.0 ) * stdlib_base_exp( stdlib_base_gammaln( (double)n ) - ( (double)n * stdlib_base_ln( x ) ) ); + } + xsquared = x * x; + + // Start by setting `partTerm` to `(n-1)! / x^(n+1)`, which is common to both the first term of the series (with k = 1) and to the leading part. We can then get to the leading term by: `partTerm * (n + 2 * x) / 2` and to the first term in the series (excluding the Bernoulli number) by: `partTerm n * (n + 1) / (2x)`. If either the factorial would over- or the power term underflow, set `partTerm` to 0 and then we know that we have to use logs for the initial terms: + if ( n > MAX_FACTORIAL && ( (double)n * (double)n ) > STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + partTerm = 0.0; + } else { + partTerm = stdlib_base_factorial( (double)n - 1.0 ) * stdlib_base_pow( x, -(double)n - 1.0 ); + } + if ( partTerm == 0.0 ) { + // Either `n` is very large, or the power term underflows. Set the initial values of `partTerm`, `term`, and `sum` via logs: + partTerm = stdlib_base_gammaln( (double)n ) - ( ( (double)n + 1.0 ) * stdlib_base_ln( x ) ); + sum = stdlib_base_exp( partTerm + stdlib_base_ln( (double)n + ( 2.0 * x ) ) - STDLIB_CONSTANT_FLOAT64_LN2 ); + partTerm += stdlib_base_ln( (double)n * ( (double)n + 1.0 ) ) - STDLIB_CONSTANT_FLOAT64_LN2 - stdlib_base_ln( x ); + partTerm = stdlib_base_exp( partTerm ); + } else { + sum = partTerm * ( (double)n + ( 2.0 * x ) ) / 2.0; + partTerm *= ( (double)n * ( (double)n + 1.0 ) ) / 2.0; + partTerm /= x; + } + // If the leading term is 0, so is the result: + if ( sum == 0.0 ) { + return sum; // NOTE: the early return intentionally precedes the final parity adjustment and may thus return positive zero + } + for ( k = 1; ; ) { + term = partTerm * stdlib_base_bernoulli( (double)( 2 * k ) ); + sum += term; + + // Normal termination condition: + if ( stdlib_base_abs( term / sum ) < STDLIB_CONSTANT_FLOAT64_EPS ) { + break; + } + // Increment our counter, and move `partTerm` on to the next value (note: the recurrence products must be computed in double-precision arithmetic, as the equivalent 32-bit integer products may overflow): + k += 1; + k2 = 2.0 * (double)k; + partTerm *= ( (double)n + k2 - 2.0 ) * ( (double)n - 1.0 + k2 ); + partTerm /= ( k2 - 1.0 ) * k2; + partTerm /= xsquared; + if ( k > MAX_SERIES_ITERATIONS ) { + return 0.0 / 0.0; // NaN + } + } + if ( ( n - 1 ) & 1 ) { + sum = -sum; + } + return sum; +} + +/** +* Evaluates the polygamma function for intermediate values of `x` by performing a forward recursion into the asymptotic (large `x`) regime. +* +* ## References +* +* - [Wolfram Research: Recurrence identities][@wolfram:polygamma2:recurrence] +* +* [@wolfram:polygamma2:recurrence]: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0017/ +* +* @param n derivative to evaluate +* @param x input value +* @return (n+1)'th derivative +*/ +static double attransitionplus( const int32_t n, const double x ) { + double minusMminus1; + double lnterm; + double zpows; + double iter; + double sum0; + double d4d; + double N; + double z; + int32_t k; + + // Use N = (0.4 * digits) + (4 * n) for target value for x: + d4d = 0.4 * DIGITS_BASE10; + N = d4d + ( 4.0 * (double)n ); + + // NOTE: `iter` is a double having a fractional part and may be NaN (when `x` is NaN); the loops below run exactly `floor(iter)` times (zero times for NaN), matching JavaScript semantics: + iter = N - stdlib_base_trunc( x ); + + if ( iter > (double)MAX_SERIES_ITERATIONS ) { + return 0.0 / 0.0; // NaN + } + minusMminus1 = -(double)n - 1.0; + z = x; + sum0 = 0.0; + + // cppcheck-suppress unreadVariable + zpows = 0.0; // NOTE: dead initialization, kept to mirror the JavaScript implementation + + // Forward recursion to larger `x`, need to check for overflow first though: + if ( stdlib_base_ln( z + iter ) * minusMminus1 > -STDLIB_CONSTANT_FLOAT64_MAX_LN ) { + for ( k = 1; (double)k <= iter; k++ ) { + zpows = stdlib_base_pow( z, minusMminus1 ); + sum0 += zpows; + z += 1.0; + } + sum0 *= stdlib_base_factorial( (double)n ); + } else { + for ( k = 1; (double)k <= iter; k++ ) { + lnterm = ( stdlib_base_ln( z ) * minusMminus1 ) + stdlib_base_gammaln( (double)n + 1.0 ); + sum0 += stdlib_base_exp( lnterm ); + z += 1.0; + } + } + if ( ( n - 1 ) & 1 ) { + sum0 = -sum0; + } + return sum0 + atinfinityplus( n, z ); +} + +/** +* Evaluates the polygamma function near zero. +* +* ## Notes +* +* - If we take this [expansion][1] for `polygamma` and substitute in this [expression][2] for `polygamma(n, 1)`, we get an alternating series for polygamma when `x` is small in terms of zeta functions of integer arguments (which are easy to evaluate, at least when the integer is even). +* +* [1]: http://functions.wolfram.com/06.15.06.0003.02 +* [2]: http://functions.wolfram.com/06.15.03.0009.01 +* +* @param n derivative to evaluate +* @param x input value +* @return (n+1)'th derivative +*/ +static double nearzero( const int32_t n, const double x ) { + double factorialPart; + double prefix; + double scale; + double term; + double sum; + int32_t k; + + // In order to avoid spurious overflow, save the `n!` term for later, and rescale at the end: + scale = stdlib_base_factorial( (double)n ); + + // "factorialPart" contains everything except the zeta function evaluations in each term: + factorialPart = 1.0; + + // "prefix" is what we'll be adding the accumulated sum to, it will be `n! / z^(n+1)`, but since we're scaling by `n!` it is just `1 / z^(n+1)` for now: + prefix = stdlib_base_pow( x, (double)n + 1.0 ); + if ( prefix == 0.0 ) { + // NOTE: positive infinity is returned for all `n`, irrespective of parity (e.g., when `x` is `±0` or the power underflows), matching the JavaScript implementation: + return STDLIB_CONSTANT_FLOAT64_PINF; + } + prefix = 1.0 / prefix; + + // First term in the series is necessarily `< zeta(2) < 2`, so ignore the sum if it will have no effect on the result: + if ( prefix > 2.0 / STDLIB_CONSTANT_FLOAT64_EPS ) { + // NOTE: the JavaScript implementation wraps the following returns in overflow guards of the form `AX/prefix < scale ? ±infinity : ...`; however, as `AX` is an undefined (never assigned) variable there, the comparison is always false and the guards never fire. We reproduce the effective behavior (unguarded products), as natural IEEE 754 overflow to infinity coincides with the intended outcome: + if ( n & 1 ) { + return prefix * scale; + } + return -( prefix * scale ); + } + sum = prefix; + for ( k = 0; ; ) { + // Get the k'th term (note: the argument sums must be computed in double-precision arithmetic, as the equivalent 32-bit integer sums may overflow): + term = factorialPart * stdlib_base_zeta( (double)k + (double)n + 1.0 ); + sum += term; + + // Termination condition: + if ( stdlib_base_abs( term ) < stdlib_base_abs( sum * STDLIB_CONSTANT_FLOAT64_EPS ) ) { + break; + } + // Move on `k` and `factorialPart`: + k += 1; + factorialPart *= ( -x * ( (double)n + (double)k ) ) / (double)k; + + // Last chance exit: + if ( k > MAX_SERIES_ITERATIONS ) { + return 0.0 / 0.0; // NaN + } + } + // We need to multiply by the scale, at each stage checking for overflow: + if ( STDLIB_CONSTANT_FLOAT64_MAX / scale < sum ) { + // NOTE: positive infinity is returned irrespective of parity (and also when `scale` is infinite, as `MAX` divided by infinity is `+0`), matching the JavaScript implementation: + return STDLIB_CONSTANT_FLOAT64_PINF; + } + sum *= scale; + return ( n & 1 ) ? sum : -sum; +} + +/** +* Evaluates the polygamma function. +* +* @param n order of derivative +* @param x input value +* @return (n+1)'th derivative +* +* @example +* double out = stdlib_base_polygamma( 3, 1.2 ); +* // returns ~3.245 +* +* @example +* double out = stdlib_base_polygamma( 5, 1.2 ); +* // returns ~41.39 +* +* @example +* double out = stdlib_base_polygamma( 3, -4.9 ); +* // returns ~60014.239 +* +* @example +* double out = stdlib_base_polygamma( -1, 5.3 ); +* // returns NaN +* +* @example +* double out = stdlib_base_polygamma( 2, -2.0 ); +* // returns NaN +* +* @example +* double out = stdlib_base_polygamma( 1, 0.0 / 0.0 ); +* // returns NaN +*/ +double stdlib_base_polygamma( const int32_t n, const double x ) { + double result; + double limit; + double z; + double v; + + if ( n < 0 ) { + // NOTE: the JavaScript implementation returns NaN for any `n` which is not a nonnegative integer; non-integer and NaN `n` are unrepresentable here, as, at the C/N-API boundary, `n` is an `int32_t` (the add-on truncates non-integer values and maps NaN to zero): + return 0.0 / 0.0; // NaN + } + if ( n == 0 ) { + return stdlib_base_digamma( x ); + } + if ( n == 1 ) { + return stdlib_base_trigamma( x ); + } + if ( x < 0.0 ) { + if ( stdlib_base_floor( x ) == x ) { + // Result is infinity if `x` is an odd negative integer, and NaN (pole error) if `x` is an even negative integer. Parity is determined by halving the truncated value and comparing against its floor, which is equivalent to the JavaScript implementation's `trunc( x ) & 1` (ToInt32 wrapping preserves parity for every integer-valued double) for all inputs, including `x` equal to negative infinity, which is considered "even": + v = stdlib_base_trunc( x ) / 2.0; + if ( v != stdlib_base_floor( v ) ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + return 0.0 / 0.0; // NaN + } + z = 1.0 - x; + result = stdlib_base_polygamma( n, z ) + ( STDLIB_CONSTANT_FLOAT64_PI * polycotpi( n, z, x ) ); + return ( n & 1 ) ? -result : result; + } + // Limit for use of small-x series is chosen so that the series doesn't go too divergent in the first few terms. Ordinarily, this would mean setting the limit to `~1/n`, but we can tolerate a small amount of divergence: + limit = stdlib_base_min( 5.0 / (double)n, 0.25 ); + if ( x < limit ) { + return nearzero( n, x ); + } + if ( x > ( 0.4 * DIGITS_BASE10 ) + ( 4.0 * (double)n ) ) { + return atinfinityplus( n, x ); + } + if ( x == 1.0 ) { + return ( ( n & 1 ) ? 1.0 : -1.0 ) * stdlib_base_factorial( (double)n ) * stdlib_base_zeta( (double)n + 1.0 ); + } + if ( x == 0.5 ) { + result = ( ( n & 1 ) ? 1.0 : -1.0 ) * stdlib_base_factorial( (double)n ) * stdlib_base_zeta( (double)n + 1.0 ); + if ( stdlib_base_abs( result ) >= stdlib_base_ldexp( STDLIB_CONSTANT_FLOAT64_MAX, -n - 1 ) ) { + return ( stdlib_base_signum( result ) == 1.0 ) ? STDLIB_CONSTANT_FLOAT64_PINF : STDLIB_CONSTANT_FLOAT64_NINF; + } + // NOTE: the following line only executes when the previous guard failed, implying a finite factorial and, in turn, `n <= 170`, so `n+1` cannot overflow: + result *= stdlib_base_ldexp( 1.0, n + 1 ) - 1.0; + return result; + } + return attransitionplus( n, x ); +} diff --git a/lib/node_modules/@stdlib/math/base/special/polygamma/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/polygamma/test/test.native.js new file mode 100644 index 000000000000..c3cef78f0bd7 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/polygamma/test/test.native.js @@ -0,0 +1,267 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +'use strict'; + +// MODULES // + +var resolve = require( 'path' ).resolve; +var tape = require( 'tape' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var factorial = require( '@stdlib/math/base/special/factorial' ); +var incrspace = require( '@stdlib/array/base/incrspace' ); +var ldexp = require( '@stdlib/math/base/special/ldexp' ); +var zeta = require( '@stdlib/math/base/special/riemann-zeta' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); +var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' ); +var tryRequire = require( '@stdlib/utils/try-require' ); + + +// VARIABLES // + +var polygamma = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( polygamma instanceof Error ) +}; + + +// FIXTURES // + +var largePositive = require( './fixtures/cpp/large_positive.json' ); +var largeNegative = require( './fixtures/cpp/large_negative.json' ); +var smallPositive = require( './fixtures/cpp/small_positive.json' ); +var hugePositive = require( './fixtures/cpp/huge_positive.json' ); +var positive = require( './fixtures/cpp/positive.json' ); +var negative = require( './fixtures/cpp/negative.json' ); + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof polygamma, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'if provided `NaN` for `x`, the function returns `NaN`', opts, function test( t ) { + var val = polygamma( 1, NaN ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = polygamma( 2, NaN ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided a negative `n`, the function returns `NaN`', opts, function test( t ) { + var val = polygamma( -1, 5.3 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = polygamma( -3, 2.5 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided an even negative integer for `x`, the function returns `NaN`', opts, function test( t ) { + var val = polygamma( 2, -2.0 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = polygamma( 2, -4.0 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'if provided an odd negative integer for `x`, the function returns positive infinity', opts, function test( t ) { + var val = polygamma( 4, -1.0 ); + t.strictEqual( val, PINF, 'returns expected value' ); + + val = polygamma( 4, -3.0 ); + t.strictEqual( val, PINF, 'returns expected value' ); + t.end(); +}); + +tape( 'the function evaluates the polygamma function for small positive numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = smallPositive.n; + x = smallPositive.x; + expected = smallPositive.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function evaluates the polygamma function for positive numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = positive.n; + x = positive.x; + expected = positive.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 16 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function evaluates the polygamma function for negative numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = negative.n; + x = negative.x; + expected = negative.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 2900 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function evaluates the polygamma function for large positive numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = largePositive.n; + x = largePositive.x; + expected = largePositive.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 20 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function evaluates the polygamma function for large negative numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = largeNegative.n; + x = largeNegative.x; + expected = largeNegative.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 2800 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'the function evaluates the polygamma function for huge positive numbers', opts, function test( t ) { + var expected; + var n; + var x; + var y; + var i; + + n = hugePositive.n; + x = hugePositive.x; + expected = hugePositive.y; + for ( i = 0; i < x.length; i++ ) { + y = polygamma( n[ i ], x[ i ] ); + t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'for positive odd `n` and `x` equal to `1.0`, the function computes the factorial of `n` multiplied by the Riemann zeta function evaluated at `n+1`', opts, function test( t ) { + var expected; + var actual; + var values; + var i; + var n; + + values = incrspace( 1, 99, 2 ); + for ( i = 0; i < values.length; i++ ) { + n = values[ i ]; + expected = factorial( n ) * zeta( n+1 ); + actual = polygamma( n, 1.0 ); + t.strictEqual( isAlmostSameValue( actual, expected, 0 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'for positive even `n` and `x` equal to `1.0`, the function computes the negative factorial of `n` multiplied by the Riemann zeta function evaluated at `n+1`', opts, function test( t ) { + var expected; + var actual; + var values; + var i; + var n; + + values = incrspace( 2, 100, 2 ); + for ( i = 0; i < values.length; i++ ) { + n = values[ i ]; + expected = (-1.0) * factorial( n ) * zeta( n+1 ); + actual = polygamma( n, 1.0 ); + t.strictEqual( isAlmostSameValue( actual, expected, 0 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'for positive odd `n` and `x` equal to `0.5`, the function computes the factorial of `n` multiplied by the Riemann zeta function evaluated at `n+1`', opts, function test( t ) { + var expected; + var actual; + var values; + var i; + var n; + + values = incrspace( 1, 199, 2 ); + for ( i = 0; i < values.length; i++ ) { + n = values[ i ]; + expected = factorial( n ) * zeta( n+1 ); + expected *= ldexp( 1.0, n + 1 ) - 1.0; + actual = polygamma( n, 0.5 ); + t.strictEqual( isAlmostSameValue( actual, expected, 0 ), true, 'returns expected value' ); + } + t.end(); +}); + +tape( 'for positive even `n` and `x` equal to `0.5`, the function computes the negative factorial of `n` multiplied by the Riemann zeta function evaluated at `n+1`', opts, function test( t ) { + var expected; + var actual; + var values; + var i; + var n; + + values = incrspace( 2, 200, 2 ); + for ( i = 0; i < values.length; i++ ) { + n = values[ i ]; + expected = (-1.0) * factorial( n ) * zeta( n+1 ); + expected *= ldexp( 1.0, n + 1 ) - 1.0; + actual = polygamma( n, 0.5 ); + t.strictEqual( isAlmostSameValue( actual, expected, 0 ), true, 'returns expected value' ); + } + t.end(); +});