#include <cmath>
#include <cstdio>
#include <limits>

static float f(float x)
{
    return x * x - 2.f;
}

static float bisect(float (*f)(float), float a, float b, float xtol, float ftol)
{
    float fa = f(a), fb = f(b);
    for (unsigned nmax = 100; nmax; --nmax) { // avoid getting stuck
        const float mid = 0.5f * (a + b);
        const float fmid = f(mid);
        //std::printf("a %e fa %e b %e fb %e mid %e fmid %e\n", a, fa, b, fb, mid, fmid);
        // look for interval in which f changes sign
        if (fa * fmid <= 0.f) b = mid, fb = fmid;
        else if (fb * fmid <= 0.f) a = mid, fa = fmid;
        if (std::abs(b - a) <= xtol || std::abs(fmid) <= ftol) return mid;
    }
    return std::numeric_limits<float>::quiet_NaN();
}

int main()
{
    std::printf("sqrt(2) from bisection: %e\n", bisect(f, 1.f, 2.f, 2.f * std::numeric_limits<float>::epsilon(), 2.f * std::numeric_limits<float>::epsilon()));
    return 0;
}