fix: avoid cancellation in math/base/special/vercos and havercos

lib/node_modules/@stdlib/math/base/special/havercos/lib/main.js

@@ -20,6 +20,7 @@
 
 // MODULES //
 
+var abs2 = require( '@stdlib/math/base/special/abs2' );
 var cos = require( '@stdlib/math/base/special/cos' );
 
 
@@ -28,6 +29,22 @@ var cos = require( '@stdlib/math/base/special/cos' );
 /**
 * Computes the half-value versed cosine.
 *
+* ## Method
+*
+* Evaluating the textbook formula
+*
+* ```tex
+* \operatorname{havercos}(x) = \frac{1 + \cos(x)}{2}
+* ```
+*
+* loses significant digits when \\(\cos(x)\\) is close to negative one (i.e., for `x` near odd multiples of π). Accordingly, we use the equivalent half-angle identity
+*
+* ```tex
+* \operatorname{havercos}(x) = \cos^2\left(\tfrac{x}{2}\right)
+* ```
+*
+* which does not suffer from cancellation.
+*
 * @param {number} x - input value (in radians)
 * @returns {number} half-value versed cosine
 *
@@ -48,7 +65,7 @@ var cos = require( '@stdlib/math/base/special/cos' );
 * // returns NaN
 */
 function havercos( x ) {
-	return (1.0 + cos( x )) / 2.0;
+	return abs2( cos( x / 2.0 ) );
 }
 
 

lib/node_modules/@stdlib/math/base/special/havercos/manifest.json

@@ -37,6 +37,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/napi/unary",
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     },
@@ -51,6 +52,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     },
@@ -65,6 +67,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     }

lib/node_modules/@stdlib/math/base/special/havercos/src/main.c

@@ -17,11 +17,28 @@
 */
 
 #include "stdlib/math/base/special/havercos.h"
+#include "stdlib/math/base/special/abs2.h"
 #include "stdlib/math/base/special/cos.h"
 
 /**
 * Computes the half-value versed cosine.
 *
+* ## Method
+*
+* Evaluating the textbook formula
+*
+* ```tex
+* \operatorname{havercos}(x) = \frac{1 + \cos(x)}{2}
+* ```
+*
+* loses significant digits when `cos(x)` is close to negative one (i.e., for `x` near odd multiples of π). Accordingly, we use the equivalent half-angle identity
+*
+* ```tex
+* \operatorname{havercos}(x) = \cos^2\left(\tfrac{x}{2}\right)
+* ```
+*
+* which does not suffer from cancellation.
+*
 * @param x    input value (in radians)
 * @return     half-value versed cosine
 *
@@ -30,5 +47,5 @@
 * // returns 0.5
 */
 double stdlib_base_havercos( const double x ) {
-	return ( 1.0 + stdlib_base_cos( x ) ) / 2.0;
+	return stdlib_base_abs2( stdlib_base_cos( x / 2.0 ) );
 }

lib/node_modules/@stdlib/math/base/special/havercos/test/fixtures/julia/runner.jl

@@ -19,7 +19,7 @@
 import JSON
 
 """
-    gen( domain, filepath )
+	gen( domain, filepath )
 
 Generate fixture data and write to file.
 
@@ -37,7 +37,12 @@ julia> gen( x, \"./data.json\" );
 """
 function gen( domain, filepath )
 	x = collect( domain );
