iii. Root Finding and Integration

Solving $f(x)=0$ and computing ${\int }_a^{b}f$ numerically. For contest purposes, binary search with a fixed iteration count is often the right root finder (slower per step, unconditionally correct); Newton is the tool when speed or precision doubling matters, and it reappears everywhere as a doubling scheme.

Method	Convergence	Needs $f^{\prime }$	Note
bisection	linear, 1 bit per step	no	always works given a sign change
Newton	quadratic near a simple root	yes	correct digits double each step
secant	superlinear, order $\phi \approx 1.618$	no	finite difference instead of $f^{\prime }$
regula falsi	linear to superlinear	no
Brent / safeguarded Newton	superlinear, guaranteed	no / yes	hybrid with a bracket
fixed-point iteration $x\leftarrow g(x)$	linear if $|g^{\prime }|<1$	no

Newton’s method

$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}.$

Failure	Symptom	Fix
$f^{\prime }(x_n)\approx 0$	huge jump, divergence	bisection fallback
bad starting point	another root, or divergence	bracket first
multiple root ($f^{\prime }=0$ at the root)	only linear convergence	$x-m\frac{f}{f^{\prime }}$ for multiplicity $m$
oscillation	cycles between two points	damping or a hybrid

Optimization via Newton: applying it to $f^{\prime }$, $x\leftarrow x-\frac{f^{\prime }(x)}{f^{\prime \prime }(x)}$, finds a stationary point quadratically, but any critical point, not necessarily a minimum; ternary search is safer for unimodal functions in contests.

Newton beyond real numbers. The same iteration works wherever there is a “derivative”, and each step doubles the precision:

Setting	Iteration	Doubles
reals	$x-f/f^{\prime }$	correct digits
power series inverse	$B\leftarrow B(2-AB)$	correct coefficients: Power Series
power series  ​, $\log$, $\exp$	$B\leftarrow \frac{1}{2}(B+A/B)$, $B\leftarrow B(1-\log B+A)$	coefficients
$p$-adic (Hensel lifting)	$r-f(r)f^{\prime }(r{)}^{-1}$	the power of $p$: Modular Square Roots
matrix inverse (Newton–Schulz)	$X(2I-AX)$	correct digits
big-integer reciprocal	$x(2-ax)$	bits: how big-number division is implemented

Because the last doubling step dominates, the total cost is $T(n)=T(n/2)+O(M(n))=O(M(n))$, the cost of one multiplication: this turns “invert a power series in $O(n\log n)$” into a solved problem.

template<class F, class D> double newton(F f, D df, double x, int iters = 60) {
    for (int i = 0; i < iters; i++) {
        double d = df(x); if (d == 0) break;
        double nx = x - f(x) / d;
        if (fabs(nx - x) < 1e-15 * max(1.0, fabs(x))) return nx;
        x = nx;
    }
    return x;
}
 
template<class F, class D> double safeNewton(F f, D df, double lo, double hi) {   // f(lo), f(hi) of opposite signs
    double x = (lo + hi) / 2; bool negLo = f(lo) < 0;
    for (int i = 0; i < 200; i++) {
        double fx = f(x); if (fx == 0) return x;
        if ((fx < 0) == negLo) lo = x; else hi = x;
        double d = df(x), nx = d != 0 ? x - fx / d : (lo + hi) / 2;
        x = (nx > lo && nx < hi) ? nx : (lo + hi) / 2;     // fall back to bisection
    }
    return x;
}
 
ll isqrtFloor(ll n) {                                      // floor(sqrt(n)), n >= 0: correct the double guess
    ll x = (ll)sqrtl((long double)n);
    while (x > 0 && (i128)x * x > n) x--;
    while ((i128)(x + 1) * (x + 1) <= n) x++;
    return x;
}
 
ull isqrtNewton(ull n) {                                   // floor(sqrt(n)) with integer Newton steps
    if (n < 2) return n;
    ull x = n, y = x / 2 + (x & 1);                        // (x + n / x) / 2 for x = n, without overflow
    while (y < x) x = y, y = (x + n / x) / 2;
    return x;
}
 
