Skip to content

fix: avoid cancellation near |x| = 1 in math/base/special/acsc #16113

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Abhist17 wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft

fix: avoid cancellation near |x| = 1 in math/base/special/acsc #16113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/acsc/lib/main.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@
// MODULES //

var asin = require( '@stdlib/math/base/special/asin' );
var atan = require( '@stdlib/math/base/special/atan' );
var sqrt = require( '@stdlib/math/base/special/sqrt' );


// MAIN //
Expand Down Expand Up @@ -54,6 +56,13 @@ var asin = require( '@stdlib/math/base/special/asin' );
* // returns NaN
*/
function acsc( x ) {
// Near `|x| = 1`, `asin( 1/x )` amplifies the rounding of `1/x` (the slope of `asin` is unbounded at one). Use `acsc(x) = atan( 1/sqrt(x^2-1) ) = atan( 1/sqrt((x-1)*(x+1)) )` instead, where `x-1` is exact for `1 <= x <= 2`:
if ( x >= 1.0 && x <= 2.0 ) {
return atan( 1.0 / sqrt( ( x - 1.0 ) * ( x + 1.0 ) ) );
}
if ( x <= -1.0 && x >= -2.0 ) {
return -atan( 1.0 / sqrt( ( -x - 1.0 ) * ( 1.0 - x ) ) );
}
return asin( 1.0 / x );
}

Expand Down
12 changes: 9 additions & 3 deletions lib/node_modules/@stdlib/math/base/special/acsc/manifest.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,9 @@
"libpath": [],
"dependencies": [
"@stdlib/math/base/napi/unary",
"@stdlib/math/base/special/asin"
"@stdlib/math/base/special/asin",
"@stdlib/math/base/special/atan",
"@stdlib/math/base/special/sqrt"
]
},
{
Expand All @@ -51,7 +53,9 @@
"libraries": [],
"libpath": [],
"dependencies": [
"@stdlib/math/base/special/asin"
"@stdlib/math/base/special/asin",
"@stdlib/math/base/special/atan",
"@stdlib/math/base/special/sqrt"
]
},
{
Expand All @@ -65,7 +69,9 @@
"libraries": [],
"libpath": [],
"dependencies": [
"@stdlib/math/base/special/asin"
"@stdlib/math/base/special/asin",
"@stdlib/math/base/special/atan",
"@stdlib/math/base/special/sqrt"
]
}
]
Expand Down
9 changes: 9 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/acsc/src/main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@

#include "stdlib/math/base/special/acsc.h"
#include "stdlib/math/base/special/asin.h"
#include "stdlib/math/base/special/atan.h"
#include "stdlib/math/base/special/sqrt.h"

/**
* Computes the arccosecant of a double-precision floating-point number.
Expand All @@ -30,5 +32,12 @@
* // returns ~1.57
*/
double stdlib_base_acsc( const double x ) {
// Near `|x| = 1`, `asin( 1/x )` amplifies the rounding of `1/x` (the slope of `asin` is unbounded at one). Use `acsc(x) = atan( 1/sqrt(x^2-1) ) = atan( 1/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( 1.0 / stdlib_base_sqrt( ( x - 1.0 ) * ( x + 1.0 ) ) );
}
if ( x <= -1.0 && x >= -2.0 ) {
return -stdlib_base_atan( 1.0 / stdlib_base_sqrt( ( -x - 1.0 ) * ( 1.0 - x ) ) );
}
return stdlib_base_asin( 1.0 / x );
}
1 change: 1 addition & 0 deletions lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/near_negative_one.json
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/near_positive_one.json
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/negative.json
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/positive.json
View file
Open in desktop

Large diffs are not rendered by default.

12 changes: 11 additions & 1 deletion lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/runner.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
import JSON

