Skip to content

Add FRAC_1_SQRT_2PI constant #383

Closed
Closed
@sunsided

Description

@sunsided

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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    ACP-acceptedAPI Change Proposal is accepted (seconded with no objections)I-libs-api-nominatedIndicates that an issue has been nominated for discussion during a team meeting.T-libs-apiapi-change-proposalA proposal to add or alter unstable APIs in the standard libraries

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions