fix: avoid cancellation near sin(x) = -1 in math/base/special/covercos and hacovercos

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

@@ -21,6 +21,7 @@
 // MODULES //
 
 var sin = require( '@stdlib/math/base/special/sin' );
+var cos = require( '@stdlib/math/base/special/cos' );
 
 
 // MAIN //
@@ -48,7 +49,16 @@ var sin = require( '@stdlib/math/base/special/sin' );
 * // returns NaN
 */
 function covercos( x ) {
-	return 1.0 + sin( x );
+	var s;
+	var c;
+
+	s = sin( x );
+	if ( s < 0.0 ) {
+		// Where `sin(x)` is close to minus one, `1 + sin(x)` cancels. Use `1 + sin(x) = cos(x)^2 / (1 - sin(x))` where `sin(x) < 0`, since `cos(x)` is computed with full relative accuracy near its zeros and the denominator does not cancel:
+		c = cos( x );
+		return ( c * c ) / ( 1.0 - s );
+	}
+	return 1.0 + s;
 }
 
 

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

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

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

@@ -18,6 +18,7 @@
 
 #include "stdlib/math/base/special/covercos.h"
 #include "stdlib/math/base/special/sin.h"
+#include "stdlib/math/base/special/cos.h"
 
 /**
 * Computes the coversed cosine.
@@ -30,5 +31,14 @@
 * // returns 2.0
 */
 double stdlib_base_covercos( const double x ) {
-	return 1.0 + stdlib_base_sin( x );
+	double s;
+	double c;
+
+	s = stdlib_base_sin( x );
+	if ( s < 0.0 ) {
+		// Where `sin(x)` is close to minus one, `1 + sin(x)` cancels. Use `1 + sin(x) = cos(x)^2 / (1 - sin(x))` where `sin(x) < 0`, since `cos(x)` is computed with full relative accuracy near its zeros and the denominator does not cancel:
+		c = stdlib_base_cos( x );
+		return ( c * c ) / ( 1.0 - s );
+	}
+	return 1.0 + s;
 }

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

@@ -19,7 +19,7 @@
 import JSON
 
 """
     gen( domain, filepath )
 
 Generate fixture data and write to file.
 
@@ -37,7 +37,9 @@
 """
 function gen( domain, filepath )
 	x = collect( domain );
