Function nikisas::exp

pub fn exp(x: f32) -> f32

Computes exponentiation function of a number.

Notes

The input domain is limited to approximately [ln(min(positive f32)), ln(max(f32))] ≈ [-87.3, 88.7] due to limits of machine representation.

Example

use nikisas::{exp, consts::E};
assert_eq!(exp(1.0), E);

Implementation details

First, special cases are handled. If x is 1, then the result is simply Euler's number. If x is near zero, then the result is simply 1. Otherwise, input x is reduced to an integer k and real z such that

  x = k * ln(2) + z and |z| ≤ ln(2) / 2

Exponentiation of z is done using polynomial in the form:

  exp(z) ≈ 1 + z + 1/2 * z^2 + z^3 * P(z)

The "prefix" corresponds to coefficients of low-degree Taylor polynomial of exp(z) for z = 0 and P is found using special minimax algorithm in Sollya.

The reconstruction follows this identity:

  exp(x) = exp(k * ln(2) + z) = exp(ln(2))^k * exp(z) = 2^k * exp(z)

Computation of exp(z) is done using aforementioned polynomial approximation and multiplying by 2^k can be implemented exactly using bit manipulation of floating point number representation.