Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add FRAC_1_SQRT_2PI constant #383

Open
sunsided opened this issue May 26, 2024 · 8 comments
Open

Add FRAC_1_SQRT_2PI constant #383

sunsided opened this issue May 26, 2024 · 8 comments
Labels
api-change-proposal A proposal to add or alter unstable APIs in the standard libraries T-libs-api

Comments

@sunsided
Copy link

sunsided commented May 26, 2024

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 rust-lang/rust#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 rust-lang/rust#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:

Related proposals:

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.
@sunsided sunsided added api-change-proposal A proposal to add or alter unstable APIs in the standard libraries T-libs-api labels May 26, 2024
@the8472
Copy link
Member

the8472 commented May 26, 2024

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
Copy link
Member

kennytm commented May 26, 2024

@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
Copy link

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
Copy link
Author

sunsided commented May 27, 2024

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

@kennytm
Copy link
Member

kennytm commented May 27, 2024

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
Copy link
Member

joshtriplett commented May 28, 2024

(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
Copy link
Member

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
Copy link

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
api-change-proposal A proposal to add or alter unstable APIs in the standard libraries T-libs-api
Projects
None yet
Development

No branches or pull requests

6 participants