[][src]Function nikisas::pow

pub fn pow(x: f32, p: f32) -> f32

Computes a number raised to a power.

Notes

For negative bases, the power must be an integer. If it's not, then NaN is returned.

Approximation of power function is hard. A straightforward way would be to use the following identity:

  x^p = exp(ln(x^p)) = exp(p * ln(x))

but that would produce very inaccurate results in general, since even a tiny error in calculation of natural logarithm is exponentiated afterwards (not to say that calculation exponentiation has also some error).

We fallback to partially using square-and-multiply algorithm for doing power on integral parts of the exponent, which might be up to 32 iterations of a small loop. This is a bit more cycles than doing some fancy approximation using a fixed-degree polynomial or rational function, but in most cases should be fine.

Examples

use nikisas::pow;
assert_eq!(pow(2.0, -1.0), 0.5);

Implementation details

First, special cases are handled:

  • if x is near 1, then the result is simply 1,
  • if p is near 1, then the result is simply x,
  • if p is near 0, then the result is simply 1,
  • if x is near 2, then specialized pow2 is used, and
  • if x is near 10, then specialized pow10 is used.

If x is non-negative, the procedure goes like this. First, x is decomposed to real y and integer n, such that

  x = y * 2^n, where 1 ≤ y < 2

Second, p and q = p * n are decomposed as follows:

  p = pi + pf, such that pi is integer and 0 ≤ pf < 1
  q = qi + qf, such that qi is integer and |qf| ≤ 1/2

With all this, we can then use the following identity:

  x^p = (y * 2^n)^p
      = y^p * 2^(pn)
      = y^(pi + pf) * 2^(qi + qf)
      = y^pi * y^pf * 2^qi * 2^qf

For y^pi we use square-and-multiply loop algorithm, for y^pf the aforementioned identity exp(pf * ln(y)) is used with hope that it does not introduce too big error as it is only one term in the whole computation, for 2^qf we use pow2 routine, and multiplying by 2^qi can be implemented exactly using bit manipulation of floating point number representation.

If x is negative, the p must be an integer. This is true when z is zero, where z is the fractional part of p = k + z. If this is a case, we again decompose x into x = y * 2^n. Then the same procedure as before is used, only that the fractional exponents don't exist, thus

  y^pf = 2^qf = 1, y^pi = y^k, 2^qi = 2^(k * n)