-	y = (1.0 .+ cos.( x )) ./ 2.0;
+
+	# Evaluate the textbook formula in extended precision (4096 bits resolve the cancellation even for `x` closest to an odd multiple of pi), and round the result to double precision:
+	y = setprecision( BigFloat, 4096 ) do
+		Float64.( ( 1 .+ cos.( big.( x ) ) ) ./ 2 )
+	end
+
 	data = Dict([
 		("x", x),
 		("expected", y)
@@ -54,6 +59,11 @@ file = @__FILE__;
 # Extract the directory in which this file resides:
 dir = dirname( file );
 
+# Values near odd multiples of pi, where `1 + cos(x)` is close to zero:
+x = vec( [ (2*k+1)*pi + s*d for d in exp10.( range( -12.0, stop = -1.0, length = 100 ) ), k in -10:9, s in (-1, 1) ] )
+out = joinpath( dir, "near_odd_pi.json" );
+gen( x, out );
+
 # Negative medium sized values:
 x = range( -256.0*pi, stop = 0.0, length = 4000 )
 out = joinpath( dir, "medium_negative.json" );

lib/node_modules/@stdlib/math/base/special/havercos/test/test.js

@@ -21,6 +21,7 @@
 // MODULES //
 
 var tape = require( 'tape' );
+var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
@@ -29,6 +30,7 @@ var havercos = require( './../lib' );
 
 // FIXTURES //
 
+var nearOddPi = require( './fixtures/julia/near_odd_pi.json' );
 var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var mediumPositive = require( './fixtures/julia/medium_positive.json' );
 var largeNegative = require( './fixtures/julia/large_negative.json' );
@@ -45,6 +47,23 @@ tape( 'main export is a function', function test( t ) {
 	t.end();
 });
 
+tape( 'the function computes the half-value versed cosine (for x near odd multiples of pi)', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearOddPi.x;
+	expected = nearOddPi.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = havercos( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
 tape( 'the function computes the half-value versed cosine (for -256*pi < x < 0 )', function test( t ) {
 	var expected;
 	var x;
@@ -56,7 +75,7 @@ tape( 'the function computes the half-value versed cosine (for -256*pi < x < 0 )
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();
@@ -73,7 +92,7 @@ tape( 'the function computes the half-value versed cosine (for 0 < x < 256*pi )'
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();
@@ -90,7 +109,7 @@ tape( 'the function computes the half-value versed cosine (for -2**60 (pi/2) < x
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();
@@ -107,7 +126,7 @@ tape( 'the function computes the half-value versed cosine (for 2**20 (pi/2) < x
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();
@@ -124,7 +143,7 @@ tape( 'the function computes the half-value versed cosine (for x <= -2**60 (PI/2
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();
@@ -141,7 +160,7 @@ tape( 'the function computes the half-value versed cosine (for x >= 2**60 (PI/2)
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = havercos( x[i] );
-		t.strictEqual( y, expected[ i ], 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 
 	t.end();

lib/node_modules/@stdlib/math/base/special/havercos/test/test.native.js

@@ -39,6 +39,7 @@ var opts = {
 
 // FIXTURES //
 
+var nearOddPi = require( './fixtures/julia/near_odd_pi.json' );
 var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var mediumPositive = require( './fixtures/julia/medium_positive.json' );
 var largeNegative = require( './fixtures/julia/large_negative.json' );
@@ -55,6 +56,23 @@ tape( 'main export is a function', opts, function test( t ) {
 	t.end();
 });
 
+tape( 'the function computes the half-value versed cosine (for x near odd multiples of pi)', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearOddPi.x;
+	expected = nearOddPi.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = havercos( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
 tape( 'the function computes the half-value versed cosine (for -256*pi < x < 0 )', opts, function test( t ) {
 	var expected;
 	var x;

lib/node_modules/@stdlib/math/base/special/vercos/lib/main.js

@@ -20,6 +20,7 @@
 
 // MODULES //
 
+var abs2 = require( '@stdlib/math/base/special/abs2' );
 var cos = require( '@stdlib/math/base/special/cos' );
 
 
@@ -28,6 +29,22 @@ var cos = require( '@stdlib/math/base/special/cos' );
 /**
 * Computes the versed cosine.
 *
+* ## Method
+*
+* Evaluating the textbook formula
+*
+* ```tex
+* \operatorname{vercos}(x) = 1 + \cos(x)
+* ```
+*
+* loses significant digits when \\(\cos(x)\\) is close to negative one (i.e., for `x` near odd multiples of π). Accordingly, we use the equivalent half-angle identity
+*
+* ```tex
+* \operatorname{vercos}(x) = 2 \cos^2\left(\tfrac{x}{2}\right)
+* ```
+*
+* which does not suffer from cancellation.
+*
 * @param {number} x - input value (in radians)
 * @returns {number} versed cosine
 *
@@ -48,7 +65,7 @@ var cos = require( '@stdlib/math/base/special/cos' );
 * // returns NaN
 */
 function vercos( x ) {
-	return 1.0 + cos( x );
+	return 2.0 * abs2( cos( x / 2.0 ) );
 }
 
 

lib/node_modules/@stdlib/math/base/special/vercos/manifest.json

@@ -37,6 +37,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/napi/unary",
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     },
@@ -51,6 +52,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     },
@@ -65,6 +67,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
+        "@stdlib/math/base/special/abs2",
         "@stdlib/math/base/special/cos"
       ]
     }

lib/node_modules/@stdlib/math/base/special/vercos/src/main.c

@@ -17,11 +17,28 @@
 */
 
 #include "stdlib/math/base/special/vercos.h"
+#include "stdlib/math/base/special/abs2.h"
 #include "stdlib/math/base/special/cos.h"
 
 /**
 * Computes the versed cosine.
 *
+* ## Method
+*
+* Evaluating the textbook formula
+*
+* ```tex
+* \operatorname{vercos}(x) = 1 + \cos(x)
+* ```
+*
+* loses significant digits when `cos(x)` is close to negative one (i.e., for `x` near odd multiples of π). Accordingly, we use the equivalent half-angle identity
+*
+* ```tex
+* \operatorname{vercos}(x) = 2 \cos^2\left(\tfrac{x}{2}\right)
+* ```
+*
+* which does not suffer from cancellation.
+*
 * @param x    input value (in radians)
 * @return     versed cosine
 *
@@ -30,5 +47,5 @@
 * // returns 1.0
 */
 double stdlib_base_vercos( const double x ) {
-	return 1.0 + stdlib_base_cos( x );
+	return 2.0 * stdlib_base_abs2( stdlib_base_cos( x / 2.0 ) );
 }

lib/node_modules/@stdlib/math/base/special/vercos/test/fixtures/julia/runner.jl

@@ -19,7 +19,7 @@
 import JSON
 
 """
-    gen( domain, filepath )
+	gen( domain, filepath )
 
 Generate fixture data and write to file.
 
@@ -37,7 +37,12 @@ julia> gen( x, \"./data.json\" );
 """
 function gen( domain, filepath )
 	x = collect( domain );
-	y = 1.0 .+ cos.( x );
+
+	# Evaluate the textbook formula in extended precision (4096 bits resolve the cancellation even for `x` closest to an odd multiple of pi), and round the result to double precision:
+	y = setprecision( BigFloat, 4096 ) do
+		Float64.( 1 .+ cos.( big.( x ) ) )
+	end
+
 	data = Dict([
 		("x", x),
 		("expected", y)
@@ -54,6 +59,11 @@ file = @__FILE__;
 # Extract the directory in which this file resides:
 dir = dirname( file );
 
+# Values near odd multiples of pi, where `1 + cos(x)` is close to zero:
+x = vec( [ (2*k+1)*pi + s*d for d in exp10.( range( -12.0, stop = -1.0, length = 100 ) ), k in -10:9, s in (-1, 1) ] )
+out = joinpath( dir, "near_odd_pi.json" );
+gen( x, out );
+
 # Negative medium sized values:
 x = range( -256.0*pi, stop = 0.0, length = 4000 )
 out = joinpath( dir, "medium_negative.json" );