"""
gen( domain, name )

Check failure on line 22 in lib/node_modules/@stdlib/math/base/special/acsc/test/fixtures/julia/runner.jl

View workflow job for this annotation

GitHub Actions / Lint Changed Files 

wrong indentation type (spaces instead of tabs)

Generate fixture data and write to file.

Expand All @@ -37,7 +37,9 @@
"""
function gen( domain, name )
x = collect( domain );
y = acsc.( x );
y = setprecision( BigFloat, 256 ) do
Float64.( acsc.( BigFloat.( x ) ) );
end

# Store data to be written to file as a collection:
data = Dict([
Expand Down Expand Up @@ -68,3 +70,11 @@
# Generate fixture data for positive values:
x = range( 1.0, stop = 1000.0, length = 2003 );
gen( x, "positive.json" );

# Values close to one (in `Float64`, `acsc` is `asin(1/x)`, which amplifies the rounding of `1/x`):
x = 1.0 .+ exp10.( range( -15.0, stop = -1.0, length = 1000 ) );
gen( x, "near_positive_one.json" );

x = -1.0 .- exp10.( range( -15.0, stop = -1.0, length = 1000 ) );
gen( x, "near_negative_one.json" );

39 changes: 37 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/acsc/test/test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@

var tape = require( 'tape' );
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' );
var randu = require( '@stdlib/random/base/randu' );
var acsc = require( './../lib' );

Expand All @@ -30,6 +31,8 @@ var acsc = require( './../lib' );

var negative = require( './fixtures/julia/negative.json' );
var positive = require( './fixtures/julia/positive.json' );
var nearPositiveOne = require( './fixtures/julia/near_positive_one.json' );
var nearNegativeOne = require( './fixtures/julia/near_negative_one.json' );


// TESTS //
Expand All @@ -51,7 +54,7 @@ tape( 'the function computes the arccosecant (negative values)', function test(

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( y, expected[ i ], 'returns expected value' );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});
Expand All @@ -67,7 +70,39 @@ tape( 'the function computes the arccosecant (positive values)', function test(

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( y, expected[ i ], 'returns expected value' );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});

tape( 'the function computes the arccosecant (values close to `1`)', function test( t ) {
var expected;
var x;
var y;
var i;

x = nearPositiveOne.x;
expected = nearPositiveOne.expected;

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});

tape( 'the function computes the arccosecant (values close to `-1`)', function test( t ) {
var expected;
var x;
var y;
var i;

x = nearNegativeOne.x;
expected = nearNegativeOne.expected;

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});
Expand Down
39 changes: 37 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/acsc/test/test.native.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
var resolve = require( 'path' ).resolve;
var tape = require( 'tape' );
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' );
var randu = require( '@stdlib/random/base/randu' );
var tryRequire = require( '@stdlib/utils/try-require' );

Expand All @@ -31,6 +32,8 @@ var tryRequire = require( '@stdlib/utils/try-require' );

var negative = require( './fixtures/julia/negative.json' );
var positive = require( './fixtures/julia/positive.json' );
var nearPositiveOne = require( './fixtures/julia/near_positive_one.json' );
var nearNegativeOne = require( './fixtures/julia/near_negative_one.json' );


// VARIABLES //
Expand Down Expand Up @@ -60,7 +63,7 @@ tape( 'the function computes the arccosecant (negative values)', opts, function

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( y, expected[ i ], 'returns expected value' );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});
Expand All @@ -76,7 +79,39 @@ tape( 'the function computes the arccosecant (positive values)', opts, function

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( y, expected[ i ], 'returns expected value' );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});

tape( 'the function computes the arccosecant (values close to `1`)', opts, function test( t ) {
var expected;
var x;
var y;
var i;

x = nearPositiveOne.x;
expected = nearPositiveOne.expected;

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});

tape( 'the function computes the arccosecant (values close to `-1`)', opts, function test( t ) {
var expected;
var x;
var y;
var i;

x = nearNegativeOne.x;
expected = nearNegativeOne.expected;

for ( i = 0; i < x.length; i++ ) {
y = acsc( x[i] );
t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
}
t.end();
});
Expand Down
Loading