Report a bug
If you spot a problem with this page, click here to create a Github issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

mir.numeric

Base numeric algorithms.
Reworked part of std.numeric.
Authors:
Ilya Yaroshenko (API, findLocalMin, findRoot extension), Don Clugston (findRoot)
enum mir_find_root_status: int;

alias FindRootStatus = mir_find_root_status;
success
Success
badBounds
nanX
nanY
struct mir_find_root_result(T);

alias FindRootResult = mir_find_root_result(T);
T ax;
Left bound
T bx;
Rifht bound
T ay;
f(ax) or f(ax).fabs.fmin(T.max / 2).copysign(f(ax)).
T by;
f(bx) or f(bx).fabs.fmin(T.max / 2).copysign(f(bx)).
uint iterations;
Amount of target function calls.
inout ref auto validate() return;
Returns:
self Required_versions:D_Exceptions
FindRootStatus status();
T x();
A bound that corresponds to the minimal absolute function value.
Returns:
!(ay.fabs > by.fabs) ? ax : bx
T y();
The minimal of absolute function values.
Returns:
!(ay.fabs > by.fabs) ? ay : by
FindRootResult!T findRoot(alias f, alias tolerance = null, T)(const T ax, const T bx, const T fax = T.nan, const T fbx = T.nan, const T lowerBound = T.nan, const T upperBound = T.nan, uint maxIterations = T.sizeof * 16, uint steps = 0)
if (isFloatingPoint!T && __traits(compiles, T(f(T.init))) && (is(typeof(tolerance) == typeof(null)) || __traits(compiles, () { auto _ = bool(tolerance(T.init, T.init)); } )));
Find root of a real function f(x) by bracketing, allowing the termination condition to be specified.
Given a function f and a range [a .. b] such that f(a) and f(b) have opposite signs or at least one of them equals ±0, returns the value of x in the range which is closest to a root of f(x). If f(x) has more than one root in the range, one will be chosen arbitrarily. If f(x) returns NaN, NaN will be returned; otherwise, this algorithm is guaranteed to succeed.
Uses an algorithm based on TOMS748, which uses inverse cubic interpolation whenever possible, otherwise reverting to parabolic or secant interpolation. Compared to TOMS748, this implementation improves worst-case performance by a factor of more than 100, and typical performance by a factor of 2. For 80-bit reals, most problems require 8 to 15 calls to f(x) to achieve full machine precision. The worst-case performance (pathological cases) is approximately twice the number of bits.

References "On Enclosing Simple Roots of Nonlinear Equations", G. Alefeld, F.A. Potra, Yixun Shi, Mathematics of Computation 61, pp733-744 (1993). Fortran code available from www.netlib.org as algorithm TOMS478.

Parameters:
f Function to be analyzed. f(ax) and f(bx) should have opposite signs, or lowerBound and/or upperBound should be defined to perform initial interval extension.
tolerance Defines an early termination condition. Receives the current upper and lower bounds on the root. The delegate must return true when these bounds are acceptable. If this function always returns false or it is null, full machine precision will be achieved.
T ax Left inner bound of initial range of f known to contain the root.
T bx Right inner bound of initial range of f known to contain the root. Can be equal to ax.
T fax Value of f(ax) (optional).
T fbx Value of f(bx) (optional).
T lowerBound Lower outer bound for interval extension (optional)
T upperBound Upper outer bound for interval extension (optional)
uint maxIterations Appr. maximum allowed number of function calls (inluciding function calls for grid).
uint steps Number of nodes in the left side and right side regular grids (optional). If it equals to 0 (default), then the algorithm uses ieeeMean for searching. The algorithm performs searching if sgn(fax) equals to sgn(fbx) and at least one outer bound allows to extend the inner initial range.
Examples:
import mir.math.common: log, exp;

auto logRoot = findRoot!log(0, double.infinity).validate.x;
assert(logRoot == 1);

auto shift = 1;
auto expm1Root = findRoot!(x => exp(x) - shift)
    (-double.infinity, double.infinity).validate.x;
assert(expm1Root == 0);

auto approxLogRoot = findRoot!(log, (a, b) => fabs(a - b) < 1e-5)(0, double.infinity).validate.x;
assert(fabs(approxLogRoot - 1) < 1e-5);
Examples:
With adaptive bounds
import mir.math.common: log, exp, fabs;

auto logRoot = findRoot!log(
        10, 10, // assume we have one initial point
        double.nan, double.nan, // fa, fb aren't provided by user
        -double.infinity, double.infinity, // all space is available for the bounds extension.
    ).validate.x;
assert(logRoot == 1);

auto shift = 1;
alias expm1Fun = (double x) => exp(x) - shift;
auto expm1RootRet = findRoot!expm1Fun
    (
        11, 10, // reversed order for interval is always OK
        expm1Fun(11), expm1Fun(10), // the order must be the same as bounds
        0, double.infinity, // space for the bounds extension.
    ).validate;
assert(expm1Fun(expm1RootRet.x) == 0);

auto approxLogRoot = findRoot!(log, (a, b) => fabs(a - b) < 1e-5)(
        -1e10, +1e10,
        double.nan, double.nan,
        0, double.infinity,
    ).validate.x;
