Numerical Integration Calculator (cpp, written by Codex)

Codex wrote this cpp project from the task below. It does not run on a clean Ubuntu 22.04 machine as written.

Task given to the agent:

TASK: Numerical Integration Calculator

Write a program that computes definite integrals of mathematical functions using multiple numerical methods, comparing accuracy and convergence across methods.

FUNCTIONAL REQUIREMENTS:
- Accept a mathematical expression as a command-line argument via --function flag (e.g., --function "sin(x)*exp(-x)")
- Accept integration bounds via --lower and --upper flags
- Support multiple numerical integration methods selectable via --method flag: trapezoidal rule, Simpson's rule, Simpson's 3/8 rule, Gaussian quadrature, and Romberg integration
- Support a configurable number of subintervals via --intervals flag (default: 1000) for methods that use subdivision
- Run all methods and compare results via --compare flag, showing each method's result, estimated error, and computation time
- Support adaptive integration: automatically refine the interval count until the result converges within a specified tolerance via --tolerance flag (default: 1e-10)
- Parse mathematical expressions supporting: basic operators (+, -, *, /, ^), standard functions (sin, cos, tan, exp, log, sqrt, abs), constants (pi, e), and nested parentheses
- Support improper integrals with infinite bounds via --infinite flag using appropriate limit-based techniques
- Support tabulated data integration: read (x, y) pairs from a CSV file via --data flag and integrate using the available methods
- Print results to console: integral value, estimated error, method used, intervals used, and computation time
- Save results as JSON with --output flag (default: integration_result.json)
- If no input is given, compute several well-known integrals (integral of sin(x) from 0 to pi = 2, integral of exp(-x^2) from 0 to infinity = sqrt(pi)/2, integral of 1/x from 1 to e = 1) using all methods and display a comparison table with exact vs computed values and relative errors
- Handle errors: division by zero within the integration range, non-convergent integrals, invalid mathematical expressions, and bounds where lower > upper

Create a complete C++ project for a clean Ubuntu 22.04 machine with only G++ 12+ and CMake 3.22+ installed. Include:
- Source code
- CMakeLists.txt with all dependencies (direct and transitive) pinned to exact versions
- README.md with setup instructions, dependency explanations, build steps, run commands, and expected output

[]

