Skip to content

Instantly share code, notes, and snippets.

@teryror
Last active May 5, 2021 11:55
Show Gist options
  • Save teryror/3992cb34ca6531c2e9d2c5812197debb to your computer and use it in GitHub Desktop.
Save teryror/3992cb34ca6531c2e9d2c5812197debb to your computer and use it in GitHub Desktop.
// 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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment