fix: evaluate large arguments accurately in math/base/special/fresnel, fresnelc and fresnels

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

@@ -36,6 +36,7 @@
 
 var sincos = require( '@stdlib/math/base/special/sincos' ).assign;
 var abs = require( '@stdlib/math/base/special/abs' );
+var fmod = require( '@stdlib/math/base/special/fmod' );
 var HALF_PI = require( '@stdlib/constants/float64/half-pi' );
 var PI = require( '@stdlib/constants/float64/pi' );
 var polyS = require( './rational_psqs.js' );
@@ -113,6 +114,8 @@ var sc = [ 0.0, 0.0 ];
 */
 function fresnel( x, out, stride, offset ) {
 	var x2;
+	var hi;
+	var lo;
 	var xa;
 	var f;
 	var g;
@@ -125,7 +128,8 @@ function fresnel( x, out, stride, offset ) {
 		t = x2 * x2;
 		out[ offset ] = xa * x2 * polyS( t );
 		out[ offset+stride ] = xa * polyC( t );
-	} else if ( xa > 36974.0 ) {
+	} else if ( xa > 1.0e16 ) {
+		// Above `1e16`, the oscillating terms, which are at most `1/(πx)`, are below half an ulp of `0.5`:
 		out[ offset+stride ] = 0.5;
 		out[ offset ] = 0.5;
 	} else {
@@ -136,8 +140,13 @@ function fresnel( x, out, stride, offset ) {
 		t = 1.0 / t;
 		f = 1.0 - ( u * polyF( u ) );
 		g = t * polyG( u );
-		t = HALF_PI * x2;
-		sincos( 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 `134217729 = 2^27+1`) so that `hi*hi`, `2*hi*lo` and `lo*lo` are exact, and reduce each term:
+		t = 134217729.0 * xa;
+		hi = t - ( t - xa );
+		lo = xa - hi;
+		t = fmod( hi*hi, 4.0 ) + fmod( 2.0*hi*lo, 4.0 ) + ( lo*lo );
+		sincos( HALF_PI * t, sc, 1, 0 );
 		t = PI * xa;
 		out[ offset+stride ] = 0.5 + ( ( (f*sc[0]) - (g*sc[1]) ) / t );
 		out[ offset ] = 0.5 - ( ( (f*sc[1]) + (g*sc[0]) ) / t );

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

@@ -41,6 +41,7 @@
         "@stdlib/napi/argv-float64array",
         "@stdlib/napi/export",
         "@stdlib/math/base/special/abs",
+        "@stdlib/math/base/special/fmod",
         "@stdlib/math/base/special/sincos",
         "@stdlib/constants/float64/pi",
         "@stdlib/constants/float64/half-pi"
@@ -58,6 +59,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/special/abs",
+        "@stdlib/math/base/special/fmod",
         "@stdlib/math/base/special/sincos",
         "@stdlib/constants/float64/pi",
         "@stdlib/constants/float64/half-pi"
@@ -75,6 +77,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/special/abs",
+        "@stdlib/math/base/special/fmod",
         "@stdlib/math/base/special/sincos",
         "@stdlib/constants/float64/pi",
         "@stdlib/constants/float64/half-pi"

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

@@ -33,6 +33,7 @@
 #include "stdlib/math/base/special/fresnel.h"
 #include "stdlib/math/base/special/sincos.h"
 #include "stdlib/math/base/special/abs.h"
+#include "stdlib/math/base/special/fmod.h"
 #include "stdlib/constants/float64/half_pi.h"
 #include "stdlib/constants/float64/pi.h"
 
@@ -249,6 +250,8 @@ static double rational_pgqg( const double x ) {
 */
 void stdlib_base_fresnel( const double x, double* S, double* C ) {
 	double x2;
+	double hi;
+	double lo;
 	double xa;
 	double f;
 	double g;
@@ -263,7 +266,8 @@ void stdlib_base_fresnel( const double x, double* S, double* C ) {
 		t = x2 * x2;
 		*S = xa * x2 * rational_psqs( t );
 		*C = xa * rational_pcqc( t );
-	} else if ( xa > 36974.0 ) {
+	} else if ( xa > 1.0e16 ) {
+		// Above `1e16`, the oscillating terms, which are at most `1/(πx)`, are below half an ulp of `0.5`:
 		*C = 0.5;
 		*S = 0.5;
 	} else {
@@ -274,8 +278,13 @@ void stdlib_base_fresnel( const double x, double* S, double* C ) {
 		t = 1.0 / t;
 		f = 1.0 - ( u * rational_pfqf( u ) );
 		g = t * rational_pgqg( u );
-		t = STDLIB_CONSTANT_FLOAT64_HALF_PI * x2;
-		stdlib_base_sincos( 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 `134217729 = 2^27+1`) so that `hi*hi`, `2*hi*lo` and `lo*lo` are exact, and reduce each term:
+		t = 134217729.0 * xa;
+		hi = t - ( t - xa );
+		lo = xa - hi;
+		t = stdlib_base_fmod( hi * hi, 4.0 ) + stdlib_base_fmod( 2.0 * hi * lo, 4.0 ) + ( lo * lo );
+		stdlib_base_sincos( STDLIB_CONSTANT_FLOAT64_HALF_PI * t, &s, &c );
 		t = STDLIB_CONSTANT_FLOAT64_PI * xa;
 		*C = 0.5 + ( ( ( f * s ) - ( g * c ) ) / t );
 		*S = 0.5 - ( ( ( f * c ) + ( g * s ) ) / t );