fix: avoid cancellation in math/base/special/versin and haversin

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

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

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

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

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

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

lib/node_modules/@stdlib/math/base/special/haversin/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 the smallest `x`), 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,16 @@ file = @__FILE__;
 # Extract the directory in which this file resides:
 dir = dirname( file );
 
+# Negative small values:
+x = -exp10.( range( -300.0, stop = 0.0, length = 4000 ) )
+out = joinpath( dir, "small_negative.json" );
+gen( x, out );
+
+# Positive small values:
+x = exp10.( range( -300.0, stop = 0.0, length = 4000 ) )
+out = joinpath( dir, "small_positive.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/haversin/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,8 @@ var haversin = require( './../lib' );
 
 // FIXTURES //
 
+var smallNegative = require( './fixtures/julia/small_negative.json' );
+var smallPositive = require( './fixtures/julia/small_positive.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 +48,40 @@ tape( 'main export is a function', function test( t ) {
 	t.end();
 });
 
+tape( 'the function computes the half-value versed sine (for -1 <= x <= -1e-300)', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = smallNegative.x;
+	expected = smallNegative.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = haversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
+tape( 'the function computes the half-value versed sine (for 1e-300 <= x <= 1)', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = smallPositive.x;
+	expected = smallPositive.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = haversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
 tape( 'the function computes the half-value versed sine (for -256*pi < x < 0 )', function test( t ) {
 	var expected;
 	var x;
@@ -56,7 +93,7 @@ tape( 'the function computes the half-value versed sine (for -256*pi < x < 0 )',
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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 +110,7 @@ tape( 'the function computes the half-value versed sine (for 0 < x < 256*pi )',
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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 +127,7 @@ tape( 'the function computes the half-value versed sine (for -2**60 (pi/2) < x <
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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 +144,7 @@ tape( 'the function computes the half-value versed sine (for 2**20 (pi/2) < x <
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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 +161,7 @@ tape( 'the function computes the half-value versed sine (for x <= -2**60 (PI/2)
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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 +178,7 @@ tape( 'the function computes the half-value versed sine (for x >= 2**60 (PI/2) )
 
 	for ( i = 0; i < x.length; i++ ) {
 		y = haversin( 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/haversin/test/test.native.js

@@ -39,6 +39,8 @@ var opts = {
 
 // FIXTURES //
 
+var smallNegative = require( './fixtures/julia/small_negative.json' );
+var smallPositive = require( './fixtures/julia/small_positive.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 +57,40 @@ tape( 'main export is a function', opts, function test( t ) {
 	t.end();
 });
 
+tape( 'the function computes the half-value versed sine (for -1 <= x <= -1e-300)', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = smallNegative.x;
+	expected = smallNegative.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = haversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
+tape( 'the function computes the half-value versed sine (for 1e-300 <= x <= 1)', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = smallPositive.x;
+	expected = smallPositive.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = haversin( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+
+	t.end();
+});
+
 tape( 'the function computes the half-value versed sine (for -256*pi < x < 0 )', opts, function test( t ) {
 	var expected;
 	var x;

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

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

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

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