From 3d0b53b97c98fb7e3f3b5b6add68da48d6316369 Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Sat, 7 Dec 2024 21:50:36 +0530 Subject: [PATCH 1/7] feat: add C implementation for math/base/special/falling-factorial --- .../falling-factorial/benchmark/benchmark.js | 11 +- .../benchmark/benchmark.native.js | 55 +++++++ .../benchmark/c/native/Makefile | 146 ++++++++++++++++++ .../benchmark/c/native/benchmark.c | 138 +++++++++++++++++ 4 files changed, 345 insertions(+), 5 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.js index f50f7840a0bf..4ceb8fe532a5 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.js @@ -21,8 +21,8 @@ // MODULES // var bench = require( '@stdlib/bench' ); -var randu = require( '@stdlib/random/base/randu' ); -var ceil = require( '@stdlib/math/base/special/ceil' ); +var randu = require( '@stdlib/random/array/uniform' ); +var discreteUniform = require( '@stdlib/random/array/discrete-uniform' ); var isnan = require( '@stdlib/math/base/assert/is-nan' ); var pkg = require( './../package.json' ).name; var fallingFactorial = require( './../lib' ); @@ -36,11 +36,12 @@ bench( pkg, function benchmark( b ) { var y; var i; + x = randu( 100, 1.0, 41.0 ); + n = discreteUniform( 100, 0, 40 ); + b.tic(); for ( i = 0; i < b.iterations; i++ ) { - x = ( randu()*40.0 ) + 1.0; - n = ceil( randu()*40.0 ); - y = fallingFactorial( x, n ); + y = fallingFactorial( x[ i % x.length ], n[ i % n.length ] ); if ( isnan( y ) ) { b.fail( 'should not return NaN' ); } diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js new file mode 100644 index 000000000000..4ceb8fe532a5 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js @@ -0,0 +1,55 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2018 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +'use strict'; + +// MODULES // + +var bench = require( '@stdlib/bench' ); +var randu = require( '@stdlib/random/array/uniform' ); +var discreteUniform = require( '@stdlib/random/array/discrete-uniform' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var pkg = require( './../package.json' ).name; +var fallingFactorial = require( './../lib' ); + + +// MAIN // + +bench( pkg, function benchmark( b ) { + var x; + var n; + var y; + var i; + + x = randu( 100, 1.0, 41.0 ); + n = discreteUniform( 100, 0, 40 ); + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + y = fallingFactorial( x[ i % x.length ], n[ i % n.length ] ); + 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/falling-factorial/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile new file mode 100644 index 000000000000..f69e9da2b4d3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2024 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/falling-factorial/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c new file mode 100644 index 000000000000..d3ce002e4264 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c @@ -0,0 +1,138 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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/falling_factorial.h" +#include +#include +#include +#include +#include + +#define NAME "falling-factorial" +#define ITERATIONS 100 +#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 [0,1). +* +* @return random number +*/ +static double rand_double( void ) { + int r = rand(); + return (double)r / ( (double)RAND_MAX + 1.0 ); +} + +/** +* Runs a benchmark. +* +* @return elapsed time in seconds +*/ +static double benchmark( void ) { + double elapsed; + int32_t n[ 100 ]; + double x[ 100 ]; + double t; + double y; + int i; + + for ( i = 0; i < 100; i++ ) { + x[ i ] = ( 40.0 * rand_double() ) + 1.0; + n[ i ] = (int32_t)( 41.0 * rand_double() ); + } + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + y = stdlib_base_falling_factorial( x[ i % 100 ], n[ i % 100 ] ); + if ( y != y ) { + printf( "should return a nonnegative integer\n" ); + break; + } + } + elapsed = tic() - t; + if ( y != y ) { + printf( "should return a nonnegative integer\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 ); +} From 28af1a95aae4dc2746723344935586a385404341 Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Sat, 7 Dec 2024 23:01:08 +0530 Subject: [PATCH 2/7] feat: add C implementation for math/base/special/falling-factorial --- .../benchmark/benchmark.native.js | 15 +++++++-- .../falling-factorial/examples/c/Makefile | 0 .../falling-factorial/examples/c/example.c | 31 +++++++++++++++++++ 3 files changed, 43 insertions(+), 3 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js index 4ceb8fe532a5..eed82c096450 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2018 The Stdlib Authors. +* Copyright (c) 2024 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. @@ -20,17 +20,26 @@ // MODULES // +var resolve = require( 'path' ).resolve; var bench = require( '@stdlib/bench' ); var randu = require( '@stdlib/random/array/uniform' ); var discreteUniform = require( '@stdlib/random/array/discrete-uniform' ); var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var tryRequire = require( '@stdlib/utils/try-require' ); var pkg = require( './../package.json' ).name; -var fallingFactorial = require( './../lib' ); + + +// VARIABLES // + +var fallingFactorial = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( fallingFactorial instanceof Error ) +}; // MAIN // -bench( pkg, function benchmark( b ) { +bench( pkg+'::native', opts, function benchmark( b ) { var x; var n; var y; diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c new file mode 100644 index 000000000000..3d12d5086ada --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c @@ -0,0 +1,31 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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/falling_factorial.h" +#include + +int main( void ) { + const double x[] = { -10.0, -7.78, -5.56, -3.33, -1.11, 1.11, 3.33, 5.56, 7.78, 10.0 }; + + double v; + int i; + for ( i = 0; i < 10; i++ ) { + v = stdlib_base_falling_factorial( x[ i ] ); + printf( "x: %lf, falling_factorial(x): %lf\n", x[ i ], v ); + } +} From 7c04bd708097200cf19fda42f186c810b7d9888b Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Sun, 8 Dec 2024 11:31:35 +0530 Subject: [PATCH 3/7] feat: add C implementation for math/base/special/falling-factorial --- .../falling-factorial/examples/c/Makefile | 146 ++++ .../falling-factorial/examples/c/example.c | 6 +- .../falling-factorial/examples/index.js | 13 +- .../math/base/special/falling_factorial.h | 40 + .../special/falling-factorial/lib/native.js | 59 ++ .../special/falling-factorial/src/Makefile | 70 ++ .../special/falling-factorial/src/addon.c | 22 + .../base/special/falling-factorial/src/main.c | 730 ++++++++++++++++++ 8 files changed, 1078 insertions(+), 8 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile index e69de29bb2d1..6aed70daf167 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2024 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/falling-factorial/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c index 3d12d5086ada..800863ab6941 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c @@ -18,14 +18,16 @@ #include "stdlib/math/base/special/falling_factorial.h" #include +#include int main( void ) { const double x[] = { -10.0, -7.78, -5.56, -3.33, -1.11, 1.11, 3.33, 5.56, 7.78, 10.0 }; + const int32_t n[] = { 5, 4, 3, 2, 1, 0, -1, -2, -3, -4 }; double v; int i; for ( i = 0; i < 10; i++ ) { - v = stdlib_base_falling_factorial( x[ i ] ); - printf( "x: %lf, falling_factorial(x): %lf\n", x[ i ], v ); + v = stdlib_base_falling_factorial( x[ i ], n[ i ] ); + printf( "x: %lf, n: %d, falling_factorial(x,n): %lf\n", x[ i ], n[ i ], v ); } } diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js index caf9ada7144c..828cb7aeff64 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2018 The Stdlib Authors. +* Copyright (c) 2024 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. @@ -18,16 +18,17 @@ 'use strict'; -var randu = require( '@stdlib/random/base/randu' ); -var ceil = require( '@stdlib/math/base/special/ceil' ); +var randu = require( '@stdlib/random/array/uniform' ); +var discreteUniform = require( '@stdlib/random/array/discrete-uniform' ); var fallingFactorial = require( './../lib' ); var n; var x; var i; +x = randu( 100, -20.0, 20.0 ); +n = discreteUniform( 100, 0, 20 ); + for ( i = 0; i < 100; i++ ) { - x = ( randu()*40.0 ) - 20.0; - n = ceil( randu()*20.0 ); - console.log( 'fallingFactorial(%d,%d) = %d', x, n, fallingFactorial( x, n ) ); + console.log( 'fallingFactorial(%d,%d) = %d', x[ i ], n[ i ], fallingFactorial( x[ i ], n[ i ] ) ); } diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h new file mode 100644 index 000000000000..14d77f9f0fbf --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h @@ -0,0 +1,40 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +#ifndef STDLIB_MATH_BASE_SPECIAL_FALLING_FACTORIAL_H +#define STDLIB_MATH_BASE_SPECIAL_FALLING_FACTORIAL_H + +#include + +/* +* If C++, prevent name mangling so that the compiler emits a binary file having undecorated names, thus mirroring the behavior of a C compiler. +*/ +#ifdef __cplusplus +extern "C" { +#endif + +/** +* Computes the falling factorial of `x` and `n`. +*/ +double stdlib_base_falling_factorial( const double x, const int32_t n ); + +#ifdef __cplusplus +} +#endif + +#endif // !STDLIB_MATH_BASE_SPECIAL_FALLING_FACTORIAL_H diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js new file mode 100644 index 000000000000..3a7b0100d2cc --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js @@ -0,0 +1,59 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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 // + +/** +* Computes the falling factorial of `x` and `n`. +* +* @private +* @param {number} x - first function parameter +* @param {NonNegativeInteger} n - second function parameter +* @returns {number} function value +* +* @example +* var v = fallingFactorial( 0.9, 5 ); +* // returns ~0.644 +* +* @example +* var v = fallingFactorial( -9.0, 3 ); +* // returns -990.0 +* +* @example +* var v = fallingFactorial( 0.0, 2 ); +* // returns 0.0 +* +* @example +* var v = fallingFactorial( 3.0, -2 ); +* // returns NaN +*/ +function fallingFactorial( x, n ) { + return addon( x, n ); +} + + +// EXPORTS // + +module.exports = fallingFactorial; diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile new file mode 100644 index 000000000000..bcf18aa46655 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile @@ -0,0 +1,70 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2024 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/falling-factorial/src/addon.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c new file mode 100644 index 000000000000..e8beff97c10d --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c @@ -0,0 +1,22 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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/falling_factorial.h" +#include "stdlib/math/base/napi/binary.h" + +STDLIB_MATH_BASE_NAPI_MODULE_DI_D( stdlib_base_falling_factorial ) diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c new file mode 100644 index 000000000000..7aa85e46e096 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c @@ -0,0 +1,730 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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 following copyright and license were part of the original implementation available as part of [Boost Library]{@link https://www.boost.org/doc/libs/1_86_0/boost/math/special_functions/factorials.hpp}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. +* +* Developed at SunPro, a Sun Microsystems, Inc. business. +* Permission to use, copy, modify, and distribute this +* software is freely granted, provided that this notice +* is preserved. +* ``` +*/ + +#include "stdlib/math/base/special/pow.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/constants/float64/pinf.h" +#include "stdlib/constants/float64/high_word_abs_mask.h" +#include "stdlib/number/float64/base/get_high_word.h" +#include "stdlib/math/base/assert/is_odd.h" +#include "stdlib/math/base/special/copysign.h" +#include "stdlib/constants/float64/ninf.h" +#include "stdlib/number/float64/base/set_high_word.h" +#include "stdlib/number/float64/base/set_low_word.h" +#include "stdlib/math/base/special/ldexp.h" +#include "stdlib/constants/float64/ln_two.h" +#include "stdlib/constants/float64/exponent_bias.h" +#include "stdlib/constants/float64/high_word_significand_mask.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/assert/is_infinite.h" +#include "stdlib/math/base/assert/is_integer.h" +#include "stdlib/math/base/special/sqrt.h" +#include "stdlib/number/float64/base/to_words.h" +#include "stdlib/constants/float64/num_high_word_significand_bits.h" + +// 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1 +static const int32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; + +static const double HUGE_VALUE = 1.0e300; + +static const double TINY_VALUE = 1.0e-300; + +// 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 +static const int32_t HIGH_MIN_NORMAL_EXP = 0x00100000; + +// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 +static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; + +// High: LN2 +static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 + +// Low: LN2 +static const double LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39 + +// 1/LN2 +static const double INV_LN2 = 1.44269504088896338700e+00; // 0x3FF71547, 0x652B82FE + +// High (24 bits): 1/LN2 +static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60000000 + +// Low: 1/LN2 +static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 + +// 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 +static const int32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; + +// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 +static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; + +// 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 +static const int32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; + +// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 +static const int32_t HIGH_SIGNIFICAND_HALF = 0x00080000; + +static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 + +// 2/(3*LN2) +static const double CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD + +// (float)CP +static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 + +// Low: CP_HI +static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 + +static const double BP [2] = { 1.0, 1.5 }; + +static const double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 + +static const double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 + +// 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 +static const int32_t HIGH_BIASED_EXP_31 = 0x41e00000; + +// 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 +static const int32_t HIGH_BIASED_EXP_64 = 0x43f00000; + +// 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 +static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; + +// 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 +static const int32_t HIGH_1075 = 0x4090cc00; + +// 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 +static const uint32_t HIGH_NEG_1075 = 0xc090cc00; + +static const int32_t HIGH_NUM_NONSIGN_BITS = 31; + +// -(1024-log2(ovfl+.5ulp)) +static const double OVT = 8.0085662595372944372e-17; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_l + +/** +* 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_l( const double x ) { + return 0.5999999999999946 + (x * (0.4285714285785502 + (x * (0.33333332981837743 + (x * (0.272728123808534 + (x * (0.23066074577556175 + (x * 0.20697501780033842))))))))); +} + +// END: polyval_l + +// BEGIN: polyval_w + +/** +* 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_w( const double x ) { + return 0.5 + (x * (-0.3333333333333333 + (x * 0.25))); +} + +// END: polyval_w + +// BEGIN: polyval_p + +/** +* 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_p( const double x ) { + return 0.16666666666666602 + (x * (-0.0027777777777015593 + (x * (0.00006613756321437934 + (x * (-0.0000016533902205465252 + (x * 4.1381367970572385e-8))))))); +} + +// END: polyval_p + +/* End auto-generated functions. */ + +/** +* Evaluates the exponential function when \\( y = \pm \infty\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_infinite( -1.0, 1.0/0.0 ); +* // returns NaN +* +* @example +* double out = y_is_infinite( 1.5, -1.0/0.0 ); +* // returns 0.0 +*/ +static double y_is_infinite( const double x, const double y ) { + if ( x == -1.0 ) { + return 0.0 / 0.0; // NaN + } + if ( x == 1.0 ) { + return 1.0; + } + // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) + if ( ( stdlib_base_abs( x ) < 1.0 ) == ( y == STDLIB_CONSTANT_FLOAT64_PINF ) ) { + return 0.0; + } + // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Evaluates the exponential function when \\(|y| > 2^64\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_huge( 9.0, 3.6893488147419103e19 ); +* // returns Infinity +* +* @example +* double out = y_is_huge( -3.14, -3.6893488147419103e19 ); +* // returns 0.0 +*/ +static double y_is_huge( const double x, const double y ) { + uint32_t ahx; + uint32_t hx; + + stdlib_base_float64_get_high_word( x, &hx ); + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + if ( ahx <= HIGH_MAX_NEAR_UNITY ) { + if ( y < 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; + } + // `x` has a biased exponent greater than or equal to `0`... + + if ( y > 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; +} + +/** +* Evaluates the exponential function when \\(|x| = 0\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = x_is_zero( 0.0, 2 ); +* // returns 0.0 +* +* @example +* double out = x_is_zero( 0.0, -9 ); +* // returns Infinity +*/ +static double x_is_zero( const double x, const double y ) { + if ( y == STDLIB_CONSTANT_FLOAT64_NINF ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( y == STDLIB_CONSTANT_FLOAT64_PINF ) { + return 0.0; + } + if ( y > 0.0 ) { + if ( stdlib_base_is_odd( y ) ) { + return x; // handles +-0 + } + return 0.0; + } + // y < 0.0 + if ( stdlib_base_is_odd( y ) ) { + return stdlib_base_copysign( STDLIB_CONSTANT_FLOAT64_PINF, x ); // handles +-0 + } + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Computes \\(2^{\mathrm{hp} + \mathrm{lp}\\). +* +* @param j high word of `hp + lp` +* @param hp first power summand +* @param lp second power summand +* @return function value +* +* @example +* double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); +* // returns ~0.79 +*/ +static double pow2( uint32_t j, const double hp, const double lp ) { + uint32_t tmp; + double hpc; + int32_t jc; + int32_t nc; + int32_t i; + int32_t k; + int32_t n; + double t1; + double t; + double r; + double u; + double v; + double w; + double z; + + hpc = hp; + jc = (int32_t)j; + i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + k = ( ( i >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); + n = 0; + nc = (int32_t)n; + + // `|z| > 0.5`, set `n = z+0.5` + if ( i > HIGH_BIASED_EXP_NEG_1 ) { + n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ); + k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n + tmp = ( ( n & ~( HIGH_SIGNIFICAND_MASK >> k ) ) ); + t = 0.0; + stdlib_base_float64_set_high_word( tmp, &t ); + n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS - k ) ); + nc = (int32_t)n; + if ( jc < 0 ) { + nc = -nc; + } + hpc -= t; + } + t = lp + hpc; + stdlib_base_float64_set_low_word( 0, &t ); + u = t * LN2_HI; + v = ( ( lp - ( t - hpc ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); + z = u + v; + w = v - ( z - u ); + t = z * z; + t1 = z - ( t * polyval_p( t ) ); + r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); + z = 1.0 - ( r - z ); + stdlib_base_float64_get_high_word( z, &j ); + j = j + ( nc << STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ); + jc = (int32_t)j; + + // Check for subnormal output... + if ( ( jc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) <= 0 ) { + z = stdlib_base_ldexp( z, nc ); + } else { + stdlib_base_float64_set_high_word( j, &z ); + } + return z; +} + +/** +* Computes \\(\operatorname{log}(x)\\) assuming \\(|1-x|\\) is small and using the approximation \\(x - x^2/2 + x^3/3 - x^4/4\\). +* +* @param ax absolute value of x +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* logx( 9.0, &o1, &o2 ); +*/ +void logx( const double ax, double *o1, double *o2 ) { + double t2; + double t1; + double t; + double w; + double u; + double v; + + t = ax - 1.0; // `t` has `20` trailing zeros + w = t * t * polyval_w( t ); + u = INV_LN2_HI * t; // `INV_LN2_HI` has `21` significant bits + v = ( t * INV_LN2_LO ) - ( w * INV_LN2 ); + t1 = u + v; + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = v - (t1 - u); + + *o1 = t1; + *o2 = t2; +} + +/** +* Computes \\(\operatorname{log2}(ax)\\). +* +* @param ax absolute value of x +* @param ahx high word of ax +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* log2ax( 9.0, 1075970048, &o1, &o2 ); +*/ +void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { + uint32_t ahxcc; + uint32_t tmp; + int32_t ahxc; + double axc; + int32_t n; + int32_t j; + int32_t k; + double ss; // hs + ls + double s2; // ss squared + double hs; + double ls; + double ht; + double lt; + double bp; // BP constant + double dp; // DP constant + double hp; + double lp; + double hz; + double lz; + double t1; + double t2; + double t; + double r; + double u; + double v; + + n = 0; + axc = ax; + ahxc = ahx; + ahxcc = (uint32_t)ahx; + // Check if `x` is subnormal... + if ( ahxc < HIGH_MIN_NORMAL_EXP ) { + axc *= TWO53; + n -= 53; + stdlib_base_float64_get_high_word( axc, &ahxcc ); + ahxc = (int32_t)ahxcc; + } + // Extract the unbiased exponent of `x`: + n += ( ( ahxc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); + + // Isolate the significand bits of `x`: + j = ( ahxc & HIGH_SIGNIFICAND_MASK ); + + // Normalize `ahx` by setting the (biased) exponent to `1023`: + ahxc = ( j|HIGH_BIASED_EXP_0 ); + + // Determine the interval of `|x|` by comparing significand bits... + + // |x| < sqrt(3/2) + if ( j <= 0x3988E ) { // 0 00000000000 00111001100010001110 + k = 0; + } + // |x| < sqrt(3) + else if ( j < 0xBB67A ) { // 0 00000000000 10111011011001111010 + k = 1; + } + // |x| >= sqrt(3) + else { + k = 0; + n += 1; + ahxc -= HIGH_MIN_NORMAL_EXP; + } + ahxcc = (uint32_t)ahxc; + + // Load the normalized high word into `|x|`: + stdlib_base_float64_set_high_word( ahxcc, &axc ); + + // Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: + bp = BP[ k ]; // BP[0] = 1.0, BP[1] = 1.5 + u = axc - bp; // (x-1) || (x-1.5) + v = 1.0 / (axc + bp); // 1/(x+1) || 1/(x+1.5) + ss = u * v; + hs = ss; + stdlib_base_float64_set_low_word( 0, &hs ); // set all low word (less significant significand) bits to 0s + + // Compute `ht = ax + bp` (via manipulation, i.e., bit flipping, of the high word): + tmp = ( ( ahxc >> 1 ) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; + tmp += ( k << 18 ); // `(k<<18)` can be considered the word equivalent of `1.0` or `1.5` + ht = 0.0; + stdlib_base_float64_set_high_word( tmp, &ht ); + lt = axc - ( ht - bp ); + ls = v * ( ( u - ( hs * ht ) ) - ( hs * lt ) ); + + // Compute `log(ax)`... + + s2 = ss * ss; + r = s2 * s2 * polyval_l( s2 ); + r += ls * (hs + ss); + s2 = hs * hs; + ht = 3.0 + s2 + r; + stdlib_base_float64_set_low_word( 0, &ht ); + lt = r - ( ( ht - 3.0 ) - s2 ); + + // u+v = ss*(1+...): + u = hs * ht; + v = ( ls * ht ) + ( lt * ss ); + + // 2/(3LN2) * (ss+...): + hp = u + v; + stdlib_base_float64_set_low_word( 0, &hp ); + lp = v - ( hp - u ); + hz = CP_HI * hp; // CP_HI+CP_LO = 2/(3*LN2) + lz = ( CP_LO * hp ) + ( lp*CP ) + DP_LO[ k ]; + + // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz + dp = DP_HI[ k ]; + t = n; + t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); + + *o1 = t1; + *o2 = t2; +} + +/** +* Evaluates the exponential function. +* +* @param x base +* @param y exponent +* @return function value +* +* @example +* double out = stdlib_base_pow( 2.0, 3.0 ); +* // returns 8.0 +* +* @example +* double out = stdlib_base_pow( 4.0, 0.5 ); +* // returns 2.0 +*/ +double stdlib_base_pow( const double x, const double y ) { + uint32_t hx; // high word `x` + uint32_t lx; // low word `x` + uint32_t hy; // high word `y` + uint32_t ly; // low word `y` + int32_t ahx; // absolute value high word `x` + int32_t ahy; // absolute value high word `y` + uint32_t i; + uint32_t j; + int32_t ic; + int32_t jc; + int32_t sx; // sign `x` + int32_t sy; // sign `y` + double ax; // absolute value `x` + double y1; + double hp; + double lp; + double t1; + double t2; + double xc; + double z; // y prime + + xc = x; + if ( stdlib_base_is_nan( xc ) || stdlib_base_is_nan( y ) ) { + return 0.0/0.0; // NaN + } + + // Split `y` into high and low words: + stdlib_base_float64_to_words( y, &hy, &ly ); + + // Special cases `y`... + if ( ly == 0 ) { + if ( y == 0.0 ) { + return 1.0; + } + if ( y == 1.0 ) { + return xc; + } + if ( y == -1.0 ) { + return 1.0 / xc; + } + if ( y == 0.5 ) { + return stdlib_base_sqrt( xc ); + } + if ( y == -0.5 ) { + return 1.0 / stdlib_base_sqrt( xc ); + } + if ( y == 2.0 ) { + return xc * xc; + } + if ( y == 3.0 ) { + return xc * xc * xc; + } + if ( y == 4.0 ) { + xc *= xc; + return xc * xc; + } + if ( stdlib_base_is_infinite( y ) ) { + return y_is_infinite( xc, y ); + } + } + // Split `x` into high and low words: + stdlib_base_float64_to_words( xc, &hx, &lx ); + + // Special cases `x`... + if ( lx == 0 ) { + if ( hx == 0 ) { + return x_is_zero( xc, y ); + } + if ( xc == 1.0 ) { + return 1.0; + } + if ( + xc == -1.0 && + stdlib_base_is_odd( y ) + ) { + return -1.0; + } + if ( stdlib_base_is_infinite( xc ) ) { + if ( xc == STDLIB_CONSTANT_FLOAT64_NINF ) { + // `stdlib_base_pow( 1/x, -y )` + return stdlib_base_pow( -0.0, -y ); + } + if ( y < 0.0 ) { + return 0.0; + } + return STDLIB_CONSTANT_FLOAT64_PINF; + } + } + if ( + xc < 0.0 && + stdlib_base_is_integer( y ) == false + ) { + // Signal NaN... + return 0.0/0.0; + } + ax = stdlib_base_abs( x ); + + // Remove the sign bits (i.e., get absolute values): + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + ahy = ( hy & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + // Extract the sign bits: + sx = ( hx >> HIGH_NUM_NONSIGN_BITS ); + sy = ( hy >> HIGH_NUM_NONSIGN_BITS ); + + // Determine the sign of the result... + if ( sx && stdlib_base_is_odd( y ) ) { + sx = -1; + } else { + sx = 1; + } + // Case 1: `|y|` is huge... + + // |y| > 2^31 + if ( ahy > HIGH_BIASED_EXP_31 ) { + // `|y| > 2^64`, then must over- or underflow... + if ( ahy > HIGH_BIASED_EXP_64 ) { + return y_is_huge( xc, y ); + } + // Over- or underflow if `x` is not close to unity... + if ( ahx < HIGH_MAX_NEAR_UNITY ) { + // y < 0 + if ( sy == 1 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( ahx > HIGH_BIASED_EXP_0 ) { + // y > 0 + if ( sy == 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + // At this point, `|1-x|` is tiny (`<= 2^-20`). Suffice to compute `log(x)` by `x - x^2/2 + x^3/3 - x^4/4`. + logx( ax, &t1, &t2 ); + } + // Case 2: `|y|` is not huge... + else { + log2ax( ax, ahx, &t1, &t2 ); + } + + // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... + y1 = y; + stdlib_base_float64_set_low_word( 0, &y1 ); + lp = ( ( y - y1 ) * t1 ) + ( y * t2 ); + hp = y1 * t1; + z = lp + hp; + + // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner. + stdlib_base_float64_to_words( z, &j, &i ); + jc = (int32_t)j; + ic = (int32_t)i; + + // z >= 1024 + if ( jc >= HIGH_BIASED_EXP_10 ) { + // z > 1024 + if ( ( ( jc - HIGH_BIASED_EXP_10 )|ic ) != 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + if ( ( lp + OVT ) > ( z - hp ) ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + } + // z <= -1075 + else if ( ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { + // z < -1075 + if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( lp <= ( z - hp ) ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + } + // Compute `2^(hp+lp)`... + z = pow2( j, hp, lp ); + return sx * z; +} From 2efc1cf0ab9a8b63c011c194bbf4432b67ffe553 Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Thu, 12 Dec 2024 14:07:52 +0530 Subject: [PATCH 4/7] docs: add include path in package.json --- .../base/special/falling-factorial/README.md | 91 +++ .../benchmark/c/native/benchmark.c | 2 +- .../special/falling-factorial/binding.gyp | 170 ++++ .../special/falling-factorial/include.gypi | 53 ++ .../special/falling-factorial/manifest.json | 96 +++ .../special/falling-factorial/package.json | 3 + .../special/falling-factorial/src/addon.c | 1 + .../base/special/falling-factorial/src/main.c | 771 +++--------------- .../special/falling-factorial/test/test.js | 24 +- .../falling-factorial/test/test.native.js | 162 ++++ 10 files changed, 706 insertions(+), 667 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/manifest.json create mode 100644 lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/README.md b/lib/node_modules/@stdlib/math/base/special/falling-factorial/README.md index 2d0fe36c6667..40a00b3bbd49 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/README.md +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/README.md @@ -100,6 +100,97 @@ for ( i = 0; i < 100; i++ ) { + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/falling_factorial.h" +``` + +#### stdlib_base_falling_factorial( x, n ) + +Evaluates the [falling factorial][falling-and-rising-factorials] of `x` and `n`. + +```c +double out = stdlib_base_falling_factorial( 0.9, 5 ); +// returns ~0.644 + +out = stdlib_base_falling_factorial( -9.0, 3 ); +// returns -990.0 +``` + +The function accepts the following arguments: + +- **x**: `[in] double` first function parameter. +- **n**: `[in] int32_t` second function parameter. + +```c +double stdlib_base_falling_factorial( const double x, const int32_t n ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/falling_factorial.h" +#include +#include + +int main( void ) { + const double x[] = { -10.0, -7.78, -5.56, -3.33, -1.11, 1.11, 3.33, 5.56, 7.78, 10.0 }; + const int32_t n[] = { 5, 4, 3, 2, 1, 0, -1, -2, -3, -4 }; + + double v; + int i; + for ( i = 0; i < 10; i++ ) { + v = stdlib_base_falling_factorial( x[ i ], n[ i ] ); + printf( "x: %lf, n: %d, falling_factorial(x,n): %lf\n", x[ i ], n[ i ], v ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c index d3ce002e4264..37e91447fece 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c @@ -90,9 +90,9 @@ static double rand_double( void ) { * @return elapsed time in seconds */ static double benchmark( void ) { - double elapsed; int32_t n[ 100 ]; double x[ 100 ]; + double elapsed; double t; double y; int i; diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp b/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp new file mode 100644 index 000000000000..ec3992233442 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp @@ -0,0 +1,170 @@ +# @license Apache-2.0 +# +# Copyright (c) 2024 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/falling-factorial/include.gypi b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi new file mode 100644 index 000000000000..575cb043c0bf --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi @@ -0,0 +1,53 @@ +# @license Apache-2.0 +# +# Copyright (c) 2024 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 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1 -static const int32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; - -static const double HUGE_VALUE = 1.0e300; - -static const double TINY_VALUE = 1.0e-300; - -// 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 -static const int32_t HIGH_MIN_NORMAL_EXP = 0x00100000; - -// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 -static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; - -// High: LN2 -static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 - -// Low: LN2 -static const double LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39 - -// 1/LN2 -static const double INV_LN2 = 1.44269504088896338700e+00; // 0x3FF71547, 0x652B82FE - -// High (24 bits): 1/LN2 -static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60000000 - -// Low: 1/LN2 -static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 - -// 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 -static const int32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; - -// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 -static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; - -// 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 -static const int32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; - -// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 -static const int32_t HIGH_SIGNIFICAND_HALF = 0x00080000; - -static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 - -// 2/(3*LN2) -static const double CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD - -// (float)CP -static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 - -// Low: CP_HI -static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 - -static const double BP [2] = { 1.0, 1.5 }; - -static const double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 - -static const double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 - -// 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 -static const int32_t HIGH_BIASED_EXP_31 = 0x41e00000; - -// 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 -static const int32_t HIGH_BIASED_EXP_64 = 0x43f00000; - -// 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 -static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; - -// 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 -static const int32_t HIGH_1075 = 0x4090cc00; - -// 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 -static const uint32_t HIGH_NEG_1075 = 0xc090cc00; - -static const int32_t HIGH_NUM_NONSIGN_BITS = 31; - -// -(1024-log2(ovfl+.5ulp)) -static const double OVT = 8.0085662595372944372e-17; - -/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ - -// BEGIN: polyval_l +#include "stdlib/constants/float64/max_safe_nth_factorial.h" +#include +#include /** -* Evaluates a polynomial. +* Computes the rising factorial of `x` and `n`. * * ## 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_l( const double x ) { - return 0.5999999999999946 + (x * (0.4285714285785502 + (x * (0.33333332981837743 + (x * (0.272728123808534 + (x * (0.23066074577556175 + (x * 0.20697501780033842))))))))); -} - -// END: polyval_l - -// BEGIN: polyval_w - -/** -* Evaluates a polynomial. +* - The rising factorial is defined as * -* ## Notes +* ```tex +* \operatorname{risingFactorial}(x, n) = x (x-1) (x-2) (x-3) \ldots (x-n+1) +* ``` * -* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* or equivalently * -* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* ```tex +* \operatorname{risingFactorial}(x, n) = \frac{ \Gamma(x + n) }{ \Gamma(x) }; +* ``` * -* @param x value at which to evaluate the polynomial -* @return evaluated polynomial -*/ -static double polyval_w( const double x ) { - return 0.5 + (x * (-0.3333333333333333 + (x * 0.25))); -} - -// END: polyval_w - -// BEGIN: polyval_p - -/** -* 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_p( const double x ) { - return 0.16666666666666602 + (x * (-0.0027777777777015593 + (x * (0.00006613756321437934 + (x * (-0.0000016533902205465252 + (x * 4.1381367970572385e-8))))))); -} - -// END: polyval_p - -/* End auto-generated functions. */ - -/** -* Evaluates the exponential function when \\( y = \pm \infty\\). -* -* @param x base -* @param y exponent -* @return output value +* @param x first function parameter +* @param n second function parameter +* @return function value * * @example -* double out = y_is_infinite( -1.0, 1.0/0.0 ); -* // returns NaN +* double out = risingFactorial( 0.9, 5 ); +* // returns ~94.766 * * @example -* double out = y_is_infinite( 1.5, -1.0/0.0 ); -* // returns 0.0 -*/ -static double y_is_infinite( const double x, const double y ) { - if ( x == -1.0 ) { +* double out = risingFactorial( -9.0, 3 ); +* // returns -504.0 + */ +static double risingFactorial( const double x, const int32_t n ) { + double result; + int32_t nc; + double xc; + bool inv; + + if ( stdlib_base_is_nan( x ) ) { return 0.0 / 0.0; // NaN } - if ( x == 1.0 ) { - return 1.0; - } - // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) - if ( ( stdlib_base_abs( x ) < 1.0 ) == ( y == STDLIB_CONSTANT_FLOAT64_PINF ) ) { - return 0.0; - } - // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) - return STDLIB_CONSTANT_FLOAT64_PINF; -} - -/** -* Evaluates the exponential function when \\(|y| > 2^64\\). -* -* @param x base -* @param y exponent -* @return output value -* -* @example -* double out = y_is_huge( 9.0, 3.6893488147419103e19 ); -* // returns Infinity -* -* @example -* double out = y_is_huge( -3.14, -3.6893488147419103e19 ); -* // returns 0.0 -*/ -static double y_is_huge( const double x, const double y ) { - uint32_t ahx; - uint32_t hx; - - stdlib_base_float64_get_high_word( x, &hx ); - ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); - - if ( ahx <= HIGH_MAX_NEAR_UNITY ) { - if ( y < 0 ) { - // Signal overflow... - return HUGE_VALUE * HUGE_VALUE; + nc = n; + xc = x; + inv = false; + if ( xc < 0.0 ) { + // For `xc < 0`, we really have a falling factorial, modulo a possible change of sign. Note that the falling factorial isn't defined for negative `nc`, so we'll get rid of that case first: + if ( nc < 0 ) { + xc += nc; + nc = -nc; + inv = true; } - // Signal underflow... - return TINY_VALUE * TINY_VALUE; - } - // `x` has a biased exponent greater than or equal to `0`... - - if ( y > 0 ) { - // Signal overflow... - return HUGE_VALUE * HUGE_VALUE; - } - // Signal underflow... - return TINY_VALUE * TINY_VALUE; -} - -/** -* Evaluates the exponential function when \\(|x| = 0\\). -* -* @param x base -* @param y exponent -* @return output value -* -* @example -* double out = x_is_zero( 0.0, 2 ); -* // returns 0.0 -* -* @example -* double out = x_is_zero( 0.0, -9 ); -* // returns Infinity -*/ -static double x_is_zero( const double x, const double y ) { - if ( y == STDLIB_CONSTANT_FLOAT64_NINF ) { - return STDLIB_CONSTANT_FLOAT64_PINF; + result = ( ( nc & 1 ) ? -1.0 : 1.0 ) * stdlib_base_falling_factorial( -xc, nc ); + if ( inv ) { + result = 1.0 / result; + } + return result; } - if ( y == STDLIB_CONSTANT_FLOAT64_PINF ) { - return 0.0; + if ( nc == 0 ) { + return 1.0; } - if ( y > 0.0 ) { - if ( stdlib_base_is_odd( y ) ) { - return x; // handles +-0 + if ( xc == 0.0 ) { + if ( nc < 0 ) { + return -stdlib_base_gamma_delta_ratio( xc + 1.0, -nc ); } return 0.0; } - // y < 0.0 - if ( stdlib_base_is_odd( y ) ) { - return stdlib_base_copysign( STDLIB_CONSTANT_FLOAT64_PINF, x ); // handles +-0 - } - return STDLIB_CONSTANT_FLOAT64_PINF; -} - -/** -* Computes \\(2^{\mathrm{hp} + \mathrm{lp}\\). -* -* @param j high word of `hp + lp` -* @param hp first power summand -* @param lp second power summand -* @return function value -* -* @example -* double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); -* // returns ~0.79 -*/ -static double pow2( uint32_t j, const double hp, const double lp ) { - uint32_t tmp; - double hpc; - int32_t jc; - int32_t nc; - int32_t i; - int32_t k; - int32_t n; - double t1; - double t; - double r; - double u; - double v; - double w; - double z; - - hpc = hp; - jc = (int32_t)j; - i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); - k = ( ( i >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); - n = 0; - nc = (int32_t)n; - - // `|z| > 0.5`, set `n = z+0.5` - if ( i > HIGH_BIASED_EXP_NEG_1 ) { - n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ); - k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n - tmp = ( ( n & ~( HIGH_SIGNIFICAND_MASK >> k ) ) ); - t = 0.0; - stdlib_base_float64_set_high_word( tmp, &t ); - n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS - k ) ); - nc = (int32_t)n; - if ( jc < 0 ) { - nc = -nc; - } - hpc -= t; + if ( ( xc < 1.0 ) && ( xc + nc < 0.0 ) ) { + result = stdlib_base_gamma_delta_ratio( 1.0 - xc, -nc ); + return ( nc & 1 ) ? -result : result; } - t = lp + hpc; - stdlib_base_float64_set_low_word( 0, &t ); - u = t * LN2_HI; - v = ( ( lp - ( t - hpc ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); - z = u + v; - w = v - ( z - u ); - t = z * z; - t1 = z - ( t * polyval_p( t ) ); - r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); - z = 1.0 - ( r - z ); - stdlib_base_float64_get_high_word( z, &j ); - j = j + ( nc << STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ); - jc = (int32_t)j; - // Check for subnormal output... - if ( ( jc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) <= 0 ) { - z = stdlib_base_ldexp( z, nc ); - } else { - stdlib_base_float64_set_high_word( j, &z ); - } - return z; + // We don't optimize this for small `nc`, because `stdlib_base_gamma_delta_ratio` is already optimized for that use case: + return 1.0 / stdlib_base_gamma_delta_ratio( xc, nc ); } /** -* Computes \\(\operatorname{log}(x)\\) assuming \\(|1-x|\\) is small and using the approximation \\(x - x^2/2 + x^3/3 - x^4/4\\). +* Computes the falling factorial of `x` and `n`. * -* @param ax absolute value of x -* @param o1 destination for output1 -* @param o2 destination for output2 +* ## Notes * -* @example -* logx( 9.0, &o1, &o2 ); -*/ -void logx( const double ax, double *o1, double *o2 ) { - double t2; - double t1; - double t; - double w; - double u; - double v; - - t = ax - 1.0; // `t` has `20` trailing zeros - w = t * t * polyval_w( t ); - u = INV_LN2_HI * t; // `INV_LN2_HI` has `21` significant bits - v = ( t * INV_LN2_LO ) - ( w * INV_LN2 ); - t1 = u + v; - stdlib_base_float64_set_low_word( 0, &t1 ); - t2 = v - (t1 - u); - - *o1 = t1; - *o2 = t2; -} - -/** -* Computes \\(\operatorname{log2}(ax)\\). +* - The falling factorial is defined as * -* @param ax absolute value of x -* @param ahx high word of ax -* @param o1 destination for output1 -* @param o2 destination for output2 +* ```tex +* \operatorname{fallingFactorial}(x, n) = x (x-1) (x-2) (x-3) \ldots (x-n+1) +* ``` * -* @example -* log2ax( 9.0, 1075970048, &o1, &o2 ); -*/ -void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { - uint32_t ahxcc; - uint32_t tmp; - int32_t ahxc; - double axc; - int32_t n; - int32_t j; - int32_t k; - double ss; // hs + ls - double s2; // ss squared - double hs; - double ls; - double ht; - double lt; - double bp; // BP constant - double dp; // DP constant - double hp; - double lp; - double hz; - double lz; - double t1; - double t2; - double t; - double r; - double u; - double v; - - n = 0; - axc = ax; - ahxc = ahx; - ahxcc = (uint32_t)ahx; - // Check if `x` is subnormal... - if ( ahxc < HIGH_MIN_NORMAL_EXP ) { - axc *= TWO53; - n -= 53; - stdlib_base_float64_get_high_word( axc, &ahxcc ); - ahxc = (int32_t)ahxcc; - } - // Extract the unbiased exponent of `x`: - n += ( ( ahxc >> STDLIB_CONSTANT_FLOAT64_NUM_HIGH_WORD_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); - - // Isolate the significand bits of `x`: - j = ( ahxc & HIGH_SIGNIFICAND_MASK ); - - // Normalize `ahx` by setting the (biased) exponent to `1023`: - ahxc = ( j|HIGH_BIASED_EXP_0 ); - - // Determine the interval of `|x|` by comparing significand bits... - - // |x| < sqrt(3/2) - if ( j <= 0x3988E ) { // 0 00000000000 00111001100010001110 - k = 0; - } - // |x| < sqrt(3) - else if ( j < 0xBB67A ) { // 0 00000000000 10111011011001111010 - k = 1; - } - // |x| >= sqrt(3) - else { - k = 0; - n += 1; - ahxc -= HIGH_MIN_NORMAL_EXP; - } - ahxcc = (uint32_t)ahxc; - - // Load the normalized high word into `|x|`: - stdlib_base_float64_set_high_word( ahxcc, &axc ); - - // Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: - bp = BP[ k ]; // BP[0] = 1.0, BP[1] = 1.5 - u = axc - bp; // (x-1) || (x-1.5) - v = 1.0 / (axc + bp); // 1/(x+1) || 1/(x+1.5) - ss = u * v; - hs = ss; - stdlib_base_float64_set_low_word( 0, &hs ); // set all low word (less significant significand) bits to 0s - - // Compute `ht = ax + bp` (via manipulation, i.e., bit flipping, of the high word): - tmp = ( ( ahxc >> 1 ) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; - tmp += ( k << 18 ); // `(k<<18)` can be considered the word equivalent of `1.0` or `1.5` - ht = 0.0; - stdlib_base_float64_set_high_word( tmp, &ht ); - lt = axc - ( ht - bp ); - ls = v * ( ( u - ( hs * ht ) ) - ( hs * lt ) ); - - // Compute `log(ax)`... - - s2 = ss * ss; - r = s2 * s2 * polyval_l( s2 ); - r += ls * (hs + ss); - s2 = hs * hs; - ht = 3.0 + s2 + r; - stdlib_base_float64_set_low_word( 0, &ht ); - lt = r - ( ( ht - 3.0 ) - s2 ); - - // u+v = ss*(1+...): - u = hs * ht; - v = ( ls * ht ) + ( lt * ss ); - - // 2/(3LN2) * (ss+...): - hp = u + v; - stdlib_base_float64_set_low_word( 0, &hp ); - lp = v - ( hp - u ); - hz = CP_HI * hp; // CP_HI+CP_LO = 2/(3*LN2) - lz = ( CP_LO * hp ) + ( lp*CP ) + DP_LO[ k ]; - - // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz - dp = DP_HI[ k ]; - t = n; - t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) - stdlib_base_float64_set_low_word( 0, &t1 ); - t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); - - *o1 = t1; - *o2 = t2; -} - -/** -* Evaluates the exponential function. -* -* @param x base -* @param y exponent +* @param x first function parameter +* @param n second function parameter * @return function value * * @example -* double out = stdlib_base_pow( 2.0, 3.0 ); -* // returns 8.0 +* double out = stdlib_base_falling_factorial( 0.9, 5 ); +* // returns ~0.644 * * @example -* double out = stdlib_base_pow( 4.0, 0.5 ); -* // returns 2.0 +* double out = stdlib_base_falling_factorial( -9.0, 3 ); +* // returns -990.0 */ -double stdlib_base_pow( const double x, const double y ) { - uint32_t hx; // high word `x` - uint32_t lx; // low word `x` - uint32_t hy; // high word `y` - uint32_t ly; // low word `y` - int32_t ahx; // absolute value high word `x` - int32_t ahy; // absolute value high word `y` - uint32_t i; - uint32_t j; - int32_t ic; - int32_t jc; - int32_t sx; // sign `x` - int32_t sy; // sign `y` - double ax; // absolute value `x` - double y1; - double hp; - double lp; +double stdlib_base_falling_factorial( const double x, const int32_t n ) { + double result; + double xp1; + double n2; double t1; double t2; double xc; - double z; // y prime - xc = x; - if ( stdlib_base_is_nan( xc ) || stdlib_base_is_nan( y ) ) { - return 0.0/0.0; // NaN - } - - // Split `y` into high and low words: - stdlib_base_float64_to_words( y, &hy, &ly ); - - // Special cases `y`... - if ( ly == 0 ) { - if ( y == 0.0 ) { - return 1.0; - } - if ( y == 1.0 ) { - return xc; - } - if ( y == -1.0 ) { - return 1.0 / xc; - } - if ( y == 0.5 ) { - return stdlib_base_sqrt( xc ); - } - if ( y == -0.5 ) { - return 1.0 / stdlib_base_sqrt( xc ); - } - if ( y == 2.0 ) { - return xc * xc; - } - if ( y == 3.0 ) { - return xc * xc * xc; - } - if ( y == 4.0 ) { - xc *= xc; - return xc * xc; - } - if ( stdlib_base_is_infinite( y ) ) { - return y_is_infinite( xc, y ); - } + if ( stdlib_base_is_nan( x ) || !stdlib_base_is_nonnegative_integer( n ) ) { + return 0.0 / 0.0; // NaN } - // Split `x` into high and low words: - stdlib_base_float64_to_words( xc, &hx, &lx ); - - // Special cases `x`... - if ( lx == 0 ) { - if ( hx == 0 ) { - return x_is_zero( xc, y ); - } - if ( xc == 1.0 ) { - return 1.0; - } - if ( - xc == -1.0 && - stdlib_base_is_odd( y ) - ) { - return -1.0; - } - if ( stdlib_base_is_infinite( xc ) ) { - if ( xc == STDLIB_CONSTANT_FLOAT64_NINF ) { - // `stdlib_base_pow( 1/x, -y )` - return stdlib_base_pow( -0.0, -y ); - } - if ( y < 0.0 ) { - return 0.0; - } - return STDLIB_CONSTANT_FLOAT64_PINF; - } + if ( x == 0.0 ) { + return 0.0; } - if ( - xc < 0.0 && - stdlib_base_is_integer( y ) == false - ) { - // Signal NaN... - return 0.0/0.0; + if ( x < 0.0 ) { + // For `x < 0`, we really have a rising factorial modulo a possible change of sign: + return ( ( n & 1 ) ? -1.0 : 1.0 ) * risingFactorial( -x, n ); } - ax = stdlib_base_abs( x ); - - // Remove the sign bits (i.e., get absolute values): - ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); - ahy = ( hy & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); - - // Extract the sign bits: - sx = ( hx >> HIGH_NUM_NONSIGN_BITS ); - sy = ( hy >> HIGH_NUM_NONSIGN_BITS ); - - // Determine the sign of the result... - if ( sx && stdlib_base_is_odd( y ) ) { - sx = -1; - } else { - sx = 1; + if ( n == 0 ) { + return 1.0; } - // Case 1: `|y|` is huge... - - // |y| > 2^31 - if ( ahy > HIGH_BIASED_EXP_31 ) { - // `|y| > 2^64`, then must over- or underflow... - if ( ahy > HIGH_BIASED_EXP_64 ) { - return y_is_huge( xc, y ); - } - // Over- or underflow if `x` is not close to unity... - if ( ahx < HIGH_MAX_NEAR_UNITY ) { - // y < 0 - if ( sy == 1 ) { - // Signal overflow... - return sx * HUGE_VALUE * HUGE_VALUE; - } - // Signal underflow... - return sx * TINY_VALUE * TINY_VALUE; - } - if ( ahx > HIGH_BIASED_EXP_0 ) { - // y > 0 - if ( sy == 0 ) { - // Signal overflow... - return sx * HUGE_VALUE * HUGE_VALUE; + if ( x < 0.5 ) { + // Computing `1 + x` will throw away digits, so split up calculation... + if ( n > STDLIB_CONSTANT_FLOAT64_MAX_SAFE_NTH_FACTORIAL - 2 ) { + // Given a ratio of two very large numbers, we need to split the calculation up into two blocks: + t1 = x * stdlib_base_falling_factorial( x - 1.0, STDLIB_CONSTANT_FLOAT64_MAX_SAFE_NTH_FACTORIAL - 2 ); + t2 = stdlib_base_falling_factorial( x - STDLIB_CONSTANT_FLOAT64_MAX_SAFE_NTH_FACTORIAL + 1.0, n - STDLIB_CONSTANT_FLOAT64_MAX_SAFE_NTH_FACTORIAL + 1 ); + if ( ( STDLIB_CONSTANT_FLOAT64_MAX / stdlib_base_abs( t1 ) ) < stdlib_base_abs( t2 ) ) { + return STDLIB_CONSTANT_FLOAT64_PINF; } - // Signal underflow... - return sx * TINY_VALUE * TINY_VALUE; + return t1 * t2; } - // At this point, `|1-x|` is tiny (`<= 2^-20`). Suffice to compute `log(x)` by `x - x^2/2 + x^3/3 - x^4/4`. - logx( ax, &t1, &t2 ); + return x * stdlib_base_falling_factorial( x - 1.0, n - 1 ); } - // Case 2: `|y|` is not huge... - else { - log2ax( ax, ahx, &t1, &t2 ); - } - - // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... - y1 = y; - stdlib_base_float64_set_low_word( 0, &y1 ); - lp = ( ( y - y1 ) * t1 ) + ( y * t2 ); - hp = y1 * t1; - z = lp + hp; - - // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner. - stdlib_base_float64_to_words( z, &j, &i ); - jc = (int32_t)j; - ic = (int32_t)i; - - // z >= 1024 - if ( jc >= HIGH_BIASED_EXP_10 ) { - // z > 1024 - if ( ( ( jc - HIGH_BIASED_EXP_10 )|ic ) != 0 ) { - // Signal overflow... - return sx * HUGE_VALUE * HUGE_VALUE; - } - if ( ( lp + OVT ) > ( z - hp ) ) { - // Signal overflow... - return sx * HUGE_VALUE * HUGE_VALUE; - } - } - // z <= -1075 - else if ( ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { - // z < -1075 - if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { - // Signal underflow... - return sx * TINY_VALUE * TINY_VALUE; - } - if ( lp <= ( z - hp ) ) { - // Signal underflow... - return sx * TINY_VALUE * TINY_VALUE; - } - } - // Compute `2^(hp+lp)`... - z = pow2( j, hp, lp ); - return sx * z; + xc = x; + if ( xc <= n - 1.0 ) { + // `xc + 1 - n` will be negative and computing the ratio of two gammas will not work, so split the product up into three parts: + xp1 = xc + 1.0; + n2 = stdlib_base_abs( stdlib_base_floor( xp1 ) ); + if ( n2 == xp1 ) { + return 0.0; + } + result = stdlib_base_gamma_delta_ratio( xp1, -(int32_t)n2 ); + xc -= n2; + result *= xc; + n2 += 1.0; + if ( (int32_t)n2 < n ) { + result *= stdlib_base_falling_factorial( xc - 1.0, n - (int32_t)n2 ); + } + return result; + } + // Simple case: just the ratio of two (positive argument) gamma functions. Note that we don't optimize this for small `n`, because `gammaDeltaRatio` is already optimized for that use case: + return stdlib_base_gamma_delta_ratio( xc + 1.0, -n ); } diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.js index 926b23acd385..84b19b1618bf 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.js @@ -45,47 +45,47 @@ tape( 'main export is a function', function test( t ) { tape( 'the function returns `NaN` if provided a negative `n`', function test( t ) { var val = fallingFactorial( 0.5, -1.0 ); - t.strictEqual( isnan( val ), true, 'returns NaN' ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); val = fallingFactorial( 0.5, -2.0 ); - t.strictEqual( isnan( val ), true, 'returns NaN' ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); t.end(); }); tape( 'the function returns `NaN` if provided a non-integer `n`', function test( t ) { var val = fallingFactorial( 0.5, 1.5 ); - t.strictEqual( isnan( val ), true, 'returns NaN' ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); val = fallingFactorial( 0.5, 0.3 ); - t.strictEqual( isnan( val ), true, 'returns NaN' ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); val = fallingFactorial( 2.0, -0.3 ); - t.strictEqual( isnan( val ), true, 'returns NaN' ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); t.end(); }); tape( 'the function returns `0` for `x = 0`', function test( t ) { var val = fallingFactorial( 0.0, 2 ); - t.strictEqual( val, 0, 'returns 0' ); + t.strictEqual( val, 0, 'returns expected value' ); val = fallingFactorial( 0.0, 1 ); - t.strictEqual( val, 0, 'returns 0' ); + t.strictEqual( val, 0, 'returns expected value' ); t.end(); }); tape( 'the function returns `1` for `n = 0`', function test( t ) { var val = fallingFactorial( 2.0, 0 ); - t.strictEqual( val, 1, 'returns 1' ); + t.strictEqual( val, 1, 'returns expected value' ); val = fallingFactorial( -1.7, 0 ); - t.strictEqual( val, 1, 'returns 1' ); + t.strictEqual( val, 1, 'returns expected value' ); t.end(); }); tape( 'the function returns `+infinity` for small `x` and large `n`', function test( t ) { - t.equal( fallingFactorial( 0.2, 200 ), PINF, 'equals +infinity' ); - t.equal( fallingFactorial( 0.4, 175 ), PINF, 'equals +infinity' ); - t.equal( fallingFactorial( 0.1, 250 ), PINF, 'equals +infinity' ); + t.equal( fallingFactorial( 0.2, 200 ), PINF, 'returns expected value' ); + t.equal( fallingFactorial( 0.4, 175 ), PINF, 'returns expected value' ); + t.equal( fallingFactorial( 0.1, 250 ), PINF, 'returns expected value' ); t.end(); }); diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js new file mode 100644 index 000000000000..d698824ff7d5 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js @@ -0,0 +1,162 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2024 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 abs = require( '@stdlib/math/base/special/abs' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var tryRequire = require( '@stdlib/utils/try-require' ); + + +// FIXTURES // + +var negative = require( './fixtures/cpp/negative.json' ); +var large = require( './fixtures/cpp/large.json' ); +var small = require( './fixtures/cpp/small.json' ); + + +// VARIABLES // + +var fallingFactorial = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( fallingFactorial instanceof Error ) +}; + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof fallingFactorial, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'the function returns `NaN` if provided a negative `n`', opts, function test( t ) { + var val = fallingFactorial( 0.5, -1.0 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + + val = fallingFactorial( 0.5, -2.0 ); + t.strictEqual( isnan( val ), true, 'returns expected value' ); + t.end(); +}); + +tape( 'the function returns `0` for `x = 0`', opts, function test( t ) { + var val = fallingFactorial( 0.0, 2 ); + t.strictEqual( val, 0, 'returns expected value' ); + + val = fallingFactorial( 0.0, 1 ); + t.strictEqual( val, 0, 'returns expected value' ); + t.end(); +}); + +tape( 'the function returns `1` for `n = 0`', opts, function test( t ) { + var val = fallingFactorial( 2.0, 0 ); + t.strictEqual( val, 1, 'returns expected value' ); + + val = fallingFactorial( -1.7, 0 ); + t.strictEqual( val, 1, 'returns expected value' ); + t.end(); +}); + +tape( 'the function returns `+infinity` for small `x` and large `n`', opts, function test( t ) { + t.equal( fallingFactorial( 0.2, 200 ), PINF, 'returns expected value' ); + t.equal( fallingFactorial( 0.4, 175 ), PINF, 'returns expected value' ); + t.equal( fallingFactorial( 0.1, 250 ), PINF, 'returns expected value' ); + t.end(); +}); + +tape( 'the function evaluates the falling factorial for large `x`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var n; + var x; + var y; + + expected = large.expected; + x = large.x; + n = large.n; + for ( i = 0; i < x.length; i++ ) { + y = fallingFactorial( x[i], n[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+'. n: '+n[i]+', y: '+y+'. expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 80.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. n: '+n[i]+'. y: '+y+'. E: '+expected[i]+'. Δ: '+delta+'. tol: '+tol ); + } + } + t.end(); +}); + +tape( 'the function evaluates the falling factorial for small `x`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var n; + var x; + var y; + + expected = small.expected; + x = small.x; + n = small.n; + for ( i = 0; i < x.length; i++ ) { + y = fallingFactorial( x[i], n[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+'. n: '+n[i]+', y: '+y+'. expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 80.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. n: '+n[i]+'. y: '+y+'. E: '+expected[i]+'. Δ: '+delta+'. tol: '+tol ); + } + } + t.end(); +}); + +tape( 'the function evaluates the falling factorial for negative `x`', opts, function test( t ) { + var expected; + var delta; + var tol; + var i; + var n; + var x; + var y; + + expected = negative.expected; + x = negative.x; + n = negative.n; + for ( i = 0; i < x.length; i++ ) { + y = fallingFactorial( x[i], n[i] ); + if ( y === expected[i] ) { + t.strictEqual( y, expected[i], 'x: '+x[i]+'. n: '+n[i]+', y: '+y+'. expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 100.0 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. n: '+n[i]+'. y: '+y+'. E: '+expected[i]+'. Δ: '+delta+'. tol: '+tol ); + } + } + t.end(); +}); From 98eba77156b224d117855adf3e207fdd87e1d7ad Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Thu, 12 Dec 2024 14:20:02 +0530 Subject: [PATCH 5/7] docs: reorder directory paths --- .../@stdlib/math/base/special/falling-factorial/package.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/package.json b/lib/node_modules/@stdlib/math/base/special/falling-factorial/package.json index 65372cd3dabb..ede81d490ea9 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/package.json +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/package.json @@ -20,8 +20,8 @@ "doc": "./docs", "example": "./examples", "include": "./include", - "src": "./src", "lib": "./lib", + "src": "./src", "test": "./test" }, "types": "./docs/types", From 8fbad2207766a342af0d376f139cd984aca5b4d9 Mon Sep 17 00:00:00 2001 From: stdlib-bot <82920195+stdlib-bot@users.noreply.github.com> Date: Wed, 1 Jan 2025 20:09:16 +0000 Subject: [PATCH 6/7] chore: update copyright years --- .../special/falling-factorial/benchmark/benchmark.native.js | 2 +- .../base/special/falling-factorial/benchmark/c/native/Makefile | 2 +- .../special/falling-factorial/benchmark/c/native/benchmark.c | 2 +- .../@stdlib/math/base/special/falling-factorial/binding.gyp | 2 +- .../math/base/special/falling-factorial/examples/c/Makefile | 2 +- .../math/base/special/falling-factorial/examples/c/example.c | 2 +- .../@stdlib/math/base/special/falling-factorial/include.gypi | 2 +- .../include/stdlib/math/base/special/falling_factorial.h | 2 +- .../@stdlib/math/base/special/falling-factorial/lib/native.js | 2 +- .../@stdlib/math/base/special/falling-factorial/src/Makefile | 2 +- .../@stdlib/math/base/special/falling-factorial/src/addon.c | 2 +- .../@stdlib/math/base/special/falling-factorial/src/main.c | 2 +- .../math/base/special/falling-factorial/test/test.native.js | 2 +- 13 files changed, 13 insertions(+), 13 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js index eed82c096450..3f5f0271e7a7 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/benchmark.native.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile index f69e9da2b4d3..a4bd7b38fd74 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/Makefile @@ -1,7 +1,7 @@ #/ # @license Apache-2.0 # -# Copyright (c) 2024 The Stdlib Authors. +# Copyright (c) 2025 The Stdlib Authors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c index 37e91447fece..1c28a3ac773b 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/benchmark/c/native/benchmark.c @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp b/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp index ec3992233442..68a1ca11d160 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/binding.gyp @@ -1,6 +1,6 @@ # @license Apache-2.0 # -# Copyright (c) 2024 The Stdlib Authors. +# Copyright (c) 2025 The Stdlib Authors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile index 6aed70daf167..25ced822f96a 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/Makefile @@ -1,7 +1,7 @@ #/ # @license Apache-2.0 # -# Copyright (c) 2024 The Stdlib Authors. +# Copyright (c) 2025 The Stdlib Authors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c index 800863ab6941..80a7af13c5ec 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/c/example.c @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi index 575cb043c0bf..ecfaf82a3279 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include.gypi @@ -1,6 +1,6 @@ # @license Apache-2.0 # -# Copyright (c) 2024 The Stdlib Authors. +# Copyright (c) 2025 The Stdlib Authors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h index 14d77f9f0fbf..b0990c2dfdca 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/include/stdlib/math/base/special/falling_factorial.h @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js index 3a7b0100d2cc..dfc7d7b341ba 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/lib/native.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile index bcf18aa46655..7733b6180cb4 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/Makefile @@ -1,7 +1,7 @@ #/ # @license Apache-2.0 # -# Copyright (c) 2024 The Stdlib Authors. +# Copyright (c) 2025 The Stdlib Authors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c index cfdc47ba5d9e..5321f79f0a13 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/addon.c @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c index ce11bb4bc6ab..207466016c71 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/src/main.c @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js index d698824ff7d5..9897234f431a 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/test/test.native.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. From 31f0f1f61659e600178cd2b54a264992854dc395 Mon Sep 17 00:00:00 2001 From: Philipp Burckhardt Date: Wed, 1 Jan 2025 20:52:02 -0500 Subject: [PATCH 7/7] chore: revert year change in license comment Signed-off-by: Philipp Burckhardt --- .../math/base/special/falling-factorial/examples/index.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js index 828cb7aeff64..3eb35464d32a 100644 --- a/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js +++ b/lib/node_modules/@stdlib/math/base/special/falling-factorial/examples/index.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2024 The Stdlib Authors. +* Copyright (c) 2018 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.