Passing callable in lambda capture

Question

While reading SICP's section about "Formulating Abstractions with Higher-OrderFunctions" I was curious to implement a simple, root solving procedure using Newtons' Method in C++.

However my implementation is faulty but I cannot figure out why. The code compiles but I got a runtime crash due to stack overflow because the loop recursion inside fixedPoint function cannot reach MAX_IT. I run it in Debug mode so tail recursion does not apply.

The problem is that the function g inside fixedPoint is never evaluated correctly (I always get inf. result). I think that I miss something about callable passing and lambda captures (I try to pass std::function in lambda). Below I have the sample code which runs a square root problem.

#include <iostream>
#include <functional>
#include <cmath>

static constexpr double DELTA = 0.0000001;
static constexpr int MAX_IT = 50000;
using my_function = std::function<double(double)>;

double derivative(const my_function& f, const double x) {
    return (f(x + DELTA) - f(x)) / DELTA;
}

bool isConverge(const double p1, const double p2) {
    return std::abs(p1 - p2) < DELTA;
}

double fixedPoint(const my_function& g, double x = 1) {

    int itCnt = 0;
    std::function<double(double, double)> loop;
    loop = [&g, &itCnt, &loop](double x0, double x1) -> double {
        if (isConverge(x0, x1) || (++itCnt > MAX_IT)) return x1;
        return loop(x1, g(x1));
    };

    return loop(x, g(x));
}

my_function newtonsMethod(const my_function& f) {
    return [&f](double x) -> double {return x - f(x) / derivative(f, x); };
}

double findRoot(const my_function& f) {
    return fixedPoint(newtonsMethod(f));
}

double sqRoot(double x) { 
    return findRoot([x](int y) {return y * y - x; });
}

int main() {
    const double num = 133;
    std::cout << "Square root " << num << ": " << sqRoot(num) << '\n';
    return 0;
}

Answer 1

Your issue is int truncation from double with

return findRoot([x](int y) {return y * y - x; });

So when you compute derivative, you obtain 0.

Change to

return findRoot([x](double y) {return y * y - x; });

Demo