fix: avoid amplifying the rounding of x^2 for negative x in math/base/special/erfcx

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

@@ -49,6 +49,7 @@
 
 var exp = require( '@stdlib/math/base/special/exp' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
+var floor = require( '@stdlib/math/base/special/floor' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var erfcxY100 = require( './erfcx_y100.js' );
 
@@ -100,7 +101,9 @@ var INV_SQRT_PI = 0.56418958354775628694807945156; // 1 / sqrt(pi)
 * // returns NaN
 */
 function erfcx( x ) {
+	var ex2;
 	var x2;
+	var xh;
 
 	if ( isnan( x ) ) {
 		return x;
@@ -120,11 +123,13 @@ function erfcx( x ) {
 	if ( x < -26.7 ) {
 		return PINF;
 	}
-	x2 = x * x;
+	// Forming `x*x` rounds away bits that `exp` then amplifies by a factor of `x^2` (several hundred ULPs near `x = -26`). Split `x = xh + xl`, with `xh` holding the leading bits of `x` so that `xh*xh` is exact, and use `exp(x^2) = exp(xh^2) * exp(xl*(x+xh))`, whose second factor is close to one:
+	xh = floor( x * 65536.0 ) / 65536.0;
+	ex2 = exp( xh * xh ) * exp( ( x - xh ) * ( x + xh ) );
 	if ( x < -6.1 ) {
-		return 2.0 * exp( x2 );
+		return 2.0 * ex2;
 	}
-	return ( 2.0*exp( x2 ) ) - erfcxY100( 400.0/(4.0-x) );
+	return ( 2.0*ex2 ) - erfcxY100( 400.0/(4.0-x) );
 }
 
 

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

@@ -40,7 +40,8 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/floor",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/math/base/special/exp"
+        "@stdlib/math/base/special/exp",
+        "@stdlib/math/base/special/floor"
       ]
     },
     {
@@ -57,7 +58,8 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/floor",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/math/base/special/exp"
+        "@stdlib/math/base/special/exp",
+        "@stdlib/math/base/special/floor"
       ]
     },
     {
@@ -74,7 +76,8 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/floor",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/math/base/special/exp"
+        "@stdlib/math/base/special/exp",
+        "@stdlib/math/base/special/floor"
       ]
     }
   ]

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

@@ -45,6 +45,7 @@
 
 #include "stdlib/math/base/special/erfcx.h"
 #include "stdlib/math/base/special/exp.h"
+#include "stdlib/math/base/special/floor.h"
 #include "stdlib/constants/float64/pinf.h"
 #include "stdlib/math/base/special/floor.h"
 #include "stdlib/math/base/assert/is_nan.h"
