ACP: Adding `exp2_m1`

[Apersoma]

Proposal

Adding a function for floats that computes 2^x -1.

Problem statement

This solves the same problem as exp_m1 but for powers of 2.

Motivating examples or use cases

I personally was looking for this to use to compare a polynomial approximation of exp2_m1 for 0 <= x < 1 against it to optimize the coefficients, and it is generally use full for the same reasons as exp_m1.

Solution sketch

/// Returns 2^(self) - 1 in a way that is accurate even if the number is close to zero.
///
  /// # Unspecified precision
  ///
  /// The precision of this function is non-deterministic. This means it varies by platform, Rust version, and
  /// can even differ within the same execution from one invocation to the next.
  /// This function currently corresponds to the `expm1` from libc on Unix
  /// and Windows. Note that this might change in the future.
  ///
  /// # Examples
  ///
  /// ```
  /// let x = 1e-16_f64;
  ///
  /// // for very small x, 2^x is approximately 1 + x*ln(2) + (x*ln(2))^2 / 2
  /// let x_ln2 = x * core::f64::const::LN_2;
  /// let approx = x_ln2  + x_ln2  * x_ln2 / 2.0;
  /// let abs_difference = (x.exp_m1() - approx).abs();
  ///
  /// assert!(abs_difference < 1e-20);
  /// ```
pub fn exp2_m1(self) -> f64 {..}

Alternatives

User implementation or importing a crate that implements it.

Links and related work

exp_m1 docs
exp2 docs

What happens now?

This issue contains an API change proposal (or ACP) and is part of the libs-api team feature lifecycle. Once this issue is filed, the libs-api team will review open proposals as capability becomes available. Current response times do not have a clear estimate, but may be up to several months.

Possible responses

The libs team may respond in various different ways. First, the team will consider the problem (this doesn't require any concrete solution or alternatives to have been proposed):

• We think this problem seems worth solving, and the standard library might be the right place to solve it.
• We think that this probably doesn't belong in the standard library.

Second, if there's a concrete solution:

• We think this specific solution looks roughly right, approved, you or someone else should implement this. (Further review will still happen on the subsequent implementation PR.)
• We're not sure this is the right solution, and the alternatives or other materials don't give us enough information to be sure about that. Here are some questions we have that aren't answered, or rough ideas about alternatives we'd want to see discussed.

[tgross35]

Note that this is a C23 math function so support is pretty slim (e.g. no musl yet), so we would need an implementation in our libm first. Which I'm happy to add at any time, regardless of whether it's public in std.

Others in the same boat that we could also consider: compound, exp10, exp10m1, exp2m1, log10p1, log2p1, powr, rootn, rsqrtn, sinpi and similar.

[scottmcm]

Definitely feels like if we add exp2_m1 we should add log2_1p at the same time. I don't see any reason to have one without the other.

That said, if it's really new I think we'd probably be fine waiting for it for a while. Getting it from a non-std libm crate for now seems plenty tolerable, and I don't really know how to pick a subset of C23 stuff meaningfully. Why would exp2_m1 be worth adding but not exp10_m1, for example? But at the same time it's not obvious that we should just "meh, add everything right away" either.

Certainly the motivation of "I want to use this to optimize my polynomial" isn't sufficient justification to rush it in.

Unnatural bases like this don't even have the nice exp_m1(ε) ≈ ε ≈ log_1p(ε) property, sadly.

(TBH, I do individually wish we had sinpi and friends to avoid the f64::PI.sin() problem, but don't know that I can really make a great case that rust should add them right now either...)

[kennytm]

ISTM among all the new C23 math functions, exp2_m1 would be one of those that need least help from the compiler/std.

Function	Formula
exp2m1	$2^x - 1$	(x * LN_2).exp_m1()
exp10	$10^x$	(Should use the llvm.exp10 intrinsic)
exp10m1	$10^x - 1$	(x * LN_10).exp_m1()
log2p1	$\log_2(x + 1)$	x.ln_1p() / LN_2
log10p1	$\log_{10}(x + 1)$	x.ln_1p() / LN_10
rsqrt	$x^{-1/2}$	x.sqrt().recip()
compound	$(1+x)^y$	(x.ln_1p() * y).exp()
rootn	$x^{1/y}$
powr	$e^{y \ln x}$	(x.ln() * y).exp()
acospi	$(\cos^{-1} x) / \pi$
asinpi	$(\sin^{-1} x) / \pi$
atanpi	$(\tan^{-1} x) / \pi$
atan2pi	$(\mathrm{atan2}(x, y)) / \pi$
cospi	$\cos(\pi x)$
sinpi	$\sin(\pi x)$
tanpi	$\tan(\pi x)$

The substitute code for exp2_m1 seems to be precise enough

use std::f64::consts::LN_2;
fn main() {
    let x = 1e-5;
    dbg!((
        6.931495828305653209089800561681e-6, // exactly computed
        (x * LN_2).exp_m1(), // substitute code
        x.exp2() - 1.0, // numerically unstable formula
        x * LN_2 + x * x * LN_2 * LN_2 * 0.5, // quadratic approximation
        x * LN_2, // linear approximation
    ));
}
/* output = (
    6.931495828305653e-6,
    6.931495828305654e-6,
    6.93149582819963e-6,
    6.931495828250149e-6,
    6.931471805599453e-6,
) */

[Apersoma]

The substitute code for exp2_m1 seems to be precise enough

I'd argue your test isn't great because it should be testing for x < EPSILON, but its good there too.

[scottmcm]

(x * LN_2).exp_m1()

Is that ±½ULP, though? With two roundings I'd assume it's not (it'd only be ±1ULP), and because it's so simple people can write that if they don't want it more precise.

For monotonic unary things like exp2_m1, I'd really rather Rust move to committing to ±½ULP, and not move further away from it unless we really have to.

(If it can't practically for non-monotonic non-unary things like atan2 that's one thing, but the "simpler" ones we ought to.)

[Apersoma]

(x * LN_2).exp_m1()

Is that ±½ULP, though? With two roundings I'd assume it's not (it'd only be ±1ULP), and because it's so simple people can write that if they don't want it more precise.

For monotonic unary things like exp2_m1, I'd really rather Rust move to committing to ±½ULP, and not move further away from it unless we really have to.

(If it can't practically for non-monotonic non-unary things like atan2 that's one thing, but the "simpler" ones we ought to.)

It would be far greater than ±1ULP because the second op is raising a number > 1 to a rounded value.

[kennytm]

I'd argue your test isn't great because it should be testing for x < EPSILON, but its good there too.

The test x < EPSILON isn't great because in that region even the linear approximation x * LN2 is precise enough.

Is that ±½ULP, though? With two roundings I'd assume it's not (it'd only be ±1ULP), and because it's so simple people can write that if they don't want it more precise.

Yes, all the substitution will be less precise than the native implementation. But note also that almost every floating-point functions have no ULP guarantee, e.g. https://doc.rust-lang.org/std/primitive.f64.html#method.exp_m1:

• Unspecified precision

  The precision of this function is non-deterministic. This means it varies by platform, Rust version, and can even differ within the same execution from one invocation to the next. [...]

So I don't see how exp2_m1 can enforce "±1 ULP", much less half. If the Rust function just forwards to libc's exp2m1 like everything else, this depends entirely on the system implementation. For instance glibc only guarantees ±9 ULP.

[nia-e]

We're in favour of this, and log2p1 to provide the inverse also. In general the sentiment is in favour of adding more C23 maths functions if there aren't any blockers, so further ACPs are welcome. Thanks :D