From 3f43dc14331be1b6638cac664d264f824aae5779 Mon Sep 17 00:00:00 2001 From: aman-095 Date: Sun, 3 Mar 2024 18:00:40 +0530 Subject: [PATCH 01/18] Add C implementation for pow --- .../@stdlib/math/base/special/pow/README.md | 95 +- .../special/pow/benchmark/benchmark.native.js | 62 + .../special/pow/benchmark/c/native/Makefile | 146 +++ .../pow/benchmark/c/native/benchmark.c | 138 +++ .../@stdlib/math/base/special/pow/binding.gyp | 170 +++ .../math/base/special/pow/examples/c/Makefile | 146 +++ .../base/special/pow/examples/c/example.c | 35 + .../math/base/special/pow/include.gypi | 53 + .../include/stdlib/math/base/special/pow.h | 38 + .../math/base/special/pow/lib/native.js | 62 + .../math/base/special/pow/lib/polyval_l.js | 2 +- .../math/base/special/pow/lib/polyval_p.js | 2 +- .../math/base/special/pow/lib/polyval_w.js | 2 +- .../math/base/special/pow/manifest.json | 123 ++ .../math/base/special/pow/scripts/evalpoly.js | 58 +- .../math/base/special/pow/src/Makefile | 70 ++ .../@stdlib/math/base/special/pow/src/addon.c | 23 + .../@stdlib/math/base/special/pow/src/main.c | 702 +++++++++++ .../math/base/special/pow/test/test.native.js | 1030 +++++++++++++++++ 19 files changed, 2952 insertions(+), 5 deletions(-) create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/binding.gyp create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/include.gypi create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/include/stdlib/math/base/special/pow.h create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/lib/native.js create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/manifest.json create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/src/Makefile create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/src/addon.c create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/src/main.c create mode 100644 lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js diff --git a/lib/node_modules/@stdlib/math/base/special/pow/README.md b/lib/node_modules/@stdlib/math/base/special/pow/README.md index 68c771136232..b9794e9b6f7f 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/README.md +++ b/lib/node_modules/@stdlib/math/base/special/pow/README.md @@ -2,7 +2,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. @@ -113,6 +113,99 @@ for ( i = 0; i < 100; i++ ) { + + +* * * + +
+ +## C APIs + + + +
+ +
+ + + + + +
+ +### Usage + +```c +#include "stdlib/math/base/special/pow.h" +``` + +#### stdlib_base_pow( base, exponent ) + +Evaluates the [exponential function][exponential-function]. + +```c +double out = stdlib_base_pow( 3.141592653589793, 5.0 ); +// returns ~306.0197 + +out = stdlib_base_pow( 4.0, 0.5 ); +// returns 2.0 +``` + +The function accepts the following arguments: + +- **b**: `[in] double` input value. +- **x**: `[in] double` input value. + +```c +double stdlib_base_pow( const double b, const double x ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/pow.h" +#include +#include + +int main( void ) { + double b; + double x; + double out; + int i; + + for ( i = 0; i < 100; i++ ) { + b = round( randu()*10.0 ); + x = round( randu()*10.0 ) - 5.0; + out = stdlib_base_pow ( b, x ); + printf( "pow(%lf, %lf) = %lf\n", b, x, v ); + } +} +``` + +
+ + + +
+ + +
diff --git a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js new file mode 100644 index 000000000000..007e3b2fabc4 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js @@ -0,0 +1,62 @@ +/** +* @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( '@stdlib/math/base/assert/is-nan' ); +var tryRequire = require( '@stdlib/utils/try-require' ); +var pkg = require( './../package.json' ).name; + + +// VARIABLES // + +var pow = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( pow instanceof Error ) +}; + + +// MAIN // + +bench( pkg+'::native', opts, function benchmark( b ) { + var v; + var x; + var y; + var i; + + b.tic(); + for ( i = 0; i < b.iterations; i++ ) { + v = round( randu()*10.0 ); + x = round( randu()*10.0 ) - 5.0; + y = pow( v, x ); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + } + b.toc(); + if ( isnan( y ) ) { + b.fail( 'should not return NaN' ); + } + b.pass( 'benchmark finished' ); + b.end(); +}); diff --git a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/Makefile b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/Makefile new file mode 100644 index 000000000000..f69e9da2b4d3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/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/pow/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c new file mode 100644 index 000000000000..64d35e2f73b0 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/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. +*/ + +/** +* Benchmark `ellipk`. +*/ +#include "stdlib/math/base/special/pow.h" +#include +#include +#include +#include +#include + +#define NAME "pow" +#define ITERATIONS 1000000 +#define REPEATS 3 + +/** +* Prints the TAP version. +*/ +void print_version() { + printf( "TAP version 13\n" ); +} + +/** +* Prints the TAP summary. +* +* @param total total number of tests +* @param passing total number of passing tests +*/ +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 +*/ +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 +*/ +double tic() { + 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 +*/ +double rand_double() { + int r = rand(); + return (double)r / ( (double)RAND_MAX + 1.0 ); +} + +/** +* Runs a benchmark. +* +* @return elapsed time in seconds +*/ +double benchmark() { + double elapsed; + double b; + double x; + double y; + double t; + int i; + + t = tic(); + for ( i = 0; i < ITERATIONS; i++ ) { + b = rand_double() * 10.0; + x = rand_double() * 5.0; + y = stdlib_base_pow( b, 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/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp b/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp new file mode 100644 index 000000000000..507cb00291e7 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/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 +} \ No newline at end of file diff --git a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile new file mode 100644 index 000000000000..91d364d19fc3 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2022 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/pow/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c new file mode 100644 index 000000000000..04ab94c1850b --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c @@ -0,0 +1,35 @@ +/** +* @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/pow.h" +#include +#include + +int main( void ) { + double b; + double x; + double out; + int i; + + for ( i = 0; i < 100; i++ ) { + b = round( randu()*10.0 ); + x = round( randu()*10.0 ) - 5.0; + out = stdlib_base_pow ( b, x ); + printf( "pow(%lf, %lf) = %lf\n", b, x, v ); + } +} diff --git a/lib/node_modules/@stdlib/math/base/special/pow/include.gypi b/lib/node_modules/@stdlib/math/base/special/pow/include.gypi new file mode 100644 index 000000000000..c6495fc1da3f --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/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 +#define HIGH_MAX_NEAR_UNITY 0x3fefffff|0 + +#define HUGE 1.0e300 + +#define TINY 1.0e-300 + +// 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 +#define HIGH_MIN_NORMAL_EXP 0x00100000|0 // asm type annotation + +// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 +#define HIGH_BIASED_EXP_NEG_1 0x3fe00000|0 // asm type annotation + +// TODO: consider making into an external constant +#define HIGH_NUM_SIGNIFICAND_BITS 20|0 // asm type annotation + +// High: LN2 +#define LN2_HI 6.93147182464599609375e-01 // 0x3FE62E43, 0x00000000 + +// Low: LN2 +#define LN2_LO -1.90465429995776804525e-09 // 0xBE205C61, 0x0CA86C39 + +// 1/LN2 +#define INV_LN2 1.44269504088896338700e+00 // 0x3FF71547, 0x652B82FE + +// High (24 bits): 1/LN2 +#define INV_LN2_HI 1.44269502162933349609e+00 // 0x3FF71547, 0x60000000 + +// Low: 1/LN2 +#define INV_LN2_LO 1.92596299112661746887e-08 // 0x3E54AE0B, 0xF85DDF44 + +// 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 +#define HIGH_SIGNIFICAND_MASK 0x000fffff|0 // asm type annotation + +// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 +#define HIGH_BIASED_EXP_0 0x3ff00000|0 // asm type annotation + +// 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 +#define HIGH_BIASED_EXP_NEG_512 0x20000000|0 // asm type annotation + +// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 +#define HIGH_SIGNIFICAND_HALF 0x00080000|0 // asm type annotation + +#define TWO53 9007199254740992.0 // 0x43400000, 0x00000000 + +// 2/(3*LN2) +#define CP 9.61796693925975554329e-01 // 0x3FEEC709, 0xDC3A03FD + +// (float)CP +#define CP_HI 9.61796700954437255859e-01 // 0x3FEEC709, 0xE0000000 + +// Low: CP_HI +#define CP_LO -7.02846165095275826516e-09 // 0xBE3E2FE0, 0x145B01F5 + +double BP [2] = { 1.0, 1.5 }; + +double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 + +double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 + +// 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 +#define HIGH_BIASED_EXP_31 0x41e00000|0 // asm type annotation + +// 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 +#define HIGH_BIASED_EXP_64 0x43f00000|0 // asm type annotation + +// 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 +#define HIGH_BIASED_EXP_10 0x40900000|0 // asm type annotation + +// 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 +#define HIGH_1075 0x4090cc00|0 // asm type annotation + +// 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 +#define HIGH_NEG_1075 0xc090cc00>>0 // asm type annotation + +#define HIGH_NUM_NONSIGN_BITS 31|0 // asm type annotation + +// -(1024-log2(ovfl+.5ulp)) +#define OVT = 8.0085662595372944372e-17; + +// High/low words workspace: +double WORDS[2] = { 0|0, 0|0 }; + +// Log workspace: +double LOG_WORKSPACE[2] = { 0.0, 0.0 }; + + +// 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 + + +// MAIN // + +/** +* Evaluates the exponential function when \\( y = \pm \infty\\). +* +* @param x input value +* @param y input value +* @return output value +* +* @example +* double out = y_is_infinite( -1.0, Infinity ); +* // returns NaN +* +* @example +* double out = y_is_infinite( 1.5, -Infinity ); +* // returns 0.0 +*/ +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 input value +* @param y input value +* @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 +*/ +double y_is_huge( const double x, const double y ) { + double ahx; + double hx; + + hx = stdlib_base_float64_get_high_word ( x ); + ahx = ( hx && STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + if ( ahx <= HIGH_MAX_NEAR_UNITY ) { + if ( y < 0 ) { + // Signal overflow... + return HUGE * HUGE; + } + // Signal underflow... + return TINY * TINY; + } + // `x` has a biased exponent greater than or equal to `0`... + + if ( y > 0 ) { + // Signal overflow... + return HUGE * HUGE; + } + // Signal underflow... + return TINY * TINY; +} + +/** +* Evaluates the exponential function when \\(|x| = 0\\). +* +* @param x input value +* @param y input value +* @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 +*/ +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 +*/ +double pow2( const double j, const double hp, const double lp ) { + double tmp; + double t1; + double t; + double r; + double u; + double v; + double w; + double z; + double n; + double i; + double k; + + i = ( j & stdlib_base_high_base_abs_mask )|0; // asm type annotation + k = ( ( i >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS )|0; // asm type annotation + n = 0; + + // `|z| > 0.5`, set `n = z+0.5` + if ( i > HIGH_BIASED_EXP_NEG_1 ) { + n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ) >> 0; // asm type annotation + k = ( ( ( n & stdlib_base_high_base_abs_mask ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS )|0; // new k for n + tmp = ( ( n & !( HIGH_SIGNIFICAND_MASK >> k ) ) ) >> 0; // asm type annotation + t = stdlib_base_float64_set_high_word( 0.0, tmp ); + n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( HIGH_NUM_SIGNIFICAND_BITS - k ) ) >> 0; // eslint-disable-line max-len + if ( j < 0 ) { + n = -n; + } + hp -= t; + } + t = lp + hp; + t = stdlib_base_float64_set_low_word ( t, 0 ); + u = t * LN2_HI; + v = ( ( lp - ( t - hp ) ) * 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 ); + j = stdlib_base_float64_get_high_word ( z ); + j = j|0; + j += ( n << HIGH_NUM_SIGNIFICAND_BITS ) >> 0; // asm type annotation + + // Check for subnormal output... + if ( ( j >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { + z = stdlib_base_ldexp ( z, n ); + } else { + z = stdlib_base_float64_set_high_word ( z, j ); + } + 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 out output array +* @param ax absolute value of x +* @return output array containing a tuple comprised of high and low parts +* +* @example +* double out[2] = logx( [ 0.0, 0.0 ], 9.0 ); // => [ t1, t2 ] +* // returns [ -1265.7236328125, -0.0008163940840404393 ] +*/ +void logx( const double *out, const double ax ) { + 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; + t1 = stdlib_base_float64_set_low_word ( t1, 0 ); + t2 = v - (t1 - u); + + out[ 0 ] = t1; + out[ 1 ] = t2; +} + +/** +* Computes \\(\operatorname{log2}(ax)\\). +* +* @param out output array +* @param ax absolute value of x +* @param ahx high word of ax +* @return output array containing a tuple comprised of high and low parts +* +* @example +* double out[2] = logax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] +* // returns [ 3.169923782348633, 0.0000012190936795504075 ] +*/ +void logax( const double *out, const double ax, const double ahx ) { + double tmp; + 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; + double n; + double j; + double k; + + n = 0|0; // asm type annotation + + // Check if `x` is subnormal... + if ( ahx < HIGH_MIN_NORMAL_EXP ) { + ax *= TWO53; + n -= 53|0; // asm type annotation + ahx = stdlib_base_float64_get_high_word ( ax ); + } + // Extract the unbiased exponent of `x`: + n += ( ( ahx >> HIGH_NUM_SIGNIFICAND_BITS ) - BIAS )|0; // asm type annotation + + // Isolate the significand bits of `x`: + j = ( ahx & HIGH_SIGNIFICAND_MASK )|0; // asm type annotation + + // Normalize `ahx` by setting the (biased) exponent to `1023`: + ahx = ( j|HIGH_BIASED_EXP_0 )|0; // asm type annotation + + // 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|0; // asm type annotation + ahx -= HIGH_MIN_NORMAL_EXP; + } + // Load the normalized high word into `|x|`: + ax = stdlib_base_float64_set_high_word ( ax, ahx ); + + // 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 = ax - bp; // (x-1) || (x-1.5) + v = 1.0 / (ax + bp); // 1/(x+1) || 1/(x+1.5) + ss = u * v; + hs = stdlib_base_float64_set_low_word ( ss, 0 ); // 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 = ( ( ahx >> 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 = stdlib_base_float64_set_high_word ( 0.0, tmp ); + lt = ax - ( 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; + ht = stdlib_base_float64_set_low_word ( ht, 0 ); + lt = r - ( ( ht - 3.0 ) - s2 ); + + // u+v = ss*(1+...): + u = hs * ht; + v = ( ls * ht ) + ( lt * ss ); + + // 2/(3LN2) * (ss+...): + hp = u + v; + hp = stdlib_base_float64_set_low_word ( hp, 0 ); + 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) + t1 = stdlib_base_float64_set_low_word ( t1, 0 ); + t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); + + out[ 0 ] = t1; + out[ 1 ] = t2; +} + +/** +* Evaluates the exponential function. +* +* @param x base +* @param y exponent +* @return function value +* +* @example +* double out = pow( 2.0, 3.0 ); +* // returns 8.0 +* +* @example +* double out = pow( 4.0, 0.5 ); +* // returns 2.0 +*/ +double pow( const double x, const double y ) { + double ahx; // absolute value high word `x` + double ahy; // absolute value high word `y` + double ax; // absolute value `x` + uint32_t hx; // high word `x` + uint32_t lx; // low word `x` + uint32_t hy; // high word `y` + uint32_t ly; // low word `y` + double sx; // sign `x` + double sy; // sign `y` + double y1; + double hp; + double lp; + double t; + double z; // y prime + uint32_t j; + uint32_t i; + if ( stdlib_base_is_nan ( x ) || stdlib_base_is_nan ( y ) ) { + return 0.0/0.0; + } + // 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 x; + } + if ( y == -1.0 ) { + return 1.0 / x; + } + if ( y == 0.5 ) { + return stdlib_base_sqrt( x ); + } + if ( y == -0.5 ) { + return 1.0 / stdlib_base_sqrt( x ); + } + if ( y == 2.0 ) { + return x * x; + } + if ( y == 3.0 ) { + return x * x * x; + } + if ( y == 4.0 ) { + x *= x; + return x * x; + } + if ( stdlib_base_is_infinite( y ) ) { + return y_is_infinite( x, y ); + } + } + // Split `x` into high and low words: + stdlib_base_float64_to_words ( x, &hx, &lx ); + + // Special cases `x`... + if ( lx == 0 ) { + if ( hx == 0 ) { + return x_is_zero( x, y ); + } + if ( x == 1.0 ) { + return 1.0; + } + if ( + x == -1.0 && + stdlib_base_is_odd( y ) + ) { + return -1.0; + } + if ( stdlib_base_is_infinite( x ) ) { + if ( x == STDLIB_CONSTANT_FLOAT64_NINF ) { + // `pow( 1/x, -y )` + return pow( -0.0, -y ); + } + if ( y < 0.0 ) { + return 0.0; + } + return STDLIB_CONSTANT_FLOAT64_PINF; + } + } + if ( + x < 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 & ABS_MASK )|0; // asm type annotation + ahy = ( hy & ABS_MASK )|0; // asm type annotation + + // Extract the sign bits: + sx = ( hx >> HIGH_NUM_NONSIGN_BITS )|0; // asm type annotation + sy = ( hy >> HIGH_NUM_NONSIGN_BITS )|0; // asm type annotation + + // Determine the sign of the result... + if ( sx && stdlib_base_is_odd( y ) ) { + sx = -1.0; + } else { + sx = 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( x, 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 * HUGE; + } + // Signal underflow... + return sx * TINY * TINY; + } + if ( ahx > HIGH_BIASED_EXP_0 ) { + // y > 0 + if ( sy == 0 ) { + // Signal overflow... + return sx * HUGE * HUGE; + } + // Signal underflow... + return sx * TINY * TINY; + } + // 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`. + t = logx ( LOG_WORKSPACE, ax ); + } + // Case 2: `|y|` is not huge... + else { + t = log2ax ( LOG_WORKSPACE, ax, ahx ); + } + // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... + y1 = stdlib_base_float64_set_low_word ( y, 0 ); + lp = ( ( y - y1 ) * t[0] ) + ( y * t[1] ); + hp = y1 * t[0]; + z = lp + hp; + + // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner. + toWords.assign( z, &j, &i ); + j = j|0; + i = i|0; + + // z >= 1024 + if ( j >= HIGH_BIASED_EXP_10 ) { + // z > 1024 + if ( ( ( j - HIGH_BIASED_EXP_10 )|i ) != 0 ) { + // Signal overflow... + return sx * HUGE * HUGE; + } + if ( ( lp + OVT ) > ( z - hp ) ) { + // Signal overflow... + return sx * HUGE * HUGE; + } + } + // z <= -1075 + else if ( ( j & stdlib_base_high_base_abs_mask ) >= HIGH_1075 ) { + // z < -1075 + if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { + // Signal underflow... + return sx * TINY * TINY; + } + if ( lp <= ( z - hp ) ) { + // Signal underflow... + return sx * TINY * TINY; + } + } + // Compute `2^(hp+lp)`... + z = pow2( j, hp, lp ); + + return sx * z; +} \ No newline at end of file diff --git a/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js new file mode 100644 index 000000000000..639a00b35600 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js @@ -0,0 +1,1030 @@ +/** +* @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 abs = require( '@stdlib/math/base/special/abs' ); +var sqrt = require( '@stdlib/math/base/special/sqrt' ); +var randu = require( '@stdlib/random/base/randu' ); +var round = require( '@stdlib/math/base/special/round' ); +var PINF = require( '@stdlib/constants/float64/pinf' ); +var NINF = require( '@stdlib/constants/float64/ninf' ); +var EPS = require( '@stdlib/constants/float64/eps' ); +var MAX_SAFE_INTEGER = require( '@stdlib/constants/float64/max-safe-integer' ); +var isnan = require( '@stdlib/math/base/assert/is-nan' ); +var isPositiveZero = require( '@stdlib/math/base/assert/is-positive-zero' ); +var isNegativeZero = require( '@stdlib/math/base/assert/is-negative-zero' ); +var tryRequire = require( '@stdlib/utils/try-require' ); + + +// FIXTURES // + +var squaredSmall = require( './fixtures/julia/squared_small.json' ); +var squaredLarge = require( './fixtures/julia/squared_large.json' ); +var cubedSmall = require( './fixtures/julia/cubed_small.json' ); +var cubedLarge = require( './fixtures/julia/cubed_large.json' ); +var toTheFourthSmall = require( './fixtures/julia/to_the_fourth_small.json' ); +var toTheFourthLarge = require( './fixtures/julia/to_the_fourth_large.json' ); +var baseNearUnitySmall = require( './fixtures/julia/base_near_unity_small.json' ); +var baseNearUnityLarge = require( './fixtures/julia/base_near_unity_large.json' ); +var baseNearUnityHuge = require( './fixtures/julia/base_near_unity_huge.json' ); +var subnormalResults = require( './fixtures/julia/subnormal_results.json' ); +var negativeExpSmallSmall = require( './fixtures/julia/negative_exp_small_small.json' ); +var negativeExpSmallLarge = require( './fixtures/julia/negative_exp_small_large.json' ); +var negativeExpLargeSmall = require( './fixtures/julia/negative_exp_large_small.json' ); +var negativeExpLargeLarge = require( './fixtures/julia/negative_exp_large_large.json' ); +var smallSmall = require( './fixtures/julia/small_small.json' ); +var smallLarge = require( './fixtures/julia/small_large.json' ); +var largeSmall = require( './fixtures/julia/large_small.json' ); +var decimalDecimal = require( './fixtures/julia/decimal_decimal.json' ); +var decimalInteger = require( './fixtures/julia/decimal_integer.json' ); +var integerDecimal = require( './fixtures/julia/integer_decimal.json' ); +var integerInteger = require( './fixtures/julia/integer_integer.json' ); + + +// VARIABLES // + +var pow = tryRequire( resolve( __dirname, './../lib/native.js' ) ); +var opts = { + 'skip': ( pow instanceof Error ) +}; + + +// TESTS // + +tape( 'main export is a function', opts, function test( t ) { + t.ok( true, __filename ); + t.strictEqual( typeof pow, 'function', 'main export is a function' ); + t.end(); +}); + +tape( 'the function accepts two parameters: a base and an exponent', opts, function test( t ) { + t.equal( pow.length, 2.0, 'arity is 2' ); + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 2` and `x` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = squaredSmall.x; + y = squaredSmall.y; + expected = squaredSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 2` and `x` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = squaredLarge.x; + y = squaredLarge.y; + expected = squaredLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 3` and `x` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = cubedSmall.x; + y = cubedSmall.y; + expected = cubedSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 3` and `x` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = cubedLarge.x; + y = cubedLarge.y; + expected = cubedLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 4` and `x` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = toTheFourthSmall.x; + y = toTheFourthSmall.y; + expected = toTheFourthSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`y = 4` and `x` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = toTheFourthLarge.x; + y = toTheFourthLarge.y; + expected = toTheFourthLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`x ~ 1`, `y` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = baseNearUnitySmall.x; + y = baseNearUnitySmall.y; + expected = baseNearUnitySmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`x ~ 1`, `y` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = baseNearUnityLarge.x; + y = baseNearUnityLarge.y; + expected = baseNearUnityLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (`x ~ 1`, `y` huge)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = baseNearUnityHuge.x; + y = baseNearUnityHuge.y; + expected = baseNearUnityHuge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (subnormal results)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = subnormalResults.x; + y = subnormalResults.y; + expected = subnormalResults.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (negative exponents; `x` small, `y` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = negativeExpSmallSmall.x; + y = negativeExpSmallSmall.y; + expected = negativeExpSmallSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + if ( expected[i] === null ) { + t.equal( actual, PINF, 'pow('+x[i]+','+y[i]+') returns +infinity' ); + } else { + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (negative exponents; `x` small, `y` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = negativeExpSmallLarge.x; + y = negativeExpSmallLarge.y; + expected = negativeExpSmallLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + if ( expected[i] === null ) { + t.equal( actual, PINF, 'pow('+x[i]+','+y[i]+') returns +infinity' ); + } else { + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (negative exponents; `x` large, `y` small)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = negativeExpLargeSmall.x; + y = negativeExpLargeSmall.y; + expected = negativeExpLargeSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + if ( expected[i] === null ) { + t.equal( actual, PINF, 'pow('+x[i]+','+y[i]+') returns +infinity' ); + } + else if ( expected[i] === 5.0e-324 ) { + // Julia, Python, and R occasionally return the min subnormal value, when, in fact, the result should underflow; e.g., `pow(571413458040000, -21.92803603935001)` returns `0` for this implementation, and `5.0e-324` on other platforms. The answer checked on Wolfram Alpha is `2.564828778186e-324`, which is less than the min subnormal. This discrepancy arises, at least when compared to Julia's Openlibm, due to differences in our `ldexp` and their `scalbn`. + t.equal( actual === expected[i] || actual === 0.0, true, 'pow('+x[i]+','+y[i]+') returns '+expected[i]+' or 0' ); + } + else { + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (negative exponents; `x` large, `y` large)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = negativeExpLargeLarge.x; + y = negativeExpLargeLarge.y; + expected = negativeExpLargeLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + if ( expected[i] === null ) { + t.equal( actual, PINF, 'pow('+x[i]+','+y[i]+') returns +infinity' ); + } else { + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (near overflow)', opts, function test( t ) { + var x; + var y; + var v; + + x = 2.0; + y = 1024; + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + x = 2.0; + y = 1023; + v = pow( x, y ); + t.equal( v, 8.98846567431158e307, 'pow('+x+','+y+') returns 8.98846567431158e307' ); + + x = 2.001355585648339; + y = 1023; + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + x = 2.0; + y = 1025; + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + x = 1.998647974895499; + y = 1025; + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + x = 2.0 + ( EPS*2.0 ); + y = 1024; + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + t.end(); +}); + +tape( 'the function evaluates the exponential function (near underflow)', opts, function test( t ) { + var x; + var y; + var v; + + x = 2.0; + y = -1075; + v = pow( x, y ); + t.equal( v, 0.0, 'pow('+x+','+y+') returns 0' ); + + x = 1.9987108395124398; + y = -1075; + v = pow( x, y ); + t.equal( v, 5.0e-324, 'pow('+x+','+y+') returns 5e-324' ); + + // Wolfram Alpha: 2.47032822920649828116748364864574508116622553173015e-324 + x = 1.9999999999999998; + y = -1075; + v = pow( x, y ); + t.equal( v, 0.0, 'pow('+x+','+y+') returns 0' ); + + x = 1.9999999999999999; + y = -1075; + v = pow( x, y ); + t.equal( v, 0.0, 'pow('+x+','+y+') returns 0' ); + + t.end(); +}); + +tape( 'the function evaluates the exponential function (small `x`, large `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = smallLarge.x; + y = smallLarge.y; + expected = smallLarge.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + if ( expected[i] === 5.0e-324 ) { + t.equal( actual === expected[i] || actual === 0.0, true, 'pow('+x[i]+','+y[i]+') returns 5e-324 or 0' ); + } else { + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (large `x`, small `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = largeSmall.x; + y = largeSmall.y; + expected = largeSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (small `x`, small `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = smallSmall.x; + y = smallSmall.y; + expected = smallSmall.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (decimal `x`, decimal `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = decimalDecimal.x; + y = decimalDecimal.y; + expected = decimalDecimal.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (decimal `x`, integer `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = decimalInteger.x; + y = decimalInteger.y; + expected = decimalInteger.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (integer `x`, decimal `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = integerDecimal.x; + y = integerDecimal.y; + expected = integerDecimal.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function evaluates the exponential function (integer `x`, integer `y`)', opts, function test( t ) { + var expected; + var actual; + var x; + var y; + var i; + + x = integerInteger.x; + y = integerInteger.y; + expected = integerInteger.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[i] ); + t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } + t.end(); +}); + +tape( 'the function returns `NaN` if provided `NaN` for the exponent', opts, function test( t ) { + var v; + + v = pow( -3.0, NaN ); + t.equal( isnan( v ), true, 'returns NaN' ); + + v = pow( 0.0, NaN ); + t.equal( isnan( v ), true, 'returns NaN' ); + + t.end(); +}); + +tape( 'the function returns `NaN` if provided `NaN` for the base', opts, function test( t ) { + var v; + + v = pow( NaN, 5.0 ); + t.equal( isnan( v ), true, 'returns NaN' ); + + v = pow( NaN, 1.0 ); + t.equal( isnan( v ), true, 'returns NaN' ); + + v = pow( NaN, 0.0 ); + t.equal( isnan( v ), true, 'returns NaN' ); + + t.end(); +}); + +tape( 'any number raised to the zero power is `1`', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = ( randu()*1.0e6 ) - 5.0e5; + v = pow( x, 0.0 ); + t.equal( v, 1.0, 'pow('+x+',0.0) returns 1' ); + } + t.equal( pow( PINF, 0.0 ), 1.0, 'inf^0 = 1' ); + t.equal( pow( NINF, 0.0 ), 1.0, '(-inf)^0 = 1' ); + t.equal( pow( 0.0, 0.0 ), 1.0, '0^0 = 1' ); + t.end(); +}); + +tape( 'any number raised to the `-1` power is the reciprocal of that number', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = ( randu()*1.0e6 ) - 5.0e5; + v = pow( x, -1.0 ); + t.equal( v, 1/x, 'pow('+x+',0.0) returns '+(1/x) ); + } + t.equal( isPositiveZero( pow( PINF, -1.0 ) ), true, 'inf^-1 = 0' ); + t.equal( isNegativeZero( pow( NINF, -1.0 ) ), true, '(-inf)^0 = -0' ); + t.equal( pow( 0.0, -1.0 ), PINF, '0^-1 = +infinity' ); + t.equal( pow( -0.0, -1.0 ), NINF, '(-0)^-1 = -infinity' ); + t.end(); +}); + +tape( 'unity raised to any power is `1`', opts, function test( t ) { + var y; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + y = ( randu()*1.0e6 ) - 5.0e5; + v = pow( 1.0, y ); + t.equal( v, 1.0, 'pow(1.0,'+y+') returns 1' ); + } + t.equal( pow( 1.0, PINF ), 1.0, '1^inf = 1' ); + t.equal( pow( 1.0, NINF ), 1.0, '1^(-inf) = 1' ); + t.equal( pow( 1.0, 0.0 ), 1.0, '1^0 = 1' ); + t.end(); +}); + +tape( '`-1` raised to any odd integer is `-1`', opts, function test( t ) { + var v; + var i; + for ( i = -51; i < 53; i += 2 ) { + v = pow( -1.0, i ); + t.equal( v, -1.0, 'pow(1.0,'+i+') returns -1' ); + } + t.end(); +}); + +tape( '`+-0` raised to `-infinity` is `+infinity`', opts, function test( t ) { + var v; + + v = pow( +0.0, NINF ); + t.equal( v, PINF, 'returns +infinity' ); + + v = pow( -0.0, NINF ); + t.equal( v, PINF, 'returns +infinity' ); + + t.end(); +}); + +tape( '`+-0` raised to `+infinity` is `0`', opts, function test( t ) { + var v; + + v = pow( +0.0, PINF ); + t.equal( isPositiveZero( v ), true, 'returns +0' ); + + v = pow( -0.0, PINF ); + t.equal( isPositiveZero( v ), true, 'returns +0' ); + + t.end(); +}); + +tape( '`+-0` raised to any positive finite number which is not an odd integer is `0`', opts, function test( t ) { + var y; + var i; + var v; + + y = [ + 2, + 3.14, + 1.0e5, + 5.1, + 9.999999 + ]; + + for ( i = 0; i < y.length; i++ ) { + v = pow( +0.0, y[ i ] ); + t.equal( isPositiveZero( v ), true, 'pow(0,'+y[i]+') returns +0' ); + + v = pow( -0.0, y[ i ] ); + t.equal( isPositiveZero( v ), true, 'pow(-0,'+y[i]+') returns +0' ); + } + t.end(); +}); + +tape( '`+-0` raised to any positive odd integer is `+-0`', opts, function test( t ) { + var i; + var v; + + for ( i = 1; i < 101; i += 2 ) { + v = pow( +0.0, i ); + t.equal( isPositiveZero( v ), true, 'pow(0,'+i+') returns +0' ); + + v = pow( -0.0, i ); + t.equal( isNegativeZero( v ), true, 'pow(-0,'+i+') returns -0' ); + } + t.end(); +}); + +tape( '`+-0` raised to any negative finite number which is not an odd integer is `+infinity`', opts, function test( t ) { + var y; + var i; + var v; + + y = [ + -2, + -3.14, + -1.0e5, + -5.1, + -9.999999 + ]; + + for ( i = 0; i < y.length; i++ ) { + v = pow( +0.0, y[ i ] ); + t.equal( v, PINF, 'pow(0,'+y[i]+') returns +infinity' ); + + v = pow( -0.0, y[ i ] ); + t.equal( v, PINF, 'pow(-0,'+y[i]+') returns +infinity' ); + } + t.end(); +}); + +tape( '`+-0` raised to any negative odd integer is `+-infinity`', opts, function test( t ) { + var i; + var v; + + for ( i = -101; i < 0; i += 2 ) { + v = pow( +0.0, i ); + t.equal( v, PINF, 'pow(0,'+i+') returns +infinity' ); + + v = pow( -0.0, i ); + t.equal( v, NINF, 'pow(-0,'+i+') returns -infinity' ); + } + t.end(); +}); + +tape( '`-1` raised to `+-infinity` is indeterminate and the function returns `NaN`', opts, function test( t ) { + var v; + + v = pow( -1.0, PINF ); + t.equal( isnan( v ), true, 'returns NaN' ); + + v = pow( -1.0, NINF ); + t.equal( isnan( v ), true, 'returns NaN' ); + + t.end(); +}); + +tape( 'raising any number to the `1/2` power returns the same results as evaluating the principal square root', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 500; i++ ) { + x = randu()*1.0e5; + v = pow( x, 0.5 ); + t.equal( v, sqrt( x ), 'returns sqrt('+x+')' ); + } + t.end(); +}); + +tape( 'raising any number to the `-1/2` power returns the same results as evaluating the inverse of the principal square root', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 500; i++ ) { + x = randu()*1.0e5; + v = pow( x, -0.5 ); + t.equal( v, 1.0/sqrt( x ), 'returns sqrt('+x+')' ); + } + t.end(); +}); + +tape( 'the function returns `+infinity` if a positive number greater than `1` is raised to `+infinity`', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = 1.0 + ( randu()*1.0e5 ); + v = pow( x, PINF ); + t.equal( v, PINF, 'pow('+x+',+infinity) returns +infinity' ); + } + t.end(); +}); + +tape( 'the function returns `+infinity` if a negative number less than `-1` is raised to `+infinity`', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = ( -randu()*1.0e5 ) - 1.0; + v = pow( x, PINF ); + t.equal( v, PINF, 'pow('+x+',+infinity) returns +infinity' ); + } + t.end(); +}); + +tape( 'the function returns `+infinity` if `-infinity` is raised to `+infinity`', opts, function test( t ) { + var v = pow( NINF, PINF ); + t.equal( v, PINF, 'returns +infinity' ); + t.end(); +}); + +tape( 'the function returns `(-0)^y` if `-infinity` is raised to a numeric `y` power', opts, function test( t ) { + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + y = ( randu()*100.0 ) - 50.0; + v = pow( NINF, y ); + t.equal( v, pow( -0.0, -y ), 'pow(-infinity,'+y+') returns pow(-0,-'+y+')' ); + } + v = pow( NINF, NINF ); + t.equal( v, pow( -0.0, PINF ), 'pow(-infinity,-infinity) returns pow(-0,+infinity)' ); + + v = pow( NINF, PINF ); + t.equal( v, pow( -0.0, NINF ), 'pow(-infinity,+infinity) returns pow(-0,-infinity)' ); + + v = pow( NINF, 5.0 ); + t.equal( v, pow( -0.0, -5.0 ), 'pow(-infinity,5) returns pow(-0,-5)' ); + + v = pow( NINF, 2.0 ); + t.equal( v, pow( -0.0, -2.0 ), 'pow(-infinity,2) returns pow(-0,-2)' ); + + v = pow( NINF, -5.0 ); + t.equal( v, pow( -0.0, 5.0 ), 'pow(-infinity,-5) returns pow(-0,5)' ); + + v = pow( NINF, -2.0 ); + t.equal( v, pow( -0.0, 2.0 ), 'pow(-infinity,-2) returns pow(-0,2)' ); + + t.end(); +}); + +tape( 'the function returns `0` if a positive number greater than `1` is raised to `-infinity`', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = 1.0 + ( randu()*1.0e5 ); + v = pow( x, NINF ); + t.equal( v, 0.0, 'pow('+x+',-infinity) returns 0' ); + } + t.end(); +}); + +tape( 'the function returns `0` if `+infinity` is raised to `-infinity`', opts, function test( t ) { + var v = pow( PINF, NINF ); + t.equal( v, 0.0, 'returns 0' ); + t.end(); +}); + +tape( 'the function returns `0` if `+infinity` is raised to any negative finite number', opts, function test( t ) { + var y; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + y = -(randu() * 1.0e5) - EPS; + v = pow( PINF, y ); + t.equal( isPositiveZero( v ), true, 'returns 0' ); + } + t.end(); +}); + +tape( 'the function returns `+infinity` if `+infinity` is raised to any positive finite number', opts, function test( t ) { + var y; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + y = ( randu()*1.0e5 ) + EPS; + v = pow( PINF, y ); + t.equal( v, PINF, 'returns +infinity' ); + } + t.end(); +}); + +tape( 'the function returns `0` if a negative number less than `-1` is raised to `-infinity`', opts, function test( t ) { + var x; + var v; + var i; + for ( i = 0; i < 100; i++ ) { + x = ( -randu()*1.0e5 ) - 1.0; + v = pow( x, NINF ); + t.equal( v, 0.0, 'pow('+x+',-infinity) returns 0' ); + } + t.end(); +}); + +tape( 'the function returns `NaN` if asked to raise a negative number to an exponent whose absolute value is greater than 0 but less than 1', opts, function test( t ) { + var x; + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + x = ( -randu()*1.0e5 ) - EPS; + y = randu(); + if ( y === 0.0 || y === 1.0 ) { + continue; + } + v = pow( x, y ); + t.equal( isnan( v ), true, 'pow('+x+','+y+') returns NaN' ); + } + t.end(); +}); + +tape( 'if `y > 2^64` and `x` is greater than `1`, the function overflows', opts, function test( t ) { + var x; + var y; + var v; + + x = 1.0 + EPS; // smallest value greater than 1.0 + y = 3.6893488147419103e19; // 2.0^65 + v = pow( x, y ); + t.equal( v, PINF, 'returns +infinity' ); + t.end(); +}); + +tape( 'if `y < -(2^64)` and `x` is greater than `1`, the function underflows', opts, function test( t ) { + var x; + var y; + var v; + + x = 1.0 + EPS; // smallest value greater than 1.0 + y = -3.6893488147419103e19; // -(2.0)^65 + v = pow( x, y ); + t.equal( v, 0.0, 'returns 0' ); + t.end(); +}); + +tape( 'if `y > 2^64` and `x` is less than `1`, the function underflows', opts, function test( t ) { + var x; + var y; + var v; + + x = 1.0 - EPS; // greatest value less than 1.0 + y = 3.6893488147419103e19; // 2.0^65 + v = pow( x, y ); + t.equal( v, 0.0, 'returns 0' ); + t.end(); +}); + +tape( 'if `y < -(2^64)` and `x` is less than `1`, the function overflows', opts, function test( t ) { + var x; + var y; + var v; + + x = 1.0 - EPS; // greatest value less than 1.0 + y = -3.6893488147419103e19; // -(2.0^65) + v = pow( x, y ); + t.equal( v, PINF, 'returns +infinity' ); + t.end(); +}); + +tape( 'if `y > 2^31` and `|x|` significantly less than unity, the function underflows', opts, function test( t ) { + var exp; + var x; + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + x = 0.5 * randu(); + + exp = round( randu()*32.0 ) + 32; + y = pow( 2.0, exp ); + + v = pow( x, y ); + t.equal( isPositiveZero( v ), true, 'pow('+x+','+y+') returns +0' ); + + if ( y < MAX_SAFE_INTEGER ) { + y += 1; + v = pow( -x, y ); + t.equal( isNegativeZero( v ), true, 'pow(-'+x+','+y+') returns -0' ); + } + } + t.end(); +}); + +tape( 'if `y < -(2^31)` and `|x|` significantly less than unity, the function overflows', opts, function test( t ) { + var exp; + var x; + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + x = 0.5 * randu(); + + exp = round( randu()*32.0 ) + 32; + y = -pow( 2.0, exp ); + + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + if ( abs(y) < MAX_SAFE_INTEGER ) { + y -= 1; + v = pow( -x, y ); + t.equal( v, NINF, 'pow(-'+x+','+y+') returns -infinity' ); + } + } + t.end(); +}); + +tape( 'if `y < -(2^31)` and `|x|` significantly greater than unity, the function underflows', opts, function test( t ) { + var exp; + var x; + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + x = 2.0 + ( randu()*1.0e5 ); + + exp = round( randu()*32.0 ) + 32; + y = -pow( 2.0, exp ); + + v = pow( x, y ); + t.equal( isPositiveZero( v ), true, 'pow('+x+','+y+') returns +0' ); + + if ( abs(y) < MAX_SAFE_INTEGER ) { + y -= 1; + v = pow( -x, y ); + t.equal( isNegativeZero( v ), true, 'pow(-'+x+','+y+') returns -0' ); + } + } + t.end(); +}); + +tape( 'if `y > 2^31` and `|x|` significantly greater than unity, the function overflows', opts, function test( t ) { + var exp; + var x; + var y; + var v; + var i; + + for ( i = 0; i < 100; i++ ) { + x = 2.0 + ( randu()*1.0e5 ); + + exp = round( randu()*32.0 ) + 32; + y = pow( 2.0, exp ); + + v = pow( x, y ); + t.equal( v, PINF, 'pow('+x+','+y+') returns +infinity' ); + + if ( y < MAX_SAFE_INTEGER ) { + y += 1; + v = pow( -x, y ); + t.equal( v, NINF, 'pow(-'+x+','+y+') returns -infinity' ); + } + } + t.end(); +}); From 3a78f0c811334174d1831320036f4cdd93b818d6 Mon Sep 17 00:00:00 2001 From: aman-095 Date: Wed, 6 Mar 2024 21:55:33 +0530 Subject: [PATCH 02/18] chore: changes --- .../@stdlib/math/base/special/pow/README.md | 2 +- .../@stdlib/math/base/special/pow/src/main.c | 158 +++++++++--------- 2 files changed, 81 insertions(+), 79 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/README.md b/lib/node_modules/@stdlib/math/base/special/pow/README.md index b9794e9b6f7f..c4174cdff03b 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/README.md +++ b/lib/node_modules/@stdlib/math/base/special/pow/README.md @@ -141,7 +141,7 @@ for ( i = 0; i < 100; i++ ) { #### stdlib_base_pow( base, exponent ) -Evaluates the [exponential function][exponential-function]. +Evaluates the exponential function. ```c double out = stdlib_base_pow( 3.141592653589793, 5.0 ); diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index f7a5b7556709..0538155125ac 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -38,62 +38,59 @@ - -// VARIABLES // - // 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1 -#define HIGH_MAX_NEAR_UNITY 0x3fefffff|0 +static const uint32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; -#define HUGE 1.0e300 +static const double HUGE = 1.0e300; -#define TINY 1.0e-300 +static const double TINY = 1.0e-300; // 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 -#define HIGH_MIN_NORMAL_EXP 0x00100000|0 // asm type annotation +static const uint32_t HIGH_MIN_NORMAL_EXP = 0x00100000; // 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 -#define HIGH_BIASED_EXP_NEG_1 0x3fe00000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; // TODO: consider making into an external constant -#define HIGH_NUM_SIGNIFICAND_BITS 20|0 // asm type annotation +static const uint32_t HIGH_NUM_SIGNIFICAND_BITS = 20; // High: LN2 -#define LN2_HI 6.93147182464599609375e-01 // 0x3FE62E43, 0x00000000 +static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 // Low: LN2 -#define LN2_LO -1.90465429995776804525e-09 // 0xBE205C61, 0x0CA86C39 +static const double LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39 // 1/LN2 -#define INV_LN2 1.44269504088896338700e+00 // 0x3FF71547, 0x652B82FE +static const double INV_LN2 = 1.44269504088896338700e+00; // 0x3FF71547, 0x652B82FE // High (24 bits): 1/LN2 -#define INV_LN2_HI 1.44269502162933349609e+00 // 0x3FF71547, 0x60000000 +static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60000000 // Low: 1/LN2 -#define INV_LN2_LO 1.92596299112661746887e-08 // 0x3E54AE0B, 0xF85DDF44 +static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 // 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 -#define HIGH_SIGNIFICAND_MASK 0x000fffff|0 // asm type annotation +static const uint32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; // 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 -#define HIGH_BIASED_EXP_0 0x3ff00000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_0 = 0x3ff00000; // 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 -#define HIGH_BIASED_EXP_NEG_512 0x20000000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; // 0x00080000 = 524288 => 0 00000000000 10000000000000000000 -#define HIGH_SIGNIFICAND_HALF 0x00080000|0 // asm type annotation +static const uint32_t HIGH_SIGNIFICAND_HALF = 0x00080000; -#define TWO53 9007199254740992.0 // 0x43400000, 0x00000000 +static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 // 2/(3*LN2) -#define CP 9.61796693925975554329e-01 // 0x3FEEC709, 0xDC3A03FD +static const double CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD // (float)CP -#define CP_HI 9.61796700954437255859e-01 // 0x3FEEC709, 0xE0000000 +static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 // Low: CP_HI -#define CP_LO -7.02846165095275826516e-09 // 0xBE3E2FE0, 0x145B01F5 +static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 double BP [2] = { 1.0, 1.5 }; @@ -102,31 +99,32 @@ double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x4000000 double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 // 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 -#define HIGH_BIASED_EXP_31 0x41e00000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_31 = 0x41e00000; // 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 -#define HIGH_BIASED_EXP_64 0x43f00000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_64 = 0x43f00000; // 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 -#define HIGH_BIASED_EXP_10 0x40900000|0 // asm type annotation +static const uint32_t HIGH_BIASED_EXP_10 = 0x40900000; // 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 -#define HIGH_1075 0x4090cc00|0 // asm type annotation +static const uint32_t HIGH_1075 = 0x4090cc00; // 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 -#define HIGH_NEG_1075 0xc090cc00>>0 // asm type annotation +static const uint32_t HIGH_NEG_1075 = 0xc090cc00; -#define HIGH_NUM_NONSIGN_BITS 31|0 // asm type annotation +static const uint32_t HIGH_NUM_NONSIGN_BITS = 31; // -(1024-log2(ovfl+.5ulp)) -#define OVT = 8.0085662595372944372e-17; +static const double OVT = 8.0085662595372944372e-17; // High/low words workspace: -double WORDS[2] = { 0|0, 0|0 }; +uint32_t WORDS[2] = { 0, 0 }; // Log workspace: double LOG_WORKSPACE[2] = { 0.0, 0.0 }; +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ // BEGIN: polyval_l @@ -146,7 +144,9 @@ 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 +// END: polyval_l + +// BEGIN: polyval_w /** * Evaluates a polynomial. @@ -164,7 +164,9 @@ static double polyval_w( const double x ) { return 0.5 + (x * (-0.3333333333333333 + (x * 0.25))); } -// END: polyval_w// BEGIN: polyval_p +// END: polyval_w + +// BEGIN: polyval_p /** * Evaluates a polynomial. @@ -184,8 +186,7 @@ static double polyval_p( const double x ) { // END: polyval_p - -// MAIN // +/* End auto-generated functions. */ /** * Evaluates the exponential function when \\( y = \pm \infty\\). @@ -233,10 +234,12 @@ double y_is_infinite( const double x, const double y ) { * // returns 0.0 */ double y_is_huge( const double x, const double y ) { - double ahx; - double hx; + uint32_t hx; + double ahx; + double xc; - hx = stdlib_base_float64_get_high_word ( x ); + xc = x; + stdlib_base_float64_get_high_word ( xc, &hx ); ahx = ( hx && STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); if ( ahx <= HIGH_MAX_NEAR_UNITY ) { @@ -245,7 +248,7 @@ double y_is_huge( const double x, const double y ) { return HUGE * HUGE; } // Signal underflow... - return TINY * TINY; + return TINY * TINY; } // `x` has a biased exponent greater than or equal to `0`... @@ -254,7 +257,7 @@ double y_is_huge( const double x, const double y ) { return HUGE * HUGE; } // Signal underflow... - return TINY * TINY; + return TINY * TINY; } /** @@ -304,7 +307,7 @@ double x_is_zero( const double x, const double y ) { * double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); * // returns ~0.79 */ -double pow2( const double j, const double hp, const double lp ) { +double pow2( const uint32_t j, const double hp, const double lp ) { double tmp; double t1; double t; @@ -313,18 +316,18 @@ double pow2( const double j, const double hp, const double lp ) { double v; double w; double z; - double n; - double i; - double k; + int32_t n; + int32_t i; + int32_t k; - i = ( j & stdlib_base_high_base_abs_mask )|0; // asm type annotation - k = ( ( i >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS )|0; // asm type annotation + i = ( j & stdlib_base_high_base_abs_mask ); + k = ( ( i >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); n = 0; // `|z| > 0.5`, set `n = z+0.5` if ( i > HIGH_BIASED_EXP_NEG_1 ) { n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ) >> 0; // asm type annotation - k = ( ( ( n & stdlib_base_high_base_abs_mask ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS )|0; // new k for n + k = ( ( ( n & stdlib_base_high_base_abs_mask ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n tmp = ( ( n & !( HIGH_SIGNIFICAND_MASK >> k ) ) ) >> 0; // asm type annotation t = stdlib_base_float64_set_high_word( 0.0, tmp ); n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( HIGH_NUM_SIGNIFICAND_BITS - k ) ) >> 0; // eslint-disable-line max-len @@ -334,7 +337,7 @@ double pow2( const double j, const double hp, const double lp ) { hp -= t; } t = lp + hp; - t = stdlib_base_float64_set_low_word ( t, 0 ); + stdlib_base_float64_set_low_word ( 0, &t ); u = t * LN2_HI; v = ( ( lp - ( t - hp ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); z = u + v; @@ -343,7 +346,7 @@ double pow2( const double j, const double hp, const double lp ) { t1 = z - ( t * polyval_p ( t ) ); r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); z = 1.0 - ( r - z ); - j = stdlib_base_float64_get_high_word ( z ); + stdlib_base_float64_get_high_word ( z, &j ); j = j|0; j += ( n << HIGH_NUM_SIGNIFICAND_BITS ) >> 0; // asm type annotation @@ -351,7 +354,7 @@ double pow2( const double j, const double hp, const double lp ) { if ( ( j >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { z = stdlib_base_ldexp ( z, n ); } else { - z = stdlib_base_float64_set_high_word ( z, j ); + stdlib_base_float64_set_high_word ( j, &z ); } return z; } @@ -380,7 +383,7 @@ void logx( const double *out, const double ax ) { u = INV_LN2_HI * t; // `INV_LN2_HI` has `21` significant bits v = ( t * INV_LN2_LO ) - ( w * INV_LN2 ); t1 = u + v; - t1 = stdlib_base_float64_set_low_word ( t1, 0 ); + stdlib_base_float64_set_low_word ( 0, &t1 ); t2 = v - (t1 - u); out[ 0 ] = t1; @@ -399,7 +402,7 @@ void logx( const double *out, const double ax ) { * double out[2] = logax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] * // returns [ 3.169923782348633, 0.0000012190936795504075 ] */ -void logax( const double *out, const double ax, const double ahx ) { +void log2ax( const double *out, const double ax, const uint32_t ahx ) { double tmp; double ss; // hs + ls double s2; // ss squared @@ -419,26 +422,26 @@ void logax( const double *out, const double ax, const double ahx ) { double r; double u; double v; - double n; - double j; - double k; + uint32_t n; + uint32_t j; + int32_t k; - n = 0|0; // asm type annotation + n = 0; // Check if `x` is subnormal... if ( ahx < HIGH_MIN_NORMAL_EXP ) { ax *= TWO53; - n -= 53|0; // asm type annotation - ahx = stdlib_base_float64_get_high_word ( ax ); + n -= 53; + stdlib_base_float64_get_high_word ( ax, &ahx ); } // Extract the unbiased exponent of `x`: - n += ( ( ahx >> HIGH_NUM_SIGNIFICAND_BITS ) - BIAS )|0; // asm type annotation + n += ( ( ahx >> HIGH_NUM_SIGNIFICAND_BITS ) - BIAS ); // Isolate the significand bits of `x`: - j = ( ahx & HIGH_SIGNIFICAND_MASK )|0; // asm type annotation + j = ( ahx & HIGH_SIGNIFICAND_MASK ); // Normalize `ahx` by setting the (biased) exponent to `1023`: - ahx = ( j|HIGH_BIASED_EXP_0 )|0; // asm type annotation + ahx = ( j|HIGH_BIASED_EXP_0 ); // Determine the interval of `|x|` by comparing significand bits... @@ -453,11 +456,11 @@ void logax( const double *out, const double ax, const double ahx ) { // |x| >= sqrt(3) else { k = 0; - n += 1|0; // asm type annotation + n += 1; ahx -= HIGH_MIN_NORMAL_EXP; } // Load the normalized high word into `|x|`: - ax = stdlib_base_float64_set_high_word ( ax, ahx ); + stdlib_base_float64_set_high_word ( ahx, &ax ); // 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 @@ -480,7 +483,7 @@ void logax( const double *out, const double ax, const double ahx ) { r += ls * (hs + ss); s2 = hs * hs; ht = 3.0 + s2 + r; - ht = stdlib_base_float64_set_low_word ( ht, 0 ); + stdlib_base_float64_set_low_word ( 0, &ht ); lt = r - ( ( ht - 3.0 ) - s2 ); // u+v = ss*(1+...): @@ -489,7 +492,7 @@ void logax( const double *out, const double ax, const double ahx ) { // 2/(3LN2) * (ss+...): hp = u + v; - hp = stdlib_base_float64_set_low_word ( hp, 0 ); + 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 ]; @@ -498,7 +501,7 @@ void logax( const double *out, const double ax, const double ahx ) { dp = DP_HI[ k ]; t = n; t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) - t1 = stdlib_base_float64_set_low_word ( t1, 0 ); + stdlib_base_float64_set_low_word ( 0, &t1 ); t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); out[ 0 ] = t1; @@ -521,15 +524,15 @@ void logax( const double *out, const double ax, const double ahx ) { * // returns 2.0 */ double pow( const double x, const double y ) { - double ahx; // absolute value high word `x` - double ahy; // absolute value high word `y` + uint32_t ahx; // absolute value high word `x` + uint32_t ahy; // absolute value high word `y` double ax; // absolute value `x` uint32_t hx; // high word `x` uint32_t lx; // low word `x` uint32_t hy; // high word `y` uint32_t ly; // low word `y` - double sx; // sign `x` - double sy; // sign `y` + uint32_t sx; // sign `x` + uint32_t sy; // sign `y` double y1; double hp; double lp; @@ -612,12 +615,12 @@ double pow( const double x, const double y ) { ax = stdlib_base_abs ( x ); // Remove the sign bits (i.e., get absolute values): - ahx = ( hx & ABS_MASK )|0; // asm type annotation - ahy = ( hy & ABS_MASK )|0; // asm type annotation + ahx = ( hx & ABS_MASK ); + ahy = ( hy & ABS_MASK ); // Extract the sign bits: - sx = ( hx >> HIGH_NUM_NONSIGN_BITS )|0; // asm type annotation - sy = ( hy >> HIGH_NUM_NONSIGN_BITS )|0; // asm type annotation + 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 ) ) { @@ -661,15 +664,14 @@ double pow( const double x, const double y ) { t = log2ax ( LOG_WORKSPACE, ax, ahx ); } // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... - y1 = stdlib_base_float64_set_low_word ( y, 0 ); + y1 = y; + stdlib_base_float64_set_low_word ( 0, &y1 ); lp = ( ( y - y1 ) * t[0] ) + ( y * t[1] ); hp = y1 * t[0]; z = lp + hp; // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner. - toWords.assign( z, &j, &i ); - j = j|0; - i = i|0; + stdlib_base_float64_to_words ( z, &j, &i ); // z >= 1024 if ( j >= HIGH_BIASED_EXP_10 ) { @@ -699,4 +701,4 @@ double pow( const double x, const double y ) { z = pow2( j, hp, lp ); return sx * z; -} \ No newline at end of file +} From d4904b3c1fc04f91950c80d6ca3e34301218a92b Mon Sep 17 00:00:00 2001 From: aman-095 Date: Thu, 7 Mar 2024 00:46:02 +0530 Subject: [PATCH 03/18] feat: add c implementation of pow --- .../@stdlib/math/base/special/pow/README.md | 10 +- .../math/base/special/pow/examples/c/Makefile | 2 +- .../base/special/pow/examples/c/example.c | 12 +- .../math/base/special/pow/manifest.json | 12 +- .../@stdlib/math/base/special/pow/src/main.c | 223 ++++++++++-------- 5 files changed, 140 insertions(+), 119 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/README.md b/lib/node_modules/@stdlib/math/base/special/pow/README.md index c4174cdff03b..5c7451c5bb1e 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/README.md +++ b/lib/node_modules/@stdlib/math/base/special/pow/README.md @@ -188,12 +188,12 @@ int main( void ) { double x; double out; int i; - + for ( i = 0; i < 100; i++ ) { - b = round( randu()*10.0 ); - x = round( randu()*10.0 ) - 5.0; - out = stdlib_base_pow ( b, x ); - printf( "pow(%lf, %lf) = %lf\n", b, x, v ); + b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); + x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; + out = stdlib_base_pow( b, x ); + printf( "pow(%lf, %lf) = %lf\n", b, x, out ); } } ``` diff --git a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile index 91d364d19fc3..6aed70daf167 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile @@ -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/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c index 04ab94c1850b..f4880a23de60 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c @@ -21,15 +21,15 @@ #include int main( void ) { + double out; double b; double x; - double out; int i; - + for ( i = 0; i < 100; i++ ) { - b = round( randu()*10.0 ); - x = round( randu()*10.0 ) - 5.0; - out = stdlib_base_pow ( b, x ); - printf( "pow(%lf, %lf) = %lf\n", b, x, v ); + b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); + x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; + out = stdlib_base_pow( b, x ); + printf( "pow(%lf, %lf) = %lf\n", b, x, out ); } } diff --git a/lib/node_modules/@stdlib/math/base/special/pow/manifest.json b/lib/node_modules/@stdlib/math/base/special/pow/manifest.json index 724bc3468bd8..e466fe9a877a 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/manifest.json +++ b/lib/node_modules/@stdlib/math/base/special/pow/manifest.json @@ -47,14 +47,14 @@ "@stdlib/number/float64/base/set-high-word", "@stdlib/number/float64/base/set-low-word", "@stdlib/math/base/special/ldexp", - "@stdlib/contants/float64/ln-two", + "@stdlib/constants/float64/ln-two", "@stdlib/constants/float64/exponent-bias", "@stdlib/constants/float64/high-word-significand-mask", "@stdlib/math/base/assert/is-nan", "@stdlib/math/base/assert/is-infinite", "@stdlib/math/base/assert/is-integer", "@stdlib/math/base/special/sqrt", - "@stdlib/math/base/assert/to-words" + "@stdlib/number/float64/base/to-words" ] }, { @@ -78,14 +78,14 @@ "@stdlib/number/float64/base/set-high-word", "@stdlib/number/float64/base/set-low-word", "@stdlib/math/base/special/ldexp", - "@stdlib/contants/float64/ln-two", + "@stdlib/constants/float64/ln-two", "@stdlib/constants/float64/exponent-bias", "@stdlib/constants/float64/high-word-significand-mask", "@stdlib/math/base/assert/is-nan", "@stdlib/math/base/assert/is-infinite", "@stdlib/math/base/assert/is-integer", "@stdlib/math/base/special/sqrt", - "@stdlib/math/base/assert/to-words" + "@stdlib/number/float64/base/to-words" ] }, { @@ -109,14 +109,14 @@ "@stdlib/number/float64/base/set-high-word", "@stdlib/number/float64/base/set-low-word", "@stdlib/math/base/special/ldexp", - "@stdlib/contants/float64/ln-two", + "@stdlib/constants/float64/ln-two", "@stdlib/constants/float64/exponent-bias", "@stdlib/constants/float64/high-word-significand-mask", "@stdlib/math/base/assert/is-nan", "@stdlib/math/base/assert/is-infinite", "@stdlib/math/base/assert/is-integer", "@stdlib/math/base/special/sqrt", - "@stdlib/math/base/assert/to-words" + "@stdlib/number/float64/base/to-words" ] } ] diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index 0538155125ac..5c75bbad7660 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -27,14 +27,14 @@ #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/contants/float64/ln_two.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" -#include "stdlib/math/base/assert/to_words.h" +#include "stdlib/math/base/special/sqrt.h" +#include "stdlib/number/float64/base/to_words.h" @@ -105,22 +105,19 @@ static const uint32_t HIGH_BIASED_EXP_31 = 0x41e00000; static const uint32_t HIGH_BIASED_EXP_64 = 0x43f00000; // 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 -static const uint32_t HIGH_BIASED_EXP_10 = 0x40900000; +static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; // 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 -static const uint32_t HIGH_1075 = 0x4090cc00; +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_NEG_1075 = 0xc090cc00; static const uint32_t HIGH_NUM_NONSIGN_BITS = 31; // -(1024-log2(ovfl+.5ulp)) static const double OVT = 8.0085662595372944372e-17; -// High/low words workspace: -uint32_t WORDS[2] = { 0, 0 }; - // Log workspace: double LOG_WORKSPACE[2] = { 0.0, 0.0 }; @@ -211,7 +208,7 @@ double y_is_infinite( const double x, const double y ) { return 1.0; } // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) - if ( (stdlib_base_abs ( x ) < 1.0) == (y == STDLIB_CONSTANT_FLOAT64_PINF) ) { + if ( (stdlib_base_abs( x ) < 1.0) == (y == STDLIB_CONSTANT_FLOAT64_PINF) ) { return 0.0; } // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) @@ -239,8 +236,8 @@ double y_is_huge( const double x, const double y ) { double xc; xc = x; - stdlib_base_float64_get_high_word ( xc, &hx ); - ahx = ( hx && STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + stdlib_base_float64_get_high_word( xc, &hx ); + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); if ( ahx <= HIGH_MAX_NEAR_UNITY ) { if ( y < 0 ) { @@ -283,13 +280,13 @@ double x_is_zero( const double x, const double y ) { return 0.0; } if ( y > 0.0 ) { - if ( stdlib_base_is_odd ( y ) ) { + if ( stdlib_base_is_odd( y ) ) { return x; // handles +-0 } return 0.0; } // y < 0.0 - if ( stdlib_base_is_odd ( y ) ) { + if ( stdlib_base_is_odd( y ) ) { return stdlib_base_copysign( STDLIB_CONSTANT_FLOAT64_PINF, x ); // handles +-0 } return STDLIB_CONSTANT_FLOAT64_PINF; @@ -308,7 +305,7 @@ double x_is_zero( const double x, const double y ) { * // returns ~0.79 */ double pow2( const uint32_t j, const double hp, const double lp ) { - double tmp; + int32_t tmp; double t1; double t; double r; @@ -317,44 +314,51 @@ double pow2( const uint32_t j, const double hp, const double lp ) { double w; double z; int32_t n; - int32_t i; + uint32_t i; int32_t k; - - i = ( j & stdlib_base_high_base_abs_mask ); + int32_t jc; + uint32_t jcc; + double hpc; + + hpc = hp; + jc = j; + jcc = j; + i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); k = ( ( i >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); n = 0; // `|z| > 0.5`, set `n = z+0.5` if ( i > HIGH_BIASED_EXP_NEG_1 ) { n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ) >> 0; // asm type annotation - k = ( ( ( n & stdlib_base_high_base_abs_mask ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n + k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n tmp = ( ( n & !( HIGH_SIGNIFICAND_MASK >> k ) ) ) >> 0; // asm type annotation - t = stdlib_base_float64_set_high_word( 0.0, tmp ); + t = 0.0; + stdlib_base_float64_set_high_word( tmp, &t ); n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( HIGH_NUM_SIGNIFICAND_BITS - k ) ) >> 0; // eslint-disable-line max-len - if ( j < 0 ) { + if ( jc < 0 ) { n = -n; } - hp -= t; + hpc -= t; } - t = lp + hp; - stdlib_base_float64_set_low_word ( 0, &t ); + t = lp + hpc; + stdlib_base_float64_set_low_word( 0, &t ); u = t * LN2_HI; - v = ( ( lp - ( t - hp ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); + 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 ) ); + 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|0; - j += ( n << HIGH_NUM_SIGNIFICAND_BITS ) >> 0; // asm type annotation + stdlib_base_float64_get_high_word( z, &jcc ); + jc = jcc|0; + jc += ( n << HIGH_NUM_SIGNIFICAND_BITS ) >> 0; // asm type annotation // Check for subnormal output... - if ( ( j >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { - z = stdlib_base_ldexp ( z, n ); + if ( ( jc >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { + z = stdlib_base_ldexp( z, n ); } else { - stdlib_base_float64_set_high_word ( j, &z ); + stdlib_base_float64_set_high_word( jcc, &z ); } return z; } @@ -370,7 +374,7 @@ double pow2( const uint32_t j, const double hp, const double lp ) { * double out[2] = logx( [ 0.0, 0.0 ], 9.0 ); // => [ t1, t2 ] * // returns [ -1265.7236328125, -0.0008163940840404393 ] */ -void logx( const double *out, const double ax ) { +void logx( double *out, const double ax ) { double t2; double t1; double t; @@ -379,11 +383,11 @@ void logx( const double *out, const double ax ) { double v; t = ax - 1.0; // `t` has `20` trailing zeros - w = t * t * polyval_w ( t ); + 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 ); + stdlib_base_float64_set_low_word( 0, &t1 ); t2 = v - (t1 - u); out[ 0 ] = t1; @@ -402,8 +406,9 @@ void logx( const double *out, const double ax ) { * double out[2] = logax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] * // returns [ 3.169923782348633, 0.0000012190936795504075 ] */ -void log2ax( const double *out, const double ax, const uint32_t ahx ) { - double tmp; +void log2ax( double *out, const double ax, const uint32_t ahx ) { + uint32_t ahxc; + uint32_t tmp; double ss; // hs + ls double s2; // ss squared double hs; @@ -422,26 +427,29 @@ void log2ax( const double *out, const double ax, const uint32_t ahx ) { double r; double u; double v; - uint32_t n; - uint32_t j; - int32_t k; + int32_t n; + int32_t j; + uint32_t k; + double axc; n = 0; + axc = ax; + ahxc = ahx; // Check if `x` is subnormal... - if ( ahx < HIGH_MIN_NORMAL_EXP ) { - ax *= TWO53; + if ( ahxc < HIGH_MIN_NORMAL_EXP ) { + axc *= TWO53; n -= 53; - stdlib_base_float64_get_high_word ( ax, &ahx ); + stdlib_base_float64_get_high_word( axc, &ahxc ); } // Extract the unbiased exponent of `x`: - n += ( ( ahx >> HIGH_NUM_SIGNIFICAND_BITS ) - BIAS ); + n += ( ( ahxc >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // Isolate the significand bits of `x`: - j = ( ahx & HIGH_SIGNIFICAND_MASK ); + j = ( ahxc & HIGH_SIGNIFICAND_MASK ); // Normalize `ahx` by setting the (biased) exponent to `1023`: - ahx = ( j|HIGH_BIASED_EXP_0 ); + ahxc = ( j|HIGH_BIASED_EXP_0 ); // Determine the interval of `|x|` by comparing significand bits... @@ -457,33 +465,35 @@ void log2ax( const double *out, const double ax, const uint32_t ahx ) { else { k = 0; n += 1; - ahx -= HIGH_MIN_NORMAL_EXP; + ahxc -= HIGH_MIN_NORMAL_EXP; } // Load the normalized high word into `|x|`: - stdlib_base_float64_set_high_word ( ahx, &ax ); + stdlib_base_float64_set_high_word( ahxc, &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 = ax - bp; // (x-1) || (x-1.5) - v = 1.0 / (ax + bp); // 1/(x+1) || 1/(x+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 = stdlib_base_float64_set_low_word ( ss, 0 ); // set all low word (less significant significand) bits to 0s + 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 = ( ( ahx >> 1 ) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; + 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 = stdlib_base_float64_set_high_word ( 0.0, tmp ); - lt = ax - ( ht - bp ); + 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 = 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 ); + stdlib_base_float64_set_low_word( 0, &ht ); lt = r - ( ( ht - 3.0 ) - s2 ); // u+v = ss*(1+...): @@ -492,7 +502,7 @@ void log2ax( const double *out, const double ax, const uint32_t ahx ) { // 2/(3LN2) * (ss+...): hp = u + v; - stdlib_base_float64_set_low_word ( 0, &hp ); + 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 ]; @@ -501,7 +511,7 @@ void log2ax( const double *out, const double ax, const uint32_t ahx ) { dp = DP_HI[ k ]; t = n; t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) - stdlib_base_float64_set_low_word ( 0, &t1 ); + stdlib_base_float64_set_low_word( 0, &t1 ); t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); out[ 0 ] = t1; @@ -516,14 +526,14 @@ void log2ax( const double *out, const double ax, const uint32_t ahx ) { * @return function value * * @example -* double out = pow( 2.0, 3.0 ); +* double out = stdlib_base_pow( 2.0, 3.0 ); * // returns 8.0 * * @example -* double out = pow( 4.0, 0.5 ); +* double out = stdlib_base_pow( 4.0, 0.5 ); * // returns 2.0 */ -double pow( const double x, const double y ) { +double stdlib_base_pow( const double x, const double y ) { uint32_t ahx; // absolute value high word `x` uint32_t ahy; // absolute value high word `y` double ax; // absolute value `x` @@ -531,20 +541,25 @@ double pow( const double x, const double y ) { uint32_t lx; // low word `x` uint32_t hy; // high word `y` uint32_t ly; // low word `y` - uint32_t sx; // sign `x` - uint32_t sy; // sign `y` + int32_t sx; // sign `x` + int32_t sy; // sign `y` double y1; double hp; double lp; - double t; + double t[ 2 ] = {0.0, 0.0}; double z; // y prime uint32_t j; uint32_t i; - if ( stdlib_base_is_nan ( x ) || stdlib_base_is_nan ( y ) ) { - return 0.0/0.0; + double xc; + int32_t jc; + int32_t ic; + + 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 ); + stdlib_base_float64_to_words( y, &hy, &ly ); // Special cases `y`... if ( ly == 0 ) { @@ -552,52 +567,52 @@ double pow( const double x, const double y ) { return 1.0; } if ( y == 1.0 ) { - return x; + return xc; } if ( y == -1.0 ) { - return 1.0 / x; + return 1.0 / xc; } if ( y == 0.5 ) { - return stdlib_base_sqrt( x ); + return stdlib_base_sqrt( xc ); } if ( y == -0.5 ) { - return 1.0 / stdlib_base_sqrt( x ); + return 1.0 / stdlib_base_sqrt( xc ); } if ( y == 2.0 ) { - return x * x; + return xc * xc; } if ( y == 3.0 ) { - return x * x * x; + return xc * xc * xc; } if ( y == 4.0 ) { - x *= x; - return x * x; + xc *= xc; + return xc * xc; } if ( stdlib_base_is_infinite( y ) ) { - return y_is_infinite( x, y ); + return y_is_infinite( xc, y ); } } // Split `x` into high and low words: - stdlib_base_float64_to_words ( x, &hx, &lx ); + stdlib_base_float64_to_words( xc, &hx, &lx ); // Special cases `x`... if ( lx == 0 ) { if ( hx == 0 ) { - return x_is_zero( x, y ); + return x_is_zero( xc, y ); } - if ( x == 1.0 ) { + if ( xc == 1.0 ) { return 1.0; } if ( - x == -1.0 && + xc == -1.0 && stdlib_base_is_odd( y ) ) { return -1.0; } - if ( stdlib_base_is_infinite( x ) ) { - if ( x == STDLIB_CONSTANT_FLOAT64_NINF ) { - // `pow( 1/x, -y )` - return pow( -0.0, -y ); + 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; @@ -606,17 +621,17 @@ double pow( const double x, const double y ) { } } if ( - x < 0.0 && - stdlib_base_is_integer ( y ) == false + xc < 0.0 && + stdlib_base_is_integer( y ) == false ) { // Signal NaN... return 0.0/0.0; } - ax = stdlib_base_abs ( x ); + ax = stdlib_base_abs( x ); // Remove the sign bits (i.e., get absolute values): - ahx = ( hx & ABS_MASK ); - ahy = ( hy & ABS_MASK ); + 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 ); @@ -624,9 +639,9 @@ double pow( const double x, const double y ) { // Determine the sign of the result... if ( sx && stdlib_base_is_odd( y ) ) { - sx = -1.0; + sx = -1; } else { - sx = 1.0; + sx = 1; } // Case 1: `|y|` is huge... @@ -634,7 +649,7 @@ double pow( const double x, const double y ) { if ( ahy > HIGH_BIASED_EXP_31 ) { // `|y| > 2^64`, then must over- or underflow... if ( ahy > HIGH_BIASED_EXP_64 ) { - return y_is_huge( x, y ); + return y_is_huge( xc, y ); } // Over- or underflow if `x` is not close to unity... @@ -657,26 +672,32 @@ double pow( const double x, const double y ) { return sx * TINY * TINY; } // 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`. - t = logx ( LOG_WORKSPACE, ax ); + t[0] = LOG_WORKSPACE[0]; + t[1] = LOG_WORKSPACE[1]; + logx( t, ax ); } // Case 2: `|y|` is not huge... else { - t = log2ax ( LOG_WORKSPACE, ax, ahx ); + t[0] = LOG_WORKSPACE[0]; + t[1] = LOG_WORKSPACE[1]; + log2ax( t, ax, ahx ); } // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... y1 = y; - stdlib_base_float64_set_low_word ( 0, &y1 ); + stdlib_base_float64_set_low_word( 0, &y1 ); lp = ( ( y - y1 ) * t[0] ) + ( y * t[1] ); hp = y1 * t[0]; 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 ); + stdlib_base_float64_to_words( z, &j, &i ); + jc = j; + ic = i; // z >= 1024 - if ( j >= HIGH_BIASED_EXP_10 ) { + if ( jc >= HIGH_BIASED_EXP_10 ) { // z > 1024 - if ( ( ( j - HIGH_BIASED_EXP_10 )|i ) != 0 ) { + if ( ( ( jc - HIGH_BIASED_EXP_10 )|ic ) != 0 ) { // Signal overflow... return sx * HUGE * HUGE; } @@ -686,9 +707,9 @@ double pow( const double x, const double y ) { } } // z <= -1075 - else if ( ( j & stdlib_base_high_base_abs_mask ) >= HIGH_1075 ) { + else if ( ( jc & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { // z < -1075 - if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { + if ( ( ( jc - HIGH_NEG_1075 )|ic ) != 0 ) { // Signal underflow... return sx * TINY * TINY; } @@ -698,7 +719,7 @@ double pow( const double x, const double y ) { } } // Compute `2^(hp+lp)`... - z = pow2( j, hp, lp ); + z = pow2( jc, hp, lp ); return sx * z; } From c1090fb02841a4cf23c6889dde795f9f175101c3 Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 17:59:54 -0800 Subject: [PATCH 04/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/README.md b/lib/node_modules/@stdlib/math/base/special/pow/README.md index 5c7451c5bb1e..28d0d3c55a93 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/README.md +++ b/lib/node_modules/@stdlib/math/base/special/pow/README.md @@ -153,11 +153,11 @@ out = stdlib_base_pow( 4.0, 0.5 ); The function accepts the following arguments: -- **b**: `[in] double` input value. -- **x**: `[in] double` input value. +- **base**: `[in] double` base. +- **exponent**: `[in] double` exponent. ```c -double stdlib_base_pow( const double b, const double x ); +double stdlib_base_pow( const double base, const double exponent ); ```
@@ -184,9 +184,9 @@ double stdlib_base_pow( const double b, const double x ); #include int main( void ) { + double out; double b; double x; - double out; int i; for ( i = 0; i < 100; i++ ) { From 57f8a165426e59cf53911b4e1ac45c6067df6db7 Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:01:13 -0800 Subject: [PATCH 05/18] Apply suggestions from code review Signed-off-by: Athan --- .../math/base/special/pow/benchmark/c/native/benchmark.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c index 64d35e2f73b0..880a68893fa1 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/c/native/benchmark.c @@ -17,7 +17,7 @@ */ /** -* Benchmark `ellipk`. +* Benchmark `pow`. */ #include "stdlib/math/base/special/pow.h" #include From 763b5c49dba5e2c2b3e9a31c4ce1abf335bbb986 Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:01:59 -0800 Subject: [PATCH 06/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/binding.gyp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp b/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp index 507cb00291e7..ec3992233442 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp +++ b/lib/node_modules/@stdlib/math/base/special/pow/binding.gyp @@ -167,4 +167,4 @@ ], # end actions }, # end target copy_addon ], # end targets -} \ No newline at end of file +} From 348099a5e99f3f8c91025e76a7884ecf7c54d1c9 Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:02:58 -0800 Subject: [PATCH 07/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/include.gypi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/include.gypi b/lib/node_modules/@stdlib/math/base/special/pow/include.gypi index c6495fc1da3f..575cb043c0bf 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/include.gypi +++ b/lib/node_modules/@stdlib/math/base/special/pow/include.gypi @@ -50,4 +50,4 @@ ' Date: Wed, 6 Mar 2024 18:04:04 -0800 Subject: [PATCH 08/18] Apply suggestions from code review Signed-off-by: Athan --- .../base/special/pow/include/stdlib/math/base/special/pow.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/include/stdlib/math/base/special/pow.h b/lib/node_modules/@stdlib/math/base/special/pow/include/stdlib/math/base/special/pow.h index 4e56dd164877..5fddb2ea091f 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/include/stdlib/math/base/special/pow.h +++ b/lib/node_modules/@stdlib/math/base/special/pow/include/stdlib/math/base/special/pow.h @@ -29,10 +29,10 @@ extern "C" { /** * Evaluates the exponential function. */ -double stdlib_base_pow( const double b, const double x ); +double stdlib_base_pow( const double base, const double exponent ); #ifdef __cplusplus } #endif -#endif // !STDLIB_MATH_BASE_SPECIAL_POW_H \ No newline at end of file +#endif // !STDLIB_MATH_BASE_SPECIAL_POW_H From eabee3a47bc23d94249862fe9ec747af8f2882cb Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:07:35 -0800 Subject: [PATCH 09/18] Apply suggestions from code review Signed-off-by: Athan --- .../math/base/special/pow/lib/native.js | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js index ff804513f5e9..16b7078f8e08 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js @@ -29,31 +29,31 @@ var addon = require( './../src/addon.node' ); Evaluates the exponential function. * * @private -* @param {number} b - input value -* @param {number} x - input value +* @param {number} x - base +* @param {number} y - exponent * @returns {number} exponential value * * @example * var v = pow( 2.0, 3.0 ); * // returns 8.0 * -* v = ellipk( 4.0, 0.5 ); +* v = pow( 4.0, 0.5 ); * // returns 2.0 * -* v = ellipk( 100.0, 0.0 ); +* v = pow( 100.0, 0.0 ); * // returns 1.0 * -* v = ellipk( 3.141592653589793, -0.2 ); +* v = pow( 3.141592653589793, -0.2 ); * // returns ~0.7954 * -* v = ellipk( NaN, 3.0 ); +* v = pow( NaN, 3.0 ); * // returns NaN * -* v = ellipk( 5.0, NaN ); +* v = pow( 5.0, NaN ); * // returns NaN */ -function pow( b, x ) { - return addon( b, x ); +function pow( x, y ) { + return addon( x, y ); } From ad9750aab405e5bdf467b8ab8ff8b244cbe0e40d Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:07:53 -0800 Subject: [PATCH 10/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/lib/native.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js index 16b7078f8e08..9a527c6001a0 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js @@ -26,7 +26,7 @@ var addon = require( './../src/addon.node' ); // MAIN // /** -Evaluates the exponential function. +* Evaluates the exponential function. * * @private * @param {number} x - base From a27095aeb7e868fa53c3d0652a6c003af124124c Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:08:40 -0800 Subject: [PATCH 11/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/lib/native.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js index 9a527c6001a0..8c7b26489b92 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js @@ -31,7 +31,7 @@ var addon = require( './../src/addon.node' ); * @private * @param {number} x - base * @param {number} y - exponent -* @returns {number} exponential value +* @returns {number} function value * * @example * var v = pow( 2.0, 3.0 ); From aa62c7139e663a998152a2b0ffc4ae063bbde6c3 Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:09:46 -0800 Subject: [PATCH 12/18] Apply suggestions from code review Signed-off-by: Athan --- .../math/base/special/pow/lib/native.js | 22 ++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js index 8c7b26489b92..e37dc181675c 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/lib/native.js @@ -37,20 +37,32 @@ var addon = require( './../src/addon.node' ); * var v = pow( 2.0, 3.0 ); * // returns 8.0 * -* v = pow( 4.0, 0.5 ); +* @example +* var v = pow( 4.0, 0.5 ); * // returns 2.0 * -* v = pow( 100.0, 0.0 ); +* @example +* var v = pow( 100.0, 0.0 ); * // returns 1.0 * -* v = pow( 3.141592653589793, -0.2 ); +* @example +* var v = pow( 3.141592653589793, 5.0 ); +* // returns ~306.0197 +* +* @example +* var v = pow( 3.141592653589793, -0.2 ); * // returns ~0.7954 * -* v = pow( NaN, 3.0 ); +* @example +* var v = pow( NaN, 3.0 ); * // returns NaN * -* v = pow( 5.0, NaN ); +* @example +* var v = pow( 5.0, NaN ); * // returns NaN +* +* @example +* var v = pow( NaN, NaN ); */ function pow( x, y ) { return addon( x, y ); From 019b21e6f7d857b4f86cf6e4e1b4b26b2b882abc Mon Sep 17 00:00:00 2001 From: Athan Date: Wed, 6 Mar 2024 18:12:07 -0800 Subject: [PATCH 13/18] Apply suggestions from code review Signed-off-by: Athan --- lib/node_modules/@stdlib/math/base/special/pow/src/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/Makefile b/lib/node_modules/@stdlib/math/base/special/pow/src/Makefile index 81bb164286c3..bcf18aa46655 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/Makefile +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/Makefile @@ -67,4 +67,4 @@ clean-addon: #/ clean: clean-addon -.PHONY: clean \ No newline at end of file +.PHONY: clean From 86115e5a63bedd267d8c6a58655f2f77af4b8285 Mon Sep 17 00:00:00 2001 From: aman-095 Date: Fri, 8 Mar 2024 15:20:19 +0530 Subject: [PATCH 14/18] fix: incorrect implementation of pow --- .../special/pow/benchmark/benchmark.native.js | 4 +- .../base/special/pow/examples/c/example.c | 6 +- .../@stdlib/math/base/special/pow/src/addon.c | 1 - .../@stdlib/math/base/special/pow/src/main.c | 228 +++++++++--------- 4 files changed, 115 insertions(+), 124 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js index 007e3b2fabc4..1713caa183dd 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/benchmark/benchmark.native.js @@ -46,8 +46,8 @@ bench( pkg+'::native', opts, function benchmark( b ) { b.tic(); for ( i = 0; i < b.iterations; i++ ) { - v = round( randu()*10.0 ); - x = round( randu()*10.0 ) - 5.0; + v = ( randu()*10.0 ); + x = ( randu()*10.0 ) - 5.0; y = pow( v, x ); if ( isnan( y ) ) { b.fail( 'should not return NaN' ); diff --git a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c index f4880a23de60..6eb35dea0850 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c @@ -21,14 +21,14 @@ #include int main( void ) { - double out; + double out; double b; double x; int i; for ( i = 0; i < 100; i++ ) { - b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); - x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; + b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); + x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; out = stdlib_base_pow( b, x ); printf( "pow(%lf, %lf) = %lf\n", b, x, out ); } diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/addon.c b/lib/node_modules/@stdlib/math/base/special/pow/src/addon.c index 89a949de9196..ff8767b9e48c 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/addon.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/addon.c @@ -19,5 +19,4 @@ #include "stdlib/math/base/special/pow.h" #include "stdlib/math/base/napi/binary.h" -// cppcheck-suppress shadowFunction STDLIB_MATH_BASE_NAPI_MODULE_DD_D( stdlib_base_pow ) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index 5c75bbad7660..11f85dc28827 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -36,23 +36,21 @@ #include "stdlib/math/base/special/sqrt.h" #include "stdlib/number/float64/base/to_words.h" - - // 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1 -static const uint32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; +static const int32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; -static const double HUGE = 1.0e300; +static const double HUGE_VALUE = 1.0e300; -static const double TINY = 1.0e-300; +static const double TINY_VALUE = 1.0e-300; // 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 -static const uint32_t HIGH_MIN_NORMAL_EXP = 0x00100000; +static const int32_t HIGH_MIN_NORMAL_EXP = 0x00100000; // 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 -static const uint32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; +static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; // TODO: consider making into an external constant -static const uint32_t HIGH_NUM_SIGNIFICAND_BITS = 20; +static const int32_t HIGH_NUM_SIGNIFICAND_BITS = 20; // High: LN2 static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 @@ -70,16 +68,16 @@ static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60 static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 // 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 -static const uint32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; +static const int32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; // 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 -static const uint32_t HIGH_BIASED_EXP_0 = 0x3ff00000; +static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; // 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 -static const uint32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; +static const int32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; // 0x00080000 = 524288 => 0 00000000000 10000000000000000000 -static const uint32_t HIGH_SIGNIFICAND_HALF = 0x00080000; +static const int32_t HIGH_SIGNIFICAND_HALF = 0x00080000; static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 @@ -92,17 +90,17 @@ static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE000000 // Low: CP_HI static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 -double BP [2] = { 1.0, 1.5 }; +static const double BP [2] = { 1.0, 1.5 }; -double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 +static const double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 -double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 +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 uint32_t HIGH_BIASED_EXP_31 = 0x41e00000; +static const int32_t HIGH_BIASED_EXP_31 = 0x41e00000; // 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 -static const uint32_t HIGH_BIASED_EXP_64 = 0x43f00000; +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; @@ -111,16 +109,13 @@ static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; static const int32_t HIGH_1075 = 0x4090cc00; // 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 -static const int32_t HIGH_NEG_1075 = 0xc090cc00; +static const uint32_t HIGH_NEG_1075 = 0xc090cc00; -static const uint32_t HIGH_NUM_NONSIGN_BITS = 31; +static const int32_t HIGH_NUM_NONSIGN_BITS = 31; // -(1024-log2(ovfl+.5ulp)) static const double OVT = 8.0085662595372944372e-17; -// Log workspace: -double LOG_WORKSPACE[2] = { 0.0, 0.0 }; - /* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ // BEGIN: polyval_l @@ -188,19 +183,19 @@ static double polyval_p( const double x ) { /** * Evaluates the exponential function when \\( y = \pm \infty\\). * -* @param x input value -* @param y input value +* @param x base +* @param y exponent * @return output value * * @example -* double out = y_is_infinite( -1.0, Infinity ); +* double out = y_is_infinite( -1.0, 1.0/0.0 ); * // returns NaN * * @example -* double out = y_is_infinite( 1.5, -Infinity ); +* double out = y_is_infinite( 1.5, -1.0/0.0 ); * // returns 0.0 */ -double y_is_infinite( const double x, const double y ) { +static double y_is_infinite( const double x, const double y ) { if ( x == -1.0 ) { return 0.0 / 0.0; // NaN } @@ -208,7 +203,7 @@ double y_is_infinite( const double x, const double y ) { return 1.0; } // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) - if ( (stdlib_base_abs( x ) < 1.0) == (y == STDLIB_CONSTANT_FLOAT64_PINF) ) { + if ( ( stdlib_base_abs( x ) < 1.0 ) == ( y == STDLIB_CONSTANT_FLOAT64_PINF ) ) { return 0.0; } // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) @@ -218,8 +213,8 @@ double y_is_infinite( const double x, const double y ) { /** * Evaluates the exponential function when \\(|y| > 2^64\\). * -* @param x input value -* @param y input value +* @param x base +* @param y exponent * @return output value * * @example @@ -230,39 +225,37 @@ double y_is_infinite( const double x, const double y ) { * double out = y_is_huge( -3.14, -3.6893488147419103e19 ); * // returns 0.0 */ -double y_is_huge( const double x, const double y ) { +static double y_is_huge( const double x, const double y ) { + uint32_t ahx; uint32_t hx; - double ahx; - double xc; - xc = x; - stdlib_base_float64_get_high_word( xc, &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 * HUGE; + return HUGE_VALUE * HUGE_VALUE; } // Signal underflow... - return TINY * TINY; + return TINY_VALUE * TINY_VALUE; } // `x` has a biased exponent greater than or equal to `0`... if ( y > 0 ) { // Signal overflow... - return HUGE * HUGE; + return HUGE_VALUE * HUGE_VALUE; } // Signal underflow... - return TINY * TINY; + return TINY_VALUE * TINY_VALUE; } /** * Evaluates the exponential function when \\(|x| = 0\\). * -* @param x input value -* @param y input value -* @return output value +* @param x base +* @param y exponent +* @return output value * * @example * double out = x_is_zero( 0.0, 2 ); @@ -272,7 +265,7 @@ double y_is_huge( const double x, const double y ) { * double out = x_is_zero( 0.0, -9 ); * // returns Infinity */ -double x_is_zero( const double x, const double y ) { +static double x_is_zero( const double x, const double y ) { if ( y == STDLIB_CONSTANT_FLOAT64_NINF ) { return STDLIB_CONSTANT_FLOAT64_PINF; } @@ -295,17 +288,24 @@ double x_is_zero( const double x, const double y ) { /** * 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 +* @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 */ -double pow2( const uint32_t j, const double hp, const double lp ) { - int32_t tmp; +static double pow2( const int32_t j, const double hp, const double lp ) { + uint32_t jcc; + 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; @@ -313,30 +313,26 @@ double pow2( const uint32_t j, const double hp, const double lp ) { double v; double w; double z; - int32_t n; - uint32_t i; - int32_t k; - int32_t jc; - uint32_t jcc; - double hpc; hpc = hp; jc = j; - jcc = j; + jcc = (uint32_t)j; i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); k = ( ( i >> HIGH_NUM_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 ) ) ) >> 0; // asm type annotation + n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ); k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n - tmp = ( ( n & !( HIGH_SIGNIFICAND_MASK >> k ) ) ) >> 0; // asm type annotation + 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 ) >> ( HIGH_NUM_SIGNIFICAND_BITS - k ) ) >> 0; // eslint-disable-line max-len + n = ( ( ( n & HIGH_SIGNIFICAND_MASK ) | HIGH_MIN_NORMAL_EXP ) >> ( HIGH_NUM_SIGNIFICAND_BITS - k ) ); + nc = (int32_t)n; if ( jc < 0 ) { - n = -n; + nc = -nc; } hpc -= t; } @@ -351,12 +347,13 @@ double pow2( const uint32_t j, const double hp, const double lp ) { r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); z = 1.0 - ( r - z ); stdlib_base_float64_get_high_word( z, &jcc ); - jc = jcc|0; - jc += ( n << HIGH_NUM_SIGNIFICAND_BITS ) >> 0; // asm type annotation + jc = (int32_t)jcc; + jcc = jc + ( nc << HIGH_NUM_SIGNIFICAND_BITS ); + jc = (int32_t)jcc; // Check for subnormal output... if ( ( jc >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { - z = stdlib_base_ldexp( z, n ); + z = stdlib_base_ldexp( z, nc ); } else { stdlib_base_float64_set_high_word( jcc, &z ); } @@ -366,15 +363,14 @@ double pow2( const uint32_t j, const double hp, const double lp ) { /** * Computes \\(\operatorname{log}(x)\\) assuming \\(|1-x|\\) is small and using the approximation \\(x - x^2/2 + x^3/3 - x^4/4\\). * -* @param out output array * @param ax absolute value of x -* @return output array containing a tuple comprised of high and low parts +* @param o1 destination for output1 +* @param o2 destination for output2 * * @example -* double out[2] = logx( [ 0.0, 0.0 ], 9.0 ); // => [ t1, t2 ] -* // returns [ -1265.7236328125, -0.0008163940840404393 ] +* logx( 9.0, &o1, &o2 ); */ -void logx( double *out, const double ax ) { +void logx( const double ax, double *o1, double *o2 ) { double t2; double t1; double t; @@ -390,25 +386,29 @@ void logx( double *out, const double ax ) { stdlib_base_float64_set_low_word( 0, &t1 ); t2 = v - (t1 - u); - out[ 0 ] = t1; - out[ 1 ] = t2; + *o1 = t1; + *o2 = t2; } /** * Computes \\(\operatorname{log2}(ax)\\). * -* @param out output array * @param ax absolute value of x * @param ahx high word of ax -* @return output array containing a tuple comprised of high and low parts +* @param o1 destination for output1 +* @param o2 destination for output2 * * @example -* double out[2] = logax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] -* // returns [ 3.169923782348633, 0.0000012190936795504075 ] +* log2ax( 9.0, 1075970048, &o1, &o2 ); */ -void log2ax( double *out, const double ax, const uint32_t ahx ) { - uint32_t ahxc; +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; @@ -427,20 +427,17 @@ void log2ax( double *out, const double ax, const uint32_t ahx ) { double r; double u; double v; - int32_t n; - int32_t j; - uint32_t k; - double axc; 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, &ahxc ); + stdlib_base_float64_get_high_word( axc, &ahxcc ); + ahxc = (int32_t)ahxcc; } // Extract the unbiased exponent of `x`: n += ( ( ahxc >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); @@ -467,8 +464,9 @@ void log2ax( double *out, const double ax, const uint32_t ahx ) { n += 1; ahxc -= HIGH_MIN_NORMAL_EXP; } + ahxcc = (uint32_t)ahxc; // Load the normalized high word into `|x|`: - stdlib_base_float64_set_high_word( ahxc, &axc ); + 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 @@ -514,8 +512,8 @@ void log2ax( double *out, const double ax, const uint32_t ahx ) { stdlib_base_float64_set_low_word( 0, &t1 ); t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); - out[ 0 ] = t1; - out[ 1 ] = t2; + *o1 = t1; + *o2 = t2; } /** @@ -534,30 +532,32 @@ void log2ax( double *out, const double ax, const uint32_t ahx ) { * // returns 2.0 */ double stdlib_base_pow( const double x, const double y ) { - uint32_t ahx; // absolute value high word `x` - uint32_t ahy; // absolute value high word `y` - double ax; // absolute value `x` 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 t[ 2 ] = {0.0, 0.0}; - double z; // y prime - uint32_t j; - uint32_t i; + double t1; + double t2; double xc; - int32_t jc; - int32_t ic; + 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 ); @@ -594,7 +594,6 @@ double stdlib_base_pow( const double x, const double 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 ) { @@ -632,11 +631,9 @@ double stdlib_base_pow( const double x, const double y ) { // 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; @@ -657,69 +654,64 @@ double stdlib_base_pow( const double x, const double y ) { // y < 0 if ( sy == 1 ) { // Signal overflow... - return sx * HUGE * HUGE; + return sx * HUGE_VALUE * HUGE_VALUE; } // Signal underflow... - return sx * TINY * TINY; + return sx * TINY_VALUE * TINY_VALUE; } if ( ahx > HIGH_BIASED_EXP_0 ) { // y > 0 if ( sy == 0 ) { // Signal overflow... - return sx * HUGE * HUGE; + return sx * HUGE_VALUE * HUGE_VALUE; } // Signal underflow... - return sx * TINY * TINY; + 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`. - t[0] = LOG_WORKSPACE[0]; - t[1] = LOG_WORKSPACE[1]; - logx( t, ax ); + logx( ax, &t1, &t2 ); } // Case 2: `|y|` is not huge... else { - t[0] = LOG_WORKSPACE[0]; - t[1] = LOG_WORKSPACE[1]; - log2ax( t, ax, ahx ); + 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 ) * t[0] ) + ( y * t[1] ); - hp = y1 * t[0]; + 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 = j; - ic = 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 * HUGE; + return sx * HUGE_VALUE * HUGE_VALUE; } if ( ( lp + OVT ) > ( z - hp ) ) { // Signal overflow... - return sx * HUGE * HUGE; + return sx * HUGE_VALUE * HUGE_VALUE; } } // z <= -1075 - else if ( ( jc & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { + else if ( ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { // z < -1075 - if ( ( ( jc - HIGH_NEG_1075 )|ic ) != 0 ) { + if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { // Signal underflow... - return sx * TINY * TINY; + return sx * TINY_VALUE * TINY_VALUE; } if ( lp <= ( z - hp ) ) { // Signal underflow... - return sx * TINY * TINY; + return sx * TINY_VALUE * TINY_VALUE; } } // Compute `2^(hp+lp)`... z = pow2( jc, hp, lp ); - return sx * z; } From 300f9d3604dbcd38200bef08c71fdb236163acc2 Mon Sep 17 00:00:00 2001 From: aman-095 Date: Wed, 13 Mar 2024 14:47:09 +0530 Subject: [PATCH 15/18] chore: apply review changes --- .../base/special/pow/examples/c/example.c | 20 ++++---- .../@stdlib/math/base/special/pow/src/main.c | 46 +++++++++---------- 2 files changed, 33 insertions(+), 33 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c index 6eb35dea0850..ee26d99bea68 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c @@ -21,15 +21,15 @@ #include int main( void ) { - double out; - double b; - double x; - int i; + double out; + double b; + double x; + int i; - for ( i = 0; i < 100; i++ ) { - b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); - x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; - out = stdlib_base_pow( b, x ); - printf( "pow(%lf, %lf) = %lf\n", b, x, out ); - } + for ( i = 0; i < 100; i++ ) { + b = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ); + x = ( ( (double)rand() / (double)RAND_MAX ) * 10.0 ) - 5.0; + out = stdlib_base_pow( b, x ); + printf( "pow(%lf, %lf) = %lf\n", b, x, out ); + } } diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index 11f85dc28827..1ba8af724bda 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -196,18 +196,18 @@ static double polyval_p( const double x ) { * // 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; - } + 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; + return STDLIB_CONSTANT_FLOAT64_PINF; } /** @@ -229,23 +229,23 @@ 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 ); + 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 ) { + if ( ahx <= HIGH_MAX_NEAR_UNITY ) { + if ( y < 0 ) { // Signal overflow... - return HUGE_VALUE * HUGE_VALUE; - } + return HUGE_VALUE * HUGE_VALUE; + } // Signal underflow... - return TINY_VALUE * TINY_VALUE; - } + return TINY_VALUE * TINY_VALUE; + } // `x` has a biased exponent greater than or equal to `0`... - if ( y > 0 ) { + if ( y > 0 ) { // Signal overflow... - return HUGE_VALUE * HUGE_VALUE; - } + return HUGE_VALUE * HUGE_VALUE; + } // Signal underflow... return TINY_VALUE * TINY_VALUE; } @@ -266,9 +266,9 @@ static double y_is_huge( const double x, const double y ) { * // 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_NINF ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } if ( y == STDLIB_CONSTANT_FLOAT64_PINF ) { return 0.0; } @@ -323,7 +323,7 @@ static double pow2( const int32_t j, const double hp, const double lp ) { nc = (int32_t)n; // `|z| > 0.5`, set `n = z+0.5` - if ( i > HIGH_BIASED_EXP_NEG_1 ) { + if ( i > HIGH_BIASED_EXP_NEG_1 ) { n = ( j + ( HIGH_MIN_NORMAL_EXP >> ( k + 1 ) ) ); k = ( ( ( n & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); // new k for n tmp = ( ( n & ~( HIGH_SIGNIFICAND_MASK >> k ) ) ); @@ -335,7 +335,7 @@ static double pow2( const int32_t j, const double hp, const double lp ) { nc = -nc; } hpc -= t; - } + } t = lp + hpc; stdlib_base_float64_set_low_word( 0, &t ); u = t * LN2_HI; From 922a55e23537bffe6e9322b2cac9c413cdc5f695 Mon Sep 17 00:00:00 2001 From: Philipp Burckhardt Date: Fri, 7 Jun 2024 10:14:11 -0400 Subject: [PATCH 16/18] style: fix indentation, add empty lines before comments where applicable --- .../@stdlib/math/base/special/pow/src/main.c | 44 ++++++++++--------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index 1ba8af724bda..6a653e2c3b82 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -185,7 +185,7 @@ static double polyval_p( const double x ) { * * @param x base * @param y exponent -* @return output value +* @return output value * * @example * double out = y_is_infinite( -1.0, 1.0/0.0 ); @@ -215,7 +215,7 @@ static double y_is_infinite( const double x, const double y ) { * * @param x base * @param y exponent -* @return output value +* @return output value * * @example * double out = y_is_huge( 9.0, 3.6893488147419103e19 ); @@ -253,9 +253,9 @@ static double y_is_huge( const double x, const double y ) { /** * Evaluates the exponential function when \\(|x| = 0\\). * -* @param x base -* @param y exponent -* @return output value +* @param x base +* @param y exponent +* @return output value * * @example * double out = x_is_zero( 0.0, 2 ); @@ -288,10 +288,10 @@ static double x_is_zero( const double x, const double y ) { /** * 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 +* @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 ); @@ -393,10 +393,10 @@ void logx( const double ax, double *o1, double *o2 ) { /** * Computes \\(\operatorname{log2}(ax)\\). * -* @param ax absolute value of x +* @param ax absolute value of x * @param ahx high word of ax -* @param o1 destination for output1 -* @param o2 destination for output2 +* @param o1 destination for output1 +* @param o2 destination for output2 * * @example * log2ax( 9.0, 1075970048, &o1, &o2 ); @@ -465,6 +465,7 @@ void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { ahxc -= HIGH_MIN_NORMAL_EXP; } ahxcc = (uint32_t)ahxc; + // Load the normalized high word into `|x|`: stdlib_base_float64_set_high_word( ahxcc, &axc ); @@ -521,7 +522,7 @@ void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { * * @param x base * @param y exponent -* @return function value +* @return function value * * @example * double out = stdlib_base_pow( 2.0, 3.0 ); @@ -532,10 +533,10 @@ void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { * // 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` + 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; @@ -544,14 +545,14 @@ double stdlib_base_pow( const double x, const double y ) { int32_t jc; int32_t sx; // sign `x` int32_t sy; // sign `y` - double ax; // absolute value `x` + double ax; // absolute value `x` double y1; double hp; double lp; double t1; double t2; double xc; - double z; // y prime + double z; // y prime xc = x; if ( stdlib_base_is_nan( xc ) || stdlib_base_is_nan( y ) ) { @@ -594,6 +595,7 @@ double stdlib_base_pow( const double x, const double 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 ) { @@ -631,9 +633,11 @@ double stdlib_base_pow( const double x, const double y ) { // 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; @@ -649,7 +653,6 @@ double stdlib_base_pow( const double x, const double y ) { 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 ) { @@ -675,6 +678,7 @@ double stdlib_base_pow( const double x, const double y ) { 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 ); From a20a1fb3f0623fcff28cf360bbd33f02d0fe03e4 Mon Sep 17 00:00:00 2001 From: Philipp Burckhardt Date: Wed, 19 Jun 2024 20:37:34 +0000 Subject: [PATCH 17/18] test: increase tolerance due to divergence from compiler optimizations --- .../math/base/special/pow/test/test.native.js | 75 +++++++++++++++++-- 1 file changed, 68 insertions(+), 7 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js index 639a00b35600..9fafa645e19d 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js +++ b/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js @@ -221,6 +221,8 @@ tape( 'the function evaluates the exponential function (`x ~ 1`, `y` large)', op tape( 'the function evaluates the exponential function (`x ~ 1`, `y` huge)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -230,7 +232,14 @@ tape( 'the function evaluates the exponential function (`x ~ 1`, `y` huge)', opt expected = baseNearUnityHuge.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); @@ -411,6 +420,8 @@ tape( 'the function evaluates the exponential function (near underflow)', opts, tape( 'the function evaluates the exponential function (small `x`, large `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -422,8 +433,13 @@ tape( 'the function evaluates the exponential function (small `x`, large `y`)', actual = pow( x[i], y[i] ); if ( expected[i] === 5.0e-324 ) { t.equal( actual === expected[i] || actual === 0.0, true, 'pow('+x[i]+','+y[i]+') returns 5e-324 or 0' ); + } else if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); } else { - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); } } t.end(); @@ -449,6 +465,8 @@ tape( 'the function evaluates the exponential function (large `x`, small `y`)', tape( 'the function evaluates the exponential function (small `x`, small `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -458,7 +476,14 @@ tape( 'the function evaluates the exponential function (small `x`, small `y`)', expected = smallSmall.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); @@ -466,6 +491,8 @@ tape( 'the function evaluates the exponential function (small `x`, small `y`)', tape( 'the function evaluates the exponential function (decimal `x`, decimal `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -475,7 +502,14 @@ tape( 'the function evaluates the exponential function (decimal `x`, decimal `y` expected = decimalDecimal.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); @@ -483,6 +517,8 @@ tape( 'the function evaluates the exponential function (decimal `x`, decimal `y` tape( 'the function evaluates the exponential function (decimal `x`, integer `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -492,7 +528,14 @@ tape( 'the function evaluates the exponential function (decimal `x`, integer `y` expected = decimalInteger.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); @@ -500,6 +543,8 @@ tape( 'the function evaluates the exponential function (decimal `x`, integer `y` tape( 'the function evaluates the exponential function (integer `x`, decimal `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -509,7 +554,14 @@ tape( 'the function evaluates the exponential function (integer `x`, decimal `y` expected = integerDecimal.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); @@ -517,6 +569,8 @@ tape( 'the function evaluates the exponential function (integer `x`, decimal `y` tape( 'the function evaluates the exponential function (integer `x`, integer `y`)', opts, function test( t ) { var expected; var actual; + var delta; + var tol; var x; var y; var i; @@ -526,7 +580,14 @@ tape( 'the function evaluates the exponential function (integer `x`, integer `y` expected = integerInteger.expected; for ( i = 0; i < x.length; i++ ) { actual = pow( x[i], y[i] ); - t.equal( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + if ( actual === expected[i] ) { + t.strictEqual( actual, expected[i], 'pow('+x[i]+','+y[i]+') returns '+expected[i] ); + } else { + // NOTE: the results may differ slightly from the reference and JavaScript implementations due to compiler optimizations which may be performed resulting in result divergence. For discussion, see https://github.com/stdlib-js/stdlib/pull/2298#discussion_r1624765205 + delta = abs( actual - expected[i] ); + tol = EPS * abs( expected[i] ); + t.ok( delta <= tol, 'within tolerance. pow('+x[i]+','+y[i]+'). Expected: '+expected[i]+'. Actual: '+actual+'. Delta: '+delta+'. Tol: '+tol+'.' ); + } } t.end(); }); From feb89d6768d1092b7e7d1226b32cb59f95d08d8c Mon Sep 17 00:00:00 2001 From: Philipp Burckhardt Date: Wed, 19 Jun 2024 21:04:45 +0000 Subject: [PATCH 18/18] fix: address cppcheck lint errors --- .../@stdlib/math/base/special/pow/src/main.c | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c index 6a653e2c3b82..29360167b486 100644 --- a/lib/node_modules/@stdlib/math/base/special/pow/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/pow/src/main.c @@ -297,8 +297,7 @@ static double x_is_zero( const double x, const double y ) { * double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); * // returns ~0.79 */ -static double pow2( const int32_t j, const double hp, const double lp ) { - uint32_t jcc; +static double pow2( uint32_t j, const double hp, const double lp ) { uint32_t tmp; double hpc; int32_t jc; @@ -315,8 +314,7 @@ static double pow2( const int32_t j, const double hp, const double lp ) { double z; hpc = hp; - jc = j; - jcc = (uint32_t)j; + jc = (int32_t)j; i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); k = ( ( i >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); n = 0; @@ -346,16 +344,15 @@ static double pow2( const int32_t j, const double hp, const double lp ) { 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, &jcc ); - jc = (int32_t)jcc; - jcc = jc + ( nc << HIGH_NUM_SIGNIFICAND_BITS ); - jc = (int32_t)jcc; + stdlib_base_float64_get_high_word( z, &j ); + j = j + ( nc << HIGH_NUM_SIGNIFICAND_BITS ); + jc = (int32_t)j; // Check for subnormal output... if ( ( jc >> HIGH_NUM_SIGNIFICAND_BITS ) <= 0 ) { z = stdlib_base_ldexp( z, nc ); } else { - stdlib_base_float64_set_high_word( jcc, &z ); + stdlib_base_float64_set_high_word( j, &z ); } return z; } @@ -716,6 +713,6 @@ double stdlib_base_pow( const double x, const double y ) { } } // Compute `2^(hp+lp)`... - z = pow2( jc, hp, lp ); + z = pow2( j, hp, lp ); return sx * z; }