Double.ExpM1 loses sign, precision, and underflows to zero for small inputs

[LEI-Hongfaan]

Description

double.ExpM1 currently exhibits three related correctness problems for small inputs:

• Signed-zero loss: ExpM1(-0.0) returns +0.0 instead of -0.0.
• Precision loss for tiny inputs: For small positive and negative x near zero the implementation returns values with much-larger-than-expected ULP error compared to a high-quality reference.
• Underflow-to-zero for tiny negatives and positives: Tiny inputs that should produce a nonzero result are returned as 0.0, losing both magnitude and sign.

These behaviors break IEEE 754 expectations for expm1.

Reproduction Steps

using System.Diagnostics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;

public static class ImprovedMath {

    // rational interpolation of ((e^x - 1) / x) on [0, 0.5], order (5, 5)
    internal static ReadOnlySpan<Vector128<double>> ExpM1Internal_Double_A0_RI5 => ExpM1Internal_Double_A0_RI5_;

    private static readonly Vector128<double>[] ExpM1Internal_Double_A0_RI5_ = [
        Vector128.Create([3.4075547183307054E-06D, -1.6256158393388274E-05D]),
        Vector128.Create([0.0001342714090530565D, 0.00058772540449553428D]),
        Vector128.Create([0.0014218291431196866D, -0.0096368893730589423D]),
        Vector128.Create([0.030639089470444997D, 0.08840746548950737D]),
        Vector128.Create([0.051129914628542009D, -0.44887008537145801D]),
        Vector128.Create([1D, 1D]),
    ];

    // rational interpolation of ((e^x - 1) / x) on [-0.5, 0], order (5, 5)
    internal static ReadOnlySpan<Vector128<double>> ExpM1Internal_Double_A1_RI5 => ExpM1Internal_Double_A1_RI5_;

    private static readonly Vector128<double>[] ExpM1Internal_Double_A1_RI5_ = [
        Vector128.Create([2.6538062829946136E-06D, -2.002143187522963E-05D]),
        Vector128.Create([0.00011922976878158256D, 0.00067797260132220947D]),
        Vector128.Create([0.0011060597069885104D, -0.010583989992818755D]),
        Vector128.Create([0.030007938899324207D, 0.093457753733520998D]),
        Vector128.Create([0.039767036998273018D, -0.460232963001727D]),
        Vector128.Create([1D, 1D]),
    ];

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static double ExpM1Internal_Double_A(double x, ReadOnlySpan<Vector128<double>> coeffs) {
        Debug.Assert(coeffs.Length > 0);
        Debug.Assert(-0.5 <= x && x <= 0.5);
        var xx = Vector128.Create(x);
        var f = MemoryMarshal.GetReference(coeffs); // coeffs[0]
        for (int i = 1; i < coeffs.Length; i++) {
            f = Vector128.FusedMultiplyAdd(f, xx, coeffs[i]);
        }
        return f.GetElement(0) / f.GetElement(1);
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static double ExpM1Internal_Double_B(double x) {
        Debug.Assert(-0.5 <= x && x <= 0.5);
        var f = ExpM1Internal_Double_A(x, double.IsNegative(x) ? ExpM1Internal_Double_A1_RI5 : ExpM1Internal_Double_A0_RI5);
        return x * f;
    }

    public static double ExpM1(double x) {
        return double.Abs(x) <= 0.5D ? ExpM1Internal_Double_B(x) : double.Exp(x) - 1D;
    }
}

internal class Program {

    static void Main(string[] args) {

        static void TestCase(double x) {
            Console.WriteLine($"System:   {x,24:G17} -> {double.ExpM1(x),24:G17}");
            Console.WriteLine($"ThisImpl: {x,24:G17} -> {ImprovedMath.ExpM1(x),24:G17}");
        }

        TestCase(-0.0D);
        TestCase(1E-3D);
        TestCase(1E-4D);
        TestCase(-1E-5D);
        TestCase(1E-6D);
        TestCase(-1E-16D);
        TestCase(-1E-17D);
        TestCase(1E-100D);
        TestCase(-double.Epsilon);
        TestCase(0.5D);
        TestCase(-0.5D);
        TestCase(1D);
        TestCase(-1D);
    }
}

Expected behavior

.

Actual behavior

System:                         -0 ->                        0
ThisImpl:                       -0 ->                       -0
System:                      0.001 ->    0.0010005001667083846
ThisImpl:                    0.001 ->    0.0010005001667083414
System:                     0.0001 ->    0.0001000050001667141
ThisImpl:                   0.0001 ->   0.00010000500016667085
System:    -1.0000000000000001E-05 ->  -9.9999500001723973E-06
ThisImpl:  -1.0000000000000001E-05 ->   -9.999950000166668E-06
System:     9.9999999999999995E-07 ->   1.0000004999621837E-06
ThisImpl:   9.9999999999999995E-07 ->   1.0000005000001665E-06
System:    -9.9999999999999998E-17 ->  -1.1102230246251565E-16
ThisImpl:  -9.9999999999999998E-17 ->  -9.9999999999999998E-17
System:    -1.0000000000000001E-17 ->                        0
ThisImpl:  -1.0000000000000001E-17 ->  -1.0000000000000001E-17
System:                     1E-100 ->                        0
ThisImpl:                   1E-100 ->                   1E-100
System:   -4.9406564584124654E-324 ->                        0
ThisImpl: -4.9406564584124654E-324 -> -4.9406564584124654E-324
System:                        0.5 ->      0.64872127070012819
ThisImpl:                      0.5 ->      0.64872127070012808
System:                       -0.5 ->     -0.39346934028736658
ThisImpl:                     -0.5 ->     -0.39346934028736652
System:                          1 ->       1.7182818284590451
ThisImpl:                        1 ->       1.7182818284590451
System:                         -1 ->     -0.63212055882855767
ThisImpl:                       -1 ->     -0.63212055882855767

Regression?

No response

Known Workarounds

No response

Configuration

No response

Other information

No response

[MihuBot]

I'm a bot. Here are possible related and/or duplicate issues (I may be wrong):

• Implementations of Double.ExpM1, Double.LogP1 & co are not correct for arguments close to 0 #77979
• Implementations of Double.ExpM1, Double.LogP1 & co are not correct for arguments close to 0 (reminder) #115399

[dotnet-policy-service]

Tagging subscribers to this area: @dotnet/area-system-numerics
See info in area-owners.md if you want to be subscribed.