//====================================================================
// Driver program: Calculates first- and second-order derivatives
// Method: centered differences with Richardson extrapolation
//
// The derivatives are first calculated with step sizes h and h/2.
// Richardson extrapolation then removes the leading O(h^2) error.
//
// Alexander Godunov, September 2026
//====================================================================
#include <cmath>
#include <iomanip>
#include <iostream>

using ScalarFunction = double (*)(double);

void richardsonDerivative(ScalarFunction func, double x, double h, int m,
                          double& derivH, double& derivH2, double& derivR,
                          double& errorEst, int& status);

double centeredDifference(ScalarFunction func, double x, double h, int m);
double testFunction(double x);


int main()
{
    double x, h;
    double derivH, derivH2, derivR;
    double errorEst, exact, actualError;
    int status;

    //----------------------------------------------------------------
    // Step 1: specify the evaluation point and initial step size
    //
    // The example uses
    //
    //       f(x) = ln(x),
    //
    // for which f'(x) = 1/x and f''(x) = -1/x^2.
    //----------------------------------------------------------------
    x = 2.0;
    h = 0.2;

    //----------------------------------------------------------------
    // Step 2: calculate the first- and second-order derivatives
    //----------------------------------------------------------------
    std::cout << "Centered differences with Richardson extrapolation\n";
    std::cout << std::fixed << std::setprecision(3)
              << "f(x) = ln(x),  x = " << x << ",  h = " << h << "\n\n";

    std::cout << " order       exact            D(h)          D(h/2)"
              << "       Richardson     estimated error"
              << "        actual error\n";

    for (int m = 1; m <= 2; ++m) {
        richardsonDerivative(testFunction, x, h, m, derivH, derivH2,
                             derivR, errorEst, status);

        if (status != 0) {
            std::cout << "Derivative of order " << m
                      << " could not be evaluated; status = "
                      << status << '\n';
            continue;
        }

        if (m == 1) {
            exact = 1.0/x;
        } else {
            exact = -1.0/(x*x);
        }

        actualError = std::abs(derivR-exact);

        std::cout << std::scientific << std::setprecision(8)
                  << std::setw(6) << m
                  << std::setw(17) << exact
                  << std::setw(17) << derivH
                  << std::setw(17) << derivH2
                  << std::setw(17) << derivR
                  << std::setw(17) << errorEst
                  << std::setw(17) << actualError << '\n';
    }

    return 0;
}


void richardsonDerivative(ScalarFunction func, double x, double h, int m,
                          double& derivH, double& derivH2, double& derivR,
                          double& errorEst, int& status)
//====================================================================
// Calculate a first- or second-order derivative using a centered
// difference and one level of Richardson extrapolation.
//
// Input ...
// func       - function to be differentiated
// x          - point where the derivative is required
// h          - initial finite-difference step size
// m          - order of the derivative (1 or 2)
//
// Output ...
// derivH     - centered difference calculated with h
// derivH2    - centered difference calculated with h/2
// derivR     - Richardson-extrapolated derivative
// errorEst   - magnitude of the Richardson correction
// status     - 0: success
//              1: m is not 1 or 2
//              2: h is not positive
//
// A centered difference has an error proportional to h^2. For step
// refinement by a factor of 2, Richardson extrapolation gives
//
//       D_R = D(h/2) + [D(h/2)-D(h)]/(2^2-1).
//====================================================================
{
    derivH = 0.0;
    derivH2 = 0.0;
    derivR = 0.0;
    errorEst = 0.0;
    status = 0;

    // Check the order of the requested derivative
    if (m < 1 || m > 2) {
        status = 1;
        return;
    }

    // A positive finite-difference step is required
    if (h <= 0.0) {
        status = 2;
        return;
    }

    // Calculate the centered differences using h and h/2
    derivH = centeredDifference(func, x, h, m);
    derivH2 = centeredDifference(func, x, h/2.0, m);

    // Remove the leading O(h^2) error
    const double correction = (derivH2-derivH)/3.0;
    derivR = derivH2 + correction;

    // The Richardson correction provides an error estimate
    errorEst = std::abs(correction);
}


double centeredDifference(ScalarFunction func, double x, double h, int m)
//====================================================================
// Calculate a first- or second-order derivative using a second-order
// centered-difference formula.
//====================================================================
{
    if (m == 1) {
        return (func(x+h)-func(x-h))/(2.0*h);
    }

    return (func(x+h)-2.0*func(x)+func(x-h))/(h*h);
}


double testFunction(double x)
//====================================================================
// Function used in the driver program. Replace this expression with
// the function required in a practical calculation.
//====================================================================
{
    return std::log(x);
}