cmake_minimum_required(VERSION 3.22)
project(numerical_integrator VERSION 1.0.0 LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

# ---------------------------------------------------------------------------
# Dependencies
# ---------------------------------------------------------------------------
# GSL (GNU Scientific Library) — system installed
# nlohmann/json v3.11.3 — header-only JSON library
# ---------------------------------------------------------------------------

find_package(GSL REQUIRED)

include(FetchContent)

FetchContent_Declare(
    nlohmann_json
    GIT_REPOSITORY https://github.com/nlohmann/json.git
    GIT_TAG        v3.11.3
    GIT_SHALLOW    TRUE
)

set(JSON_BuildTests   OFF CACHE BOOL "" FORCE)
set(JSON_Install      OFF CACHE BOOL "" FORCE)
set(JSON_MultipleHeaders OFF CACHE BOOL "" FORCE)

FetchContent_MakeAvailable(nlohmann_json)

# ---------------------------------------------------------------------------
# Executable
# ---------------------------------------------------------------------------

add_executable(numerical_integrator src/main.cpp)

target_link_libraries(numerical_integrator PRIVATE
    GSL::gsl
    GSL::gslcblas
    nlohmann_json::nlohmann_json
)

if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
    target_compile_options(numerical_integrator PRIVATE -Wall -Wextra -Wpedantic)
endif()

# Numerical Integration Calculator (C++)

## Requirements
- G++ 12+
- CMake 3.22+
- GSL
- OpenSSL-compatible toolchain for linked packages from CMake config

## Build
```bash
cmake -S . -B build
cmake --build build --config Release
```

## Run
```bash
./build/numerical_integrator
```

## Dependencies
- GNU Scientific Library (GSL)
- `nlohmann/json` pinned via `FetchContent` to `v3.11.3`

/**
 * Numerical Integration Calculator
 * Computes definite integrals using Trapezoidal, Simpson's, Gauss-Legendre,
 * and Romberg methods with convergence comparison.
 *
 * Dependencies: GSL (system), nlohmann/json 3.11.3
 */

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
#include <chrono>
#include <algorithm>
#include <numeric>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <nlohmann/json.hpp>

using json = nlohmann::json;
using Func = std::function<double(double)>;

// ---------------------------------------------------------------------------
// Trapezoidal Rule
// ---------------------------------------------------------------------------
double trapezoidal(const Func& f, double a, double b, int n) {
    double h = (b - a) / n;
    double sum = 0.5 * (f(a) + f(b));
    for (int i = 1; i < n; i++) {
        sum += f(a + i * h);
    }
    return sum * h;
}

// ---------------------------------------------------------------------------
// Simpson's 1/3 Rule
// ---------------------------------------------------------------------------
double simpsons(const Func& f, double a, double b, int n) {
    if (n % 2 != 0) n++;
    double h = (b - a) / n;
    double sum = f(a) + f(b);
    for (int i = 1; i < n; i += 2) {
        sum += 4.0 * f(a + i * h);
    }
    for (int i = 2; i < n; i += 2) {
        sum += 2.0 * f(a + i * h);
    }
    return sum * h / 3.0;
}

// ---------------------------------------------------------------------------
// Gauss-Legendre Quadrature (using GSL)
// ---------------------------------------------------------------------------
double gauss_legendre(const Func& f, double a, double b, int n) {
    gsl_integration_glfixed_table* table = gsl_integration_glfixed_table_alloc(n);
    double result = 0.0;
    for (int i = 0; i < n; i++) {
        double xi, wi;
        gsl_integration_glfixed_point(a, b, i, &xi, &wi, table);
        result += wi * f(xi);
    }
    gsl_integration_glfixed_table_free(table);
    return result;
}

// ---------------------------------------------------------------------------
// Romberg Integration
// ---------------------------------------------------------------------------
struct RombergResult {
    double value;
    int order;
    std::vector<std::vector<double>> table;
};

RombergResult romberg(const Func& f, double a, double b, int maxOrder = 10) {
    std::vector<std::vector<double>> R(maxOrder, std::vector<double>(maxOrder, 0.0));
    double h = b - a;
    R[0][0] = 0.5 * h * (f(a) + f(b));

    for (int i = 1; i < maxOrder; i++) {
        double hi = h / std::pow(2.0, i);
        int numNew = static_cast<int>(std::pow(2.0, i - 1));
        double sum = 0.0;
        for (int k = 1; k <= numNew; k++) {
            sum += f(a + hi * (2 * k - 1));
        }
        R[i][0] = 0.5 * R[i - 1][0] + hi * sum;

        for (int j = 1; j <= i; j++) {
            double factor = std::pow(4.0, j);
            R[i][j] = (factor * R[i][j - 1] - R[i - 1][j - 1]) / (factor - 1.0);
        }

        if (std::abs(R[i][i] - R[i - 1][i - 1]) < 1e-12) {
            std::vector<std::vector<double>> trimmed(R.begin(), R.begin() + i + 1);
            return {R[i][i], i + 1, trimmed};
        }
    }
    return {R[maxOrder - 1][maxOrder - 1], maxOrder, R};
}

// ---------------------------------------------------------------------------
// GSL adaptive integration (reference)
// ---------------------------------------------------------------------------
double gsl_wrapper_fn(double x, void* params) {
    auto* fp = static_cast<Func*>(params);
    return (*fp)(x);
}

double gsl_reference(Func f, double a, double b, double& abserr) {
    gsl_integration_workspace* w = gsl_integration_workspace_alloc(10000);
    gsl_function F;
    F.function = &gsl_wrapper_fn;
    F.params = &f;
    double result;
    gsl_integration_qags(&F, a, b, 1e-14, 1e-14, 10000, w, &result, &abserr);
    gsl_integration_workspace_free(w);
    return result;
}

// ---------------------------------------------------------------------------
// Test functions
// ---------------------------------------------------------------------------
struct TestFunction {
    std::string name;
    std::string description;
    Func f;
    double a, b;
    double exactValue;
    bool hasExact;
};

std::vector<TestFunction> createTestFunctions() {
    std::vector<TestFunction> tests;

    double exactPoly = 3.0 * std::pow(2, 5) / 5.0 - 2.0 * std::pow(2, 4) / 4.0
                     + std::pow(2, 3) / 3.0 - 5.0 * 4.0 / 2.0 + 14.0;
    tests.push_back({"polynomial", "3x^4 - 2x^3 + x^2 - 5x + 7 on [0, 2]",
        [](double x) { return 3*x*x*x*x - 2*x*x*x + x*x - 5*x + 7; },
        0, 2, exactPoly, true});

    tests.push_back({"trigonometric", "sin(x)*cos(x) on [0, pi/2]",
        [](double x) { return std::sin(x) * std::cos(x); },
        0, M_PI / 2.0, 0.5, true});

    tests.push_back({"exponential", "exp(-x^2) on [0, 1] (Gaussian)",
        [](double x) { return std::exp(-x * x); },
        0, 1, 0.7468241328124271, true});

    tests.push_back({"oscillatory", "sin(10x)*exp(-x) on [0, pi]",
        [](double x) { return std::sin(10 * x) * std::exp(-x); },
        0, M_PI, 0.0, false});

    double exactSing = 2.0 * (1.0 - std::sqrt(1e-10));
    tests.push_back({"singular_endpoint", "1/sqrt(x) on [~0, 1] (near-singular)",
        [](double x) { return x > 0 ? 1.0 / std::sqrt(x) : 0.0; },
        1e-10, 1, exactSing, true});

    return tests;
}

// ---------------------------------------------------------------------------
// Convergence study
// ---------------------------------------------------------------------------
json convergenceStudy(const Func& f, double a, double b, double exact, int maxN) {
    std::vector<int> ns;
    for (int k = 1; (1 << k) <= maxN; k++) ns.push_back(1 << k);

    std::vector<double> errTrap, errSimp, errGauss, errRomb;

    for (int n : ns) {
        errTrap.push_back(std::abs(trapezoidal(f, a, b, n) - exact));
        int nSimp = (n % 2 == 0) ? n : n + 1;
        errSimp.push_back(std::abs(simpsons(f, a, b, nSimp) - exact));
        int nGauss = std::min(n, 64);
        errGauss.push_back(std::abs(gauss_legendre(f, a, b, nGauss) - exact));
        int order = std::max(2, static_cast<int>(std::log2(n)));
        errRomb.push_back(std::abs(romberg(f, a, b, order).value - exact));
    }

    return json{
        {"n_values", ns},
        {"trapezoidal_errors", errTrap},
        {"simpsons_errors", errSimp},
        {"gauss_errors", errGauss},
        {"romberg_errors", errRomb}
    };
}

// ---------------------------------------------------------------------------
// Main
// ---------------------------------------------------------------------------
int main() {
    std::cout << std::string(70, '=') << "\n";
    std::cout << "  Numerical Integration Calculator\n";
    std::cout << "  Methods: Trapezoidal, Simpson's, Gauss-Legendre, Romberg\n";
    std::cout << std::string(70, '=') << "\n";

    auto tests = createTestFunctions();
    json allResults;

    for (auto& test : tests) {
        double exact = test.exactValue;
        if (!test.hasExact) {
            double abserr;
            exact = gsl_reference(test.f, test.a, test.b, abserr);
            test.exactValue = exact;
        }

        std::cout << "\n" << std::string(70, '=') << "\n";
        std::cout << "Test: " << test.description << "\n";
        std::cout << std::fixed << std::setprecision(15);
        std::cout << "Exact value: " << exact << "\n";
        std::cout << std::string(70, '=') << "\n";

        json testResults;
        int nPoints = 100;

        // Trapezoidal
        auto t0 = std::chrono::high_resolution_clock::now();
        double val = trapezoidal(test.f, test.a, test.b, nPoints);
        auto t1 = std::chrono::high_resolution_clock::now();
        double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        testResults["Trapezoidal"] = {{"value", val}, {"error", std::abs(val - exact)}, {"time_ms", ms}};

        // Simpson's
        t0 = std::chrono::high_resolution_clock::now();
        val = simpsons(test.f, test.a, test.b, nPoints);
        t1 = std::chrono::high_resolution_clock::now();
        ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        testResults["Simpson's"] = {{"value", val}, {"error", std::abs(val - exact)}, {"time_ms", ms}};

        // Gauss-Legendre
        t0 = std::chrono::high_resolution_clock::now();
        val = gauss_legendre(test.f, test.a, test.b, 20);
        t1 = std::chrono::high_resolution_clock::now();
        ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        testResults["Gauss-Legendre"] = {{"value", val}, {"error", std::abs(val - exact)}, {"time_ms", ms}};

        // Romberg
        t0 = std::chrono::high_resolution_clock::now();
        auto rResult = romberg(test.f, test.a, test.b, 10);
        val = rResult.value;
        t1 = std::chrono::high_resolution_clock::now();
        ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        testResults["Romberg"] = {{"value", val}, {"error", std::abs(val - exact)}, {"time_ms", ms}};

        // GSL reference
        t0 = std::chrono::high_resolution_clock::now();
        double abserr;
        val = gsl_reference(test.f, test.a, test.b, abserr);
        t1 = std::chrono::high_resolution_clock::now();
        ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        testResults["GSL-QAGS"] = {{"value", val}, {"error", std::abs(val - exact)}, {"time_ms", ms}};

        // Print table
        std::cout << std::left;
        std::cout << "\n" << std::setw(20) << "Method" << std::setw(22) << "Result"
                  << std::setw(15) << "Error" << std::setw(12) << "Time (ms)" << "\n";
        std::cout << std::string(70, '-') << "\n";

        for (auto& [method, data] : testResults.items()) {
            std::cout << std::setw(20) << method
                      << std::setw(22) << std::fixed << std::setprecision(15) << data["value"].get<double>()
                      << std::setw(15) << std::scientific << std::setprecision(2) << data["error"].get<double>()
                      << std::setw(12) << std::fixed << std::setprecision(4) << data["time_ms"].get<double>()
                      << "\n";
        }

        allResults[test.name] = testResults;
    }

    // Convergence study
    std::cout << "\n\n" << std::string(70, '=') << "\n";
    std::cout << "  Convergence Study: exp(-x^2) on [0, 1]\n";
    std::cout << std::string(70, '=') << "\n";

    auto fGauss = [](double x) { return std::exp(-x * x); };
    double exactGauss = 0.7468241328124271;
    json conv = convergenceStudy(fGauss, 0, 1, exactGauss, 512);

    auto& ns = conv["n_values"];
    std::cout << "\n" << std::setw(8) << "n" << std::setw(15) << "Trapezoidal"
              << std::setw(15) << "Simpson's" << std::setw(15) << "Gauss-Legendre"
              << std::setw(15) << "Romberg" << "\n";
    std::cout << std::string(68, '-') << "\n";
    for (size_t i = 0; i < ns.size(); i++) {
        std::cout << std::left << std::setw(8) << ns[i].get<int>()
                  << std::scientific << std::setprecision(2)
                  << std::setw(15) << conv["trapezoidal_errors"][i].get<double>()
                  << std::setw(15) << conv["simpsons_errors"][i].get<double>()
                  << std::setw(15) << conv["gauss_errors"][i].get<double>()
                  << std::setw(15) << conv["romberg_errors"][i].get<double>()
                  << "\n";
    }

    // JSON output
    json output;
    output["integration_results"] = allResults;
    output["convergence_study"] = conv;
    std::cout << "\n--- JSON Output ---\n" << output.dump(2) << "\n";

    // Summary
    std::cout << "\n" << std::string(70, '=') << "\n";
    std::cout << "  Summary: Best method for each test function\n";
    std::cout << std::string(70, '=') << "\n";
    for (auto& test : tests) {
        std::string bestMethod;
        double bestError = 1e300;
        for (auto& [method, data] : allResults[test.name].items()) {
            double err = data["error"].get<double>();
            if (err < bestError) {
                bestError = err;
                bestMethod = method;
            }
        }
        std::cout << "  " << std::left << std::setw(25) << test.name
                  << " -> " << std::setw(20) << bestMethod
                  << " (error: " << std::scientific << std::setprecision(2) << bestError << ")\n";
    }

    std::cout << "\nDone.\n";
    return 0;
}