9

long time browser, first time asker here. I've written a number of scripts for doing various 1D numerical integration methods and compiled them into a library. I would like that library to be as flexible as possible regarding what it is capable of integrating.

Here I include an example: a very simple trapezoidal rule example where I pass a pointer to the function to be integrated.

// Numerically integrate (*f) from a to b
// using the trapezoidal rule.
double trap(double (*f)(double), double a, double b) {
  int N = 10000;
  double step = (b-a)/N;
  double s = 0;
  for (int i=0; i<=N; i++) {
    double xi = a + i*step;
    if (i == 0 || i == N) { s += (*f)(xi); }
    else { s += 2*(*f)(xi); }
  }
  s *= (b-a)/(2*N);
  return s;
}

This works great for simple functions that only take one argument. Example:

double a = trap(sin,0,1);

However, sometimes I may want to integrate something that has more parameters, like a quadratic polynomial. In this example, the coefficients would be defined by the user before the integration. Example code:

// arbitrary quadratic polynomial
double quad(double A, double B, double C, double x) {
  return (A*pow(x,2) + B*x + C);
}

Ideally, I would be able to do something like this to integrate it:

double b = trap(quad(1,2,3),0,1);

But clearly that doesn't work. I have gotten around this problem by defining a class that has the coefficients as members and the function of interest as a member function:

class Model {
  double A,B,C;
public:
  Model() { A = 0; B = 0; C = 0; }
  Model(double x, double y, double z) { A = x; B = y; C = z; }
  double func(double x) { return (A*pow(x,2)+B*x+C); }
};

However, then my integration function needs to change to take an object as input instead of a function pointer:

// Numerically integrate model.func from a to b
// using the trapezoidal rule.
double trap(Model poly, double a, double b) {
  int N = 10000;
  double step = (b-a)/N;
  double s = 0;
  for (int i=0; i<=N; i++) {
    double xi = a + i*step;
    if (i == 0 || i == N) { s += poly.func(xi); }
    else { s += 2*poly.func(xi); }
  }
  s *= (b-a)/(2*N);
  return s;
}

This works fine, but the resulting library is not very independent, since it needs the class Model to be defined somewhere. Also, ideally the Model should be able to change from user-to-user so I wouldn't want to fix it in a header file. I have tried to use function templates and functors to get this to work but it is not very independent since again, the template should be defined in a header file (unless you want to explicitly instantiate, which I don't).

So, to sum up: is there any way I can get my integration functions to accept arbitrary 1D functions with a variable number of input parameters while still remaining independent enough that they can be compiled into a stand-alone library? Thanks in advance for the suggestions.

t354
  • 567
  • 4
  • 12

2 Answers2

8

