math.sin-cos

// Derived from fdlibm e_rem_pio2.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.

import { type SinCos } from "./math_sin_cos_types.ts";

// pi/2 in two parts (fdlibm): PIO2_1 has 33 significant bits, so k * PIO2_1 is
// exact for every k up to 2^20, which covers the permitted range.
const PIO2_1 = 1.57079632673412561417e+00;
const PIO2_1T = 6.07710050650619224932e-11;
const INV_PIO2 = 6.36619772367581382433e-01;
const LIMIT = 1000000;

// sin and cos on |r| <= pi/4 by Taylor series, which converge there to well
// below an ulp. Coefficients are correctly rounded 1/n!.
function sinSmall(r: number): number {
  const s = r * r;
  let p = 1 / 355687428096000;
  p = -1 / 1307674368000 + s * p;
  p = 1 / 6227020800 + s * p;
  p = -1 / 39916800 + s * p;
  p = 1 / 362880 + s * p;
  p = -1 / 5040 + s * p;
  p = 1 / 120 + s * p;
  p = -1 / 6 + s * p;
  return r + r * s * p;
}

function cosSmall(r: number): number {
  const s = r * r;
  let p = -1 / 6402373705728000;
  p = 1 / 20922789888000 + s * p;
  p = -1 / 87178291200 + s * p;
  p = 1 / 479001600 + s * p;
  p = -1 / 3628800 + s * p;
  p = 1 / 40320 + s * p;
  p = -1 / 720 + s * p;
  p = 1 / 24 + s * p;
  p = -1 / 2 + s * p;
  return 1 + s * p;
}

/**
 * Sine and cosine using only +, -, * and /, which IEEE 754 defines exactly.
 *
 * Math.sin and Math.cos come from each platform's maths library and may differ
 * in the last bit between JavaScript, Python and Rust; after rounding, a value
 * next to a boundary can then go either way (an SVG arc end at x.xx5). The
 * angle is reduced to a quarter turn k and a remainder |r| <= pi/4 with a
 * two-part pi/2 (Cody and Waite), and the series above do the rest.
 */
export function sinCos(radians: number): SinCos {
  if (typeof radians !== "number" || !Number.isFinite(radians) || radians < -LIMIT || radians > LIMIT) {
    throw new RangeError(`radians must be a finite number between -1000000 and 1000000, received ${radians}`);
  }
  const k = Math.floor(radians * INV_PIO2 + 0.5);
  const r = radians - k * PIO2_1 - k * PIO2_1T;
  const s = sinSmall(r);
  const c = cosSmall(r);
  const quadrant = ((k % 4) + 4) % 4;
  if (quadrant === 0) return { sin: s + 0, cos: c + 0 };
  if (quadrant === 1) return { sin: c + 0, cos: -s + 0 };
  if (quadrant === 2) return { sin: -s + 0, cos: -c + 0 };
  return { sin: -c + 0, cos: s + 0 };
}