Function extrapolate

src/compute-engine/numerics/richardson.ts:93–161

(
  f: (x: number) => number,
  x0: number,
  options: ExtrapolateOptions = {}
)

 */
export function extrapolate(
  f: (x: number) => number,
  x0: number,
  options: ExtrapolateOptions = {}
): [val: number, err: number] {
  const {
    contract = 0.125,
    step = 1,
    // An ordinary Taylor/asymptotic series in h (see the doc comment above —
    // `power=2` is an opt-in acceleration for even functions). The default
    // was transcribed as 2, contradicting the docstring and Richardson.jl:
    // on series with odd powers (e.g. Hₙ − ln n − γ ~ 1/(2n)) the h¹ term was
    // never eliminated, the error estimate stalled above `limit()`'s 1e-6
    // acceptance threshold, and convergent limits were reported as NaN.
    power = 1,
    atol = 1e-16,
    rtol = atol > 0 ? 0 : Math.sqrt(Number.EPSILON),
    maxeval = 1e6, // Number.MAX_SAFE_INTEGER
    breaktol = 2,
    // A call reached through compiled code (`_SYS.limit`) has no deadline
    // of its own: inherit the ambient one (see interruptible.ts)
    deadline = getAmbientDeadline(),
  } = options;

  if (!isFinite(x0)) {
    // use a change of variables x = 1/u
    return extrapolate((u) => f(1 / u), 1 / x0, {
      rtol,
      atol,
      maxeval,
      contract: Math.abs(contract) > 1 ? 1 / contract : contract,
      step: 1 / step,
      power,
      deadline,
    });
  }

  return withAmbientDeadline(deadline, () => {
    let h = step;
    const invcontract = Math.pow(1 / contract, power);
    let f0 = f(x0 + h);
    const neville: number[] = [f0]; // The current diagonal of the Neville tableau
    let err = Infinity;
    let numeval = 1;

    while (numeval < maxeval) {
      // Each iteration costs a function evaluation, which may itself be
      // expensive: check the evaluation deadline between evaluations.
      checkDeadline(deadline);
      numeval += 1;
      h *= contract;
      neville.push(f(x0 + h));
      let c = invcontract;
      let minerr = Infinity;

      for (let i = neville.length - 2; i >= 0; i--) {
        const old = neville[i];
        neville[i] = neville[i + 1] + (neville[i + 1] - neville[i]) / (c - 1);