Passing a C Function Callback to the Boost Math Toolkit Root Finding Library

bisect.cxx

// compile with g++ -o bisect.cxx poly.o -I<path/to/boost>

#include <cmath>
#include <iostream>

#include <boost/math/tools/roots.hpp>

// Functor providing termination condition
// Determines convergence in x based on if the relative or absolute tolerance are satisfied
template<class T>
struct eps_relabs {
  
  T rtol{1.0e-6};  // relative tolerance
  T atol{1.0e-12}; // absolute tolerance

  bool operator()(const T& a, const T&b) const {

    T d = std::abs(b - a);
    if (d <= atol) {
      // absolue
      return true;
    } else {
      // relative
      if (a != (T) 0) {
        return (d / std::abs(a)) <= rtol;
      } else {
        return false;
      }
    }
  }
};

// === C callback prototypes ===

typedef float (*root_callback_s)(
    const float* /* x */,
    void* /* data */);

// https://stackoverflow.com/questions/18410374/using-a-template-callback-function-in-c
template<typename T>
using root_callback = T (*)(
    const T* /* x */, 
    void* /* data */
  );

// We can either inline a function signature
//
//   T (*f)(const T*, void*)
//
// or use the template alias
//
//   root_callback<T> f
//

template<typename T, typename it_t>
int bisect_gen(root_callback<T> f, T a, T b, it_t max_iter, void* data, T& x) {

  using namespace boost::math::tools;

  // Our custom termination condition
  eps_relabs<T> tol;

  auto it = static_cast<std::uintmax_t>(max_iter);

  auto [x1,x2] = bisect(
    [f,data](const T& z){ return f(&z, data); },
    a, b, tol, it);

  if (it >= static_cast<std::uintmax_t>(max_iter)) {
    // exceeded max_iter function invocations
    return -1;
  }

  // found root, tol(x1,x2) == true
  x = x1 + (x2 - x1)/2; // Midway between brackets is our result, if necessary we could
                        // return the result as an interval here.
  return 0;

}

extern "C" double poly(const double* x, void* data);

int main(int argc, char const *argv[])
{
  int max_iter = 30;
  double a = 1, b = 2;

  // check callback works
  std::cout << "poly(a) = " << poly(&a,nullptr) << '\n';
  std::cout << "poly(b) = " << poly(&b,nullptr) << '\n';

  double x;
  int res = bisect_gen(poly,a,b,max_iter, nullptr, x);

  std::cout << "x       = " << x << '\n';
  std::cout << "res     = " << res << '\n';
  std::cout << "poly(x) = " << poly(&x,nullptr) << '\n';

  return 0;
}