i128 ipowCap(ll b, int k, i128 cap) { i128 r = 1; while (k--) { r *= b; if (r > cap) return cap + 1; } return r; }
ll iroot(ll n, int k) {                                    // floor(n^(1/k)), n >= 0, k >= 1
    ll x = (ll)powl((long double)n, 1.0L / k);
    while (x > 0 && ipowCap(x, k, n) > n) x--;
    while (ipowCap(x + 1, k, n) <= n) x++;
    return x;
}

Always correct the result of sqrt / pow: on a double, sqrt(n) for $n$ near $10^{18}$ is often off by one, a classic wrong answer; the adjustment loops cost nothing.

Proof: quadratic convergence, the safeguard, and the integer roots

Quadratic. Let $r$ be a simple root and $e_k=x_k-r$. Taylor: $0=f(r)=f(x_k)-e_k{f}^{\prime }(x_k)+\frac{1}{2}{e}_k^2{f}^{\prime \prime }(\xi )$, so $e_{k+1}=e_k-\frac{f(x_k)}{f^{\prime }(x_k)}=\frac{f^{\prime \prime }(\xi )}{2f^{\prime }(x_k)}{e}_k^2$; near $r$ the factor is bounded, so $|e_{k+1}|\le C{e}_k^2$. At a root of multiplicity $m$, $f/f^{\prime }$ behaves like $e/m$, so the plain step only removes a factor $1-\frac{1}{m}$ of the error, and $x-mf/f^{\prime }$ restores quadratic convergence.

Safeguarded Newton. The bracket $[lo,hi]$ always contains a sign change and never grows; every iteration either takes a Newton step inside the bracket or bisects, and the bracket is updated with the sign at the current point, so it shrinks to the root.

Integer square root. isqrtFloor ends with $x^2\le n<(x+1{)}^2$ by its two loops (computed in __int128), whatever the initial guess. In isqrtNewton, for $x>\lfloor \sqrt{n}\rfloor$ the next value $y=\lfloor (x+\lfloor n/x\rfloor )/2\rfloor$ satisfies $y\ge \lfloor \sqrt{n}\rfloor$ (AM–GM on $x$ and $n/x$; floors only lose fractional parts) and $y<x$ (since $n/x<x$); so the sequence strictly decreases while above $\lfloor \sqrt{n}\rfloor$ and never drops below it, and at $x=\lfloor \sqrt{n}\rfloor$ we have $n/x\ge x$, so $y\ge x$ and the loop stops. The first step from $x=n$ computes $\lfloor (n+1)/2\rfloor$ as $\lfloor n/2\rfloor +(n\mathrm{mod}2)$ to avoid overflow. iroot is the same correction with powers capped at $n+1$ to avoid overflow. All checked against exact __int128 arithmetic near perfect powers.

Polynomial roots

All distinct real roots of a polynomial, by recursing on the derivative and bisecting monotone pieces. polyRoots(p) searches $|x|\le 1+{\max }_i|a_i/a_n|$ (Cauchy’s bound); polyRoots(p, lo, hi) only $[lo,hi]$. Cost $O(100n^3)$ for degree $n$: $n$ derivative levels, up to $n$ pieces each, 100 bisection steps of $O(n)$.

double evalPoly(const vector<double>& p, double x) { double r = 0; for (int i = sz(p) - 1; i >= 0; i--) r = r * x + p[i]; return r; }
 
