diff --git a/CHANGELOG.md b/CHANGELOG.md index 970099721..2478c3f5b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@
-## Unreleased (2024-08-01) +## Unreleased (2024-08-02)
@@ -3755,6 +3755,28 @@ This release closes the following issue: +
+ +#### [@stdlib/math/base/special/gammaln](https://github.com/stdlib-js/stdlib/tree/develop/lib/node_modules/%40stdlib/math/base/special/gammaln) + +
+ +
+ +##### Features + +- [`d745814`](https://github.com/stdlib-js/stdlib/commit/d74581415b30604ce6db2b3a2c4f242040eb3e8e) - add C implementation for `math/base/special/gammaln` [(#2636)](https://github.com/stdlib-js/stdlib/pull/2636) + +
+ + + +
+ +
+ + +
#### [@stdlib/math/base/special/gammasgn](https://github.com/stdlib-js/stdlib/tree/develop/lib/node_modules/%40stdlib/math/base/special/gammasgn) @@ -5730,6 +5752,7 @@ A total of 30 people contributed to this release. Thank you to the following con
+- [`d745814`](https://github.com/stdlib-js/stdlib/commit/d74581415b30604ce6db2b3a2c4f242040eb3e8e) - **feat:** add C implementation for `math/base/special/gammaln` [(#2636)](https://github.com/stdlib-js/stdlib/pull/2636) _(by Gunj Joshi)_ - [`eb7f8ed`](https://github.com/stdlib-js/stdlib/commit/eb7f8ed40d9a514a673f5152484f62042dc24e67) - **docs:** add missing `stdint` includes [(#2728)](https://github.com/stdlib-js/stdlib/pull/2728) _(by Gunj Joshi, Athan Reines)_ - [`46cda32`](https://github.com/stdlib-js/stdlib/commit/46cda325a622274755194c2dd6415d3d04ab5a5b) - **feat:** add `LL_D` in `math/base/napi/binary` [(#2726)](https://github.com/stdlib-js/stdlib/pull/2726) _(by Gunj Joshi, Athan Reines)_ - [`9816dec`](https://github.com/stdlib-js/stdlib/commit/9816dece59ddf974693cf1626f3b5823d652e3c4) - **bench:** fix description [(#2709)](https://github.com/stdlib-js/stdlib/pull/2709) _(by Gunj Joshi)_ diff --git a/base/special/gammaln/README.md b/base/special/gammaln/README.md index b645aad9f..497e0d856 100644 --- a/base/special/gammaln/README.md +++ b/base/special/gammaln/README.md @@ -83,6 +83,95 @@ for ( i = 0; i < x.length; i++ ) { + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/gammaln.h" +``` + +#### stdlib_base_gammaln( x ) + +Evaluates the [natural logarithm][@stdlib/math/base/special/ln] of the [gamma function][@stdlib/math/base/special/gamma]. + +```c +double out = stdlib_base_gammaln( 2.0 ); +// returns 0.0 + +out = stdlib_base_gammaln( 4.0 ); +// returns ~1.792 +``` + +The function accepts the following arguments: + +- **x**: `[in] double` input value. + +```c +double stdlib_base_gammaln( const double x ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/gammaln.h" +#include +#include + +int main( void ) { + const double x[] = { 4.0, -1.5, -0.5, 0.5 }; + + double y; + int i; + for ( i = 0; i < 4; i++ ) { + y = stdlib_base_gammaln( x[ i ] ); + printf( "gammaln(%lf) = %lf\n", x[ i ], y ); + } +} +``` + +
+ + + +
+ + +
diff --git a/base/special/gammaln/benchmark/benchmark.native.js b/base/special/gammaln/benchmark/benchmark.native.js new file mode 100644 index 000000000..185482b4d --- /dev/null +++ b/base/special/gammaln/benchmark/benchmark.native.js @@ -0,0 +1,60 @@ +/** +* @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 bench = require( '@stdlib/bench' ); +var randu = require( '@stdlib/random/base/randu' ); +var isnan = require( './../../../../base/assert/is-nan' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var gammaln = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( gammaln instanceof Error ) +}; + + +// MAIN // + +bench( pkg+'::native', opts, function benchmark( b ) { + var x; + var y; + var i; + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + x = ( randu() * 1000.0 ) - 500.0; + y = gammaln( x ); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); diff --git a/base/special/gammaln/benchmark/c/native/Makefile b/base/special/gammaln/benchmark/c/native/Makefile new file mode 100644 index 000000000..f69e9da2b --- /dev/null +++ b/base/special/gammaln/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/base/special/gammaln/benchmark/c/native/benchmark.c b/base/special/gammaln/benchmark/c/native/benchmark.c new file mode 100644 index 000000000..ec0342a3a --- /dev/null +++ b/base/special/gammaln/benchmark/c/native/benchmark.c @@ -0,0 +1,133 @@ +/** +* @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/gammaln.h" +#include +#include +#include +#include +#include + +#define NAME "gammaln" +#define ITERATIONS 1000000 +#define REPEATS 3 + +/** +* Prints the TAP version. +*/ +static void print_version( void ) { + printf( "TAP version 13\n" ); +} + +/** +* Prints the TAP summary. +* +* @param total total number of tests +* @param passing total number of passing tests +*/ +static void print_summary( int total, int passing ) { + printf( "#\n" ); + printf( "1..%d\n", total ); // TAP plan + printf( "# total %d\n", total ); + printf( "# pass %d\n", passing ); + printf( "#\n" ); + printf( "# ok\n" ); +} + +/** +* Prints benchmarks results. +* +* @param elapsed elapsed time in seconds +*/ +static void print_results( double elapsed ) { + double rate = (double)ITERATIONS / elapsed; + printf( " ---\n" ); + printf( " iterations: %d\n", ITERATIONS ); + printf( " elapsed: %0.9f\n", elapsed ); + printf( " rate: %0.9f\n", rate ); + printf( " ...\n" ); +} + +/** +* Returns a clock time. +* +* @return clock time +*/ +static double tic( void ) { + struct timeval now; + gettimeofday( &now, NULL ); + return (double)now.tv_sec + (double)now.tv_usec / 1.0e6; +} + +/** +* Generates a random number on the interval [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; + double x; + double y; + double t; + int i; + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + x = ( 1000.0 * rand_double() ) - 500.0; + y = stdlib_base_gammaln( x ); + if ( y != y ) { + printf( "should not return NaN\n" ); + break; + } + } + elapsed = tic() - t; + if ( y != y ) { + printf( "should not return NaN\n" ); + } + return elapsed; +} + +/** +* Main execution sequence. +*/ +int main( void ) { + double elapsed; + int i; + + // Use the current time to seed the random number generator: + srand( time( NULL ) ); + + print_version(); + for ( i = 0; i < REPEATS; i++ ) { + printf( "# c::native::%s\n", NAME ); + elapsed = benchmark(); + print_results( elapsed ); + printf( "ok %d benchmark finished\n", i+1 ); + } + print_summary( REPEATS, REPEATS ); +} diff --git a/base/special/gammaln/binding.gyp b/base/special/gammaln/binding.gyp new file mode 100644 index 000000000..ec3992233 --- /dev/null +++ b/base/special/gammaln/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/base/special/gammaln/examples/c/Makefile b/base/special/gammaln/examples/c/Makefile new file mode 100644 index 000000000..6aed70daf --- /dev/null +++ b/base/special/gammaln/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/base/special/gammaln/examples/c/example.c b/base/special/gammaln/examples/c/example.c new file mode 100644 index 000000000..e6825b98b --- /dev/null +++ b/base/special/gammaln/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/gammaln.h" +#include + +int main( void ) { + const double x[] = { 4.0, -1.5, -0.5, 0.5 }; + + double y; + int i; + for ( i = 0; i < 4; i++ ) { + y = stdlib_base_gammaln( x[ i ] ); + printf( "gammaln(%lf) = %lf\n", x[ i ], y ); + } +} diff --git a/base/special/gammaln/include.gypi b/base/special/gammaln/include.gypi new file mode 100644 index 000000000..575cb043c --- /dev/null +++ b/base/special/gammaln/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': [ + ' TT = -(tail of TF) @@ -269,7 +269,7 @@ function gammaln( x ) { } else { isNegative = false; } - // If |x| < 2**-70, return -ln(|x|) + // If |x| < 2**-56, return -ln(|x|) if ( x < TINY ) { return -ln( x ); } @@ -381,15 +381,15 @@ function gammaln( x ) { r += ln( z ); } } - // 8 <= x < 2**58 - else if ( x < TWO58 ) { + // 8 <= x < 2**56 + else if ( x < TWO56 ) { t = ln( x ); z = 1.0 / x; y = z * z; w = WC + (z*polyvalW( y )); r = ((x-0.5)*(t-1.0)) + w; } - // 2**58 <= x <= Inf + // 2**56 <= x <= Inf else { r = x * ( ln(x)-1.0 ); } diff --git a/base/special/gammaln/lib/native.js b/base/special/gammaln/lib/native.js new file mode 100644 index 000000000..962c61b10 --- /dev/null +++ b/base/special/gammaln/lib/native.js @@ -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. +*/ + +'use strict'; + +// MODULES // + +var addon = require( './../src/addon.node' ); + + +// MAIN // + +/** +* Evaluates the natural logarithm of the gamma function. +* +* @private +* @param {number} x - input value +* @returns {number} function value +* +* @example +* var v = gammaln( 1.0 ); +* // returns 0.0 +* +* @example +* var v = gammaln( 2.0 ); +* // returns 0.0 +* +* @example +* var v = gammaln( 4.0 ); +* // returns ~1.792 +* +* @example +* var v = gammaln( -0.5 ); +* // returns ~1.266 +* +* @example +* var v = gammaln( 0.5 ); +* // returns ~0.572 +* +* @example +* var v = gammaln( 0.0 ); +* // returns Infinity +* +* @example +* var v = gammaln( NaN ); +* // returns NaN +*/ +function gammaln( x ) { + return addon( x ); +} + + +// EXPORTS // + +module.exports = gammaln; diff --git a/base/special/gammaln/lib/polyval_a1.js b/base/special/gammaln/lib/polyval_a1.js index ba74cc7f3..ac2694d47 100644 --- a/base/special/gammaln/lib/polyval_a1.js +++ b/base/special/gammaln/lib/polyval_a1.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_a2.js b/base/special/gammaln/lib/polyval_a2.js index 90a5b00f1..eb2508895 100644 --- a/base/special/gammaln/lib/polyval_a2.js +++ b/base/special/gammaln/lib/polyval_a2.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_r.js b/base/special/gammaln/lib/polyval_r.js index 3767d71b6..26b94f1ee 100644 --- a/base/special/gammaln/lib/polyval_r.js +++ b/base/special/gammaln/lib/polyval_r.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_s.js b/base/special/gammaln/lib/polyval_s.js index fc05ebec1..1df93401c 100644 --- a/base/special/gammaln/lib/polyval_s.js +++ b/base/special/gammaln/lib/polyval_s.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_t1.js b/base/special/gammaln/lib/polyval_t1.js index 530e52a84..ca2c1f4a5 100644 --- a/base/special/gammaln/lib/polyval_t1.js +++ b/base/special/gammaln/lib/polyval_t1.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_t2.js b/base/special/gammaln/lib/polyval_t2.js index 5e4b1d387..b340af8be 100644 --- a/base/special/gammaln/lib/polyval_t2.js +++ b/base/special/gammaln/lib/polyval_t2.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_t3.js b/base/special/gammaln/lib/polyval_t3.js index c1a9b2cd9..002cce134 100644 --- a/base/special/gammaln/lib/polyval_t3.js +++ b/base/special/gammaln/lib/polyval_t3.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_u.js b/base/special/gammaln/lib/polyval_u.js index 667b958a3..3e63b1634 100644 --- a/base/special/gammaln/lib/polyval_u.js +++ b/base/special/gammaln/lib/polyval_u.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_v.js b/base/special/gammaln/lib/polyval_v.js index 6ba3a08cd..96fff3e69 100644 --- a/base/special/gammaln/lib/polyval_v.js +++ b/base/special/gammaln/lib/polyval_v.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/lib/polyval_w.js b/base/special/gammaln/lib/polyval_w.js index b074da9a2..52066128a 100644 --- a/base/special/gammaln/lib/polyval_w.js +++ b/base/special/gammaln/lib/polyval_w.js @@ -1,7 +1,7 @@ /** * @license Apache-2.0 * -* Copyright (c) 2022 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. diff --git a/base/special/gammaln/manifest.json b/base/special/gammaln/manifest.json new file mode 100644 index 000000000..bcc36da00 --- /dev/null +++ b/base/special/gammaln/manifest.json @@ -0,0 +1,93 @@ +{ + "options": { + "task": "build" + }, + "fields": [ + { + "field": "src", + "resolve": true, + "relative": true + }, + { + "field": "include", + "resolve": true, + "relative": true + }, + { + "field": "libraries", + "resolve": false, + "relative": false + }, + { + "field": "libpath", + "resolve": true, + "relative": false + } + ], + "confs": [ + { + "task": "build", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/napi/unary", + "@stdlib/math/base/assert/is-nan", + "@stdlib/math/base/assert/is-infinite", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/trunc", + "@stdlib/math/base/special/sinpi", + "@stdlib/constants/float64/pinf", + "@stdlib/constants/float64/pi" + ] + }, + { + "task": "benchmark", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/assert/is-nan", + "@stdlib/math/base/assert/is-infinite", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/trunc", + "@stdlib/math/base/special/sinpi", + "@stdlib/constants/float64/pinf", + "@stdlib/constants/float64/pi" + ] + }, + { + "task": "examples", + "src": [ + "./src/main.c" + ], + "include": [ + "./include" + ], + "libraries": [], + "libpath": [], + "dependencies": [ + "@stdlib/math/base/assert/is-nan", + "@stdlib/math/base/assert/is-infinite", + "@stdlib/math/base/special/abs", + "@stdlib/math/base/special/ln", + "@stdlib/math/base/special/trunc", + "@stdlib/math/base/special/sinpi", + "@stdlib/constants/float64/pinf", + "@stdlib/constants/float64/pi" + ] + } + ] +} diff --git a/base/special/gammaln/scripts/evalpoly.js b/base/special/gammaln/scripts/evalpoly.js index 127712ad5..dc709bb56 100644 --- a/base/special/gammaln/scripts/evalpoly.js +++ b/base/special/gammaln/scripts/evalpoly.js @@ -24,10 +24,15 @@ // MODULES // var resolve = require( 'path' ).resolve; +var readFileSync = require( '@stdlib/fs/read-file' ).sync; var writeFileSync = require( '@stdlib/fs/write-file' ).sync; var currentYear = require( '@stdlib/time/current-year' ); var licenseHeader = require( '@stdlib/_tools/licenses/header' ); var compile = require( './../../../../base/tools/evalpoly-compile' ); +var compileC = require( './../../../../base/tools/evalpoly-compile-c' ); +var substringBefore = require( '@stdlib/string/substring-before' ); +var substringAfter = require( '@stdlib/string/substring-after' ); +var format = require( '@stdlib/string/format' ); // VARIABLES // @@ -112,6 +117,33 @@ var header = licenseHeader( 'Apache-2.0', 'js', { header += '\n/* This is a generated file. Do not edit directly. */\n'; +// FUNCTIONS // + +/** +* Inserts a compiled function into file content. +* +* @private +* @param {string} text - source content +* @param {string} id - function identifier +* @param {string} str - function string +* @returns {string} updated content +*/ +function insert( text, id, str ) { + var before; + var after; + var begin; + var end; + + begin = '// BEGIN: '+id; + end = '// END: '+id; + + before = substringBefore( text, begin ); + after = substringAfter( text, end ); + + return format( '%s// BEGIN: %s\n\n%s\n%s%s', before, id, str, end, after ); +} + + // MAIN // /** @@ -121,7 +153,9 @@ header += '\n/* This is a generated file. Do not edit directly. */\n'; */ function main() { var fpath; + var copts; var opts; + var file; var str; opts = { @@ -167,6 +201,56 @@ function main() { fpath = resolve( __dirname, '..', 'lib', 'polyval_w.js' ); str = header + compile( W ); writeFileSync( fpath, str, opts ); + + copts = { + 'dtype': 'double', + 'name': '' + }; + + fpath = resolve( __dirname, '..', 'src', 'main.c' ); + file = readFileSync( fpath, opts ); + + copts.name = 'polyval_a1'; + str = compileC( A1, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_a2'; + str = compileC( A2, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_r'; + str = compileC( R, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_s'; + str = compileC( S, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_t1'; + str = compileC( T1, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_t2'; + str = compileC( T2, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_t3'; + str = compileC( T3, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_u'; + str = compileC( U, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_v'; + str = compileC( V, copts ); + file = insert( file, copts.name, str ); + + copts.name = 'polyval_w'; + str = compileC( W, copts ); + file = insert( file, copts.name, str ); + + writeFileSync( fpath, file, opts ); } main(); diff --git a/base/special/gammaln/src/Makefile b/base/special/gammaln/src/Makefile new file mode 100644 index 000000000..bcf18aa46 --- /dev/null +++ b/base/special/gammaln/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/base/special/gammaln/src/addon.c b/base/special/gammaln/src/addon.c new file mode 100644 index 000000000..dc14a2d5f --- /dev/null +++ b/base/special/gammaln/src/addon.c @@ -0,0 +1,23 @@ +/** +* @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/gammaln.h" +#include "stdlib/math/base/napi/unary.h" + +// cppcheck-suppress shadowFunction +STDLIB_MATH_BASE_NAPI_MODULE_D_D( stdlib_base_gammaln ) diff --git a/base/special/gammaln/src/main.c b/base/special/gammaln/src/main.c new file mode 100644 index 000000000..963808e23 --- /dev/null +++ b/base/special/gammaln/src/main.c @@ -0,0 +1,550 @@ +/** +* @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, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/12.2.0/lib/msun/src/e_lgamma_r.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (C) 1993 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/gammaln.h" +#include "stdlib/math/base/assert/is_nan.h" +#include "stdlib/math/base/assert/is_infinite.h" +#include "stdlib/math/base/special/abs.h" +#include "stdlib/math/base/special/ln.h" +#include "stdlib/math/base/special/trunc.h" +#include "stdlib/math/base/special/sinpi.h" +#include "stdlib/constants/float64/pi.h" +#include "stdlib/constants/float64/pinf.h" +#include + +static const double A1C = 7.72156649015328655494e-02; // 0x3FB3C467E37DB0C8 +static const double A2C = 3.22467033424113591611e-01; // 0x3FD4A34CC4A60FAD +static const double RC = 1.0; +static const double SC = -7.72156649015328655494e-02; // 0xBFB3C467E37DB0C8 +static const double T1C = 4.83836122723810047042e-01; // 0x3FDEF72BC8EE38A2 +static const double T2C = -1.47587722994593911752e-01; // 0xBFC2E4278DC6C509 +static const double T3C = 6.46249402391333854778e-02; // 0x3FB08B4294D5419B +static const double UC = -7.72156649015328655494e-02; // 0xBFB3C467E37DB0C8 +static const double VC = 1.0; +static const double WC = 4.18938533204672725052e-01; // 0x3FDACFE390C97D69 +static const double YMIN = 1.461632144968362245; +static const double TWO52 = 4503599627370496; // 2**52 +static const double TWO56 = 72057594037927936; // 2**56 +static const double TINY = 1.3877787807814457e-17; +static const double TC = 1.46163214496836224576e+00; // 0x3FF762D86356BE3F +static const double TF = -1.21486290535849611461e-01; // 0xBFBF19B9BCC38A42 +static const double TT = -3.63867699703950536541e-18; // 0xBC50C7CAA48A971F => TT = -(tail of TF) + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_a1 + +/** +* 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_a1( const double x ) { + return 0.06735230105312927 + (x * (0.007385550860814029 + (x * (0.0011927076318336207 + (x * (0.00022086279071390839 + (x * 0.000025214456545125733))))))); +} + +// END: polyval_a1 + +// BEGIN: polyval_a2 + +/** +* 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_a2( const double x ) { + return 0.020580808432516733 + (x * (0.0028905138367341563 + (x * (0.0005100697921535113 + (x * (0.00010801156724758394 + (x * 0.000044864094961891516))))))); +} + +// END: polyval_a2 + +// BEGIN: polyval_r + +/** +* 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_r( const double x ) { + return 1.3920053346762105 + (x * (0.7219355475671381 + (x * (0.17193386563280308 + (x * (0.01864591917156529 + (x * (0.0007779424963818936 + (x * 0.000007326684307446256))))))))); +} + +// END: polyval_r + +// BEGIN: polyval_s + +/** +* 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_s( const double x ) { + return 0.21498241596060885 + (x * (0.325778796408931 + (x * (0.14635047265246445 + (x * (0.02664227030336386 + (x * (0.0018402845140733772 + (x * 0.00003194753265841009))))))))); +} + +// END: polyval_s + +// BEGIN: polyval_t1 + +/** +* 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_t1( const double x ) { + return -0.032788541075985965 + (x * (0.006100538702462913 + (x * (-0.0014034646998923284 + (x * 0.00031563207090362595))))); +} + +// END: polyval_t1 + +// BEGIN: polyval_t2 + +/** +* 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_t2( const double x ) { + return 0.01797067508118204 + (x * (-0.0036845201678113826 + (x * (0.000881081882437654 + (x * -0.00031275416837512086))))); +} + +// END: polyval_t2 + +// BEGIN: polyval_t3 + +/** +* 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_t3( const double x ) { + return -0.010314224129834144 + (x * (0.0022596478090061247 + (x * (-0.0005385953053567405 + (x * 0.0003355291926355191))))); +} + +// END: polyval_t3 + +// BEGIN: polyval_u + +/** +* 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_u( const double x ) { + return 0.6328270640250934 + (x * (1.4549225013723477 + (x * (0.9777175279633727 + (x * (0.22896372806469245 + (x * 0.013381091853678766))))))); +} + +// END: polyval_u + +// BEGIN: polyval_v + +/** +* 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_v( const double x ) { + return 2.4559779371304113 + (x * (2.128489763798934 + (x * (0.7692851504566728 + (x * (0.10422264559336913 + (x * 0.003217092422824239))))))); +} + +// END: polyval_v + +// 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.08333333333333297 + (x * (-0.0027777777772877554 + (x * (0.0007936505586430196 + (x * (-0.00059518755745034 + (x * (0.0008363399189962821 + (x * -0.0016309293409657527))))))))); +} + +// END: polyval_w + +/* End auto-generated functions. */ + +/** +* Evaluates the natural logarithm of the gamma function. +* +* ## Method +* +* 1. Argument reduction for \\(0 < x \leq 8\\). Since \\(\Gamma(1+s) = s \Gamma(s)\\), for \\(x \in \[0,8]\\), we may reduce \\(x\\) to a number in \\(\[1.5,2.5]\\) by +* +* ```tex +* \operatorname{lgamma}(1+s) = \ln(s) + \operatorname{lgamma}(s) +* ``` +* +* For example, +* +* ```tex +* \begin{align*} +* \operatorname{lgamma}(7.3) &= \ln(6.3) + \operatorname{lgamma}(6.3) \\ +* &= \ln(6.3 \cdot 5.3) + \operatorname{lgamma}(5.3) \\ +* &= \ln(6.3 \cdot 5.3 \cdot 4.3 \cdot 3.3 \cdot2.3) + \operatorname{lgamma}(2.3) +* \end{align*} +* ``` +* +* 2. Compute a polynomial approximation of \\(\mathrm{lgamma}\\) around its minimum (\\(\mathrm{ymin} = 1.461632144968362245\\)) to maintain monotonicity. On the interval \\(\[\mathrm{ymin} - 0.23, \mathrm{ymin} + 0.27]\\) (i.e., \\(\[1.23164,1.73163]\\)), we let \\(z = x - \mathrm{ymin}\\) and use +* +* ```tex +* \operatorname{lgamma}(x) = -1.214862905358496078218 + z^2 \cdot \operatorname{poly}(z) +* ``` +* +* where \\(\operatorname{poly}(z)\\) is a \\(14\\) degree polynomial. +* +* 3. Compute a rational approximation in the primary interval \\(\[2,3]\\). Let \\( s = x - 2.0 \\). We can thus use the approximation +* +* ```tex +* \operatorname{lgamma}(x) = \frac{s}{2} + s\frac{\operatorname{P}(s)}{\operatorname{Q}(s)} +* ``` +* +* with accuracy +* +* ```tex +* \biggl|\frac{\mathrm{P}}{\mathrm{Q}} - \biggr(\operatorname{lgamma}(x)-\frac{s}{2}\biggl)\biggl| < 2^{-61.71} +* ``` +* +* The algorithms are based on the observation +* +* ```tex +* \operatorname{lgamma}(2+s) = s(1 - \gamma) + \frac{\zeta(2) - 1}{2} s^2 - \frac{\zeta(3) - 1}{3} s^3 + \ldots +* ``` +* +* where \\(\zeta\\) is the zeta function and \\(\gamma = 0.5772156649...\\) is the Euler-Mascheroni constant, which is very close to \\(0.5\\). +* +* 4. For \\(x \geq 8\\), +* +* ```tex +* \operatorname{lgamma}(x) \approx \biggl(x-\frac{1}{2}\biggr) \ln(x) - x + \frac{\ln(2\pi)}{2} + \frac{1}{12x} - \frac{1}{360x^3} + \ldots +* ``` +* +* which can be expressed +* +* ```tex +* \operatorname{lgamma}(x) \approx \biggl(x-\frac{1}{2}\biggr)(\ln(x)-1)-\frac{\ln(2\pi)-1}{2} + \ldots +* ``` +* +* Let \\(z = \frac{1}{x}\\). We can then use the approximation +* +* ```tex +* f(z) = \operatorname{lgamma}(x) - \biggl(x-\frac{1}{2}\biggr)(\ln(x)-1) +* ``` +* +* by +* +* ```tex +* w = w_0 + w_1 z + w_2 z^3 + w_3 z^5 + \ldots + w_6 z^{11} +* ``` +* +* where +* +* ```tex +* |w - f(z)| < 2^{-58.74} +* ``` +* +* 5. For negative \\(x\\), since +* +* ```tex +* -x \Gamma(-x) \Gamma(x) = \frac{\pi}{\sin(\pi x)} +* ``` +* +* where \\(\Gamma\\) is the gamma function, we have +* +* ```tex +* \Gamma(x) = \frac{\pi}{\sin(\pi x)(-x)\Gamma(-x)} +* ``` +* +* Since \\(\Gamma(-x)\\) is positive, +* +* ```tex +* \operatorname{sign}(\Gamma(x)) = \operatorname{sign}(\sin(\pi x)) +* ``` +* +* for \\(x < 0\\). Hence, for \\(x < 0\\), +* +* ```tex +* \mathrm{signgam} = \operatorname{sign}(\sin(\pi x)) +* ``` +* +* and +* +* ```tex +* \begin{align*} +* \operatorname{lgamma}(x) &= \ln(|\Gamma(x)|) \\ +* &= \ln\biggl(\frac{\pi}{|x \sin(\pi x)|}\biggr) - \operatorname{lgamma}(-x) +* \end{align*} +* ``` +* +* +* +* Note that one should avoid computing \\(\pi (-x)\\) directly in the computation of \\(\sin(\pi (-x))\\). +* +* +* +* ## Special Cases +* +* ```tex +* \begin{align*} +* \operatorname{lgamma}(2+s) &\approx s (1-\gamma) & \mathrm{for\ tiny\ s} \\ +* \operatorname{lgamma}(x) &\approx -\ln(x) & \mathrm{for\ tiny\ x} \\ +* \operatorname{lgamma}(1) &= 0 & \\ +* \operatorname{lgamma}(2) &= 0 & \\ +* \operatorname{lgamma}(0) &= \infty & \\ +* \operatorname{lgamma}(\infty) &= \infty & \\ +* \operatorname{lgamma}(-\mathrm{integer}) &= \pm \infty +* \end{align*} +* ``` +* +* @param x input value +* @returns function value +* +* @example +* double out = stdlib_base_gammaln( 1.0 ); +* // returns 0.0 +*/ +double stdlib_base_gammaln( const double x ) { + bool isNegative; + int32_t flg; + double nadj; + double xc; + double p3; + double p2; + double p1; + double p; + double q; + double t; + double w; + double y; + double z; + double r; + + // Special cases: NaN, +-infinity + if ( stdlib_base_is_nan( x ) || stdlib_base_is_infinite( x ) ) { + return x; + } + + // Special case: 0 + if ( x == 0.0 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + xc = x; + if ( xc < 0.0 ) { + isNegative = true; + xc = -xc; + } else { + isNegative = false; + } + + // If |x| < 2**-56, return -ln(|x|) + if ( xc < TINY ) { + return -stdlib_base_ln( xc ); + } + if ( isNegative ) { + // If |x| >= 2**52, must be -integer + if ( xc >= TWO52 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + t = stdlib_base_sinpi( xc ); + if ( t == 0.0 ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + nadj = stdlib_base_ln( STDLIB_CONSTANT_FLOAT64_PI / stdlib_base_abs( t * xc ) ); + } + + // If x equals 1 or 2, return 0 + if ( xc == 1.0 || xc == 2.0 ) { + return 0.0; + } + + // If x < 2, use lgamma(x) = lgamma(x+1) - log(x) + if ( xc < 2.0 ) { + if ( xc <= 0.9 ) { + r = -stdlib_base_ln( xc ); + if ( xc >= ( YMIN - 1.0 + 0.27 ) ) { // 0.7316 <= x <= 0.9 + y = 1.0 - xc; + flg = 0; + } else if ( xc >= ( YMIN - 1.0 - 0.27 ) ) { // 0.2316 <= x < 0.7316 + y = xc - ( TC - 1.0 ); + flg = 1; + } else { // 0 < x < 0.2316 + y = xc; + flg = 2; + } + } else { + r = 0.0; + if ( xc >= ( YMIN + 0.27 ) ) { // 1.7316 <= x < 2 + y = 2.0 - xc; + flg = 0; + } else if ( xc >= ( YMIN - 0.27 ) ) { // 1.2316 <= x < 1.7316 + y = xc - TC; + flg = 1; + } else { // 0.9 < x < 1.2316 + y = xc - 1.0; + flg = 2; + } + } + switch ( flg ) { // eslint-disable-line default-case + case 0: + z = y * y; + p1 = A1C + ( z * polyval_a1( z ) ); + p2 = z * ( A2C + ( z * polyval_a2( z ) ) ); + p = ( y * p1 ) + p2; + r += ( p - ( 0.5 * y ) ); + break; + case 1: + z = y * y; + w = z * y; + p1 = T1C + ( w * polyval_t1( w ) ); + p2 = T2C + ( w * polyval_t2( w ) ); + p3 = T3C + ( w * polyval_t3( w ) ); + p = ( z * p1 ) - ( TT - ( w * ( p2 + ( y * p3 ) ) ) ); + r += ( TF + p ); + break; + case 2: + p1 = y * ( UC + ( y * polyval_u( y ) ) ); + p2 = VC + ( y * polyval_v( y ) ); + r += ( -0.5 * y ) + ( p1 / p2 ); + break; + } + } else if ( xc < 8.0 ) { // 2 <= x < 8 + flg = stdlib_base_trunc( xc ); + y = xc - flg; + p = y * ( SC + ( y * polyval_s( y ) ) ); + q = RC + ( y * polyval_r( y ) ); + r = ( 0.5 * y ) + ( p / q ); + z = 1.0; // gammaln(1+s) = ln(s) + gammaln(s) + switch ( flg ) { // eslint-disable-line default-case + case 7: + z *= y + 6.0; + + /* Falls through */ + case 6: + z *= y + 5.0; + + /* Falls through */ + case 5: + z *= y + 4.0; + + /* Falls through */ + case 4: + z *= y + 3.0; + + /* Falls through */ + case 3: + z *= y + 2.0; + r += stdlib_base_ln( z ); + } + } else if ( xc < TWO56 ) { // 8 <= x < 2**56 + t = stdlib_base_ln( xc ); + z = 1.0 / xc; + y = z * z; + w = WC + ( z * polyval_w( y ) ); + r = ( ( xc - 0.5 ) * ( t - 1.0 ) ) + w; + } else { // 2**56 <= x <= Inf + r = xc * ( stdlib_base_ln( xc ) - 1.0 ); + } + if ( isNegative ) { + r = nadj - r; + } + return r; +} diff --git a/base/special/gammaln/test/test.js b/base/special/gammaln/test/test.js index 3fdf1166d..a3714ff9a 100644 --- a/base/special/gammaln/test/test.js +++ b/base/special/gammaln/test/test.js @@ -23,11 +23,11 @@ var tape = require( 'tape' ); var isnan = require( './../../../../base/assert/is-nan' ); var abs = require( './../../../../base/special/abs' ); -var max = require( './../../../../base/special/maxn' ); var ln = require( './../../../../base/special/ln' ); var pow = require( './../../../../base/special/pow' ); var PINF = require( '@stdlib/constants/float64/pinf' ); var NINF = require( '@stdlib/constants/float64/ninf' ); +var EPS = require( '@stdlib/constants/float64/eps' ); var gammaln = require( './../lib' ); @@ -87,31 +87,49 @@ tape( 'the function returns `-ln(x)` for very small x', function test( t ) { }); tape( 'the function evaluates the natural logarithm of the gamma function (positive integers)', function test( t ) { + var expected; var delta; var tol; - var v; + var x; + var y; var i; - for ( i = 0; i < data1.length; i++ ) { - v = gammaln( data1[ i ] ); - delta = abs( v - expected1[ i ] ); - tol = 4.75e-15 * max( 1.0, abs( v ), abs( expected1[ i ] ) ); - t.ok( delta <= tol, 'within tolerance. x: ' + data1[ i ] + '. Value: ' + v + '. Expected: ' + expected1[ i ] + '. Tolerance: ' + tol + '.' ); + x = data1; + expected = expected1; + + for ( i = 0; i < x.length; i++ ) { + y = gammaln( x[ i ] ); + if ( y === expected[ i ] ) { + t.equal( y, expected[ i ], 'x: '+x[i]+'. Expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 21.2 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. y: '+y+'. E: '+expected[i]+'. tol: '+tol+'. Δ: '+delta+'.' ); + } } t.end(); }); tape( 'the function evaluates the natural logarithm of the gamma function (decimal values)', function test( t ) { + var expected; var delta; var tol; - var v; + var x; + var y; var i; - for ( i = 0; i < data2.length; i++ ) { - v = gammaln( data2[ i ] ); - delta = abs( v - expected2[ i ] ); - tol = 2.5e-12 * max( 1.0, abs( v ), abs( expected2[ i ] ) ); - t.ok( delta <= tol, 'within tolerance. x: ' + data2[ i ] + '. Value: ' + v + '. Expected: ' + expected2[ i ] + '. Tolerance: ' + tol + '.' ); + x = data2; + expected = expected2; + + for ( i = 0; i < x.length; i++ ) { + y = gammaln( x[ i ] ); + if ( y === expected[ i ] ) { + t.equal( y, expected[ i ], 'x: '+x[i]+'. Expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 1.1e4 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. y: '+y+'. E: '+expected[i]+'. tol: '+tol+'. Δ: '+delta+'.' ); + } } t.end(); }); diff --git a/base/special/gammaln/test/test.native.js b/base/special/gammaln/test/test.native.js new file mode 100644 index 000000000..498976b62 --- /dev/null +++ b/base/special/gammaln/test/test.native.js @@ -0,0 +1,168 @@ +/** +* @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( './../../../../base/assert/is-nan' ); +var abs = require( './../../../../base/special/abs' ); +var ln = require( './../../../../base/special/ln' ); +var pow = require( './../../../../base/special/pow' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); +var NINF = require( '@stdlib/constants/float64/ninf' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var tryRequire = require( '@stdlib/utils/try-require' ); + + +// FIXTURES // + +var data1 = require( './fixtures/r/data1.json' ); +var expected1 = require( './fixtures/r/expected1.json' ); +var data2 = require( './fixtures/r/data2.json' ); +var expected2 = require( './fixtures/r/expected2.json' ); + + +// VARIABLES // + +var gammaln = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( gammaln instanceof Error ) +}; + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof gammaln, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'if provided `NaN`, the function returns `NaN`', opts, function test( t ) { + var v = gammaln( NaN ); + t.equal( isnan( v ), true, 'returns expected value when provided a NaN' ); + t.end(); +}); + +tape( 'the function returns `infinity` when provided `infinity`', opts, function test( t ) { + var v = gammaln( PINF ); + t.equal( v, PINF, 'returns +Inf when provided +Inf' ); + + v = gammaln( NINF ); + t.equal( v, NINF, 'returns -Inf when provided -Inf' ); + + t.end(); +}); + +tape( 'the function returns `+infinity` when provided `0`', opts, function test( t ) { + var v = gammaln( 0.0 ); + t.equal( v, PINF, 'returns +Inf when provided 0' ); + t.end(); +}); + +tape( 'the function returns `+infinity` for x smaller than `-2^52`', opts, function test( t ) { + var v = gammaln( -pow( 2.0, 53 ) ); + t.equal( v, PINF, 'returns +Inf when provided 2^53' ); + t.end(); +}); + +tape( 'the function returns `-ln(x)` for very small x', opts, function test( t ) { + var x; + var v; + + x = 2.0e-90; + v = gammaln( x ); + t.equal( v, -ln( x ), 'equals -ln(x)' ); + + t.end(); +}); + +tape( 'the function evaluates the natural logarithm of the gamma function (positive integers)', opts, function test( t ) { + var expected; + var delta; + var tol; + var x; + var y; + var i; + + x = data1; + expected = expected1; + + for ( i = 0; i < x.length; i++ ) { + y = gammaln( x[ i ] ); + if ( y === expected[ i ] ) { + t.equal( y, expected[ i ], 'x: '+x[i]+'. Expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 21.2 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. y: '+y+'. E: '+expected[i]+'. tol: '+tol+'. Δ: '+delta+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the natural logarithm of the gamma function (decimal values)', opts, function test( t ) { + var expected; + var delta; + var tol; + var x; + var y; + var i; + + x = data2; + expected = expected2; + + for ( i = 0; i < x.length; i++ ) { + y = gammaln( x[ i ] ); + if ( y === expected[ i ] ) { + t.equal( y, expected[ i ], 'x: '+x[i]+'. Expected: '+expected[i] ); + } else { + delta = abs( y - expected[ i ] ); + tol = 1.1e4 * EPS * abs( expected[ i ] ); + t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. y: '+y+'. E: '+expected[i]+'. tol: '+tol+'. Δ: '+delta+'.' ); + } + } + t.end(); +}); + +tape( 'the function evaluates the natural logarithm of the gamma function for x > 2^58', opts, function test( t ) { + var x; + var v; + + x = pow( 2.0, 59 ); + v = gammaln( x ); + t.equal( v, x * (ln(x)-1.0), 'returns x*(ln(x)-1) for x>2^58' ); + + t.end(); +}); + +tape( 'if provided a positive integer, the function returns the natural logarithm of the factorial of (n-1)', opts, function test( t ) { + t.equal( gammaln( 4.0 ), ln(6.0), 'returns ln(6)' ); + t.equal( gammaln( 5.0 ), ln(24.0), 'returns ln(24)' ); + t.equal( gammaln( 6.0 ), ln(120.0), 'returns ln(120)' ); + t.end(); +}); + +tape( 'returns `+infinity` for `x=-2^51`', opts, function test( t ) { + var v = gammaln( -pow( 2.0, 51 ) ); + t.equal( v, PINF, 'returns +Infinity when provided x=-2^51' ); + t.end(); +});