[BUG]: sincosf returning values outside of [-1, 1]

[officiallyanee]

Description

Encountered an error when writing tests for fresnelf which depends on sincosf .

For doube-precision implementation,

In [1]: base.sincos(3e8)
Out[1]: [ -0.43955205118374985, 0.8982171197990818 ]

For single-precision implementation,

In [2]: var sincosf = require( '@stdlib/math/base/special/sincosf')
Out[2]: [Function: sincosf]

In [3]: sincosf(3e8)
Out[3]: [ -24.222797393798828, 37.87095642089844 ]

Note: 3e8 can be represented exactly in single-precision

In [5]: base.float64ToFloat32(3e8)
Out[5]: 300000000

Related PR:
sincosf
rempio2f

• Read and understood the Code of Conduct.
• Searched for existing issues and pull requests.

[stdlib-bot]

👋 Hi there! 👋

And thank you for opening your first issue! We will get back to you shortly. 🏃 💨

If you have any further questions while waiting for a response, please join our Zulip community to chat with project maintainers and other community members.

[officiallyanee]

The issue seems to be from how rempio2f handles float calculations ( sincosf depends on rempio2f )
particularly during rounding:

n = stdlib_base_roundf( x * INVPIO2 );

where x*INVPIO2 is a double since INVPIO2 is a double, rounding to f32 loses important precision and cause problem here

r = x - ( n * PIO2_1 );
w = n * PIO2_1T;
*rem = r - w;
return n;

where |rem| needs <= π/4 for sincosf to use kernel_sincosf here

// Argument reduction...
n = stdlib_base_rempio2f( x, &y );

// Compute the sine and cosine together:
stdlib_base_kernel_sincosf( y, sine, cosine );

FreeBSD's rempio2f the reference implementation uses doubles for all its calculations so it doesn't face this issue.

Also a nit I noticed while looking at FreeBSD's rempio2f:
they use rnint which rounds to nearest even integer ( for .5 cases ) not towards positive infinity.