C++ Finite difference differentiation - design

for the sake of simplicity consider the function:

double foo(const std::vector<double> &a, const std::vector<double> &b, ...) {
 /* do complex stuff */
 return ...;
}

(In reality the types of a and b are more complex objects).

we want to differentiate foo() with respect to its arguments. Therefore the first order sensitivity d foo/d a is a std::vector<double> with size equal to a.size(). Same reasoning goes for d foo/d b.

A naive implementation would go as follows:

std::vector<double> a = {1, 2, 3, 4, 5};
std::vector<double> b = {1, 2, 3};

// compute d foo/d a 
std::vector<double> computeDfDa(std::vector<double> a, std::vector<double> b, ..., double da = 1.0){
  std::vector<double> dfda = {};
  for (auto i = 0; i < a.size(); ++i) {

      // bump up
      a[i] += da;
      auto up = foo(a, b);
      a[i] -= da;

      // bump down
      a[i] -= da;
      auto down = foo(a, b);
      a[i] += da;

      auto derivative = (up - down) / 2.0 / da;
      dfda.pushback(derivative);
  }
   return dfda;
}

// compute d foo/d b 
std::vector<double> computeDfDb(std::vector<double> a, std::vector<double> b, ..., double db = 1.0){
  std::vector<double> dfdb = {};
  for (auto i = 0; i < b.size(); ++i) {
      // bump up
      b[i] += db;
      auto up = foo(a, b);
      b[i] -= db;

      // bump down
      b[i] -= db;
      auto down = foo(a, b);
      b[i] += db;

      auto derivative = (up - down) / 2.0 / db;
      dfdb.pushback(derivative);
  }
  return dfdb;
}

This works well however we have basically the same code for computeDfDa() and for computeDfDb().

Is there any design pattern that would allow to have a unique (maybe templated) function that would understand automatically which input to bump?

If the complexity and the number of inputs of foo() are much greater the naive solution would generate a lot of useless code as we'd have to write a computeDfDx() function for every input x of foo().