|
8 | 8 | namespace stan { |
9 | 9 | namespace math { |
10 | 10 |
|
| 11 | +#include <cmath> |
| 12 | +#include <limits> |
| 13 | +#include <algorithm> |
| 14 | + |
| 15 | +namespace internal { |
| 16 | +template <std::size_t StencilOrder = 2, typename T> |
| 17 | +inline constexpr auto eps_root_calc() { |
| 18 | + constexpr T eps = std::numeric_limits<T>::epsilon(); |
| 19 | + if constexpr (StencilOrder == 2) { |
| 20 | + // ε^(1/3) |
| 21 | + return std::cbrt(eps); |
| 22 | + } else if constexpr (StencilOrder == 3) { |
| 23 | + // ε^(1/4) = sqrt(sqrt(ε)) |
| 24 | + return std::sqrt(std::sqrt(eps)); |
| 25 | + } else if constexpr (StencilOrder == 5) { |
| 26 | + // ε^(1/6) = sqrt(cbrt(ε)) |
| 27 | + return std::sqrt(std::cbrt(eps)); |
| 28 | + } else { |
| 29 | + // General fallback: ε^(1/(p+1)) |
| 30 | + return std::pow(eps, T(1) / T(StencilOrder + 1)); |
| 31 | + } |
| 32 | +} |
| 33 | +} |
11 | 34 | /** |
12 | | - * Return the stepsize for finite difference evaluations at the |
13 | | - * specified scalar. |
| 35 | + * @brief Compute a finite-difference step size suitable for a stencil with |
| 36 | + * leading truncation order \p StencilOrder. |
| 37 | + * |
| 38 | + * This implements the standard balance between truncation error |
| 39 | + * \f$A\,h^{p}\f$ and floating-point rounding error \f$B\,\varepsilon/h\f$, |
| 40 | + * whose minimizer is \f$h_\star \propto \varepsilon^{1/(p+1)}\f$, where |
| 41 | + * \f$p = \texttt{StencilOrder}\f$ and \f$\varepsilon\f$ is machine epsilon |
| 42 | + * for the scalar type \p T. We additionally scale by \f$\max(1, |u|)\f$ to |
| 43 | + * obtain an absolute step size near the point being perturbed. |
| 44 | + * |
| 45 | + * For the common second-order (3-point central) stencil, the function uses an |
| 46 | + * `if constexpr` branch to compute \f$\varepsilon^{1/3}\f$ via `std::cbrt` |
| 47 | + * (avoids a `pow` and is numerically tidy). For higher orders it uses |
| 48 | + * \f$\varepsilon^{1/(p+1)}\f` via `std::pow`. |
| 49 | + * |
| 50 | + * @tparam T Floating-point scalar type (e.g., float, double). |
| 51 | + * @tparam StencilOrder Leading truncation order \f$p\f$ of your stencil: |
| 52 | + * - first-derivative 3-point central → \f$p=2\f$ |
| 53 | + * - first-derivative 5-point central → \f$p=4\f$ |
| 54 | + * - first-derivative 6-point central → \f$p=6\f$ |
| 55 | + * - “derivative of Hessian” via 5-point diff → \f$p=4\f$ |
| 56 | + * |
| 57 | + * @param u The coordinate value (or local scale) at which the step will be |
| 58 | + * applied; the step is scaled by \f$\max(1, |u|)\f$. |
| 59 | + * @param c Dimensionless tuning constant multiplying the theoretically |
| 60 | + * optimal step. Default is 1.0. In practice, \f$c \in [0.5, 2]\f$ |
| 61 | + * works well: |
| 62 | + * - Increase \p c (larger h) if round-off/cancellation dominates |
| 63 | + * (noisy function values, very large |f|, many subtractions). |
| 64 | + * - Decrease \p c (smaller h) if truncation dominates (very smooth |
| 65 | + * function with large higher derivatives or low curvature scale). |
| 66 | + * If your problem has a known physical scale S (not |u|), prefer |
| 67 | + * passing \p u scaled by S or modify the scaling accordingly. |
14 | 68 | * |
15 | | - * <p>The formula used is `stepsize(u) = cbrt(epsilon) * max(1, |
16 | | - * abs(u)).` |
| 69 | + * @return A step size \f$h\f$ such that \f$u+h \neq u\f$; if the theoretical |
| 70 | + * step underflows at \p u, the function falls back to the next |
| 71 | + * representable increment. |
17 | 72 | * |
18 | | - * @param u initial value to increment |
19 | | - * @return stepsize away from u for finite differences |
| 73 | + * @note The step computed here assumes smooth (at least \(p{+}1\) times |
| 74 | + * differentiable) behavior along the perturbed coordinate. |
20 | 75 | */ |
21 | | -inline double finite_diff_stepsize(double u) { |
22 | | - using std::fabs; |
23 | | - static const double eps = std::numeric_limits<double>::epsilon(); |
24 | | - static const double eps_pow |
25 | | - = std::pow(eps, 1.0 / 7.0); // for 6th-order stencil |
26 | | - return eps_pow * std::fmax(1.0, fabs(u)); |
| 76 | +template <std::size_t StencilOrder = 2, typename T> |
| 77 | +inline T finite_diff_stepsize(T u, T c = T(1)) { |
| 78 | + const T scale = std::max(T(1), std::abs(u)); |
| 79 | + const T eps_root = internal::eps_root_calc<StencilOrder, T>(); |
| 80 | + const T h = c * eps_root * scale; |
| 81 | + // Ensure perturbation isn’t rounded away at u. |
| 82 | + if (u + h == u) { |
| 83 | + const T next = std::nextafter(u, std::numeric_limits<T>::infinity()); |
| 84 | + return std::max(h, next - u); |
| 85 | + } |
| 86 | + return h; |
27 | 87 | } |
28 | 88 |
|
29 | 89 | } // namespace math |
|
0 commit comments