fix: avoid cancellation near one in math/base/special/asec

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

@@ -22,6 +22,8 @@
 // MODULES //
 
 var acos = require( '@stdlib/math/base/special/acos' );
+var atan = require( '@stdlib/math/base/special/atan' );
+var sqrt = require( '@stdlib/math/base/special/sqrt' );
 
 
 // MAIN //
@@ -32,7 +34,7 @@
 * ## Method
 *
 * ```tex
 * \operatorname{asec}(x) = \operatorname{acos}(\frac{1}{x})
 * ```
 *
 * @param {number} x - input value
@@ -51,6 +53,10 @@
 * // returns NaN
 */
 function asec( x ) {
+	// Near `x = 1`, `acos( 1/x )` loses the low bits of `1-1/x` to the rounding of `1/x`, and those bits are most of the (small) result. Use `asec(x) = atan( sqrt(x^2-1) ) = atan( sqrt((x-1)*(x+1)) )` instead, where `x-1` is exact for `1 <= x <= 2`:
+	if ( x >= 1.0 && x <= 2.0 ) {
+		return atan( sqrt( ( x - 1.0 ) * ( x + 1.0 ) ) );
+	}
 	return acos( 1.0 / x );
 }
 

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

@@ -37,7 +37,9 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/napi/unary",
-        "@stdlib/math/base/special/acos"
+        "@stdlib/math/base/special/acos",
+        "@stdlib/math/base/special/atan",
+        "@stdlib/math/base/special/sqrt"
       ]
     },
     {
@@ -51,7 +53,9 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
-        "@stdlib/math/base/special/acos"
+        "@stdlib/math/base/special/acos",
+        "@stdlib/math/base/special/atan",
+        "@stdlib/math/base/special/sqrt"
       ]
     },
     {
@@ -65,7 +69,9 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
-        "@stdlib/math/base/special/acos"
+        "@stdlib/math/base/special/acos",
+        "@stdlib/math/base/special/atan",
+        "@stdlib/math/base/special/sqrt"
       ]
     }
   ]

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

@@ -18,6 +18,8 @@
 
 #include "stdlib/math/base/special/asec.h"
 #include "stdlib/math/base/special/acos.h"
+#include "stdlib/math/base/special/atan.h"
+#include "stdlib/math/base/special/sqrt.h"
 
 /**
 * Computes the inverse (arc) secant of a double-precision floating-point number.
@@ -30,5 +32,9 @@
 * // returns 0.0
 */
 double stdlib_base_asec( const double x ) {
+	// Near `x = 1`, `acos( 1/x )` loses the low bits of `1-1/x` to the rounding of `1/x`, and those bits are most of the (small) result. Use `asec(x) = atan( sqrt(x^2-1) ) = atan( sqrt((x-1)*(x+1)) )` instead, where `x-1` is exact for `1 <= x <= 2`:
+	if ( x >= 1.0 && x <= 2.0 ) {
+		return stdlib_base_atan( stdlib_base_sqrt( ( x - 1.0 ) * ( x + 1.0 ) ) );
+	}
 	return stdlib_base_acos( 1.0 / x );
 }

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

@@ -19,7 +19,7 @@
 import JSON
 
 """
     gen( domain, name )
 
 Generate fixture data and write to file.
 
@@ -37,7 +37,9 @@
 """
 function gen( domain, name )
 	x = collect( domain );
-	y = asec.( x );
+	y = setprecision( BigFloat, 256 ) do
+		Float64.( asec.( BigFloat.( x ) ) );
+	end
 
 	# Store data to be written to file as a collection:
 	data = Dict([
@@ -92,3 +94,8 @@
 # Huge positive values:
 x = range( 1e300, stop = 1e308, length = 503 );
 gen( x, "huge_positive.json" );
+
+# Values close to one (in `Float64`, `asec` is `acos(1/x)` and loses the low bits of `1-1/x`):
+x = 1.0 .+ exp10.( range( -15.0, stop = -1.0, length = 1000 ) );
+gen( x, "near_one.json" );
+

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

@@ -36,6 +36,7 @@ var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var mediumPositive = require( './fixtures/julia/medium_positive.json' );
 var hugeNegative = require( './fixtures/julia/huge_negative.json' );
 var hugePositive = require( './fixtures/julia/huge_positive.json' );
+var nearOne = require( './fixtures/julia/near_one.json' );
 
 
 // TESTS //
@@ -174,6 +175,22 @@ tape( 'the function computes the inverse (arc) secant on the interval `[1e300,1e
 	t.end();
 });
 
+tape( 'the function computes the inverse (arc) secant close to `1`', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearOne.x;
+	expected = nearOne.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = asec( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
 tape( 'the function returns `NaN` if provided `NaN`', function test( t ) {
 	var v = asec( NaN );
 	t.strictEqual( isnan( v ), true, 'returns expected value' );

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

@@ -37,6 +37,7 @@ var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var mediumPositive = require( './fixtures/julia/medium_positive.json' );
 var hugeNegative = require( './fixtures/julia/huge_negative.json' );
 var hugePositive = require( './fixtures/julia/huge_positive.json' );
+var nearOne = require( './fixtures/julia/near_one.json' );
 
 
 // VARIABLES //
@@ -183,6 +184,22 @@ tape( 'the function computes the inverse (arc) secant on the interval `[1e300,1e
 	t.end();
 });
 
+tape( 'the function computes the inverse (arc) secant close to `1`', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	x = nearOne.x;
+	expected = nearOne.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		y = asec( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
 tape( 'the function returns `NaN` if provided `NaN`', opts, function test( t ) {
 	var v = asec( NaN );
 	t.strictEqual( isnan( v ), true, 'returns expected value' );