fix: avoid cancellation in stats/base/dists/kumaraswamy/logcdf

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/lib/factory.js

@@ -23,8 +23,12 @@
 var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var pow = require( '@stdlib/math/base/special/pow' );
+var log1mexp = require( '@stdlib/math/base/special/log1mexp' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
+var expm1 = require( '@stdlib/math/base/special/expm1' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -68,6 +72,9 @@ function factory( a, b ) {
 	* // returns <number>
 	*/
 	function logcdf( x ) {
+		var u;
+		var v;
+		var t;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
@@ -77,7 +84,22 @@ function factory( a, b ) {
 		if ( x >= 1.0 ) {
 			return 0.0;
 		}
-		return ln( 1.0 - pow( 1.0 - pow( x, a ), b ) );
+		// Evaluate `v = ln(1-x^a)` without forming `1-x^a`, which rounds away the low-order bits of a small `x^a` and, for `x^a` close to `1`, of `1-x^a` itself:
+		u = pow( x, a );
+		if ( u < FLOAT64_SMALLEST_NORMAL ) {
+			// `x^a` is subnormal or underflows, so `ln(1-(1-x^a)^b) = ln(b) + a*ln(x)`:
+			return ln( b ) + ( a * ln( x ) );
+		}
+		if ( u < 0.5 ) {
+			v = log1p( -u );
+		} else {
+			v = ln( -expm1( a * ln( x ) ) );
+		}
+		t = b * v;
+		if ( t > -FLOAT64_SMALLEST_NORMAL ) {
+			return ln( b ) + ln( -v );
+		}
+		return log1mexp( t );
 	}
 }
 

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/lib/main.js

@@ -22,8 +22,12 @@
 
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var pow = require( '@stdlib/math/base/special/pow' );
+var log1mexp = require( '@stdlib/math/base/special/log1mexp' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
+var expm1 = require( '@stdlib/math/base/special/expm1' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -81,6 +85,9 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 * // returns NaN
 */
 function logcdf( x, a, b ) {
+	var u;
+	var v;
+	var t;
 	if (
 		isnan( x ) ||
 		isnan( a ) ||
@@ -96,7 +103,22 @@ function logcdf( x, a, b ) {
 	if ( x >= 1.0 ) {
 		return 0.0;
 	}
-	return ln( 1.0 - pow( 1.0 - pow( x, a ), b ) );
+	// Evaluate `v = ln(1-x^a)` without forming `1-x^a`, which rounds away the low-order bits of a small `x^a` and, for `x^a` close to `1`, of `1-x^a` itself:
+	u = pow( x, a );
+	if ( u < FLOAT64_SMALLEST_NORMAL ) {
+		// `x^a` is subnormal or underflows, so `ln(1-(1-x^a)^b) = ln(b) + a*ln(x)`:
+		return ln( b ) + ( a * ln( x ) );
+	}
+	if ( u < 0.5 ) {
+		v = log1p( -u );
+	} else {
+		v = ln( -expm1( a * ln( x ) ) );
+	}
+	t = b * v;
+	if ( t > -FLOAT64_SMALLEST_NORMAL ) {
+		return ln( b ) + ln( -v );
+	}
+	return log1mexp( t );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/manifest.json

@@ -41,8 +41,12 @@
         "@stdlib/math/base/napi/ternary",
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/log1mexp",
+        "@stdlib/math/base/special/log1p",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/ln",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     },
     {
@@ -59,8 +63,12 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/log1mexp",
+        "@stdlib/math/base/special/log1p",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/ln",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     },
     {
@@ -77,8 +85,12 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/log1mexp",
+        "@stdlib/math/base/special/log1p",
+        "@stdlib/math/base/special/expm1",
         "@stdlib/math/base/special/ln",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     }
   ]

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/src/main.c

@@ -19,8 +19,12 @@
 #include "stdlib/stats/base/dists/kumaraswamy/logcdf.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/pow.h"
+#include "stdlib/math/base/special/log1mexp.h"
+#include "stdlib/math/base/special/log1p.h"
+#include "stdlib/math/base/special/expm1.h"
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/constants/float64/ninf.h"
+#include "stdlib/constants/float64/smallest_normal.h"
 
 /**
 * Evaluates the natural logarithm of the cumulative distribution function (CDF) for a Kumaraswamy's double bounded distribution with first shape parameter `a` and second shape parameter `b` at a value `x`.
@@ -35,6 +39,9 @@
 * // returns ~-0.693
 */
 double stdlib_base_dists_kumaraswamy_logcdf( const double x, const double a, const double b ) {
+	double u;
+	double v;
+	double t;
 	if (
 		stdlib_base_is_nan( x ) ||
 		stdlib_base_is_nan( a ) ||
@@ -50,5 +57,20 @@ double stdlib_base_dists_kumaraswamy_logcdf( const double x, const double a, con
 	if ( x >= 1.0 ) {
 		return 0.0;
 	}
-	return stdlib_base_ln( 1.0 - stdlib_base_pow( 1.0 - stdlib_base_pow( x, a ), b ) );
+	// Evaluate `v = ln(1-x^a)` without forming `1-x^a`, which rounds away the low-order bits of a small `x^a` and, for `x^a` close to `1`, of `1-x^a` itself:
+	u = stdlib_base_pow( x, a );
+	if ( u < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		// `x^a` is subnormal or underflows, so `ln(1-(1-x^a)^b) = ln(b) + a*ln(x)`:
+		return stdlib_base_ln( b ) + ( a * stdlib_base_ln( x ) );
+	}
+	if ( u < 0.5 ) {
+		v = stdlib_base_log1p( -u );
+	} else {
+		v = stdlib_base_ln( -stdlib_base_expm1( a * stdlib_base_ln( x ) ) );
+	}
+	t = b * v;
+	if ( t > -STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		return stdlib_base_ln( b ) + stdlib_base_ln( -v );
+	}
+	return stdlib_base_log1mexp( t );
 }

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/test/fixtures/julia/REQUIRE

@@ -0,0 +1,2 @@
+julia 1.5
+JSON 0.21

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/test/fixtures/julia/runner.jl

@@ -0,0 +1,80 @@
+#!/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
+
+"""
+	gen( x, a, b, name )
+
+Generate fixture data and write to file.
+
+The logarithm of the CDF is evaluated in extended precision from the textbook formula `ln(1 - (1 - x^a)^b)`, as forming `1 - x^a` in double precision rounds away the low-order bits of `x^a` (all of them for `x^a < 2^-53`).
+
+# Arguments
+
+* `x`: input value
+* `a`: first shape parameter
+* `b`: second shape parameter
+* `name::AbstractString`: output filename
+"""
+function gen( x, a, b, name )
+	z = setprecision( BigFloat, 8192 ) do
+		Float64.( log.( 1 .- ( 1 .- BigFloat.( x ) .^ BigFloat.( a ) ) .^ BigFloat.( b ) ) )
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("x", x),
+		("a", a),
+		("b", b),
+		("expected", z)
+	]);
+
+	# 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 );
+
+# Random values:
+x = rand( 1000 );
+a = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+b = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+gen( x, a, b, "data.json" );
+
+# Small `x`:
+x = exp10.( range( -300.0, stop = -1.0, length = 1000 ) );
+a = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+b = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+gen( x, a, b, "small_x.json" );
+
+# `x` close to `1`:
+x = 1.0 .- exp10.( range( -15.0, stop = -1.0, length = 1000 ) );
+a = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+b = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+gen( x, a, b, "large_x.json" );

lib/node_modules/@stdlib/stats/base/dists/kumaraswamy/logcdf/test/test.factory.js

@@ -21,15 +21,20 @@
 // MODULES //
 
 var tape = require( 'tape' );
-var isNumber = require( '@stdlib/assert/is-number' ).isPrimitive;
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var randu = require( '@stdlib/random/base/randu' );
+var isAlmostSameValue = require( '@stdlib/assert/is-almost-same-value' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
-var EPS = require( '@stdlib/constants/float64/eps' );
 var factory = require( './../lib/factory.js' );
 
 
+// FIXTURES //
+
+var data = require( './fixtures/julia/data.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
+var largeX = require( './fixtures/julia/large_x.json' );
+
+
 // TESTS //
 
 tape( 'main export is a function', function test( t ) {
@@ -186,22 +191,65 @@ tape( 'if provided a nonpositive `a`, the created function always returns `NaN`'
 	t.end();
 });
 
-tape( 'the returned function evaluates the logcdf of a Kumaraswamy\'s double bounded distribution', function test( t ) {
+tape( 'the created function evaluates the logcdf for random `x`, `a` and `b`', function test( t ) {
+	var expected;
 	var logcdf;
 	var a;
 	var b;
+	var i;
 	var x;
 	var y;
+
+	expected = data.expected;
+	x = data.x;
+	a = data.a;
+	b = data.b;
+	for ( i = 0; i < expected.length; i++ ) {
+		logcdf = factory( a[i], b[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 25 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logcdf for small `x`', function test( t ) {
+	var expected;
+	var logcdf;
+	var a;
+	var b;
 	var i;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	a = smallX.a;
+	b = smallX.b;
+	for ( i = 0; i < expected.length; i++ ) {
+		logcdf = factory( a[i], b[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logcdf for `x` close to `1`', function test( t ) {
+	var expected;
+	var logcdf;
+	var a;
+	var b;
+	var i;
+	var x;
+	var y;
 
-	// TODO: Add fixtures
-	a = ( randu()*5.0 ) + EPS;
-	b = ( randu()*5.0 ) + EPS;
-	logcdf = factory( a, b );
-	for ( i = 0; i < 100; i++ ) {
-		x = randu();
-		y = logcdf( x, a, b );
-		t.strictEqual( isNumber( y ), true, 'returns expected value' );
+	expected = largeX.expected;
+	x = largeX.x;
+	a = largeX.a;
+	b = largeX.b;
+	for ( i = 0; i < expected.length; i++ ) {
+		logcdf = factory( a[i], b[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 148 ), true, 'returns expected value' );
 	}
 	t.end();
 });