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..28d0d3c55a93 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. + +```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: + +- **base**: `[in] double` base. +- **exponent**: `[in] double` exponent. + +```c +double stdlib_base_pow( const double base, const double exponent ); +``` + +
+ + + + + +
+ +
+ + + + + +
+ +### Examples + +```c +#include "stdlib/math/base/special/pow.h" +#include +#include + +int main( void ) { + 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 ); + } +} +``` + +
+ + + +
+ + +
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..1713caa183dd --- /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 = ( randu()*10.0 ); + x = ( 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..880a68893fa1 --- /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 `pow`. +*/ +#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..ec3992233442 --- /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 +} 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..6aed70daf167 --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/Makefile @@ -0,0 +1,146 @@ +#/ +# @license Apache-2.0 +# +# Copyright (c) 2024 The Stdlib Authors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +#/ + +# VARIABLES # + +ifndef VERBOSE + QUIET := @ +else + QUIET := +endif + +# Determine the OS ([1][1], [2][2]). +# +# [1]: https://en.wikipedia.org/wiki/Uname#Examples +# [2]: http://stackoverflow.com/a/27776822/2225624 +OS ?= $(shell uname) +ifneq (, $(findstring MINGW,$(OS))) + OS := WINNT +else +ifneq (, $(findstring MSYS,$(OS))) + OS := WINNT +else +ifneq (, $(findstring CYGWIN,$(OS))) + OS := WINNT +else +ifneq (, $(findstring Windows_NT,$(OS))) + OS := WINNT +endif +endif +endif +endif + +# Define the program used for compiling C source files: +ifdef C_COMPILER + CC := $(C_COMPILER) +else + CC := gcc +endif + +# Define the command-line options when compiling C files: +CFLAGS ?= \ + -std=c99 \ + -O3 \ + -Wall \ + -pedantic + +# Determine whether to generate position independent code ([1][1], [2][2]). +# +# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options +# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option +ifeq ($(OS), WINNT) + fPIC ?= +else + fPIC ?= -fPIC +endif + +# List of includes (e.g., `-I /foo/bar -I /beep/boop/include`): +INCLUDE ?= + +# List of source files: +SOURCE_FILES ?= + +# List of libraries (e.g., `-lopenblas -lpthread`): +LIBRARIES ?= + +# List of library paths (e.g., `-L /foo/bar -L /beep/boop`): +LIBPATH ?= + +# List of C targets: +c_targets := example.out + + +# RULES # + +#/ +# Compiles source files. +# +# @param {string} [C_COMPILER] - C compiler (e.g., `gcc`) +# @param {string} [CFLAGS] - C compiler options +# @param {(string|void)} [fPIC] - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} [INCLUDE] - list of includes (e.g., `-I /foo/bar -I /beep/boop/include`) +# @param {string} [SOURCE_FILES] - list of source files +# @param {string} [LIBPATH] - list of library paths (e.g., `-L /foo/bar -L /beep/boop`) +# @param {string} [LIBRARIES] - list of libraries (e.g., `-lopenblas -lpthread`) +# +# @example +# make +# +# @example +# make all +#/ +all: $(c_targets) + +.PHONY: all + +#/ +# Compiles C source files. +# +# @private +# @param {string} CC - C compiler (e.g., `gcc`) +# @param {string} CFLAGS - C compiler options +# @param {(string|void)} fPIC - compiler flag determining whether to generate position independent code (e.g., `-fPIC`) +# @param {string} INCLUDE - list of includes (e.g., `-I /foo/bar`) +# @param {string} SOURCE_FILES - list of source files +# @param {string} LIBPATH - list of library paths (e.g., `-L /foo/bar`) +# @param {string} LIBRARIES - list of libraries (e.g., `-lopenblas`) +#/ +$(c_targets): %.out: %.c + $(QUIET) $(CC) $(CFLAGS) $(fPIC) $(INCLUDE) -o $@ $(SOURCE_FILES) $< $(LIBPATH) -lm $(LIBRARIES) + +#/ +# Runs compiled examples. +# +# @example +# make run +#/ +run: $(c_targets) + $(QUIET) ./$< + +.PHONY: run + +#/ +# Removes generated files. +# +# @example +# make clean +#/ +clean: + $(QUIET) -rm -f *.o *.out + +.PHONY: clean diff --git a/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c b/lib/node_modules/@stdlib/math/base/special/pow/examples/c/example.c new file mode 100644 index 000000000000..ee26d99bea68 --- /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 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 ); + } +} 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..575cb043c0bf --- /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 +static const int32_t HIGH_MAX_NEAR_UNITY = 0x3fefffff; + +static const double HUGE_VALUE = 1.0e300; + +static const double TINY_VALUE = 1.0e-300; + +// 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 +static const int32_t HIGH_MIN_NORMAL_EXP = 0x00100000; + +// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 +static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; + +// TODO: consider making into an external constant +static const int32_t HIGH_NUM_SIGNIFICAND_BITS = 20; + +// High: LN2 +static const double LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000 + +// Low: LN2 +static const double LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39 + +// 1/LN2 +static const double INV_LN2 = 1.44269504088896338700e+00; // 0x3FF71547, 0x652B82FE + +// High (24 bits): 1/LN2 +static const double INV_LN2_HI = 1.44269502162933349609e+00; // 0x3FF71547, 0x60000000 + +// Low: 1/LN2 +static const double INV_LN2_LO = 1.92596299112661746887e-08; // 0x3E54AE0B, 0xF85DDF44 + +// 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 +static const int32_t HIGH_SIGNIFICAND_MASK = 0x000fffff; + +// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 +static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; + +// 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 +static const int32_t HIGH_BIASED_EXP_NEG_512 = 0x20000000; + +// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 +static const int32_t HIGH_SIGNIFICAND_HALF = 0x00080000; + +static const double TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 + +// 2/(3*LN2) +static const double CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD + +// (float)CP +static const double CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 + +// Low: CP_HI +static const double CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 + +static const double BP [2] = { 1.0, 1.5 }; + +static const double DP_HI [2] = { 0.0, 5.84962487220764160156e-01 }; // 0x3FE2B803, 0x40000000 + +static const double DP_LO [2] = { 0.0, 1.35003920212974897128e-08 }; // 0x3E4CFDEB, 0x43CFD006 + +// 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31 +static const int32_t HIGH_BIASED_EXP_31 = 0x41e00000; + +// 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64 +static const int32_t HIGH_BIASED_EXP_64 = 0x43f00000; + +// 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024 +static const int32_t HIGH_BIASED_EXP_10 = 0x40900000; + +// 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000 +static const int32_t HIGH_1075 = 0x4090cc00; + +// 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000 +static const uint32_t HIGH_NEG_1075 = 0xc090cc00; + +static const int32_t HIGH_NUM_NONSIGN_BITS = 31; + +// -(1024-log2(ovfl+.5ulp)) +static const double OVT = 8.0085662595372944372e-17; + +/* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */ + +// BEGIN: polyval_l + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_l( const double x ) { + return 0.5999999999999946 + (x * (0.4285714285785502 + (x * (0.33333332981837743 + (x * (0.272728123808534 + (x * (0.23066074577556175 + (x * 0.20697501780033842))))))))); +} + +// END: polyval_l + +// BEGIN: polyval_w + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_w( const double x ) { + return 0.5 + (x * (-0.3333333333333333 + (x * 0.25))); +} + +// END: polyval_w + +// BEGIN: polyval_p + +/** +* Evaluates a polynomial. +* +* ## Notes +* +* - The implementation uses [Horner's rule][horners-method] for efficient computation. +* +* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method +* +* @param x value at which to evaluate the polynomial +* @return evaluated polynomial +*/ +static double polyval_p( const double x ) { + return 0.16666666666666602 + (x * (-0.0027777777777015593 + (x * (0.00006613756321437934 + (x * (-0.0000016533902205465252 + (x * 4.1381367970572385e-8))))))); +} + +// END: polyval_p + +/* End auto-generated functions. */ + +/** +* Evaluates the exponential function when \\( y = \pm \infty\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_infinite( -1.0, 1.0/0.0 ); +* // returns NaN +* +* @example +* double out = y_is_infinite( 1.5, -1.0/0.0 ); +* // returns 0.0 +*/ +static double y_is_infinite( const double x, const double y ) { + if ( x == -1.0 ) { + return 0.0 / 0.0; // NaN + } + if ( x == 1.0 ) { + return 1.0; + } + // (|x| > 1 && y == NINF) || (|x| < 1 && y === PINF) + if ( ( stdlib_base_abs( x ) < 1.0 ) == ( y == STDLIB_CONSTANT_FLOAT64_PINF ) ) { + return 0.0; + } + // (|x| > 1 && y == PINF) || (|x| < 1 && y == NINF) + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Evaluates the exponential function when \\(|y| > 2^64\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = y_is_huge( 9.0, 3.6893488147419103e19 ); +* // returns Infinity +* +* @example +* double out = y_is_huge( -3.14, -3.6893488147419103e19 ); +* // returns 0.0 +*/ +static double y_is_huge( const double x, const double y ) { + uint32_t ahx; + uint32_t hx; + + stdlib_base_float64_get_high_word( x, &hx ); + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + if ( ahx <= HIGH_MAX_NEAR_UNITY ) { + if ( y < 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; + } + // `x` has a biased exponent greater than or equal to `0`... + + if ( y > 0 ) { + // Signal overflow... + return HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return TINY_VALUE * TINY_VALUE; +} + +/** +* Evaluates the exponential function when \\(|x| = 0\\). +* +* @param x base +* @param y exponent +* @return output value +* +* @example +* double out = x_is_zero( 0.0, 2 ); +* // returns 0.0 +* +* @example +* double out = x_is_zero( 0.0, -9 ); +* // returns Infinity +*/ +static double x_is_zero( const double x, const double y ) { + if ( y == STDLIB_CONSTANT_FLOAT64_NINF ) { + return STDLIB_CONSTANT_FLOAT64_PINF; + } + if ( y == STDLIB_CONSTANT_FLOAT64_PINF ) { + return 0.0; + } + if ( y > 0.0 ) { + if ( stdlib_base_is_odd( y ) ) { + return x; // handles +-0 + } + return 0.0; + } + // y < 0.0 + if ( stdlib_base_is_odd( y ) ) { + return stdlib_base_copysign( STDLIB_CONSTANT_FLOAT64_PINF, x ); // handles +-0 + } + return STDLIB_CONSTANT_FLOAT64_PINF; +} + +/** +* Computes \\(2^{\mathrm{hp} + \mathrm{lp}\\). +* +* @param j high word of `hp + lp` +* @param hp first power summand +* @param lp second power summand +* @return function value +* +* @example +* double out = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 ); +* // returns ~0.79 +*/ +static double pow2( uint32_t j, const double hp, const double lp ) { + uint32_t tmp; + double hpc; + int32_t jc; + int32_t nc; + int32_t i; + int32_t k; + int32_t n; + double t1; + double t; + double r; + double u; + double v; + double w; + double z; + + hpc = hp; + jc = (int32_t)j; + i = ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + k = ( ( i >> 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 ) ) ); + 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 ) ) ); + 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 ) ); + nc = (int32_t)n; + if ( jc < 0 ) { + nc = -nc; + } + hpc -= t; + } + t = lp + hpc; + stdlib_base_float64_set_low_word( 0, &t ); + u = t * LN2_HI; + v = ( ( lp - ( t - hpc ) ) * STDLIB_CONSTANT_FLOAT64_LN2 ) + ( t * LN2_LO ); + z = u + v; + w = v - ( z - u ); + t = z * z; + t1 = z - ( t * polyval_p( t ) ); + r = ( ( z * t1 ) / ( t1 - 2.0 ) ) - ( w + ( z * w ) ); + z = 1.0 - ( r - z ); + stdlib_base_float64_get_high_word( z, &j ); + j = j + ( nc << 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( j, &z ); + } + return z; +} + +/** +* Computes \\(\operatorname{log}(x)\\) assuming \\(|1-x|\\) is small and using the approximation \\(x - x^2/2 + x^3/3 - x^4/4\\). +* +* @param ax absolute value of x +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* logx( 9.0, &o1, &o2 ); +*/ +void logx( const double ax, double *o1, double *o2 ) { + double t2; + double t1; + double t; + double w; + double u; + double v; + + t = ax - 1.0; // `t` has `20` trailing zeros + w = t * t * polyval_w( t ); + u = INV_LN2_HI * t; // `INV_LN2_HI` has `21` significant bits + v = ( t * INV_LN2_LO ) - ( w * INV_LN2 ); + t1 = u + v; + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = v - (t1 - u); + + *o1 = t1; + *o2 = t2; +} + +/** +* Computes \\(\operatorname{log2}(ax)\\). +* +* @param ax absolute value of x +* @param ahx high word of ax +* @param o1 destination for output1 +* @param o2 destination for output2 +* +* @example +* log2ax( 9.0, 1075970048, &o1, &o2 ); +*/ +void log2ax( const double ax, const int32_t ahx, double *o1, double *o2 ) { + uint32_t ahxcc; + uint32_t tmp; + int32_t ahxc; + double axc; + int32_t n; + int32_t j; + int32_t k; + double ss; // hs + ls + double s2; // ss squared + double hs; + double ls; + double ht; + double lt; + double bp; // BP constant + double dp; // DP constant + double hp; + double lp; + double hz; + double lz; + double t1; + double t2; + double t; + double r; + double u; + double v; + + n = 0; + axc = ax; + ahxc = ahx; + ahxcc = (uint32_t)ahx; + // Check if `x` is subnormal... + if ( ahxc < HIGH_MIN_NORMAL_EXP ) { + axc *= TWO53; + n -= 53; + stdlib_base_float64_get_high_word( axc, &ahxcc ); + ahxc = (int32_t)ahxcc; + } + // Extract the unbiased exponent of `x`: + n += ( ( ahxc >> HIGH_NUM_SIGNIFICAND_BITS ) - STDLIB_CONSTANT_FLOAT64_EXPONENT_BIAS ); + + // Isolate the significand bits of `x`: + j = ( ahxc & HIGH_SIGNIFICAND_MASK ); + + // Normalize `ahx` by setting the (biased) exponent to `1023`: + ahxc = ( j|HIGH_BIASED_EXP_0 ); + + // Determine the interval of `|x|` by comparing significand bits... + + // |x| < sqrt(3/2) + if ( j <= 0x3988E ) { // 0 00000000000 00111001100010001110 + k = 0; + } + // |x| < sqrt(3) + else if ( j < 0xBB67A ) { // 0 00000000000 10111011011001111010 + k = 1; + } + // |x| >= sqrt(3) + else { + k = 0; + n += 1; + ahxc -= HIGH_MIN_NORMAL_EXP; + } + ahxcc = (uint32_t)ahxc; + + // Load the normalized high word into `|x|`: + stdlib_base_float64_set_high_word( ahxcc, &axc ); + + // Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: + bp = BP[ k ]; // BP[0] = 1.0, BP[1] = 1.5 + u = axc - bp; // (x-1) || (x-1.5) + v = 1.0 / (axc + bp); // 1/(x+1) || 1/(x+1.5) + ss = u * v; + hs = ss; + stdlib_base_float64_set_low_word( 0, &hs ); // set all low word (less significant significand) bits to 0s + + // Compute `ht = ax + bp` (via manipulation, i.e., bit flipping, of the high word): + tmp = ( ( ahxc >> 1 ) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; + tmp += ( k << 18 ); // `(k<<18)` can be considered the word equivalent of `1.0` or `1.5` + ht = 0.0; + stdlib_base_float64_set_high_word( tmp, &ht ); + lt = axc - ( ht - bp ); + ls = v * ( ( u - ( hs * ht ) ) - ( hs * lt ) ); + + // Compute `log(ax)`... + + s2 = ss * ss; + r = s2 * s2 * polyval_l( s2 ); + r += ls * (hs + ss); + s2 = hs * hs; + ht = 3.0 + s2 + r; + stdlib_base_float64_set_low_word( 0, &ht ); + lt = r - ( ( ht - 3.0 ) - s2 ); + + // u+v = ss*(1+...): + u = hs * ht; + v = ( ls * ht ) + ( lt * ss ); + + // 2/(3LN2) * (ss+...): + hp = u + v; + stdlib_base_float64_set_low_word( 0, &hp ); + lp = v - ( hp - u ); + hz = CP_HI * hp; // CP_HI+CP_LO = 2/(3*LN2) + lz = ( CP_LO * hp ) + ( lp*CP ) + DP_LO[ k ]; + + // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz + dp = DP_HI[ k ]; + t = n; + t1 = ( ( hz + lz ) + dp ) + t; // log2(ax) + stdlib_base_float64_set_low_word( 0, &t1 ); + t2 = lz - ( ( ( t1 - t ) - dp ) - hz ); + + *o1 = t1; + *o2 = t2; +} + +/** +* Evaluates the exponential function. +* +* @param x base +* @param y exponent +* @return function value +* +* @example +* double out = stdlib_base_pow( 2.0, 3.0 ); +* // returns 8.0 +* +* @example +* double out = stdlib_base_pow( 4.0, 0.5 ); +* // returns 2.0 +*/ +double stdlib_base_pow( const double x, const double y ) { + uint32_t hx; // high word `x` + uint32_t lx; // low word `x` + uint32_t hy; // high word `y` + uint32_t ly; // low word `y` + int32_t ahx; // absolute value high word `x` + int32_t ahy; // absolute value high word `y` + uint32_t i; + uint32_t j; + int32_t ic; + int32_t jc; + int32_t sx; // sign `x` + int32_t sy; // sign `y` + double ax; // absolute value `x` + double y1; + double hp; + double lp; + double t1; + double t2; + double xc; + double z; // y prime + + xc = x; + if ( stdlib_base_is_nan( xc ) || stdlib_base_is_nan( y ) ) { + return 0.0/0.0; // NaN + } + + // Split `y` into high and low words: + stdlib_base_float64_to_words( y, &hy, &ly ); + + // Special cases `y`... + if ( ly == 0 ) { + if ( y == 0.0 ) { + return 1.0; + } + if ( y == 1.0 ) { + return xc; + } + if ( y == -1.0 ) { + return 1.0 / xc; + } + if ( y == 0.5 ) { + return stdlib_base_sqrt( xc ); + } + if ( y == -0.5 ) { + return 1.0 / stdlib_base_sqrt( xc ); + } + if ( y == 2.0 ) { + return xc * xc; + } + if ( y == 3.0 ) { + return xc * xc * xc; + } + if ( y == 4.0 ) { + xc *= xc; + return xc * xc; + } + if ( stdlib_base_is_infinite( y ) ) { + return y_is_infinite( xc, y ); + } + } + // Split `x` into high and low words: + stdlib_base_float64_to_words( xc, &hx, &lx ); + + // Special cases `x`... + if ( lx == 0 ) { + if ( hx == 0 ) { + return x_is_zero( xc, y ); + } + if ( xc == 1.0 ) { + return 1.0; + } + if ( + xc == -1.0 && + stdlib_base_is_odd( y ) + ) { + return -1.0; + } + if ( stdlib_base_is_infinite( xc ) ) { + if ( xc == STDLIB_CONSTANT_FLOAT64_NINF ) { + // `stdlib_base_pow( 1/x, -y )` + return stdlib_base_pow( -0.0, -y ); + } + if ( y < 0.0 ) { + return 0.0; + } + return STDLIB_CONSTANT_FLOAT64_PINF; + } + } + if ( + xc < 0.0 && + stdlib_base_is_integer( y ) == false + ) { + // Signal NaN... + return 0.0/0.0; + } + ax = stdlib_base_abs( x ); + + // Remove the sign bits (i.e., get absolute values): + ahx = ( hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + ahy = ( hy & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ); + + // Extract the sign bits: + sx = ( hx >> HIGH_NUM_NONSIGN_BITS ); + sy = ( hy >> HIGH_NUM_NONSIGN_BITS ); + + // Determine the sign of the result... + if ( sx && stdlib_base_is_odd( y ) ) { + sx = -1; + } else { + sx = 1; + } + // Case 1: `|y|` is huge... + + // |y| > 2^31 + if ( ahy > HIGH_BIASED_EXP_31 ) { + // `|y| > 2^64`, then must over- or underflow... + if ( ahy > HIGH_BIASED_EXP_64 ) { + return y_is_huge( xc, y ); + } + // Over- or underflow if `x` is not close to unity... + if ( ahx < HIGH_MAX_NEAR_UNITY ) { + // y < 0 + if ( sy == 1 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( ahx > HIGH_BIASED_EXP_0 ) { + // y > 0 + if ( sy == 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + // At this point, `|1-x|` is tiny (`<= 2^-20`). Suffice to compute `log(x)` by `x - x^2/2 + x^3/3 - x^4/4`. + logx( ax, &t1, &t2 ); + } + // Case 2: `|y|` is not huge... + else { + log2ax( ax, ahx, &t1, &t2 ); + } + + // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`... + y1 = y; + stdlib_base_float64_set_low_word( 0, &y1 ); + lp = ( ( y - y1 ) * t1 ) + ( y * t2 ); + hp = y1 * t1; + z = lp + hp; + + // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner. + stdlib_base_float64_to_words( z, &j, &i ); + jc = (int32_t)j; + ic = (int32_t)i; + + // z >= 1024 + if ( jc >= HIGH_BIASED_EXP_10 ) { + // z > 1024 + if ( ( ( jc - HIGH_BIASED_EXP_10 )|ic ) != 0 ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + if ( ( lp + OVT ) > ( z - hp ) ) { + // Signal overflow... + return sx * HUGE_VALUE * HUGE_VALUE; + } + } + // z <= -1075 + else if ( ( j & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK ) >= HIGH_1075 ) { + // z < -1075 + if ( ( ( j - HIGH_NEG_1075 )|i ) != 0 ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + if ( lp <= ( z - hp ) ) { + // Signal underflow... + return sx * TINY_VALUE * TINY_VALUE; + } + } + // Compute `2^(hp+lp)`... + z = pow2( j, hp, lp ); + return sx * z; +} 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..9fafa645e19d --- /dev/null +++ b/lib/node_modules/@stdlib/math/base/special/pow/test/test.native.js @@ -0,0 +1,1091 @@ +/** +* @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 delta; + var tol; + 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] ); + 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(); +}); + +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 delta; + var tol; + 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 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(); +}); + +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 delta; + var tol; + 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] ); + 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(); +}); + +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; + + x = decimalDecimal.x; + y = decimalDecimal.y; + expected = decimalDecimal.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[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(); +}); + +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; + + x = decimalInteger.x; + y = decimalInteger.y; + expected = decimalInteger.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[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(); +}); + +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; + + x = integerDecimal.x; + y = integerDecimal.y; + expected = integerDecimal.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[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(); +}); + +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; + + x = integerInteger.x; + y = integerInteger.y; + expected = integerInteger.expected; + for ( i = 0; i < x.length; i++ ) { + actual = pow( x[i], y[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(); +}); + +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(); +});