d3-geo-polygon/src/newton.js

import {abs, epsilon, epsilon2} from "./math.js";

// Approximate Newton-Raphson
// Solve f(x) = y, start from x
export function solve(f, y, x) {
  var steps = 100, delta, f0, f1;
  x = x === undefined ? 0 : +x;
  y = +y;
  do {
    f0 = f(x);
    f1 = f(x + epsilon);
    if (f0 === f1) f1 = f0 + epsilon;
    x -= delta = (-1 * epsilon * (f0 - y)) / (f0 - f1);
  } while (steps-- > 0 && abs(delta) > epsilon);
  return steps < 0 ? NaN : x;
}

// Approximate Newton-Raphson in 2D
// Solve f(a,b) = [x,y]
export function solve2d(f, MAX_ITERATIONS = 40, eps = epsilon2) {
  return function(x, y, a = 0, b = 0) {
    var err2, da, db;
    for (var i = 0; i < MAX_ITERATIONS; i++) {
      var p = f(a, b),
        // diffs
        tx = p[0] - x,
        ty = p[1] - y;
      if (abs(tx) < eps && abs(ty) < eps) break; // we're there!

      // backtrack if we overshot
      var h = tx * tx + ty * ty;
      if (h > err2) {
        a -= da /= 2;
        b -= db /= 2;
        continue;
      }
      err2 = h;

      // partial derivatives
      var ea = (a > 0 ? -1 : 1) * eps,
        eb = (b > 0 ? -1 : 1) * eps,
        pa = f(a + ea, b),
        pb = f(a, b + eb),
        dxa = (pa[0] - p[0]) / ea,
        dya = (pa[1] - p[1]) / ea,
        dxb = (pb[0] - p[0]) / eb,
        dyb = (pb[1] - p[1]) / eb,
        // determinant
        D = dyb * dxa - dya * dxb,
        // newton step — or half-step for small D
        l = (abs(D) < 0.5 ? 0.5 : 1) / D;
      da = (ty * dxb - tx * dyb) * l;
      db = (tx * dya - ty * dxa) * l;
      a += da;
      b += db;
      if (abs(da) < eps && abs(db) < eps) break; // we're crawling
    }
    return [a, b];
  };
}