sunsided commented 1 month ago

Proposal

Problem statement

Statistical calculations, image and signal processing applications commonly need to calculate a scaling factor of 1/√(2π) for the Gaussian probability distribution (specifically, 1/(σ√(2π)) where sigma is application specific).

$$ f(x) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right) $$

The standard normal distribution is the same with trivially σ=1 (and x=y):

$$ f(x) = \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{x^2}{2} \right) $$

With the current constants, the calculation of the scaling term requires $\frac{1}{\sqrt{2}} \times \frac{1}{\sqrt{\pi}}$, i.e.:

FRAC_1_SQRT_PI * FRAC_1_SQRT_2 / sigma

I propose to add FRAC_1_SQRT_2PI (as 1/√(2π)) in order to simplify this to

FRAC_1_SQRT_2PI / sigma

reducing clutter and keeping mathematical precision.

This would go in line with #103883 and https://github.com/rust-lang/rust/pull/123850 which add the FRAC_1_SQRT_PI constant (1/√π).

Motivating examples or use cases

The image processing library imgproc uses this code for calculating a blur filter:

/// The gaussian probability density function with a mean of zero.
#[inline]
fn gaussian_pdf(x: f32, sigma: f32) -> f32 {
    // Equation derived from <https://en.wikipedia.org/wiki/Gaussian_function>
    (sigma * (2.0 * f32::consts::PI).sqrt()).recip() * (-x.powi(2) / (2.0 * sigma.powi(2))).exp()
}

This example is less about the choice of a functional implementation style, but it showcases the need for this particular calculation. (Do note that since this is an image processing application, the parameter sigma must be freely settable and at least parts of the scaling term must be determined ad-hoc.)

The Normal Distribution (Gaussian PDF), Bayesian Statistics as a whole, Error Function calculation and Fourier Transformations are all centered around the same fundamental calculation.

In financial applications, the Black-Scholes model comes up. As an example, the black_scholes_rust crate implements it like this.

fn inc_norm(x: f64) -> f64 {
    (-x.powi(2) / 2.0).exp() / (PI.sqrt() * SQRT_2)
}

In physics, the Maxwell-Boltzmann equation has the term as $(\bullet)^\frac{3}{2}$ (.powf(1.5)).

Solution sketch

Trivially, adding FRAC_1_SQRT_2PI as 1/√(2π) = 0.3989422804014326779399460599343818684758586311649346576659258296… (from here). I went ahead and drafted it in https://github.com/rust-lang/rust/pull/125253.

Alternatives

A combination of FRAC_1_SQRT_PI and FRAC_1_SQRT_2 can be used, however at more visual clutter. Applications can make use of const evaluation and reuse their own constant, but this may be subject to floating point processing on the individual machine (arguably; I did not analyze the error potential here).

In some cases (e.g. signal processing) the entire calculation can be done "on paper", e.g. when the standard deviation parameter sigma will never change in signal processing applications or image filtering.

Bivariate normal distributions do not have this problem, so this focuses on the most common Gaussians.

Links and related work

Implementation of the above:

https://github.com/rust-lang/rust/pull/125253

Related proposals:

352

Previous work on constants:

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.

the8472 commented 1 month ago

I'm wondering if it'd be possible have some macro that evaluates things infinite precision and spits out a float constant. E.g. fconst!(1 / sqrt(2 * PI), f32).

kennytm commented 1 month ago

@the8472 you don't need infinite precision, doing the computation in a higher precision than f32/f64/f128 should be enough.

But please don't make fconst! a compiler or std feature, there is no point bringing in e.g. internal f256 support for everyone, just to let some users be able to define arbitrary mathematical constants. They can easily be calculated properly outside the compiler.

clarfonthey commented 1 month ago

I'm fine with this constant, although the hill I will die on (for my corpse to rot away and be ignored) is that it should be FRAC_1_SQRT_TAU for consistency with the existing TAU constant.

I mean this in the sense of, I would like this change in naming, but absolutely will not block it over this.

sunsided commented 1 month ago

@clarfonthey I see your point. What's your take on (renaming and) adding #[doc(alias = "FRAC_1_SQRT_2PI")]?

kennytm commented 1 month ago

One big reason why TAU was adapted (aside from the tau-manifesto) is that 2PI is not a valid identifier, PI_2 is wrong (easily mistaken as $\pi/2$), and TWO_PI is very ugly. There is also a precedence from Python. These all do not matter for FRAC_1_SQRT_2PI.

And practically everyone used $1/\sqrt{2\pi}$ rather than $1/\sqrt\tau$. So I think making FRAC_1_SQRT_TAU the canonical name will hurt its discoverability.

joshtriplett commented 1 month ago

(Speaking for myself here.) FRAC_1_SQRT_2PI seems like the right naming here. Most people still use pi, and an alias seems suboptimal. By all means add an alias for the TAU version.

joshtriplett commented 1 month ago

We discussed this in today's @rust-lang/libs-api meeting. We observed that this constant doesn't currently exist among C's constants, but we were still sympathetic to adding this. We're inclined to add this, but we'd like a second opinion.

I think @m-ou-se has worked with the C constants before. Do you have an opinion on this?

quaternic commented 1 month ago

I propose to add FRAC_1_SQRT_2PI (as 1/√(2π)) in order to simplify this to
FRAC_1_SQRT_2PI / sigma
reducing clutter and keeping mathematical precision.

Two existing alternatives:

// You could just define the const locally
const FRAC_1_SQRT_2PI: f32 = FRAC_1_SQRT_2 * FRAC_1_SQRT_PI;

// Or write the obvious thing; LLVM will optimize it just as well
(2.0 * PI).sqrt().recip() / sigma

All of *, /, sqrt should be correctly rounded (so there should be no machine dependence), and for the above expressions (for both f32 and f64), those match the result without intermediate rounding. https://rust.godbolt.org/z/KaPGvvx9W

That's partially by chance, since for different real values 1 / sqrt(x) (as computed in f32) can be off by up to 2ULP (from the exact result rounded to f32).

In principle, I'm mildly against adding more constants that are just simple expressions, since it clutters (even if only slightly) the library and only helps the cases that need that specific constant. Once you want to further map that constant, the benefit over just writing the expression mostly disappears. We have a language for composing expressions already, so I'd rather see and work with something like math!(1 / sqrt(2 * PI)) as f32 than parse and distinguish between FRAC_1_SQRT_2PI, FRAC_2_SQRT_PI, or FRAC_1_SQRT_PI.

rust-lang / libs-team

Add `FRAC_1_SQRT_2PI` constant #383