impl/typescript.ts

// Derived from fdlibm s_atan.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_atan.c (Sun Microsystems, 1993). atan of 0.5, 1, 1.5 and infinity,
// each split into a high and a low part so the sum carries more than 53 bits.
const ATAN_HI = [
  4.63647609000806093515e-01,
  7.85398163397448278999e-01,
  9.82793723247329054082e-01,
  1.57079632679489655800e+00,
];
const ATAN_LO = [
  2.26987774529616870924e-17,
  3.06161699786838301793e-17,
  1.39033110312309984516e-17,
  6.12323399573676603587e-17,
];
const AT = [
  3.33333333333329318027e-01,
  -1.99999999998764832476e-01,
  1.42857142725034663711e-01,
  -1.11111104054623557880e-01,
  9.09088713343650656196e-02,
  -7.69187620504482999495e-02,
  6.66107313738753120669e-02,
  -5.83357013379057348645e-02,
  4.97687799461593236017e-02,
  -3.65315727442169155270e-02,
  1.62858201153657823623e-02,
];
const TWO_POW_66 = 73786976294838206464;
const TWO_POW_M29 = 1.862645149230957e-9;

/**
 * Arctangent in radians using only +, -, * and /, in the same order as the
 * Python and Rust versions, so all three return the same double. Math.atan
 * comes from each platform's maths library and may differ in the last bit.
 *
 * fdlibm's method: reduce |x| against atan(0.5), atan(1), atan(1.5) or pi/2,
 * then an odd polynomial on the remainder. fdlibm picks the interval from the
 * high word of the double; comparing |x| with the same boundaries (7/16,
 * 11/16, 19/16, 39/16, 2^66, 2^-29) is exactly equivalent.
 */
export function atan(x: number): number {
  if (typeof x !== "number" || !Number.isFinite(x)) {
    throw new RangeError(`x must be a finite number, received ${String(x)}`);
  }
  const negative = x < 0;
  const a = negative ? -x : x;
  if (a >= TWO_POW_66) {
    const r = ATAN_HI[3] + ATAN_LO[3];
    return negative ? -r : r;
  }
  let i: number;
  let t = x;
  if (a < 0.4375) {
    // atan(x) rounds to x here; + 0 turns -0 into 0.
    if (a < TWO_POW_M29) return x + 0;
    i = -1;
  } else {
    t = a;
    if (a < 1.1875) {
      if (a < 0.6875) {
        i = 0;
        t = (2.0 * t - 1.0) / (2.0 + t);
      } else {
        i = 1;
        t = (t - 1.0) / (t + 1.0);
      }
    } else if (a < 2.4375) {
      i = 2;
      t = (t - 1.5) / (1.0 + 1.5 * t);
    } else {
      i = 3;
      t = -1.0 / t;
    }
  }
  const z = t * t;
  const w = z * z;
  const s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))));
  const s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))));
  if (i < 0) return t - t * (s1 + s2);
  const r = ATAN_HI[i] - ((t * (s1 + s2) - ATAN_LO[i]) - t);
  return negative ? -r : r;
}