std.math.complex fixes

[tiehuis]

Fixes a few translation errors in the initial port and adds some fixes since changed upstream.

See #19476.

---

I have tested this predominantly using a simple fuzzing setup, comparing the results of zig vs musl (which is the primary source). Code for this can be seen at the bottom of this PR.

zig build-exe main.zig complex.c -target x86_64-linux-musl -lc --zig-lib-dir ../lib -O ReleaseSafe

For values in the range [0, 1) I'e fuzzed ~1 billion results with all matching to a small epsilon.

For larger values there are a few differences noted below, however this encompasses only ~5% of values max and generally the difference is small.

Observations

There are some known failures still present. These are not introduced by this PR and likely have been present for quite a while.

The primary reason from what I can see has to do with some minor precision differences which balloon under certain circumstances. Specifically I would investigate the difference in behaviour between where we emit builtins @log, @sqrt vs. musl which use a software implementation.

e.g. a failing test case

!input     = (8.063326472751923e18, 5.343634481833767e18)
!iteration = 0
!seed      = 11911817038683360459
!
!c         = (-1.5707963267948966e0, -6.931471805599453e0)
!z         = (0e0, inf)

Tracing the Zig and C asin implementation we deviate marginally at a certain point which causes wildly different results:

C:asin:a -5343634481833766912.000000 8063326472751923200.000000
c:sqrt:r1 -36462804330739140062412635512147804160.000000 -86174938756160439949352376809912532992.000000
c:sqrt:r2 5343634481833766912.000000 -8063326472751924224.000000
C:asin:b 5343634481833766912.000000 -8063326472751924224.000000
C:asin:c 0.000000 -1024.000000
C: log = 0.000000 -1024.000000
C:asin:r 6.931472 -1.570796
Z:asin:a -5.343634481833767e18 8.063326472751923e18
z:sqrt:r1 -3.646280433073914e37 -8.617493875616044e37
z:sqrt:r2 5.343634481833767e18 -8.063326472751923e18
Z:asin:b 5.343634481833767e18 -8.063326472751923e18
Z:asin:c 0e0 0e0
Z: log = 0e0 0e0
Z:asin:r -inf 0e0

Note however that both of the above results are incorrect. The correct value is instead 0.98553904473616791... + 44.409042651917785... i.

This looks to occur for values ~1e7 or greater for input to asin. asin is called by a few other hyperbolic functions so this can be observed in a few other places too.

---

Aside from this there are some other minor accuracy differences. Likely due to the aforementioned builtins but for which the errors do not grow exponentially.

clogf

|10000000
!input     = (2.9416764e-5, 3.3091966e-4)
!iteration = 19915237
!seed      = 857500210892449895
!
!c         = (-8.009699e0, 1.4821354e0)
!z         = (-8.0097e0, 1.4821354e0)

---

Fuzzing Code

const std = @import("std");
const Z = std.math.complex;

pub fn main() !void {
    const seed = std.crypto.random.int(u64);
    var prng = std.Random.DefaultPrng.init(seed);
    var random = prng.random();

    const max_iterations = 100_000_000;

    inline for (.{
        .{ "cacosf", f32, zig_cacosf, Z.acos },
        .{ "cacoshf", f32, zig_cacoshf, Z.acosh },
        .{ "casinf", f32, zig_casinf, Z.asin },
        .{ "casinhf", f32, zig_casinhf, Z.asinh },
        .{ "catanf", f32, zig_catanf, Z.atan },
        .{ "ccosf", f32, zig_ccosf, Z.cos },
        .{ "ccoshf", f32, zig_ccoshf, Z.cosh },
        .{ "cexpf", f32, zig_cexpf, Z.exp },
        .{ "clogf", f32, zig_clogf, Z.log },
        .{ "csinf", f32, zig_csinf, Z.sin },
        .{ "csinhf", f32, zig_csinhf, Z.sinh },
        .{ "csqrtf", f32, zig_csqrtf, Z.sqrt },
        .{ "ctanf", f32, zig_ctanf, Z.tan },
        .{ "ctanhf", f32, zig_ctanhf, Z.tanh },

        .{ "cacos", f64, zig_cacos, Z.acos },
        .{ "ccosh", f64, zig_ccosh, Z.cosh },
        .{ "cacosh", f64, zig_cacosh, Z.acosh },
        .{ "casin", f64, zig_casin, Z.asin },
        .{ "casinh", f64, zig_casinh, Z.asinh },
        .{ "catan", f64, zig_catan, Z.atan },
        .{ "ccos", f64, zig_ccos, Z.cos },
        .{ "ccosh", f64, zig_ccosh, Z.cosh },
        .{ "cexp", f64, zig_cexp, Z.exp },
        .{ "clog", f64, zig_clog, Z.log },
        .{ "csin", f64, zig_csin, Z.sin },
        .{ "csinh", f64, zig_csinh, Z.sinh },
        .{ "csqrt", f64, zig_csqrt, Z.sqrt },
        .{ "ctan", f64, zig_ctan, Z.tan },
        .{ "ctanh", f64, zig_ctanh, Z.tanh },
    }) |def| {
        const funcName = def[0];
        const T = def[1];
        const func = def[2];
        const zigFunc = def[3];
        // Our log + sqrt uses builtins which may differ from musl, so we don't expect bit-identical
        // results. This affects asin,acos,asinh,acosh,log, other functions are within one epsilon.
        const eps = 5 * std.math.floatEps(T);
        var generator = RandomZ(T).init(&random);

        std.debug.print("{s}\n\n", .{funcName});
        var iteration: usize = 0;
        while (iteration <= max_iterations) : (iteration += 1) {
            if (iteration != 0 and iteration % 10_000_000 == 0) {
                std.debug.print("|{}\n", .{iteration});
            }

            const zzz = generator.next();
            const r = zzz[0];
            const i = zzz[1];
            const z = std.math.Complex(T).init(r, i);

            const c_val = blk: {
                var rr: T = undefined;
                var ri: T = undefined;
                func(&rr, &ri, r, i);
                break :blk std.math.Complex(T).init(rr, ri);
            };
            const z_val = zigFunc(z);

            // treat nan as equal
            const re_both_nan = std.math.isNan(c_val.re) and std.math.isNan(z_val.re);
            const im_both_nan = std.math.isNan(c_val.im) and std.math.isNan(z_val.im);

            if ((!re_both_nan and !std.math.approxEqAbs(T, c_val.re, z_val.re, eps) or
                (!im_both_nan and !std.math.approxEqAbs(T, c_val.im, z_val.im, eps))))
            {
                std.debug.print(
                    \\!input     = ({}, {})
                    \\!iteration = {}
                    \\!seed      = {}
                    \\!
                    \\!c         = ({}, {})
                    \\!z         = ({}, {})
                    \\
                , .{ z.re, z.im, iteration, seed, c_val.re, c_val.im, z_val.re, z_val.im });
                break;
            }
        }
        std.debug.print("\n{s}\n-----------\n", .{if (iteration > max_iterations) "OK" else "FAIL"});
    }
}

