import std.stdio, std.math, std.algorithm;

bool nearZero(T)(in T a, in T b = T.epsilon * 4) pure nothrow {
    return abs(a) <= b;
}

T[] findroot(T)(T function(T) pure nothrow f, in T start,
                in T end, T step=cast(T)0.001L,
                T tolerance = cast(T)1e-4L) {
    T[T] result;

    if (nearZero(step))
        writefln("WARNING: step size may be too small.");

    T searchRoot(T a, T b) nothrow { // search root by simple bisection
        T root;
        int limit = 49;
        T gap = b - a;

        while (!nearZero(gap) && limit--) {
            if (nearZero(f(a))) return a;
            if (nearZero(f(b))) return b;
            root = (b + a) / 2.0L;
            if (nearZero(f(root))) return root;
            if (f(a) * f(root) < 0)
                b = root;
            else
                a = root;
            gap = b - a;
        }

        return root;
    }

    immutable T dir = cast(T)(end > start ? 1.0 : -1.0);
    step = (end > start) ? abs(step) : -abs(step);
    for (T x = start; (x * dir) <= (end * dir); x += step)
        if (f(x) * f(x + step) <= 0) {
            T r = searchRoot(x, x + step);
            result[r] = f(r);
        }

    return result.keys.sort().release(); // reduce duplacated root, if any
}

void report(T)(in T[] r, in T function(T) pure f, in T tolerance = cast(T)1e-4L) {
    if (r.length) {
        writefln("Root found (tolerance = %1.4g) :", tolerance);
        foreach (x; r) {
            immutable T y = f(x);
            if (nearZero(y))
                writefln("... EXACTLY at %+1.20f, f(x) = %+1.4g", x, y);
            else if (nearZero(y, tolerance))
                writefln(".... MAY-BE at %+1.20f, f(x) = %+1.4g", x, y);
            else
                writefln("Verify needed, f(%1.4g) = %1.4g > tolerance in magnitude", x, y);
        }
    } else
        writefln("No root found.");
}

real f(real x) pure nothrow {
    return x ^^ 3 - (3 * x ^^ 2) + 2 * x;
}

void main(){
    findroot(&f, -1.0L, 3.0L, 0.001L).report(&f);
}