Function fpsincos

src/big-decimal/utils.ts:687–824

(x: bigint, bits: number)

 * @returns  [sin(x/2^bits)*2^bits, cos(x/2^bits)*2^bits]
 */
export function fpsincos(x: bigint, bits: number): [bigint, bigint] {
  const B = BigInt(bits);
  const scale = 1n << B;
  // sin(0) = 0, cos(0) = 1
  if (x === 0n) return [0n, scale];

  const pi = fppi(bits);
  const twoPi = 2n * pi;
  const halfPi = pi / 2n;

  // Step 1: Reduce modulo 2π to [0, 2π)
  // For large arguments, x % twoPi loses precision because x has many
  // more bits than twoPi. Use extended precision for the reduction.
  let r: bigint;
  const absX = bigintAbs(x);
  if (absX > scale << 30n) {
    // Large argument: compute at extended precision.
    // Extra guard bits: bitLength(|x/scale|) + 64.
    const extraBits = bitLength(absX) - bits + 64;
    const extBits = bits + extraBits;
    const extX = x << BigInt(extraBits);
    const extPi = fppi(extBits);
    const extTwoPi = 2n * extPi;

    let extR = extX % extTwoPi;
    if (extR < 0n) extR += extTwoPi;

    // Scale back
    r = extR >> BigInt(extraBits);
  } else {
    r = x % twoPi;
  }
  if (r < 0n) r += twoPi;

  // Step 2: Quadrant reduction to [0, π/2]
  // Determine quadrant and adjust sign
  let sinSign = 1n;
  let cosSign = 1n;

  if (r > 3n * halfPi) {
    // Quadrant 4: [3π/2, 2π) → sin negative, cos positive, use 2π - r
    r = twoPi - r;
    sinSign = -1n;
  } else if (r > pi) {
    // Quadrant 3: [π, 3π/2] → sin negative, cos negative, use r - π
    r = r - pi;
    sinSign = -1n;
    cosSign = -1n;
  } else if (r > halfPi) {
    // Quadrant 2: [π/2, π] → sin positive, cos negative, use π - r
    r = pi - r;
    cosSign = -1n;
  }
  // else Quadrant 1: [0, π/2] — no change

  // Step 3: Double-angle halving
  // Each halving is cheap (integer divide by 2) but reconstruction costs O(M(p)).
  // Each Taylor term also costs O(M(p)). Balance the two for O(√p) total steps.