dotnet / dotnet/runtime

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

Open
#123,325 3 comments 0 reactions 0 assignees View on GitHub
area-System.Numerics
Dominant language
C#
Stars
18.3k
Forks
5.6k
PR merge metrics
PR metrics pending

Description

### 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

```csharp
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> ExpM1Internal_Double_A0_RI5 => ExpM1Internal_Double_A0_RI5_;

private static readonly Vector128[] 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> ExpM1Internal_Double_A1_RI5 => ExpM1Internal_Double_A1_RI5_;

private static readonly Vector128[] 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> 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_

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.