fn RandomZ(comptime T: type) type {
    return struct {
        i: usize,
        r: *std.Random,

        pub fn init(r: *std.Random) @This() {
            return .{ .i = 0, .r = r };
        }

        fn float(z: *@This()) T {
            const U = @Type(.{ .Int = .{ .signedness = .unsigned, .bits = @bitSizeOf(T) } });
            const r: T = @floatFromInt(z.r.int(U));
            return r;
        }

        pub fn next(z: *@This()) [2]T {
            const n = switch (z.i) {
                // Test special cases
                1 => .{ std.math.nan(T), std.math.nan(T) },
                0 => .{ std.math.inf(T), std.math.inf(T) },
                2 => .{ std.math.nan(T), std.math.inf(T) },
                3 => .{ std.math.inf(T), std.math.nan(T) },
                4...1000 => .{ std.math.nan(T), z.float() },
                1001...2000 => .{ std.math.inf(T), z.float() },
                2001...3000 => .{ z.float(), std.math.nan(T) },
                3001...4000 => .{ z.float(), std.math.inf(T) },
                4001...5000 => .{ 0, z.float() },
                5001...6000 => .{ 1, z.float() },
                6001...7000 => .{ -1, z.float() },
                7001...8000 => .{ z.float(), 0 },
                8001...9000 => .{ z.float(), -1 },
                9001...10000 => .{ z.float(), 1 },
                else => .{ z.float(), z.float() },
            };

            z.i += 1;
            return n;
        }
    };
}

extern fn zig_cacosf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_cacoshf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_casinf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_casinhf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_catanf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_catanhf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_ccosf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_ccoshf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_cexpf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_clogf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_csinf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_csinhf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_csqrtf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_ctanf(rr: *f32, ri: *f32, r: f32, i: f32) void;
extern fn zig_ctanhf(rr: *f32, ri: *f32, r: f32, i: f32) void;

extern fn zig_cacos(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_cacosh(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_casin(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_casinh(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_catan(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_catanh(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_ccos(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_ccosh(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_cexp(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_clog(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_csin(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_csinh(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_csqrt(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_ctan(rr: *f64, ri: *f64, r: f64, i: f64) void;
extern fn zig_ctanh(rr: *f64, ri: *f64, r: f64, i: f64) void;

#include <complex.h>

#define def32(cname)                                    \
void zig_##cname(float *rr, float *ri, float r, float i)\
{                                                       \
    float complex z = CMPLXF(r, i);                     \
    float complex rz = cname(z);                        \
    *rr = crealf(rz);                                   \
    *ri = cimagf(rz);                                   \
}

#define def64(cname)                                    \
void zig_##cname(double *rr, double *ri, double r, double i)\
{                                                       \
    double complex z = CMPLX(r, i);                     \
    double complex rz = cname(z);                       \
    *rr = creal(rz);                                    \
    *ri = cimag(rz);                                    \
}

def32(cargf)
def32(cacosf)
def32(cacoshf)
def32(casinf)
def32(casinhf)
def32(catanf)
def32(catanhf)
def32(ccosf)
def32(ccoshf)
def32(cexpf)
def32(clogf)
def32(csinf)
def32(csinhf)
def32(csqrtf)
def32(ctanf)
def32(ctanhf)

def64(carg)
def64(cacos)
def64(cacosh)
def64(casin)
def64(casinh)
def64(catan)
def64(catanh)
def64(ccos)
def64(ccosh)
def64(cexp)
def64(clog)
def64(csin)
def64(csinh)
def64(csqrt)
def64(ctan)
def64(ctanh)

[andrewrk]

Great work!