assert(fabs(approxLogRoot - 1) < 1e-5);
Examples:
ditto
import core.stdc.tgmath: atan;
import mir.math;
import std.meta: AliasSeq;

const double[2][3] boundaries = [
    [0.4, 0.6],
    [1.4, 1.6],
    [0.1, 2.1]];

enum root = 1.0;

foreach(fun; AliasSeq!(
    (double x) => x ^^ 2 - root,
    (double x) => root - x ^^ 2,
    (double x) => atan(x - root),
))
{
    foreach(ref bounds; boundaries)
    {
        auto result = findRoot!fun(
            bounds[0], bounds[1],
            double.nan, double.nan, // f(a) and f(b) not provided
            -double.max, double.max, // user provided outer bounds
        );
        assert(result.validate.x == root);
    }
}

foreach(ref bounds; boundaries)
{
    auto result = findRoot!(x => sin(x - root))(
        bounds[0], bounds[1],
        double.nan, double.nan, // f(a) and f(b) not provided
        -10, 10, // user provided outer bounds
        100, // max iterations,
        10, // steps for grid
    );
    assert(result.validate.x == root);
}
// single initial point, grid, positive direction
auto result = findRoot!((double x ) => sin(x - root))(
    0.1, 0.1, // initial point, a = b
    double.nan, double.nan, // f(a) = f(b) not provided
    -100.0, 100.0, // user provided outer bounds
    150, // max iterations,
    100, // number of steps for grid
);
assert(result.validate.x == root);

// single initial point, grid, negative direction
result = findRoot!((double x ) => sin(x - root))(
    0.1, 0.1, // initial point a = b
    double.nan, double.nan, // f(a) = f(b) not provided
    100.0, -100.0, // user provided outer bounds, reversed order
    150, // max iterations,
    100, // number of steps for grid
);
assert(result.validate.x == double(root - PI)); // Left side root!
Examples:
With adaptive bounds and single initial point. Reverse outer bound order controls first step direction in case of f(a) == f(b).
enum root = 1.0;

   // roots are +/- `root`
   alias fun = (double x) => x * x - root;

double lowerBound = -10.0;
double upperBound = 10.0;

   assert(
       findRoot!fun(
           0, 0, // initial interval
           double.nan, double.nan,
           lowerBound, upperBound,
           // positive direction has higher priority
       ).validate.x == root
   );

   assert(
       findRoot!fun(
           0, 0, // initial interval
           double.nan, double.nan,
           upperBound, lowerBound,
           // reversed order
       ).validate.x == -root // other root
   );
export pure nothrow @nogc @safe FindRootResult!float findRootImpl(float ax, float bx, float fax, float fbx, float lowerBound, float upperBound, uint maxIterations, uint steps, scope float delegate(float) pure nothrow @nogc @safe f, scope bool delegate(float, float) pure nothrow @nogc @safe tolerance);

export pure nothrow @nogc @safe FindRootResult!double findRootImpl(double ax, double bx, double fax, double fbx, double lowerBound, double upperBound, uint maxIterations, uint steps, scope double delegate(double) pure nothrow @nogc @safe f, scope bool delegate(double, double) pure nothrow @nogc @safe tolerance);

export pure nothrow @nogc @safe FindRootResult!real findRootImpl(real ax, real bx, real fax, real fbx, real lowerBound, real upperBound, uint maxIterations, uint steps, scope real delegate(real) pure nothrow @nogc @safe f, scope bool delegate(real, real) pure nothrow @nogc @safe tolerance);
findRoot implementations.
struct FindLocalMinResult(T);
T x;
T y;
T error;
ref auto validate() return;
Returns:
self Required_versions:D_Exceptions
FindRootStatus status();
FindLocalMinResult!T findLocalMin(alias f, T)(const T ax, const T bx, const T relTolerance = sqrt(T.epsilon), const T absTolerance = sqrt(T.epsilon))
if (isFloatingPoint!T && __traits(compiles, T(f(T.init))));
Find a real minimum of a real function f(x) via bracketing. Given a function f and a range (ax .. bx), returns the value of x in the range which is closest to a minimum of f(x). f is never evaluted at the endpoints of ax and bx. If f(x) has more than one minimum in the range, one will be chosen arbitrarily. If f(x) returns NaN or -Infinity, (x, f(x), NaN) will be returned; otherwise, this algorithm is guaranteed to succeed.
Parameters:
f Function to be analyzed
T ax Left bound of initial range of f known to contain the minimum.
T bx Right bound of initial range of f known to contain the minimum.
T relTolerance Relative tolerance.
T absTolerance Absolute tolerance.

Preconditions ax and bx shall be finite reals.
relTolerance shall be normal positive real.
absTolerance shall be normal positive real no less then T.epsilon*2.

Returns:
A FindLocalMinResult consisting of x, y = f(x) and error = 3 * (absTolerance * fabs(x) + relTolerance). The method used is a combination of golden section search and successive parabolic interpolation. Convergence is never much slower than that for a Fibonacci search.

References "Algorithms for Minimization without Derivatives", Richard Brent, Prentice-Hall, Inc. (1973)