vector<double> polyRoots(vector<double> p, double lo, double hi) {   // p[i] = coefficient of x^i
    while (sz(p) > 1 && p.back() == 0) p.pop_back();
    int d = sz(p) - 1;
    if (d <= 0) return {};
    if (d == 1) { double r = -p[0] / p[1]; if (lo <= r && r <= hi) return {r}; return {}; }
    vector<double> der(d);
    for (int i = 1; i <= d; i++) der[i - 1] = p[i] * i;
    vector<double> pts = polyRoots(der, lo, hi);           // critical points split [lo, hi] into monotone pieces
    pts.insert(pts.begin(), lo); pts.push_back(hi);
    auto nearZero = [&](double x) {                        // |p(x)| within Horner's rounding error bound
        double scale = 0, pw = 1;
        for (double c : p) scale += fabs(c) * pw, pw *= max(1.0, fabs(x));
        return fabs(evalPoly(p, x)) <= 4 * d * DBL_EPSILON * scale;
    };
    vector<double> roots;
    for (int i = 0; i < sz(pts); i++) {
        if (nearZero(pts[i])) { roots.push_back(pts[i]); continue; }  // even-multiplicity roots touch zero here
        if (i + 1 == sz(pts) || nearZero(pts[i + 1])) continue;
        double l = pts[i], h = pts[i + 1]; bool negLeft = evalPoly(p, l) < 0;
        if (negLeft == (evalPoly(p, h) < 0)) continue;
        for (int it = 0; it < 100; it++) { double mid = (l + h) / 2; ((evalPoly(p, mid) < 0) == negLeft ? l : h) = mid; }
        roots.push_back((l + h) / 2);
    }
    sort(all(roots));
    vector<double> res;
    for (double r : roots) if (res.empty() || r - res.back() > 1e-9 * max(1.0, fabs(r))) res.push_back(r);
    return res;
}
 
vector<double> polyRoots(vector<double> p) {
    while (sz(p) > 1 && p.back() == 0) p.pop_back();
    double bound = 1;
    for (int i = 0; i + 1 < sz(p); i++) bound = max(bound, 1 + fabs(p[i] / p.back()));
    return polyRoots(p, -bound, bound);
}

Roots of even multiplicity do not change sign: a critical point with $|p(x)|\le 4n\epsilon \sum |a_i|\max (1,|x|{)}^i$ (Horner’s error bound) is reported as a root, and roots closer than $10^{-9}$ (relative) are merged. So a multiple root is reported once, as accurately as a simple root of the derivative; two distinct roots closer than about $10^{-7}$ may merge, and a local extremum within the error bound of 0 is reported as a root. Accuracy depends on conditioning: in tests (degree $\le 8$, roots in $[-10,10]$ at least $0.05$ apart) the error stayed below $10^{-7}$. Use long double or rescale $x$ for clustered roots, high degree or huge coefficients, and strip leading coefficients that are tiny only because of rounding.

Proof: why the derivative recursion finds every root

Between two consecutive real roots of $p^{\prime }$ (or the interval ends), $p^{\prime }$ has no root, so $p$ is strictly monotone there and has at most one root, which exists iff $p$ changes sign between the endpoints; bisection then converges to it (the invariant “sign at $l$ = sign at the left end, opposite sign at $h$” is kept). A root not inside any such piece is a critical point, i.e. a multiple root, which the nearZero test catches. Cauchy’s bound: if $|x|>1+{\max }_i|a_i/a_n|$ then $|a_n{x}^n|>|{\sum }_{i<n}{a}_i{x}^i|$, so no root lies outside. Tested on random polynomials built from known roots, including double roots.

Integration

Method	Error	Exact for	Note
rectangle (midpoint)	$O(h^2)$	linear	simplest
trapezoid	$O(h^2)$	linear
Simpson	$O(h^4)$	cubics	the default
Simpson 3/8	$O(h^4)$	cubics
Romberg	$O(h^{2k})$	high degree	trapezoid + repeated Richardson extrapolation
Gauss–Legendre ($n$ points)	exact to degree $2n-1$		best per evaluation, needs tabulated nodes
Monte Carlo	$O(1/\sqrt{N})$		dimension-free: the only option in high dimensions

$\text{Simpson: }{\int }_a^{b}f\approx \frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b)).$

