Incorrect value for logarithm of small numbers (denormal floats)

[ScoreUnder]

Zig Version

0.14.1

Steps to Reproduce and Observed Behavior

calculate ln(287.26582147382784e-310) via @log or std.math.log:

test {
    // expected value sourced from Wolfram Alpha
    const expected: f64 = -711.8181974647329933991007419903441317012207784976779687026632850985;
    var value: f64 = 7.26582147382784e-310;
    const tolerance: f64 = std.math.sqrt(std.math.floatEps(f64));

    std.mem.doNotOptimizeAway(&value); // ensure the compiler doesn't cleverly replace @log at compile time (it gets it right...?)

    const other_fun_logarithms = struct {
        pub fn ln_by_log2(x: anytype) @TypeOf(x) {
            return @log2(x) * std.math.ln2;
        }
        pub fn ln_by_log10(x: anytype) @TypeOf(x) {
            return @log10(x) * std.math.ln10;
        }
    };

    try std.testing.expectApproxEqRel(expected, other_fun_logarithms.ln_by_log2(value), tolerance);
    try std.testing.expectApproxEqRel(expected, other_fun_logarithms.ln_by_log10(value), tolerance);
    try std.testing.expectApproxEqRel(expected, @log(value), tolerance);
    try std.testing.expectApproxEqRel(expected, std.math.log(f64, std.math.e, value), tolerance);
}

@log and std.math.log give -7.464873809916107e2, which is almost 5% away from the actual answer.
@log2 and @log10 give the correct answer.

Expected Behavior

ln(287.26582147382784e-310) ≈ 711.8

https://www.wolframalpha.com/input?i=ln%287.26582147382784e-310%29

[ScoreUnder]

Being compiled and run on an AMD ThreadRipper 1950x.

Ghidra decomp

anyerror test_0(void)

{
  anyerror aVar1;
  f64 *in_RDI;
  f64 tolerance_00;
  f64 fVar2;
  double actual;
  f64 expected;
  f64 value;
  u64 biased_exponent;
  f64 tolerance;
  
  tolerance_00 = sqrt__anon_39499(2.220446049250313e-16);
  doNotOptimizeAway__anon_39505(in_RDI);
  fVar2 = ln_by_log2__anon_39516(7.26582147382784e-310);
  aVar1 = expectApproxEqRelInner__anon_39523(-711.818197464733,fVar2,tolerance_00);
  if (aVar1 != 0) {
    returnError();
    return aVar1;
  }
  fVar2 = ln_by_log10__anon_39547(7.26582147382784e-310);
  aVar1 = expectApproxEqRelInner__anon_39523(-711.818197464733,fVar2,tolerance_00);
  if (aVar1 != 0) {
    returnError();
    return aVar1;
  }
  actual = log(7.26582147382784e-310);
  aVar1 = expectApproxEqRelInner__anon_39523(-711.818197464733,(f64)actual,tolerance_00);
  if (aVar1 != 0) {
    returnError();
    return aVar1;
  }
  fVar2 = log__anon_35236(2.718281828459045,7.26582147382784e-310);
  aVar1 = expectApproxEqRelInner__anon_39523(-711.818197464733,fVar2,tolerance_00);
  if (aVar1 != 0) {
    returnError();
    return aVar1;
  }
  return 0;
}

f64 log__anon_35236(f64 base,f64 x)

{
  double dVar1;
  f64 fVar2;
  double dVar3;
  f64 base_local;
  f64 x_local;
  f64 float_base;
  
  if (((double)base == 2.0) && (!NAN((double)base))) {
    fVar2 = log2__anon_35381(x);
    return fVar2;
  }
  if (((double)base == 10.0) && (!NAN((double)base))) {
    fVar2 = log10__anon_35385(x);
    return fVar2;
  }
  if (((double)base == 2.718281828459045) && (!NAN((double)base))) {
    dVar3 = log((double)x);
    return (f64)dVar3;
  }
  fVar2 = lossyCast__anon_35390(base);
  dVar3 = log((double)x);
  dVar1 = log((double)fVar2);
  return (f64)(dVar3 / dVar1);
}

