std.math.acosh of large negative numbers does not return nan #35315

Zig Version

0.16.0

Steps to Reproduce and Observed Behavior

I stumbled on this by accident, but the title describes the problem quite accurately: acosh should return nan for negative numbers, but returns finite values for some large negative values.

Steps to reproduce

Input:

 const std: type = @import("std");
 
 pub fn main() void {
 const x: f64 = -1.9e4;
 const result: f64 = std.math.acosh(x);
 std.debug.print("x: {e}, acosh(x): {e}, isNan: {}\n", .{ x, result, std.math.isNan(result) });
 }

Output:
x: -1.9e4, acosh(x): -7.288321402140781e0, isNan: false

This seems to start returning finite values for inputs smaller than -1e4. It seems to work correctly for certain values (e.g. 4e-4, 1.5e-4, 2e-4).

Reproducible with the following build modes:

• Debug Mode / Default zig run and zig build run.
• ReleaseFast, ReleaseSafe.

Observed Behaviour

• std.math.acosh(-1.9e4) returns -7.288321402140781e0
• std.math.acosh(-2e4) returns -nan

Environment Information

• OS: NixOS Yarara
• Zig version: 0.16.0
• target: x86_64-linux.6.19.8...6.19.8-gnu.2.42 (Intel Alderlake)

Expected Behavior

Expected Behaviour

• std.math.acosh(-1.9e4) should return nan
• std.math.acosh(-2e4) should return nan

The acosh implementation is a Zig port of the musl code. I downloaded a fresh cross toolchain from musl.cc (Note: it's community supported and not affiliated with the musl project, trust at your own risk).

#include <stdio.h>#include <math.h>
int main(void)
{
 printf("%f\n", acosh(-1.9e4));
 return 0;
}

$ ./x86_64-linux-musl-cross/bin/x86_64-linux-musl-cc -o acosh acosh.c -lm -static
$ ./acosh
-7.288321

So it looks like the problem originated upstream, and should probably be reported there too.

The problematic part of the code seems to be this branch in acosh64:

// |x| < 0x1p26elseif(e<0x3FF+26){return@log(2*x-1/(x+@sqrt(x*x-1)));}

The idea here is that 2x - 1/(x + sqrt(x^2 - 1)) is only ever non-negative if x >= 1, so @log will return NaN for all values that are outside the domain of arccosh. The expression isn't an approximation, and in fact:

@log(2*x-1/(x+@sqrt(x*x-1)))=@log(x+sqrt(x*x-1))=arccosh(x)

So this code isn't wrong per se, when one is working with real numbers. What's happening is that the expression 1 / (x + @sqrt(x * x - 1)) is falling victim to catastrophic cancellation. One of examples given on that page is that of arcsine(ix). Note that arcsine(ix) = i arcsinh(x), which has a fairly similar definition to arccosh(x), so this is almost a textbook example.

Let x: f128 = -1.9e4 for a moment:

1/(x+@sqrt(x*x-1))==-37999.99997368421050809155699567155>-38000

And so inner expression inside @log(...) is negative, and @log returns NaN.

But if we switch back to double precision f64:

1/(x+@sqrt(x*x-1)))==-38000.00068347437<-38000

Which gives an incorrect positive expression inside @log(...).

As a fix, I'd suggest switching to an implementation based on FreeBSD. Something like this:

conststd=@import("std");constmath=std.math;constexpect=std.testing.expect;fnacosh64(x:f64)f64{consti=@as(i64,@bitCast(x));conste:i12=@intCast(i>>52);constm=@as(u64,@bitCast(i))&0x000F_FFFF_FFFF_FFFF;if(e<0x3FF){// x < 1// I'm following the FreeBSD semantics of using a signalling NaN in this case, // but don't know if that is desired behaviour for Zigreturnstd.math.snan(f64);}elseif(e>=0x7FF){// x is inf or NaNreturnx+x;}elseif(e>=0x41B){// x > 2^28return@log(x)+math.ln2;}elseif(e==0x3FFandm==0){// x == 0return0.0;}elseif(e>0x400){// 2^28 > x > 2return@log(2*x-1/(x+@sqrt(x*x-1)));}else{returnmath.log1p(x-1+@sqrt((x-1)*(x-1)+2*(x-1)));}}testacosh64{constepsilon=0.000001;tryexpect(math.approxEqAbs(f64,acosh64(1.5),0.962424,epsilon));tryexpect(math.approxEqAbs(f64,acosh64(37.45),4.315976,epsilon));tryexpect(math.approxEqAbs(f64,acosh64(89.123),5.183133,epsilon));tryexpect(math.approxEqAbs(f64,acosh64(123123.234375),12.414088,epsilon));}test"acosh64.special"{tryexpect(math.isNan(acosh64(math.nan(f64))));tryexpect(math.isNan(acosh64(0.5)));tryexpect(math.isNan(acosh64(-1.9e-4)));}

I'd put a PR together quickly myself, but I only have access to my work machine at the moment. The above may not be exactly correct either, I just threw it together quickly, modifying it slightly to use Zig's arbitrary sized integers and reduce nesting.

I'm not too familiar with the specifics of implementing this as part of the math library, as I started using Zig very recently, but is there a particular reason all the conditions need to be handled in our case? Something like this:

conststd=@import("std");constmath=std.math;fnacosh64(x:f64)f64{consti:i64=@bitCast(x);conste:i64=@bitCast(i>>52);// I couldn't get it to compile successfully with i12if(e<0x3FF){returnmath.snan(f64);}elseif(e>=0x41B){return@log(x)+math.ln2;}else{return@log(x+@sqrt((x*x)-1));}}test"acosh64"{varxo256=std.Random.Xoshiro256.init(0);varrng=xo256.random();consttol:f64=1e-14;consttest_count:usize=1e7;constrange:f64=1e8;for(0..test_count)|_|{constx:f64=range*rng.float(f64)-0.5*range;constcurrent:f64=math.acosh(x);constproposed:f64=acosh64(x);if(x<1.0){trystd.testing.expect(math.isNan(proposed));}elseif(math.isNan(x)){trystd.testing.expect(math.isNan(proposed));}elseif(math.isInf(x)){trystd.testing.expect(math.isInf(proposed));}else{trystd.testing.expect(math.approxEqAbs(f64,current,proposed,tol));}}}

Doesn't have quite as much conditional complexity and seems to handle all cases correctly, at least from what I have been able to test so far, and should propagate infs/nans. It also seems to outperform the current implementation (~2x) based on very rudimentary testing.

Sorry, that should've been an @intCast. I've corrected it.

And I'm not sure which conditions really need to be handled here, the catastrophic cancellation just nerd sniped me early and I figured I'd go with something tried and tested.