template<class F> double simpson(F f, double a, double b, int n) {   // composite, n even
    double h = (b - a) / n, s = f(a) + f(b);
    for (int i = 1; i < n; i++) s += f(a + i * h) * (i % 2 ? 4 : 2);
    return s * h / 3;
}
 
template<class F> double adaptiveSimpson(F f, double a, double b, double eps, double whole, double fa, double fm, double fb, int depth) {
    double c = (a + b) / 2, l = (a + c) / 2, r = (c + b) / 2, fl = f(l), fr = f(r);
    double left = (c - a) / 6 * (fa + 4 * fl + fm), right = (b - c) / 6 * (fm + 4 * fr + fb);
    if (depth <= 0 || fabs(left + right - whole) <= 15 * eps) return left + right + (left + right - whole) / 15;   // Richardson
    return adaptiveSimpson(f, a, c, eps / 2, left, fa, fl, fm, depth - 1) + adaptiveSimpson(f, c, b, eps / 2, right, fm, fr, fb, depth - 1);
}
template<class F> double integrate(F f, double a, double b, double eps = 1e-10) {
    double fa = f(a), fm = f((a + b) / 2), fb = f(b);
    return adaptiveSimpson(f, a, b, eps, (b - a) / 6 * (fa + 4 * fm + fb), fa, fm, fb, 50);
}

Adaptive Simpson spends effort where the function is difficult and almost none where it is flat; reusing the already computed values halves the evaluations, and the /15 term is Richardson extrapolation, which upgrades the accuracy for free.

Integration in contests	Integrand
area under a curve, between curves	$f$, the difference
arc length	$\sqrt{1+f^{\prime }(x{)}^2}$
volume of revolution	$\pi f(x{)}^2$
volume by cross-sections (the most common)	area of the 2D slice at $x$, often itself a geometric computation
expected value, probability of a region	$xp(x)$, the density
union of circles / shapes	slice length or area as a function of $x$

Simpson assumes smoothness

Adaptive Simpson silently gives wrong answers on functions with discontinuities or kinks: the error estimate compares two smooth approximations that happen to agree. Split the interval at every known breakpoint (where shapes start or stop overlapping, where the maximum changes which curve is on top) and integrate each smooth piece. Also: remove endpoint singularities by substitution, and give oscillatory integrands a finer base subdivision (aliasing).

Exact alternatives, check first: polygon area by the shoelace formula; circle–polygon intersection by closed-form sectors and triangles; piecewise-linear functions as sums of trapezoids; polynomial integrands symbolically. Reach for numerical integration only when there is no closed form. Monte Carlo: ${\int }_{\Omega }f\approx \frac{\mathrm{vol}(\Omega )}{N}\sum f(x_i)$ with error $O(1/\sqrt{N})$ independent of the dimension (4 times the samples per extra digit); variance reduction (importance sampling, stratification, quasi-random sequences) helps a lot.

Proof: Simpson is exact for cubics, and the Richardson correction

On $[-h,h]$, Simpson gives $\frac{h}{3}(f(-h)+4f(0)+f(h))$. For $f=1,x,x^2,x^3$ this is $2h,0,\frac{2h^3}{3},0$, equal to the exact integrals; by linearity it is exact for all cubics. For $x^4$ it gives $\frac{2h^5}{3}$ instead of $\frac{2h^5}{5}$, an error of $\frac{4h^5}{15}$; with Taylor’s theorem the error on one panel of width $2h$ is $-\frac{h^5}{90}{f}^{(4)}(\xi )$, so $n$ panels give $O(h^4)$ overall. Halving the panel width divides the error by about $2^4=16$: if $S_1$ is the coarse and $S_2$ the two-half estimate, $I-S_2\approx \frac{S_2-S_1}{15}$, which is both the stopping test and the correction added to the result (this makes the answer exact for quintics on each panel). Checked: composite Simpson is exact on random cubics, and integrate matches closed forms ($\sin$, $e^x$, $\sqrt{x}$, $\frac{1}{1+x^2}$) to $10^{-9}$.