What you need is templates and std::bind() (or its boost::bind() counterpart if you can't afford C++11). For instance, this is what your trap() function would become:

template<typename F>
double trap(F&& f, double a, double b) {
  int N = 10000;
  double step = (b-a)/N;
  double s = 0;
  for (int i=0; i<=N; i++) {
    double xi = a + i*step;
    if (i == 0 || i == N) { s += f(xi); }
//                               ^
    else { s += 2* f(xi); }
//                 ^
  }
  s *= (b-a)/(2*N);
  return s;
}

Notice, that we are generalizing from function pointers and allow any type of callable objects (including a C++11 lambda, for instance) to be passed in. Therefore, the syntax for invoking the user-provided function is not *f(param) (which only works for function pointers), but just f(param).

Concerning the flexibility, let's consider two hardcoded functions (and pretend them to be meaningful):

double foo(double x)
{
    return x * 2;
}

double bar(double x, double y, double z, double t)
{
    return x + y * (z - t);
}

You can now provide both the first function directly in input to trap(), or the result of binding the last three arguments of the second function to some particular value (you have free choice on which arguments to bind):

#include <functional>

int main()
{
    trap(foo, 0, 42);
    trap(std::bind(bar, std::placeholders::_1, 42, 1729, 0), 0, 42);
}

Of course, you can get even more flexibility with lambdas:

#include <functional>
#include <iostream>

int main()
{
    trap(foo, 0, 42);
    trap(std::bind(bar, std::placeholders::_1, 42, 1729, 0), 0, 42);

    int x = 1729; // Or the result of some computation...
    int y = 42; // Or some particular state information...
    trap([&] (double d) -> double
    {
        x += 42 * d; // Or some meaningful computation...
        y = 1; // Or some meaningful operation...
        return x;
    }, 0, 42);

    std::cout << y; // Prints 1
}

And you can also pass your own stateful functors tp trap(), or some callable objects wrapped in an std::function object (or boost::function if you can't afford C++11). The choice is pretty wide.

Here is a live example.

Andy Prowl
  • 124,023
  • 23
  • 387
  • 451
  • Thanks! I do like the template example because it's very elegant, but I'm not sure it will work for me - my understanding is that a template function has to be defined and declared in the same place, while I would like to define `trap()` in a separate library and then declare it for use in other projects. I took your suggestion of using `std::function` in conjunction with `std::bind` and redefined `trap()` as `double trap(std::function f, double a, double b)`; the only annoying part is that I have to redefine simple functions like sin as `std::function`. – t354 Apr 12 '13 at 14:19
  • 2
    @t354: Yes, templates suffer from that separation problem: the definition must be in header files. Alternatively, you *could* put the definitions in a `.cpp` file and provide so-called *explicit instantiations* of your template for of all the possible arguments used in your application, but if yours is a library, then it's impossible to predict how your template will be instantiated. Putting `std::function` in the signature is OK: be aware that there is a run-time overhead, so you'll have to measure whether it is a relevant overhead for you or not. – Andy Prowl Apr 12 '13 at 14:27
  • 1
    @t354: I do not understand the last part though: what do you mean by "redefine simple functions like `sin`"? If you mean this: `std::function = sin; trap(sin, 0, 42);`, than it's not needed. You can invoke directly `trap(sin, 0, 42);`, without declaring the `std::function` object. – Andy Prowl Apr 12 '13 at 14:28
  • 1
    If I have defined `double trap(std::function f, double a, double b)` and I try to invoke `trap(sin,0,42);` as you suggested, I get an error: `error: conversion from '' to non-scalar type 'std::function' requested`. I can avoid this by re-defining sin: `std::function f_sin = (double(*)(double))&std::sin;` Regarding the run-time overhead mentioned in your first comment, would using lambdas be more efficient than my current solution? I am not very familiar with them but it might be worth the effort. – t354 Apr 12 '13 at 14:48
  • 2
    @t354: Oh, that's because [`std::sin`](http://en.cppreference.com/w/cpp/numeric/math/sin) is actually an overloaded function (there is a `float` version, a `double` version, etc., so the compiler can't tell what overload you want it to pick when you write just `sin`). So you should provide an explicit cast: `trap(static_cast(sin), 0, 42);` – Andy Prowl Apr 12 '13 at 14:50
4

What you trying to do is to make this possible

trap( quad, 1, 2, 3, 0, 1 );

With C++11 we have alias template and variadic template

template< typename... Ts >
using custom_function_t = double (*f) ( double, Ts... );

above define a custom_function_t that take a double and variable numbers of arguments.

so your trap function becomes

template< typename... Ts >
double trap( custom_function_t<Ts...> f, Ts... args, double a, double b ) {
    int N = 10000;
    double step = (b-a)/N;
    double s = 0;
    for (int i=0; i<=N; i++) {
        double xi = a + i*step;
        if (i == 0 || i == N) { s += f(xi, args...); }
        else { s += 2*f(xi, args...); }
    }
    s *= (b-a)/(2*N);
    return s;
}

Usage:

double foo ( double X ) {
  return X;
}

double quad( double X, double A, double B, double C ) {
  return(A*pow(x,2) + B*x + C);
}

int main() {
  double result_foo  = trap( foo, 0, 1 );
  double result_quad = trap( quad, 1, 2, 3, 0, 1 );  // 1, 2, 3 == A, B, C respectively
}

Tested on Apple LLVM 4.2 compiler.

yngccc
  • 5,594
  • 2
  • 23
  • 33