import mir.interpolate.spline;
import mir.math.common: sqrt;
import mir.ndslice.allocation: rcslice;
import mir.test;

alias g = (double x, uint d = 0) =>
    d == 0 ? 3 * x ^^ 3 + 5 * x ^^ 2 + 0.23 * x + 2 :
    d == 1 ? 9 * x ^^ 2 + 10 * x + 0.23 :
    double.nan;

auto s = spline!double(
    [0.1, 0.4, 0.5, 1.0].rcslice!(immutable double),
    [g(0.1), g(0.4), g(0.5), g(1.0)].rcslice!(const double)
);

auto m = s.interpolationSqrt;

m(0.7).shouldApprox == sqrt(g(0.7));
auto d = m.opCall!3(0.7);
d[0].shouldApprox == sqrt(g(0.7));
d[1].shouldApprox == 0.5 / sqrt(g(0.7)) * g(0.7, 1);
d[2].shouldApprox == 2.2292836438189414;
d[3].shouldApprox == -3.11161;