double log(double __x)

{
  double dVar1;
  double dVar2;
  uint uVar3;
  int iVar4;
  double dVar5;
  double dVar6;
  undefined1 auVar7 [16];
  undefined1 auVar8 [16];
  double dVar9;
  undefined1 auVar10 [16];
  undefined1 auVar11 [16];
  undefined1 auVar12 [16];
  
  uVar3 = (uint)((ulong)__x >> 0x20);
  if ((long)__x < 0x10000000000000) {
    dVar5 = NAN;
    if ((__x == 0.0) && (!NAN(__x))) {
      dVar5 = -INFINITY;
    }
    if (__x != 0.0 && -1 < (long)__x) {
      iVar4 = -0x36;
LAB_01145af2:
      uVar3 = uVar3 + 0x95f62;
      auVar12._8_8_ = 0x3fc39a09d078c69f;
      auVar12._0_8_ = 0x3fc39a09d078c69f;
      dVar5 = (double)((ulong)((uVar3 & 0xfffff) + 0x3fe6a09e) << 0x20 | (ulong)__x & 0xffffffff) +
              -1.0;
      dVar6 = dVar5 / (dVar5 + 2.0);
      dVar1 = dVar5 * dVar5 * 0.5;
      dVar9 = dVar6 * dVar6;
      dVar2 = dVar9 * dVar9;
      auVar10._8_8_ = dVar2;
      auVar10._0_8_ = dVar2;
      auVar7._8_8_ = 0;
      auVar7._0_8_ = dVar9;
      auVar7 = vunpcklpd_avx(auVar10,auVar7);
      auVar11._8_8_ = 0;
      auVar11._0_8_ = dVar2 * 0.14798198605116586 + 0.1818357216161805;
      auVar10 = vunpcklpd_avx(auVar12,auVar11);
      auVar8._0_8_ = auVar7._0_8_ *
                     (dVar2 * (dVar2 * auVar10._0_8_ + 0.22222198432149784) + 0.3999999999940942);
      auVar8._8_8_ = auVar7._8_8_ *
                     (dVar2 * (dVar2 * auVar10._8_8_ + 0.2857142874366239) + 0.6666666666666735);
      auVar7 = vshufpd_avx(auVar8,auVar8,1);
      dVar9 = (double)(int)(((uVar3 >> 0x14) - 0x3ff) + iVar4);
      return dVar9 * 0.6931471803691238 +
             dVar5 + ((dVar6 * (dVar1 + auVar8._0_8_ + auVar7._0_8_) +
                      dVar9 * 1.9082149292705877e-10) - dVar1);
    }
  }
  else {
    dVar5 = __x;
    if (uVar3 >> 0x14 < 0x7ff) {
      iVar4 = 0;
      dVar5 = 0.0;
      if ((ulong)__x >> 0x20 != 0x3ff00000 || ((ulong)__x & 0xffffffff) != 0) goto LAB_01145af2;
    }
  }
  return dVar5;
}

So (as expected, since AFAIK it doesn't exist in AVX) there's no hardware magic going on to calculate the logarithm.

---

@log2 of the same number is also inaccurate, but not to nearly the same extent:

expected -1.0269365834968912e3
actual   -1.0269365831820458e3
                     ^

i have worked around this in my own program:

pub fn ln(x: anytype) @TypeOf(x) {
    const frexp = std.math.frexp(x);
    return std.math.log1p(frexp.significand * 2 - 1) + std.math.ln2 * @as(@TypeOf(x), @floatFromInt(frexp.exponent - 1));
}

pub fn log2(x: anytype) @TypeOf(x) {
    const frexp = std.math.frexp(x);
    return @log2(frexp.significand) + @as(@TypeOf(x), @floatFromInt(frexp.exponent));
}

pub fn log(T: type, base: T, x: T) T {
    return log2(x) / log2(base);
}