Newton–Schulz: matrix inverse by multiplications

$X_{k+1}=X_k(2I-A{X}_k),X_0=\frac{A^{\top }}{\Vert A{\Vert }_1\Vert A{\Vert }_{\infty }}.$

The scalar analogue is Newton on $f(x)=\frac{1}{x}-a$: $x_{k+1}=x_k(2-a{x}_k)$, convergent from any start with $|1-a{x}_0|<1$. Per iteration: two matrix multiplications, $O(n^3)$ (or $O(n^{\omega })$); iterations: $O(\log \kappa )$ to enter the convergence basin plus $O(\log \log \frac{1}{\epsilon })$. Sequentially worse than Gauss–Jordan; it wins where multiplication is much cheaper than elimination: GPUs and distributed systems (embarrassingly parallel GEMMs versus sequential pivoting), settings without division (fixed-point hardware, homomorphic encryption, some rings), and refinement (one step squares the error of an approximate inverse from any source). The same iteration converges to the Moore–Penrose pseudo-inverse $A^{+}$ for rectangular $A$, and relatives compute the matrix sign function and the polar decomposition.

vector<vector<double>> newtonSchulzInverse(const vector<vector<double>>& A, int iters = 100) {
    int n = sz(A); double n1 = 0, ninf = 0;
    for (int j = 0; j < n; j++) { double s = 0; for (int i = 0; i < n; i++) s += fabs(A[i][j]); n1 = max(n1, s); }
    for (int i = 0; i < n; i++) { double s = 0; for (int j = 0; j < n; j++) s += fabs(A[i][j]); ninf = max(ninf, s); }
    auto mul = [&](const vector<vector<double>>& P, const vector<vector<double>>& Q) {
        vector<vector<double>> R(n, vector<double>(n, 0));
        for (int i = 0; i < n; i++) for (int k = 0; k < n; k++) for (int j = 0; j < n; j++) R[i][j] += P[i][k] * Q[k][j];
        return R;
    };
    vector<vector<double>> X(n, vector<double>(n));
    for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) X[i][j] = A[j][i] / (n1 * ninf);
    for (int it = 0; it < iters; it++) {
        auto AX = mul(A, X); double err = 0;
        for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) { AX[i][j] = (i == j ? 2 : 0) - AX[i][j]; err += fabs(AX[i][j] - (i == j)); }
        X = mul(X, AX);
        if (err < 1e-13 * n) break;                         // AX was already the identity
    }
    return X;
}

Proof: the residual squares, and the starting point converges

With $E_k=I-A{X}_k$: $E_{k+1}=I-A{X}_k(2I-A{X}_k)=(I-A{X}_k{)}^2=E_k^2$, so $\Vert E_{k+1}\Vert \le \Vert E_k{\Vert }^2$ and correct digits double once $\Vert E_k\Vert <1$. For the start, $A{X}_0=\frac{A{A}^{\top }}{t}$ with $t=\Vert A{\Vert }_1\Vert A{\Vert }_{\infty }$. $A{A}^{\top }$ is symmetric positive definite (for invertible $A$) with largest eigenvalue $\Vert A{\Vert }_2^2\le \Vert A^{\top }A{\Vert }_1\le \Vert A^{\top }{\Vert }_1\Vert A{\Vert }_1=\Vert A{\Vert }_{\infty }\Vert A{\Vert }_1=t$, so the eigenvalues of $E_0=I-\frac{A{A}^{\top }}{t}$ lie in $[0,1)$ and its spectral norm is below 1. The number of steps to bring $\Vert E\Vert =1-\frac{{\lambda }_{\min }}{t}$ down to $\frac{1}{2}$ is about ${\log }_2\frac{t}{{\lambda }_{\min }}=O(\log \kappa )$. Checked against Gaussian elimination on random well-conditioned matrices up to $8\times 8$.

Related notes

• Floating Point
• Power Series (Newton on series)
• Binary Search