```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 ``` ```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); } ```