fix: evaluate large arguments accurately in math/base/special/fresnelf

lib/node_modules/@stdlib/math/base/special/fresnelf/lib/assign.js

@@ -38,6 +38,7 @@
 var PI = require( '@stdlib/constants/float32/pi' );
 var isnanf = require( '@stdlib/math/base/assert/is-nanf' );
 var absf = require( '@stdlib/math/base/special/absf' );
+var fmodf = require( '@stdlib/math/base/special/fmodf' );
 var sincosf = require( '@stdlib/math/base/special/sincosf' ).assign;
 var f32 = require( '@stdlib/number/float64/base/to-float32' );
 var polyC = require( './polyval_cn.js' );
@@ -114,6 +115,8 @@
 */
 function fresnelf( x, out, stride, offset ) {
 	var x2;
+	var hi;
+	var lo;
 	var xa;
 	var f;
 	var g;
@@ -134,7 +137,8 @@
 		t = f32( x2 * x2 );
 		out[ offset ] = f32( f32( xa * x2 ) * polyS( t ) );
 		out[ offset + stride ] = f32( xa * polyC( t ) );
-	} else if ( xa > 36974.0 ) {
+	} else if ( xa > 16777216.0 ) {
+		// Above `2^24`, the oscillating terms, which are at most `1/(πx)`, are below half an ulp of `0.5`:
 		out[ offset ] = 0.5;
 		out[ offset + stride ] = 0.5;
 	} else {
@@ -144,8 +148,13 @@
 		t = f32( 1.0 / t );
 		f = f32( 1.0 - f32( u * polyF( u ) ) );
 		g = f32( t * polyG( u ) );
-		t = f32( HALF_PI * x2 );
-		sincosf( t, sc, 1, 0 );
+
+		// `sin( π x^2 / 2 )` and `cos( π x^2 / 2 )` have period 4 in `x^2`, and forming `x*x` rounds away bits that the large argument turns into whole radians. Reduce `x^2` modulo 4 exactly instead: split `x = hi + lo` (Veltkamp, with `4097 = 2^12+1`) so that `hi*hi`, `2*hi*lo` and `lo*lo` are exact, and reduce each term:
+		t = f32( 4097.0 * xa );
+		hi = f32( t - f32( t - xa ) );
+		lo = f32( xa - hi );
+		t = f32( f32( fmodf( f32( hi*hi ), 4.0 ) + fmodf( f32( 2.0*f32( hi*lo ) ), 4.0 ) ) + f32( lo*lo ) );
+		sincosf( f32( HALF_PI * t ), sc, 1, 0 );
 		t = f32( PI * xa );
 		out[ offset + stride ] = f32( 0.5 + f32( f32( f32( f * sc[0] ) - f32( g * sc[1] ) ) / t ) ); // eslint-disable-line max-len
 		out[ offset ] = f32( 0.5 - f32( f32( f32( f * sc[1] ) + f32( g * sc[0] ) ) / t ) ); // eslint-disable-line max-len

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

@@ -41,6 +41,7 @@
         "@stdlib/napi/argv-float32array",
         "@stdlib/napi/export",
         "@stdlib/math/base/special/absf",
+        "@stdlib/math/base/special/fmodf",
         "@stdlib/math/base/special/sincosf",
         "@stdlib/math/base/assert/is-nanf",
         "@stdlib/constants/float32/pi",
@@ -59,6 +60,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/special/absf",
+        "@stdlib/math/base/special/fmodf",
         "@stdlib/math/base/special/sincosf",
         "@stdlib/math/base/assert/is-nanf",
         "@stdlib/constants/float32/pi",
@@ -77,6 +79,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/special/absf",
+        "@stdlib/math/base/special/fmodf",
         "@stdlib/math/base/special/sincosf",
         "@stdlib/math/base/assert/is-nanf",
         "@stdlib/constants/float32/pi",

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

@@ -35,6 +35,7 @@
 #include "stdlib/constants/float32/pi.h"
 #include "stdlib/math/base/assert/is_nanf.h"
 #include "stdlib/math/base/special/absf.h"
+#include "stdlib/math/base/special/fmodf.h"
 #include "stdlib/math/base/special/sincosf.h"
 
 /* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */
@@ -167,6 +168,8 @@ static float polyval_gn( const float x ) {
 void stdlib_base_fresnelf( const float x, float *S, float *C ) {
 	float xa;
 	float x2;
+	float hi;
+	float lo;
 	float f;
 	float g;
 	float t;
@@ -185,7 +188,8 @@ void stdlib_base_fresnelf( const float x, float *S, float *C ) {
 		t = x2 * x2;
 		*S = xa * x2 * polyval_sn( t );
 		*C = xa * polyval_cn( t );
-	} else if ( xa > 36974.0f ) {
+	} else if ( xa > 16777216.0f ) {
+		// Above `2^24`, the oscillating terms, which are at most `1/(πx)`, are below half an ulp of `0.5`:
 		*S = 0.5f;
 		*C = 0.5f;
 	} else {
@@ -195,8 +199,13 @@ void stdlib_base_fresnelf( const float x, float *S, float *C ) {
 		t = 1.0f / t;
 		f = 1.0f - ( u * polyval_fn( u ) );
 		g = t * polyval_gn( u );
-		t = STDLIB_CONSTANT_FLOAT32_HALF_PI * x2;
-		stdlib_base_sincosf( t, &s, &c );
+
+		// `sin( π x^2 / 2 )` and `cos( π x^2 / 2 )` have period 4 in `x^2`, and forming `x*x` rounds away bits that the large argument turns into whole radians. Reduce `x^2` modulo 4 exactly instead: split `x = hi + lo` (Veltkamp, with `4097 = 2^12+1`) so that `hi*hi`, `2*hi*lo` and `lo*lo` are exact, and reduce each term:
+		t = 4097.0f * xa;
+		hi = t - ( t - xa );
+		lo = xa - hi;
+		t = ( stdlib_base_fmodf( hi * hi, 4.0f ) + stdlib_base_fmodf( 2.0f * hi * lo, 4.0f ) ) + ( lo * lo );
+		stdlib_base_sincosf( STDLIB_CONSTANT_FLOAT32_HALF_PI * t, &s, &c );
 		t = STDLIB_CONSTANT_FLOAT32_PI * xa;
 		*C = 0.5f + ( ( ( f * s ) - ( g * c ) ) / t );
 		*S = 0.5f - ( ( ( f * c ) + ( g * s ) ) / t );