C++: newton raphson and overloading functions

Question

I wrote a simple implementation of the newton raphson root finding algorithm which takes an initial guess init, a unary function f and the tolerance tol as arguments, as shown below:

bool newton_raphson(double& init,
                    double(*f)(double),
                    double tol){
    const int max_iter = 10000;
    double next_x, soln = init;
    int i = 0;

    while(++i < max_iter){
        next_x = soln - f(soln)/fp_x(f, soln);
        if(fabs(next_x - soln) < tol){
            init = next_x;
            return true;
        }
        soln = next_x;
    }
    return false;
}

double fp_x(double(*f)(double),
            double x){
    const double h = 0.000001;
    return (f(x + h) - f(x - h))/2.0/h;
}

My question is: although this works perfectly fine for unary functions, I would like to change the implementation so that it works for functions f that have more than one parameter, but all except one parameter have constant values. To clarify: if I have a function f(x) = 3x + 2 as shown below

double f(double x){
    return (3*x + 2);
}

Then my implementation works. However, I would also like it to work for any functions with any given number of arguments, but only the first argument is variable. So, if I have a function f(x,y) = 3x + 2y

double f(double x, double y){
    return (3*x + 2*y);
}

I would like to find the root of f(x,2), or f(x,3) using the same function, and so on for n arguments, not just one or two (please ignore the idea that the functions I showed in the example are simple linear functions, this is just an example). Is there any way to implement the function for a varying number of arguments or do I have to write an implementation for every case?

Thanks,

NAX

NOTE

As you could tell by now, this question isn't really about newton-raphson, but it makes it easier if I use it as an example for the actual question, which is a single implementation for functions of different numbers of arguments.

UPDATE

A few answers below use std::bind and std::function to solve the problem, which actually better address my question than the selected answer; however, they are c++11 library classes/functions, (which, don't get me wrong, is something I strongly urge every c++ programmer to go ahead and learn) and at the time of this writing, I was facing some problems using them; Eclipse Juno using g++ 4.7 (which is c++11 compliant) still somehow failed to recognize std::function, and I therefore decided to go and stick with the checked answer below, which also works nicely.

Answer 1

I think you're asking for variadic functions:

A variadic function – a function declared with a parameter list ending with ellipsis (...) – can accept a varying number of arguments of differing types. Variadic functions are flexible, but they are also hazardous. The compiler can't verify that a given call to a variadic function passes an appropriate number of arguments or that those arguments have appropriate types. Consequently, a runtime call to a variadic function that passes inappropriate arguments yields undefined behavior. Such undefined behavior could be exploited to run arbitrary code. From here:

https://www.securecoding.cert.org/confluence/display/cplusplus/DCL31-CPP.+Do+not+define+variadic+functions

However, as quoted above, there are a number of problems with them.

Most notably, it only works for compile time!

However, if you are interested in implementing one, here's an article with a nice example:

http://www.informit.com/guides/content.aspx?g=cplusplus&seqNum=138

UPDATE:

IMO, I think you're better off defining functions that take structure or object arguments (i.e. a general function object), and writing functions that work on those arguments explicitly.

The other option is to do some compile-time reflection - which would be useful, but is too much trouble to do in such an example as this. Plus, "reflection" in C++ isn't "true" reflection, but rather a bad and incomplete implementation of it.

Answer 2

For what you're trying to do here, what you're looking for is std::bind (or, if you're dealing with a C++03 compiler, std::bind1st and std::bnd2nd).

These will let you "bind" values to the other parameters, leaving you with a function (technically, a function object) that only requires a single parameter.

What you'd ideally like would be something like this:

double f(double x, double y) { 
     return 3*x + 2*y;
}

double init = 1.0;

newton_raphson(init, std::bind2nd(f, 3), 1e-4);

Unfortunately, in real use, it's not quite that simple -- to work with std::bind2nd, you can't use an actual function; you need to use a function object instead, and it has to derive from std::binary_function.

std::bind is quite a bit more flexible, so that's what you almost certainly want to use instead (if at all possible).

Answer 3

I used your question as a way to force myself to learn C++11 variadic template, here is a working example.

template< typename... Ts >
double f( Ts... Vs ) {
    double array[] = { Vs... };
    int numArg = sizeof...( Vs );
    switch (numArg) {
    case 1:
        return 3 * array[0] + 2;
    case 2:
        return 3 * array[0] + 2 * array[1];
    case 3:
        return 3 * array[0] + 2 * array[1] + 1 * array[3];
    ....
    default:
        return 0.0;
    }
}

template< typename... Ts >
double newton_raphson( double &init, double tol,
                       double (*func) ( Ts... Vs ), Ts... Vs ) {
    return func( Vs... );
}

you can call it like

newton_raphson( &init, 1.0, f, 1.0, 2.0, 3.0, 4.0, 5.0 );

Answer 4

You can use std::bind and std::function. The type std::function<double(double)> represents a functional that takes in a double and returns a double. Similarly std::function<double(int,int)> is for a functional taking 2 ints and returns a double.

#include <functional>

bool newton_raphson(double& init,
                    std::function<double(double)>& f,
                    double tol){
    const int max_iter = 10000;
    double next_x, soln = init;
    int i = 0;

    while(++i < max_iter){
        next_x = soln - f(soln)/fp_x(f, soln);
        if(fabs(next_x - soln) < tol){
            init = next_x;
            return true;
        }
        soln = next_x;
    }
    return false;
}

double myfunction(double x, double y){
    return (3*x + 2*y);
}
double fp_x(std::function<double(double)> f, double x) {
    ...
}
...
double d = 1.0;
// Here we set y=2.5 and we tell bind that 1st parameter is unbounded
// If we wanted to switch and set x=2.5 and let y be unbounded, then
// we would use  (&myfunction, 2.5, std::placeholders::_1)
newton_raphson(d, std::bind(&myfunction, std::placeholders::_1, 2.5) , 1e-6);
...