More C23 math fun today. I finished up nextup/nextdown over the weekend, which included fixing the 80-bit nextafter implementation.

Today, I wrote exp10m1 and exp2m1. I was surprised how easy those turned out to be:

/*
* Compute N ** x - 1:
*
* p = log(N) * x
* r_hi = expm1(p.hi) + 1
* r_lo = expm1(p.lo) + 1
* r = r_hi * r_lo
* m = (r - 1)
*/

The math is all done with Veltkamp-split values. The result is within 1 ulp.

Oh, I've got a couple of test cases where glibc is off by > 1ulp for exp10m1:

test exp10m1
298 -0x1.b066311e8d604p-41 got -0x1.f1d165e627b95p-40 want -0x1.f1d165e627b93p-40 ulp 2
322 -0x1.b64c6c0adc35cp-34 got -0x1.f89c1d613775ap-33 want -0x1.f89c1d6137758p-33 ulp 2
355 0x1.a5b7a35f7d50dp-24 got 0x1.e585240bca4d9p-23 want 0x1.e585240bca4d7p-23 ulp 2

@keithp worst known is 3.54 ulp: https://inria.hal.science/hal-03141101

@keithp I assume you cannot do directed rounding for the intermediate ops?

@oblomov I haven't taken things that far yet. Right now, I'm aiming to get all of the functions working within 1 ulp for round-to-nearest-even.

I'm still struggling with whether it's "worth" it to support exact operations and rounding modes for an embedded library. The appeal is strong, but the extra space and time seem hard to justify for tiny devices.

And there's plenty of work to be done to reach the current goal.

Testing exp2m1 and exp2m1f with all 4 billion 32-bit floats found a single failure -- exp2m1(1024) returned NAN instead of INF.

I'm not sure how I'd find things like that any other way.

@oblomov This code needs to run on softfp and those rarely get exceptions or rounding modes -- libgcc has everything required to support those in the source, but it's disabled; I think that's because there's no place to hold the per-CPU state.

@keithp @oblomov IIRC that's the case, plumbing it all together would require adding some ABI between GCC and glibc because the routines are split. I don't think it'd be super hard or anything, it's just that softfp tends to be used in places where people don't care that much about floating point so it's just not anyone's priority.

@oblomov That did make the C23 narrowing operations "fun" -- glibc simply performs the regular operation, and if that generates FE_INEXACT, it does the round-to-odd trick. Picolibc does the operation in double-floats, using the low float to decide when to round-to-odd. Hugely more expensive in space and time, but so cheap in source code.