numerics

A single-header C++ numerical methods library. Implements bracketed root-finding, numerical differentiation with Richardson extrapolation, and polynomial divided differences. Also includes Func<Scalar> and Poly<Scalar> — composable wrappers for building mathematical expressions that work naturally with all solvers. Inspired by Brent's method. Written from scratch.

Usage

Copy functional.hpp and numerics.hpp into your project. No dependencies beyond the standard library.

#include "functional.hpp" // Func<Scalar>, Poly<Scalar>, math factories
#include "numerics.hpp"   // zero_in, differentiation, divided_difference, richardson_extrapolation

functional.hpp

Func<Scalar>

Wraps any callable as a first-class mathematical object. Supports arithmetic between functions and scalars, and function composition.

auto f = Sin<double>() + Cos<double>();          // pointwise sum
auto g = f * 2.0 - Func<double>(1.0);           // scalar arithmetic
auto h = Sin<double>().compose(Cos<double>());   // sin(cos(x))
double val = f(M_PI / 4.0);                     // evaluate anywhere

Built-in factories: Sin, Cos, Tan, Exp, Log, Sqrt, Identity

Poly<Scalar>

Dense polynomial stored as a coefficient vector where coeffs[i] is the coefficient of x^i.

Poly<double> p({-6.0, 11.0, -6.0, 1.0}); // x^3 - 6x^2 + 11x - 6
double val = p(2.0);                       // evaluate via Horner's method

// Deflate out a known root via synthetic division
double rem;
Poly<double> q = p.deflate(1.0, &rem);    // divide out (x - 1)

Supports +, -, * between polynomials and scalars. Implicitly converts to Func<Scalar> so it works anywhere a Func is expected.

numerics.hpp

zero_in

template<typename F, typename Scalar = double>
Scalar zero_in(F func, Scalar a, Scalar b, Scalar tol = 1e-10);

Finds a root of func in the bracket [a, b]. The function must have opposite signs at the endpoints. Uses inverse quadratic interpolation (IQI) for fast superlinear convergence where possible, falling back to bisection when IQI steps are rejected (out of bounds, degenerate, or making insufficient progress).

// sqrt(2): find root of x^2 - 2 in [1, 2]
double r = zero_in([](double x){ return x*x - 2.0; }, 1.0, 2.0);
// r ≈ 1.41421356237

// pi: find root of sin(x) in [3, 3.5]
double pi = zero_in(Sin<double>(), 3.0, 3.5);

// Works directly with Poly<double> via implicit conversion
Poly<double> p({-2.0, 0.0, 1.0}); // x^2 - 2
double r2 = zero_in(static_cast<Func<double>>(p), 1.0, 2.0);

// Find all roots iteratively via deflation
double r1 = zero_in(static_cast<Func<double>>(p), 0.5, 1.5);
Poly<double> q = p.deflate(r1);
double r2 = zero_in(static_cast<Func<double>>(q), 1.5, 2.5);

Numerical Differentiation

Three finite difference formulas of increasing accuracy. All accept any callable.

// Forward difference: O(h) accuracy. h > 0 forward, h < 0 backward.
double d1 = forward_difference(func, x0, h);
double d2 = backward_difference(func, x0, h);

// Central difference: O(h^2) accuracy. More accurate than forward at the same h.
double d3 = central_difference(func, x0, h);

richardson_extrapolation

template<typename F, typename Scalar = double>
Scalar richardson_extrapolation(F func, Scalar x0, Scalar h, int order,
                                Scalar tol = 1e-10, int max = 10);

Repeatedly halves h and applies the Richardson recurrence to cancel successive error terms, achieving accuracy well beyond the base formula. order is the leading error order of the base formula (2 for central difference). Stops when successive diagonal estimates agree within tol, or after max levels.

// d/dx sin(x) at x=1, much more accurate than central_difference alone
double d = richardson_extrapolation(
    [](double x){ return std::sin(x); }, 1.0, 0.5, 2
);
// d ≈ cos(1) to near machine precision

divided_difference

template<typename Scalar = double>
std::vector<std::vector<Scalar>> divided_difference(
    const std::vector<Scalar>& points,
    const std::vector<Scalar>& vals,
    bool converge = true,
    Scalar tol = 1e-10);

Builds the Newton divided difference table for a set of interpolation points and values. The diagonal entries table[i][i] are the Newton forward difference coefficients. With converge = true, returns early when successive diagonal entries agree within tol, and the returned table is trimmed to only the rows computed.

std::vector<double> pts  = {0.0, 1.0, 2.0};
std::vector<double> vals = {0.0, 1.0, 4.0};  // f(x) = x^2
auto table = divided_difference(pts, vals);
// table[2][2] == 1.0  (second-order divided difference of x^2)

Building the Tests

Tests use GoogleTest and CMake.

mkdir build && cd build
cmake ..
make
./run_tests

Or via CTest:

ctest --output-on-failure

Requirements

• C++17
• GoogleTest (for tests.cpp only — not required to use the headers)

License

MIT