Numerical Methods

// Template: find smallest x in [lo, hi] where check(x) is true
function binarySearchReal(lo, hi, check, iterations = 100) {
    for (let i = 0; i < iterations; i++) {
        const mid = (lo + hi) / 2;
        if (check(mid)) hi = mid;
        else lo = mid;
    }
    return (lo + hi) / 2;
}

// Example: minimum radius so k circles cover all n points on a line
function minRadius(points, k) {
    points.sort((a, b) => a - b);
    const check = (r) => {
        let covered = 0, circles = 0;
        while (covered < points.length) {
            circles++;
            const start = points[covered];
            // circle covers [start, start + 2r]
            while (covered < points.length && points[covered] <= start + 2*r)
                covered++;
        }
        return circles <= k;
    };
    return binarySearchReal(0, points[points.length-1]-points[0], check);
}

// Find x in [lo, hi] minimizing f(x), where f is unimodal (convex)
function ternarySearch(lo, hi, f, iterations = 200) {
    for (let i = 0; i < iterations; i++) {
        const m1 = lo + (hi - lo) / 3;
        const m2 = hi - (hi - lo) / 3;
        if (f(m1) < f(m2)) hi = m2;
        else lo = m1;
    }
    return (lo + hi) / 2;
}

// Example: find point on line segment closest to origin
// f(t) = distance from point (t*ax + (1-t)*bx, t*ay + (1-t)*by) to (0,0)
function closestPointOnSegment(ax, ay, bx, by) {
    const f = t => {
        const x = ax + t*(bx-ax), y = ay + t*(by-ay);
        return x*x + y*y;
    };
    const t = ternarySearch(0, 1, f);
    return [ax + t*(bx-ax), ay + t*(by-ay)];
}

// Integer ternary search
function ternarySearchInt(lo, hi, f) {
    while (hi - lo > 2) {
        const m1 = Math.floor((lo*2 + hi) / 3);
        const m2 = Math.floor((lo + hi*2) / 3);
        if (f(m1) <= f(m2)) hi = m2;
        else lo = m1;
    }
    let best = lo;
    for (let x = lo+1; x <= hi; x++) if (f(x) < f(best)) best = x;
    return best;
}

function simulatedAnnealing(initialState, energy, neighbor, maxIter = 100000) {
    let state = initialState;
    let bestState = state;
    let bestEnergy = energy(state);
    let T = 1.0; // initial temperature
    const cooling = 0.9999;

    for (let i = 0; i < maxIter; i++) {
        T *= cooling;
        const next = neighbor(state);
        const delta = energy(next) - energy(state);
        // Accept if better, or with probability e^(-delta/T) if worse
        if (delta < 0 || Math.random() < Math.exp(-delta / T)) {
            state = next;
            if (energy(state) < bestEnergy) {
                bestState = state;
                bestEnergy = energy(state);
            }
        }
    }
    return bestState;
}

// Example: minimize sum of distances from point to set of points (geometric median)
// neighbor: small random perturbation of (x, y)
// energy: sum of euclidean distances

// Fast square root via Newton's method
function sqrtNewton(n) {
    if (n === 0) return 0;
    let x = n;
    while (true) {
        const next = (x + n/x) / 2;
        if (Math.abs(next - x) < 1e-9) return next;
        x = next;
    }
}

// General Newton iteration: x_{n+1} = x_n - f(x_n)/f'(x_n)
// Converges quadratically (doubles correct digits each step)
function newtonMethod(f, fPrime, x0, tolerance = 1e-9) {
    let x = x0;
    for (let i = 0; i < 100; i++) {
        const dx = f(x) / fPrime(x);
        x -= dx;
        if (Math.abs(dx) < tolerance) break;
    }
    return x;
}