Function nikisas::ln

pub fn ln(x: f32) -> f32

Computes natural logarithm of a number.

Notes

Theoretical input domain is (0, max(f32)] ≈ (0, 3.40282347e+38], but near zero the values get quite inaccurate.

Examples

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

Implementation details

First, special cases are handled. If x is 1, then the result is simply 0. If x is near Euler's number, then the result is simply 1. Otherwise, the input x is decomposed into real y and integer k such that

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

At this point, we might use

ln(x) = n * ln(2) + ln(y)

but that would lead to catastrophic cancellation for cases where n = -1 and y ≈ 2. Therefore, we adjust the decomposition if y > sqrt(2) as follows

  y <- y / 2
  n <- n + 1

This keeps the equality x = y * 2^n true, but avoids mentioned cancellation and shifts y to be symmetric around logarithm's root at 1, in range [1/sqrt(2), sqrt(2)].

For reasons unclear to me, polynomial approximation of ln(1 + z) for -a/2 ≤ z < a/2 is easier than of ln(z) for 1 - a/2 ≤ z < 1 + a/2. Thus, we set z = y - 1 and approximate ln(1 + z) using a polynomial in the form:

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

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

The reconstruction then follows already mentioned identity:

  ln(x) = n * ln(2) + ln(y) = n * ln(2) + ln(1 + z)