impl/typescript.ts

// Derived from fdlibm e_exp.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: k * LN2_HI is exact for every k this range needs.
const LN2 = 0.6931471805599453;
const LN2_HI = 6.93147180369123816490e-01;
const LN2_LO = 1.90821492927058770002e-10;

/**
 * e^x using only +, -, * and /, which IEEE 754 defines exactly.
 *
 * Math.exp's last bit may differ between platforms, which is enough to split
 * two languages at a rounding boundary. Here x = k ln 2 + r with |r| <= ln 2 / 2
 * (Cody and Waite's two-part reduction), e^r is a Taylor series to r^16, and
 * 2^k is applied by exact doubling or halving.
 */
export function exp(x: number): number {
  if (typeof x !== "number" || !Number.isFinite(x) || x < -708 || x > 709) {
    throw new RangeError(`exp needs x between -708 and 709, received ${x}`);
  }
  const k = Math.floor(x / LN2 + 0.5);
  const r = x - k * LN2_HI - k * LN2_LO;
  // 1/n! as correctly rounded divisions of exact integers.
  let p = 1 / 20922789888000;
  p = 1 / 1307674368000 + r * p;
  p = 1 / 87178291200 + r * p;
  p = 1 / 6227020800 + r * p;
  p = 1 / 479001600 + r * p;
  p = 1 / 39916800 + r * p;
  p = 1 / 3628800 + r * p;
  p = 1 / 362880 + r * p;
  p = 1 / 40320 + r * p;
  p = 1 / 5040 + r * p;
  p = 1 / 720 + r * p;
  p = 1 / 120 + r * p;
  p = 1 / 24 + r * p;
  p = 1 / 6 + r * p;
  p = 1 / 2 + r * p;
  p = 1 + r * p;
  let result = 1 + r * p;
  for (let i = 0; i < k; i++) result *= 2;
  for (let i = 0; i > k; i--) result /= 2;
  return result;
}