mir.numeric

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

Throws:

Exception if FindRootResult.status isn't mir_find_root_status.success.

FindRootStatus status();

Returns:

mir_find_root_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.

Returns:

FindRootResult

Examples:

import mir.math.common: log, exp, fabs;

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

Throws:

Exception if FindRootResult.status isn't mir_find_root_status.success.

FindRootStatus status();

Returns:

mir_find_root_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)

See Also:

findRoot, std.math.isNormal

DiffResult!T diff(alias f, T)(const T x, const T h, const T factor = T(2).sqrt, const T safe = 2);

pure nothrow @nogc @trusted DiffResult!T diffImpl(T)(scope const T delegate(T) pure nothrow @nogc @safe f, const T x, const T h, const T factor = T(2).sqrt, const T safe = 2);

Calculate the derivative of a function. This function uses Ridders' method of extrapolating the results of finite difference formulas for consecutively smaller step sizes, with an improved stopping criterion described in the Numerical Recipes books by Press et al.

This method gives a much higher degree of accuracy in the answer compared to a single finite difference calculation, but requires more function evaluations; typically 6 to 12. The maximum number of function evaluations is 24.

Parameters:

T delegate(T) pure nothrow @nogc @safe f	The function of which to take the derivative.
T x	The point at which to take the derivative.
T h	A "characteristic scale" over which the function changes.
T factor	Stepsize is decreased by factor at each iteration.
T safe	Return when error is SAFE worse than the best so far.

References

• C. J. F. Ridders, Accurate computation of F'(x) and F'(x)F''(x). Advances in Engineering Software, vol. 4 (1982), issue 2, p. 75.
• W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++ (2nd ed.). Cambridge University Press, 2003.

Examples:

import mir.math.common;

auto f(double x) { return exp(x) / (sin(x) - x ^^ 2); }
auto d(double x) { return -(exp(x) * ((x - 2) * x - sin(x) + cos(x)))/(x^^2 - sin(x))^^2; }
auto r = diff!f(1.0, 0.01);
assert (approxEqual(r.value, d(1)));

struct DiffResult(T) if (__traits(isFloating, T));

T value;

T error;

T integrate(alias f, T)(const T a, const T b, const T tolerance = T.epsilon)
if (__traits(isFloating, T));

pure nothrow @nogc @safe T integrateImpl(T)(scope const T delegate(T) pure nothrow @nogc @safe f, const T a, const T b, const T tolerance = T.epsilon)
if (__traits(isFloating, T));

Integrates function on the interval [a, b] using adaptive Gauss-Lobatto algorithm.

References W. Gander and W. Gautschi. Adaptive Quadrature — Revisited

Parameters:

T delegate(T) pure nothrow @nogc @safe f	function to integrate. f should be valid on interval [a, b] including the bounds.
T a	finite left interval bound
T b	finite right interval bound
T tolerance	(optional) relative tolerance should be greater or equal to T.epsilon

Returns:

Integral value

Examples:

import mir.math.common;
import mir.math.constant;

alias cosh = x => 0.5 * (exp(x) + exp(-x));
enum Pi = double(PI);

assert(integrate!exp(0.0, 1.0).approxEqual(double(E - 1)));
assert(integrate!(x => x >= 3)(0.0, 10.0).approxEqual(7.0));
assert(integrate!sqrt(0.0, 1.0).approxEqual(2.0 / 3));
assert(integrate!(x => 23.0 / 25 * cosh(x) - cos(x))(-1.0, 1.0).approxEqual(0.479428226688801667));
assert(integrate!(x => 1 / (x ^^ 4 + x ^^ 2 + 0.9))(-1.0, 1.0).approxEqual(1.5822329637294));
assert(integrate!(x => sqrt(x ^^ 3))(0.0, 1.0).approxEqual(0.4));
assert(integrate!(x => x ? 1 / sqrt(x) : 0)(0.0, 1.0).approxEqual(2));
assert(integrate!(x => 1 / (1 + x ^^ 4))(0.0, 1.0).approxEqual(0.866972987339911));
assert(integrate!(x => 2 / (2 + sin(10 * Pi * x)))(0.0, 1.0).approxEqual(1.1547005383793));
assert(integrate!(x => 1 / (1 + x))(0.0, 1.0).approxEqual(0.6931471805599));
assert(integrate!(x => 1 / (1 + exp(x)))(0.0, 1.0).approxEqual(0.3798854930417));
assert(integrate!(x => exp(x) - 1 ? x / (exp(x) - 1) : 0)(0.0, 1.0).approxEqual(0.777504634112248));
assert(integrate!(x => sin(100 * Pi * x) / (Pi * x))(0.1, 1.0).approxEqual(0.0090986375391668));
assert(integrate!(x => sqrt(50.0) * exp(-50 * Pi * x ^^ 2))(0.0, 10.0).approxEqual(0.5));
assert(integrate!(x => 25 * exp(-25 * x))(0.0, 10.0).approxEqual(1.0));
assert(integrate!(x => 50 / Pi * (2500 * x ^^ 2 + 1))(0.0, 10.0).approxEqual(1.3263071079268e+7));
assert(integrate!(x => 50 * (sin(50 * Pi * x) / (50 * Pi * x)) ^^ 2)(0.01, 1.0).approxEqual(0.11213930374164));
assert(integrate!(x => cos(cos(x) + 3 * sin(x) + 2 * cos(2 * x) + 3 * sin(2 * x) + 3 * cos(3 * x)))(0.0, Pi).approxEqual(0.83867634269443));
assert(integrate!(x => x > 1e-15 ? log(x) : 0)(0.0, 1.0).approxEqual(-1));
assert(integrate!(x => 1 / (x ^^ 2 + 1.005))(-1.0, 1.0).approxEqual(1.5643964440690));
assert(integrate!(x => 1 / cosh(20 * (x - 0.2)) + 1 / cosh(400 * (x - 0.04)) + 1 / cosh(8000 * (x - 0.008)))(0.0, 1.0).approxEqual(0.16349495585710));
assert(integrate!(x => 4 * Pi ^^ 2 * x * sin(20 * Pi * x) * cos(2 * Pi * x))(0.0, 1.0).approxEqual(-0.6346651825434));
assert(integrate!(x => 1 / (1 + (230 * x - 30) ^^ 2))(0.0, 1.0).approxEqual(0.013492485649468));