fix: reorder Wald mode evaluation to avoid catastrophic cancellation

Author: Planeshifter

lib/node_modules/@stdlib/stats/base/dists/wald/mode/lib/main.js

@@ -54,7 +54,7 @@ var sqrt = require( '@stdlib/math/base/special/sqrt' );
 * // returns NaN
 */
 function mode( mu, lambda ) {
-	var numerator;
+	var r;
 	if (
 		isnan( mu ) ||
 		isnan( lambda ) ||
@@ -63,8 +63,8 @@ function mode( mu, lambda ) {
 	) {
 		return NaN;
 	}
-	numerator = sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) - 3.0*mu;
-	return ( mu*numerator ) / ( 2.0*lambda );
+	r = mu / lambda;
+	return mu * ( sqrt( 1.0 + ( ( 3.0*r/2.0 ) * ( 3.0*r/2.0 ) ) ) - ( 3.0*r/2.0 ) );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/wald/mode/src/main.c

@@ -32,7 +32,7 @@
 * // returns ~0.325
 */
 double stdlib_base_dists_wald_mode( const double mu, const double lambda ) {
-	double numerator;
+	double r;
 	if (
 		stdlib_base_is_nan( mu ) ||
 		stdlib_base_is_nan( lambda ) ||
@@ -41,6 +41,6 @@ double stdlib_base_dists_wald_mode( const double mu, const double lambda ) {
 	) {
 		return 0.0/0.0; // NaN
 	}
-	numerator = stdlib_base_sqrt( 4.0*lambda*lambda + 9.0*mu*mu ) - 3.0*mu;
-	return mu * numerator / ( 2.0*lambda );
+	r = mu / lambda;
+	return mu * ( stdlib_base_sqrt( 1.0 + ( ( 3.0*r/2.0 ) * ( 3.0*r/2.0 ) ) ) - ( 3.0*r/2.0 ) );
 }

lib/node_modules/@stdlib/stats/base/dists/wald/mode/test/test.js

@@ -117,7 +117,7 @@ tape( 'the function returns the mode of a Wald distribution', function test( t )
 		if ( y === expected[i] ) {
 			t.strictEqual( y, expected[i], 'mu:'+mu[i]+', lambda: '+lambda[i]+', y: '+y+', expected: '+expected[i] );
 		} else {
-			t.ok( isAlmostSameValue( y, expected[i], 1000 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
+			t.ok( isAlmostSameValue( y, expected[i], 2 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
 		}
 	}
 	t.end();

lib/node_modules/@stdlib/stats/base/dists/wald/mode/test/test.native.js

@@ -126,7 +126,7 @@ tape( 'the function returns the mode of a Wald distribution', opts, function tes
 		if ( y === expected[i] ) {
 			t.strictEqual( y, expected[i], 'mu:'+mu[i]+', lambda: '+lambda[i]+', y: '+y+', expected: '+expected[i] );
 		} else {
-			t.ok( isAlmostSameValue( y, expected[i], 1000 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
+			t.ok( isAlmostSameValue( y, expected[i], 2 ), 'within tolerance. mu: '+mu[i]+'. lambda: '+lambda[i]+'. y: '+y+'. E: '+expected[ i ]+'.' );
 		}
 	}
 	t.end();