Tools for testing identitites

Description

Using known identities is a powerful way to test precision of our functions. Currently using identities for derivatives is a bit tedious - you need to explicitly compute derivative against each variable on both left and right hand side of the identity and compare them. The proposal is to create a function that does this for you.

Example

The goal is to be able to write something like:

auto lh = [](const var& a, const var& b) {
  return stan::math::lbeta(a, b);
}
// Alternatively: auto lh = static_cast<var(*)(const var&, const var&)>(lbeta);
auto rh = [](const var& a, const var& b) {
  return stan::math::log_sum_exp(lbeta(a + 1, b), lbeta(a, b + 1));
};

for(double x : x_to_test) {
  for(double y : y_to_test) {
    expect_identity("succesors", lh, rh, tolerance, x, y); 
  }
}

I would not expect expect_identity to check this for all possible instantiations, I believe there should be a separate test that all instantiations are numerically equal and after that we should just test with all var to make things easier.

Also note while I could test that lh - rh == 0 (with zero gradients all over) this is potentially problematic as the test tolerance could not reflect the magnitude of the actual values.

Expected output

It tests that the value and all derivates are close up to a relative tolerance using expect_near_rel

Current Version:

v3.1.0

Start with the proposed expect_identity interface, the var-based derivative example, and the existing expect_near_rel behavior described in the issue. Define completion as checking values and all derivatives within relative tolerance, while separately covering numerical equality across instantiations.