@@ -1133,7 +1134,9 @@ static double erfcxY100( const double y100 ) {
 * // returns +Infinity
 */
 double stdlib_base_erfcx( const double x ) {
+	double ex2;
 	double x2;
+	double xh;
 
 	if( stdlib_base_is_nan( x ) ) {
 		return x;
@@ -1153,9 +1156,11 @@ double stdlib_base_erfcx( const double x ) {
 	if ( x < -26.7 ) {
 		return STDLIB_CONSTANT_FLOAT64_PINF;
 	}
-	x2 = x * x;
+	// Forming `x*x` rounds away bits that `exp` then amplifies by a factor of `x^2` (several hundred ULPs near `x = -26`). Split `x = xh + xl`, with `xh` holding the leading bits of `x` so that `xh*xh` is exact, and use `exp(x^2) = exp(xh^2) * exp(xl*(x+xh))`, whose second factor is close to one:
+	xh = stdlib_base_floor( x * 65536.0 ) / 65536.0;
+	ex2 = stdlib_base_exp( xh * xh ) * stdlib_base_exp( ( x - xh ) * ( x + xh ) );
 	if ( x < -6.1 ) {
-		return 2.0 * stdlib_base_exp( x2 );
+		return 2.0 * ex2;
 	}
-	return ( 2.0 * stdlib_base_exp( x2 ) ) - erfcxY100( 400.0 / ( 4.0 - x ) );
+	return ( 2.0 * ex2 ) - erfcxY100( 400.0 / ( 4.0 - x ) );
 }

lib/node_modules/@stdlib/math/base/special/erfcx/test/fixtures/julia/REQUIRE

@@ -0,0 +1,3 @@
+julia 1.5
+JSON 0.21
+SpecialFunctions 0.10

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

@@ -0,0 +1,75 @@
+#!/usr/bin/env julia
+#
+# @license Apache-2.0
+#
+# Copyright (c) 2026 The Stdlib Authors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import JSON
+import SpecialFunctions: erfcx
+
+"""
+	gen( domain, name )
+
+Generate fixture data and write to file.
+
+The scaled complementary error function is evaluated in extended precision, as `Float64` implementations form `x^2` before exponentiating it.
+
+# Arguments
+
+* `domain`: domain
+* `name::AbstractString`: output filename
+
+# Examples
+
+``` julia
+julia> x = range( -26.0, stop = -1.0, length = 2001 );
+julia> gen( x, "data.json" );
+```
+"""
+function gen( domain, name )
+	x = collect( domain );
+	y = setprecision( BigFloat, 256 ) do
+		Float64.( erfcx.( BigFloat.( x ) ) );
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("x", x),
+		("expected", y)
+	]);
+
+	# Based on the script directory, create an output filepath:
+	filepath = joinpath( dir, name );
+
+	# Write the data to the output filepath as JSON:
+	outfile = open( filepath, "w" );
+	write( outfile, JSON.json(data) );
+	write( outfile, "\n" );
+	close( outfile );
+end
+
+# Get the filename:
+file = @__FILE__;
+
+# Extract the directory in which this file resides:
+dir = dirname( file );
+
+# Medium negative values (previously generated with SciPy, whose `erfcx` forms `x^2` in double precision as well):
+x = range( -10.0, stop = -1.0, length = 1707 );
+gen( x, "medium_negative.json" );
+
+# Negative values, down to where the function overflows:
+x = range( -26.6, stop = -1.0, length = 2001 );
+gen( x, "negative.json" );

lib/node_modules/@stdlib/math/base/special/erfcx/test/fixtures/runner.py

@@ -72,10 +72,6 @@ def main():
     x = np.linspace(2.0**-54, 1.0, 1653)
     generate(x, "small_positive.json")
 
-    # Medium negative values:
-    x = np.linspace(-10.0, -1.0, 1707)
-    generate(x, "medium_negative.json")
-
     # Small negative values:
     x = np.linspace(-1.0, -2.0**-54, 1653)
     generate(x, "small_negative.json")

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

@@ -37,9 +37,10 @@ var INV_SQRT_PI = 0.56418958354775628694807945156; // 1 / sqrt(pi)
 var mediumPositive = require( './fixtures/medium_positive.json' );
 var smallPositive = require( './fixtures/small_positive.json' );
 var smallNegative = require( './fixtures/small_negative.json' );
-var mediumNegative = require( './fixtures/medium_negative.json' );
+var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var largePositive = require( './fixtures/large_positive.json' );
 var tiny = require( './fixtures/tiny.json' );
+var negative = require( './fixtures/julia/negative.json' );
 
 
 // TESTS //
@@ -148,7 +149,22 @@ tape( 'the function computes the scaled complementary error function for negativ
 	x = mediumNegative.x;
 	for ( i = 0; i < expected.length; i++ ) {
 		y = erfcx( x[ i ] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function computes the scaled complementary error function for negative numbers (extended-precision reference)', function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	expected = negative.expected;
+	x = negative.x;
+	for ( i = 0; i < expected.length; i++ ) {
+		y = erfcx( x[ i ] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -33,9 +33,10 @@ var tryRequire = require( '@stdlib/utils/try-require' );
 var mediumPositive = require( './fixtures/medium_positive.json' );
 var smallPositive = require( './fixtures/small_positive.json' );
 var smallNegative = require( './fixtures/small_negative.json' );
-var mediumNegative = require( './fixtures/medium_negative.json' );
+var mediumNegative = require( './fixtures/julia/medium_negative.json' );
 var largePositive = require( './fixtures/large_positive.json' );
 var tiny = require( './fixtures/tiny.json' );
+var negative = require( './fixtures/julia/negative.json' );
 
 
 // VARIABLES //
@@ -157,7 +158,22 @@ tape( 'the function computes the scaled complementary error function for negativ
 	x = mediumNegative.x;
 	for ( i = 0; i < expected.length; i++ ) {
 		y = erfcx( x[ i ] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function computes the scaled complementary error function for negative numbers (extended-precision reference)', opts, function test( t ) {
+	var expected;
+	var x;
+	var y;
+	var i;
+
+	expected = negative.expected;
+	x = negative.x;
+	for ( i = 0; i < expected.length; i++ ) {
+		y = erfcx( x[ i ] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });