[libc-commits] [libc] [libc][stdfix] Implement fixed point fxdivi functions in llvm-libc (PR #210020)

[Mikhail R. Gadelha via libc-commits]

https://github.com/mikhailramalho requested changes to this pull request.

Thanks for the PR. Coincidently, I was talking with @lntue to formally verify the fixed-point calls and this PR fixes 3 of the 4 issues I've found. What remains is confined to the two instantiations where `WideFXType` *is* `FXType` and there is no headroom above the result format: `lkdivi` and `ulkdivi` (the power-of-two case below is `lkdivi` only).

1. The overflow check doesn't work at `(unsigned) long accum`: both the NR branch and the power-of-two branch fail here.
```
lkdivi (25769803776,    3)   = -2863311530.00   true  8589934592.00  -> want MAX
ulkdivi(25769803776,    3)   =  2863311530.50   true  8589934592.00  -> want MAX
lkdivi (4294967296,     1)   = -4294967296.00   true  4294967296.00  -> want MAX
lkdivi (1026497183760, -128) =   570425343.87   true -8019509248.13  -> want MIN
```

We should probably do something like:
```cpp
// before `WideFXType res;`
if constexpr (OutRep::INTEGRAL_LEN > 0) {
  constexpr UInt128 MAX_UNITS = static_cast<UInt128>(1) << OutRep::INTEGRAL_LEN;
  if (LIBC_UNLIKELY(static_cast<UInt128>(n_mag) >=
                    MAX_UNITS * static_cast<UInt128>(d_mag)))
    return result_is_negative ? OutRep::MIN() : OutRep::MAX();
}
```
I tested this locally, and I could no longer find a counterexample.

2. The `testEdgeCases` test cannot catch either bug, so it should be updated to something like:
```cpp
if constexpr (has_integral) {
  constexpr IntType over_max = static_cast<IntType>(6) *
                               (static_cast<IntType>(1) << FXRep::INTEGRAL_LEN);
  EXPECT_EQ(func(over_max, 3), fx_max);

  constexpr IntType at_max = static_cast<IntType>(1) << FXRep::INTEGRAL_LEN;
  EXPECT_EQ(func(at_max, 1), fx_max);
}
```

3. Also, the accuracy at `lkdivi`/`ulkdivi` is much weaker than the comments imply.

Measured against an exact reference (`n, d <= 256`):

| entry points | max error |
|---|---|
| `rdivi`, `urdivi`, `lrdivi`, `ulrdivi`, `kdivi`, `ukdivi` | 1 ulp |
| `lkdivi` | 43 ulp |
| `ulkdivi` | 22 ulp |

Six of eight are effectively correctly rounded. For the two no-headroom types
the error grows with the quotient, always worst at `d == 3`:

```
n,d <=   64 :  11 ulp        n,d <= 1024 : 171 ulp
n,d <=  256 :  43 ulp        n,d <= 4096 : 683 ulp
```

Absolute error scales linearly with the quotient, so the relative error is
constant at about 2.33e-10, which is half an ulp of `2^-31`. That is the
representation limit of `recip` itself: it is stored in `WideFXType` with 31
fraction bits, so it carries about half an ulp of error no matter how many
iterations refine it.

A consequence worth flagging: the `E4 <= 2.055e-20` on line 336 is
unreachable at these types. The fourth iteration guarded by `F >= 31` costs a
multiply and an add and cannot improve a reciprocal that is already at its
representable precision, the measured error matches `E3` territory and stays
there. Either the comment should say that the bound applies to the reciprocal
before rounding rather than to the result, or the iteration could be dropped
for the formats where it provably cannot help.

If a tighter result is wanted at these two types, the reciprocal needs more
fraction bits than the output format has, which is exactly the headroom
`WideFXType` provides everywhere else.

https://github.com/llvm/llvm-project/pull/210020