area_of_super_ellipse.rs

// Calculate the area of a generalized ellipse (abs(x/a)^n + abs(y/b)^n = 1).
// I'm considering this as the metric to minimize for a weird optimization
// problem, where this would be evaluated _a lot_, so I tried making this
// as fast & simple as possible.

pub fn area_of_superellipse(exponent: f64, semi_axes: (f64, f64)) -> f64 {
    let n = exponent;
    let (a, b) = semi_axes;

    // This method was derived from the following formula:
    // A(n, a, b) = 4ab * Gamma(1 + 1/n)^2/Gamma(1 + 2/n),
    // where Gamma, the real extension of the factorial function,
    // is approximated with Stirling's formula,
    // n! ~ sqrt(2pi * n) * (n/e)^n.
    //
    // Simplifying to minimize calls to sqrt() and powf(), we get:

    // Precomputed constant factor 4 * sqrt(2pi) / e
    const CONST_TERM: f64 = 3.6885480355831564;

    // n and (n + 2) appear as divisors in multiple places. Factoring them out
    // like this allows us to replace all other divisions with multiplications.
    let one_over_n_plus_two = 1.0 / (n + 2.0);
    let one_over_n = 1.0 / n;

    let pq_term = one_over_n_plus_two * one_over_n.powi(2) * (n + 1.0).powi(3);
    let sqrt_term = (one_over_n_plus_two * n).sqrt();
    let powf_term = (one_over_n_plus_two * (n + 1.0)).powf(2.0 * one_over_n);

    CONST_TERM * a * b * pq_term * sqrt_term * powf_term
}