Generating a C++ function from a list of argument

I am writing a function f to be used in a Runge Kutta integrator.

output RungeKutta(function f, initial conditions IC, etc.)

Since the function will be called many times, I am looking for a way to generate the function f at compile time.

In this case, function f depends on a fixed list of parameters vector p, where p is sparse and is fixed before the code is compiled. To be concrete,

double function f(vector<double> x) {
    return x dot p;
}

Since p is sparse, taking the dot product in f is not the most efficient. Hard-coding x dot p seems to be the way to go, but p can be very long (1000).

What are my options?

Is writing another program (taking p as input) to generate a .cpp file my only option?

Thanks for the comments. Here is a more concrete example for the differential equation.

dy/dx = f_p(x)

One example for f_p(x):

p = [0, 1, 0]; x = [x1, x2, x3]

double f_p(vector<double> x) {
     return x2;  // This is what I meant by hard-coding
}

instead of:

   double f(vector<double> p, vector<double> x) {
       double r = 0;
       for (i=0; i < p.length(); i++) {
              r += p[i]*x[i];
        }
        return r;
   }

Solution

The key problem you are trying to solve is that a "leaf" function in your calculation that will be called many times will also most often do no work given the problem domain. The hope is that the redundant work - namely multiplying a value with an element of an array known at compile time to be zero - can be collapsed as part of a compile time step.

C++ has language facilities to deal with this, namely template metaprogramming. C++ templates are very powerful (ie Turing complete) and allow for things like recursive calculations based on compile time constants.

Below is an example of how to implement your example using templates and template specialization (you can also find a runnable example I've created here http://ideone.com/BDtBt7). The basic idea behind the code is to generate a type with a static function that returns the resulting dot product of an input vector of values and a compile time constant array. The static function recursively calls instances of itself, passing a lower index value as it moves through the input/constant arrays of elements. It is also templated with whether the value in the compile time constant array p that is being evaluated is zero. If it is, we can skip calculating that and move onto the next value in the recursion. Lastly, there is a base case that stops the recursion once we have reached the first element in the array.

#include <array>
#include <iostream>
#include <vector>

constexpr std::array<double, 5> p = { 1.0, 0.0, 3.0, 5.0, 0.0 };

template<size_t index, bool isZero>
struct DotProductCalculator
{
    static double Calculate(const std::vector<double>& xArg)
    {
        return (xArg[index] * p[index]) 
            + DotProductCalculator<index - 1, p[index - 1] == 0.0>::Calculate(xArg);
    }
};

template<>
struct DotProductCalculator<0, true>
{
    static double Calculate(const std::vector<double>& xArg)
    {
        return 0.0;
    }
};

template<>
struct DotProductCalculator<0, false>
{
    static double Calculate(const std::vector<double>& xArg)
    {
        return xArg[0] * p[0];
    }
};

template<size_t index>
struct DotProductCalculator<index, true>
{
    static double Calculate(const std::vector<double>& xArg)
    {
        return 0.0 + DotProductCalculator<index - 1, p[index - 1] == 0.0>::Calculate(xArg);
    }
};

template<typename ArrayType>
double f_p_driver(const std::vector<double>& xArg, const ArrayType& pAsArgument) 
{
     return DotProductCalculator<std::tuple_size<ArrayType>::value - 1, 
                            p[std::tuple_size<ArrayType>::value -1] == 0.0>::Calculate(xArg); 
}

int main() 
{
    std::vector<double> x = { 1.0, 2.0, 3.0, 4.0, 5.0 };
    double result = f_p_driver(x, p);
    std::cout << "Result: " << result;
    return 0;
}