fix: avoid cancellation near sin(x) = 1 in math/base/special/coversin and hacoversin

lib/node_modules/@stdlib/math/base/special/coversin/README.md

@@ -60,7 +60,7 @@ var v = coversin( 0.0 );
 // returns 1.0
 
 v = coversin( 3.141592653589793/2.0 );
-// returns 0.0
+// returns ~0.0
 
 v = coversin( -3.141592653589793/6.0 );
 // returns 1.5
@@ -129,7 +129,7 @@ double out = stdlib_base_coversin( 0.0 );
 // returns 1.0
 
 out = stdlib_base_coversin( 3.141592653589793 / 2.0 );
-// returns 0.0
+// returns ~0.0
 ```
 
 The function accepts the following arguments:

lib/node_modules/@stdlib/math/base/special/coversin/docs/types/index.d.ts

@@ -30,7 +30,7 @@
 *
 * @example
 * var v = coversin( 3.141592653589793/2.0 );
-* // returns 0.0
+* // returns ~0.0
 *
 * @example
 * var v = coversin( -3.141592653589793/6.0 );

lib/node_modules/@stdlib/math/base/special/coversin/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 //
@@ -37,7 +38,7 @@ var sin = require( '@stdlib/math/base/special/sin' );
 *
 * @example
 * var v = coversin( 3.141592653589793/2.0 );
-* // returns 0.0
+* // returns ~0.0
 *
 * @example
 * var v = coversin( -3.141592653589793/6.0 );
@@ -48,7 +49,16 @@ var sin = require( '@stdlib/math/base/special/sin' );
 * // returns NaN
 */
 function coversin( x ) {
-	return 1.0 - sin( x );
+	var s;
+	var c;
+
+	s = sin( x );
+	if ( s > 0.0 ) {
+		// Where `sin(x)` is close to 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/coversin/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/coversin/src/main.c

@@ -18,6 +18,7 @@
 
 #include "stdlib/math/base/special/coversin.h"
 #include "stdlib/math/base/special/sin.h"
+#include "stdlib/math/base/special/cos.h"
 
 /**
 * Computes the coversed sine.
@@ -27,8 +28,17 @@
 *
 * @example
 * double y = stdlib_base_coversin( 3.141592653589793 / 2.0 );
-* // returns 0.0
+* // returns ~0.0
 */
 double stdlib_base_coversin( 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 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/coversin/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/coversin/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 coversine (for -256*pi < x < 0)', function test
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -73,7 +74,7 @@ tape( 'the function computes the coversine (for 0 < x < 256*pi)', function test(
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -89,7 +90,7 @@ tape( 'the function computes the coversine (for -2**60 (pi/2) < x < -2**20 (pi/2
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -105,7 +106,7 @@ tape( 'the function computes the coversine (for 2**20 (pi/2) < x < 2**60 (pi/2))
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -121,7 +122,7 @@ tape( 'the function computes the coversine (for x <= -2**60 (PI/2))', function t
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -137,7 +138,23 @@ tape( 'the function computes the coversine (for x >= 2**60 (PI/2))', function te
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
+});
+
+tape( 'the function computes the coversine 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 = coversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

lib/node_modules/@stdlib/math/base/special/coversin/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 coversine (for -256*pi < x < 0)', opts, functio
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -82,7 +83,7 @@ tape( 'the function computes the coversine (for 0 < x < 256*pi)', opts, function
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -98,7 +99,7 @@ tape( 'the function computes the coversine (for -2**60 (pi/2) < x < -2**20 (pi/2
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -114,7 +115,7 @@ tape( 'the function computes the coversine (for 2**20 (pi/2) < x < 2**60 (pi/2))
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -130,7 +131,7 @@ tape( 'the function computes the coversine (for x <= -2**60 (PI/2))', opts, func
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
 });
@@ -146,7 +147,23 @@ tape( 'the function computes the coversine (for x >= 2**60 (PI/2))', opts, funct
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = coversin( 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();
+});
+
+tape( 'the function computes the coversine 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 = coversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

lib/node_modules/@stdlib/math/base/special/hacoversin/README.md

@@ -60,7 +60,7 @@ var v = hacoversin( 0.0 );
 // returns 0.5
 
 v = hacoversin( 3.141592653589793/2.0 );
-// returns 0.0
+// returns ~0.0
 
 v = hacoversin( -3.141592653589793/6.0 );
 // returns 0.75
@@ -129,7 +129,7 @@ double out = stdlib_base_hacoversin( 0.0 );
 // returns 0.5
 
 out = stdlib_base_hacoversin( 3.141592653589793 / 2.0 );
-// returns 0.0
+// returns ~0.0
 ```
 
 The function accepts the following arguments:

lib/node_modules/@stdlib/math/base/special/hacoversin/docs/types/index.d.ts

@@ -34,7 +34,7 @@
 *
 * @example
 * var v = hacoversin( 3.141592653589793/2.0 );
-* // returns 0.0
+* // returns ~0.0
 *
 * @example
 * var v = hacoversin( -3.141592653589793/6.0 );

lib/node_modules/@stdlib/math/base/special/hacoversin/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 //
@@ -37,7 +38,7 @@ var sin = require( '@stdlib/math/base/special/sin' );
 *
 * @example
 * var v = hacoversin( 3.141592653589793/2.0 );
-* // returns 0.0
+* // returns ~0.0
 *
 * @example
 * var v = hacoversin( -3.141592653589793/6.0 );
@@ -48,7 +49,16 @@ var sin = require( '@stdlib/math/base/special/sin' );
 * // returns NaN
 */
 function hacoversin( 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 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/hacoversin/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"
       ]
     }
   ]