Taylor expansion in C++: function factory?

In some application I need to have access to the elements of some polynomial expansion (for example x^i for i = 0,1,...). My idea was to create a "function factory" like

#include <functional>
#include <iostream>

std::function<double(double)> Taylor(unsigned int i)
{
    return [i](double y) {return std::pow(y, i);};
}

int main()
{   
    unsigned int n = 3; // the wanted degree in your polynomial basis
    auto f = Taylor(n);
    double a = f(3);
    std::cout << "a = " << a << std::endl;
    return 0;
}

In the example above, we would get the output a = 27.

This would be useful in the context of function estimation (i.e., roughly, find the coefficients of some polynomial expansion given data).

Is there some better way to do this? Also, if the polynomial basis is not just x^i but something more complicated, using a lambda function could be inappropriate...

Solution 1:

C++, long before lambdas were added to it, has always had a solution for this: Just implement a class with an operator(), so that instances become callables, i.e. useful as functions. We call such objects functors.

#include <iostream>
#include <memory>
#include <vector>
// Just an interface base class – not strictly necessary
class realmap {
public:
  virtual double operator()(double parameter) = 0;
  virtual ~operator();
};

class identity : public realmap {
  double operator()(double parameter) override { return parameter; }
};
class twice : public realmap {
  double operator()(double parameter) override { return 2 * parameter; }
};

class monomial : public realmap {
private:
  const unsigned int exponent;

public:
  monomial(unsigned int exponent) : exponent(exponent) {}
  double operator()(double parameter) override {
    double result = 1;
    for (unsigned int pot = 0; pot < exponent; ++pot) {
      result *= parameter;
    }
    return result;
  }
};

int main() {
  auto cube = monomial(3);
  std::cout << cube(1.3) << std::endl;
  std::vector<std::unique_ptr<realmap>> operations;
  operations.emplace_back(new twice());
  operations.emplace_back(new monomial(2));
  operations.emplace_back(new twice());
  operations.emplace_back(new monomial(3));
  double value = 3.141;
  for (auto &operation : operations) {
    value = (*operation)(value);
  }
  // now is [(3.141 * 2)²·2]³
  std::cout << value << std::endl;
}