math.fractional-power

/** Every value is held as an integer count of 10^-18. */
export const FIXED_SCALE = 10n ** 18n;

/** Intermediate values beyond 2^512 (about 10^136 after scaling) are refused. */
const MAX_BITS = 512;

const DECIMAL = /^[0-9]+(\.[0-9]{1,18})?$/;

function tooLarge(value: bigint): boolean {
  return value.toString(2).length > MAX_BITS;
}

/** Product of two non-negative fixed-point values, floored. */
export function mulFixed(a: bigint, b: bigint): bigint {
  return (a * b) / FIXED_SCALE;
}

/**
 * x^n for a non-negative fixed-point x, by square-and-multiply from the lowest
 * bit, flooring after every product. The order of operations is part of the
 * contract: it is what makes TypeScript, Python and Rust agree to the last
 * digit. Returns bound + 1 as soon as the result must exceed `bound` (x >= 1
 * only grows), so a bisection never builds an astronomically large number.
 */
export function powFixedBounded(x: bigint, n: number, bound: bigint | null): bigint {
  let result = FIXED_SCALE;
  let base = x;
  let e = n;
  while (e > 0) {
    if (e % 2 === 1) {
      result = mulFixed(result, base);
      if (bound !== null && result > bound && x >= FIXED_SCALE) return bound + 1n;
    }
    e = Math.floor(e / 2);
    if (e > 0) {
      base = mulFixed(base, base);
      if (bound !== null && base > bound && x >= FIXED_SCALE) return bound + 1n;
      if (tooLarge(base)) throw new RangeError("fractional power result too large");
    }
  }
  if (tooLarge(result)) throw new RangeError("fractional power result too large");
  return result;
}

/** x^n in fixed point; n is a whole number of 0 or more. */
export function powFixed(x: bigint, n: number): bigint {
  return powFixedBounded(x, n, null);
}

/**
 * The q-th root of x in fixed point: the largest y with powFixed(y, q) <= x,
 * found by bisection. Because powFixed floors, this is the root as this
 * arithmetic sees it, identical in every language.
 */
export function rootFixed(x: bigint, q: number): bigint {
  if (q === 1) return x;
  let lo = 0n;
  // For x >= 1, (1 + (x - 1)/q)^q >= x (Bernoulli), so the root is below it.
  let hi = x >= FIXED_SCALE ? FIXED_SCALE + (x - FIXED_SCALE) / BigInt(q) + 2n : FIXED_SCALE + 1n;
  while (hi - lo > 1n) {
    const mid = (lo + hi) / 2n;
    if (powFixedBounded(mid, q, x) <= x) lo = mid;
    else hi = mid;
  }
  return lo;
}

/** x^(p/q) in fixed point: the q-th root first, then the power, then the reciprocal if p < 0. */
export function fractionalPowerFixed(x: bigint, p: number, q: number): bigint {
  if (!Number.isInteger(p) || p < -100000 || p > 100000) {
    throw new RangeError(`exponentNumerator must be between -100000 and 100000, received ${p}`);
  }
  if (!Number.isInteger(q) || q < 1 || q > 100000) {
    throw new RangeError(`exponentDenominator must be between 1 and 100000, received ${q}`);
  }
  if (x <= 0n) {
    throw new RangeError("base must be greater than zero");
  }
  const root = rootFixed(x, q);
  if (p >= 0) return powFixed(root, p);
  const denominator = powFixed(root, -p);
  if (denominator === 0n) throw new RangeError("fractional power result too large");
  const result = (FIXED_SCALE * FIXED_SCALE) / denominator;
  if (tooLarge(result)) throw new RangeError("fractional power result too large");
  return result;
}

/** Parse a non-negative decimal of at most 18 places into fixed point. */
export function parseFixed(text: string): bigint {
  if (typeof text !== "string" || !DECIMAL.test(text)) {
    throw new RangeError(`base must be a positive decimal with at most 18 places, received "${text}"`);
  }
  const [whole, fraction = ""] = text.split(".");
  return BigInt(whole) * FIXED_SCALE + BigInt(fraction.padEnd(18, "0"));
}

/**
 * base^(exponentNumerator / exponentDenominator), to 12 decimal places.
 *
 * Computed in 18-place fixed point with a floor after every step, then
 * rounded half-up to 12 places, so the last digit printed is right unless the
 * true value sits within about 10^-15 of a rounding boundary.
 */
export function fractionalPower(base: string, exponentNumerator: number, exponentDenominator: number): string {
  const x = parseFixed(base);
  const value = fractionalPowerFixed(x, exponentNumerator, exponentDenominator);
  const rounded = (value + 500000n) / 1000000n;
  const whole = rounded / 10n ** 12n;
  const fraction = (rounded % 10n ** 12n).toString().padStart(12, "0");
  return `${whole}.${fraction}`;
}