impl/typescript.ts

// Derived from fdlibm e_log.c.
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
// Developed at SunSoft, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.

// fdlibm's split of ln 2: the high part has trailing zero bits, so e * LN2_HI
// is exact for any exponent a double can have.
const LN2_HI = 6.93147180369123816490e-01;
const LN2_LO = 1.90821492927058770002e-10;
const SQRT2 = 1.4142135623730951;

/**
 * Natural logarithm using only +, -, * and /, which IEEE 754 defines exactly.
 *
 * Math.log is whatever the platform's maths library does, and its last bit is
 * allowed to differ between JavaScript engines, CPython builds and Rust
 * targets. Rounding the result afterwards does not hide that when a value sits
 * next to a rounding boundary, so capabilities that must agree across
 * languages (a log scale, a colour conversion) use this instead.
 *
 * x = m * 2^e with m in (sqrt(1/2), sqrt(2)], found by exact halving and
 * doubling; then ln m = 2 atanh(s) with s = (m - 1) / (m + 1), |s| < 0.172,
 * summed as an odd series in s to well below an ulp.
 */
export function ln(x: number): number {
  if (typeof x !== "number" || !Number.isFinite(x) || x <= 0) {
    throw new RangeError(`ln is only defined for finite numbers greater than zero, received ${x}`);
  }
  let m = x;
  let e = 0;
  while (m >= 2) {
    m /= 2;
    e += 1;
  }
  while (m < 1) {
    m *= 2;
    e -= 1;
  }
  if (m > SQRT2) {
    m /= 2;
    e += 1;
  }
  const f = m - 1;
  const s = f / (2 + f);
  const z = s * s;
  let p = 1 / 25;
  p = 1 / 23 + z * p;
  p = 1 / 21 + z * p;
  p = 1 / 19 + z * p;
  p = 1 / 17 + z * p;
  p = 1 / 15 + z * p;
  p = 1 / 13 + z * p;
  p = 1 / 11 + z * p;
  p = 1 / 9 + z * p;
  p = 1 / 7 + z * p;
  p = 1 / 5 + z * p;
  p = 1 / 3 + z * p;
  const lnm = 2 * s + 2 * s * (z * p);
  return e * LN2_HI + (e * LN2_LO + lnm);
}