-	y = 1.0 .+ sin.( x );
+	y = setprecision( BigFloat, 256 ) do
+		Float64.( 1 .+ sin.( BigFloat.( x ) ) );
+	end
 	data = Dict([
 		("x", x),
 		("expected", y)
@@ -83,3 +85,14 @@
 x = range( 2.0^60*(pi/2.0), stop = 2.0^1000*(pi/2.0), length = 4000 )
 out = joinpath( dir, "huge_positive.json" );
 gen( x, out );
+
+# Values close to where `sin(x)` is one or minus one, where the textbook formulas cancel:
+x = Float64[];
+for k in -10:10
+	for d in exp10.( range( -8.0, stop = -1.0, length = 50 ) )
+		push!( x, (2*k+0.5)*pi + d, (2*k+0.5)*pi - d, (2*k-0.5)*pi + d, (2*k-0.5)*pi - d );
+	end
+end
+out = joinpath( dir, "near_extrema.json" );
+gen( x, out );
+

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

@@ -36,6 +36,7 @@ var largeNegative = require( './fixtures/julia/large_negative.json' );
 var largePositive = require( './fixtures/julia/large_positive.json' );
 var hugeNegative = require( './fixtures/julia/huge_negative.json' );
 var hugePositive = require( './fixtures/julia/huge_positive.json' );
+var nearExtrema = require( './fixtures/julia/near_extrema.json' );
 
 
 // TESTS //
@@ -57,7 +58,7 @@ tape( 'the function computes the coversed cosine (for -256*pi < x < 0)', functio
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -73,7 +74,7 @@ tape( 'the function computes the coversed cosine (for 0 < x < 256*pi)', function
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -89,7 +90,7 @@ tape( 'the function computes the coversed cosine (for -2**60 (pi/2) < x < -2**20
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -105,7 +106,7 @@ tape( 'the function computes the coversed cosine (for 2**20 (pi/2) < x < 2**60 (
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -121,7 +122,7 @@ tape( 'the function computes the coversed cosine (for x <= -2**60 (PI/2))', func
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -137,7 +138,23 @@ tape( 'the function computes the coversed cosine (for x >= 2**60 (PI/2))', funct
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function computes the coversed cosine close to where `sin(x)` is one or minus one', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearExtrema.x;
+	expected = nearExtrema.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = covercos( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -45,6 +45,7 @@ var largeNegative = require( './fixtures/julia/large_negative.json' );
 var largePositive = require( './fixtures/julia/large_positive.json' );
 var hugeNegative = require( './fixtures/julia/huge_negative.json' );
 var hugePositive = require( './fixtures/julia/huge_positive.json' );
+var nearExtrema = require( './fixtures/julia/near_extrema.json' );
 
 
 // TESTS //
@@ -66,7 +67,7 @@ tape( 'the function computes the coversed cosine (for -256*pi < x < 0)', opts, f
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -82,7 +83,7 @@ tape( 'the function computes the coversed cosine (for 0 < x < 256*pi)', opts, fu
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -98,7 +99,7 @@ tape( 'the function computes the coversed cosine (for -2**60 (pi/2) < x < -2**20
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -114,7 +115,7 @@ tape( 'the function computes the coversed cosine (for 2**20 (pi/2) < x < 2**60 (
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -130,7 +131,7 @@ tape( 'the function computes the coversed cosine (for x <= -2**60 (PI/2))', opts
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -146,7 +147,23 @@ tape( 'the function computes the coversed cosine (for x >= 2**60 (PI/2))', opts,
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = covercos( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function computes the coversed cosine close to where `sin(x)` is one or minus one', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearExtrema.x;
+	expected = nearExtrema.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = covercos( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -21,6 +21,7 @@
 // MODULES //
 
 var sin = require( '@stdlib/math/base/special/sin' );
+var cos = require( '@stdlib/math/base/special/cos' );
 
 
 // MAIN //
@@ -48,7 +49,16 @@ var sin = require( '@stdlib/math/base/special/sin' );
 * // returns NaN
 */
 function hacovercos( x ) {
-	return (1.0 + sin( x )) / 2.0;
+	var s;
+	var c;
+
+	s = sin( x );
+	if ( s < 0.0 ) {
+		// Where `sin(x)` is close to minus one, `1 + sin(x)` cancels. Use `1 + sin(x) = cos(x)^2 / (1 - sin(x))` where `sin(x) < 0`, since `cos(x)` is computed with full relative accuracy near its zeros and the denominator does not cancel:
+		c = cos( x );
+		return ( c * c ) / ( 2.0 * ( 1.0 - s ) );
+	}
+	return ( 1.0 + s ) / 2.0;
 }
 
 

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

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

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

@@ -18,6 +18,7 @@
 
 #include "stdlib/math/base/special/hacovercos.h"
 #include "stdlib/math/base/special/sin.h"
+#include "stdlib/math/base/special/cos.h"
 
 /**
 * Computes the half-value coversed cosine.
@@ -30,5 +31,14 @@
 * // returns 1.0
 */
 double stdlib_base_hacovercos( const double x ) {
-	return ( 1.0 + stdlib_base_sin( x ) ) / 2.0;
+	double s;
+	double c;
+
+	s = stdlib_base_sin( x );
+	if ( s < 0.0 ) {
+		// Where `sin(x)` is close to minus one, `1 + sin(x)` cancels. Use `1 + sin(x) = cos(x)^2 / (1 - sin(x))` where `sin(x) < 0`, since `cos(x)` is computed with full relative accuracy near its zeros and the denominator does not cancel:
+		c = stdlib_base_cos( x );
+		return ( c * c ) / ( 2.0 * ( 1.0 - s ) );
+	}
+	return ( 1.0 + s ) / 2.0;
 }

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

@@ -19,7 +19,7 @@
 import JSON
 
 """
     gen( x, filepath )
 
 Generate fixture data and write to file.
 
@@ -36,7 +36,9 @@
 ```
 """
 function gen( x, filepath )
-	y = (1.0 .+ sin.( x )) ./ 2.0;
+	y = setprecision( BigFloat, 256 ) do
+		Float64.( (1 .+ sin.( BigFloat.( x ) )) ./ 2 );
+	end
 	data = Dict([
 		("x", x),
 		("expected", y)
@@ -82,3 +84,14 @@
 x = range( 2.0^60*(pi/2.0), stop = 2.0^1000*(pi/2.0), length = 4000 )
 out = joinpath( dir, "huge_positive.json" );
 gen( x, out );
+
+# Values close to where `sin(x)` is one or minus one, where the textbook formulas cancel:
+x = Float64[];
+for k in -10:10
+	for d in exp10.( range( -8.0, stop = -1.0, length = 50 ) )
+		push!( x, (2*k+0.5)*pi + d, (2*k+0.5)*pi - d, (2*k-0.5)*pi + d, (2*k-0.5)*pi - d );
+	end
+end
+out = joinpath( dir, "near_extrema.json" );
+gen( x, out );
+