Numerical Integration Calculator (cpp, written by Claude Code)

--- /dev/null
+++ b/setup.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+# System packages this project needs on a clean Ubuntu machine.
+set -e
+export DEBIAN_FRONTEND=noninteractive
+apt-get update -qq
+apt-get install -y -qq --no-install-recommends libboost-dev

Claude Code 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

[
  {
    "category": "underspecification",
    "rule": "signature.missing_system_requirement",
    "source": "failure_signature",
    "excerpt": "Could NOT find Boost (missing: Boost_INCLUDE_DIR)"
  },
  {
    "category": "underspecification",
    "rule": "diff.adds_external_environment_requirement",
    "source": "manifest_diff:setup.sh",
    "excerpt": "export DEBIAN_FRONTEND=noninteractive"
  },
  {
    "category": "underspecification",
    "rule": "diff.adds_external_environment_requirement",
    "source": "manifest_diff:setup.sh",
    "excerpt": "apt-get install -y -qq --no-install-recommends libboost-dev"
  }
]

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

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(Boost REQUIRED)

# CLI11 - fetch via FetchContent or find locally
include(FetchContent)
FetchContent_Declare(
    cli11
    GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
    GIT_TAG v2.4.1
)
FetchContent_MakeAvailable(cli11)

add_executable(numerical_integrator main.cpp)

target_include_directories(numerical_integrator PRIVATE ${Boost_INCLUDE_DIRS})
target_link_libraries(numerical_integrator PRIVATE CLI11::CLI11)
target_compile_options(numerical_integrator PRIVATE -Wall -Wextra -O2)

/**
 * Numerical Integration Calculator (C++ - Trial 2)
 *
 * Computes definite integrals via trapezoidal, Simpson, Gauss-Legendre,
 * and Romberg methods with convergence comparison.
 *
 * Dependencies: boost-math, CLI11
 */

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
#include <map>

#include <boost/math/quadrature/trapezoidal.hpp>
#include <boost/math/quadrature/gauss.hpp>
#include <boost/math/constants/constants.hpp>
#include <CLI/CLI.hpp>

using Function = std::function<double(double)>;
constexpr double PI = boost::math::constants::pi<double>();

// ---------------------------------------------------------------------------
// Test functions
// ---------------------------------------------------------------------------

struct TestFunction {
    std::string name;
    Function fn;
};

std::map<std::string, TestFunction> get_test_functions() {
    return {
        {"sin", {"sin(x)", [](double x) { return std::sin(x); }}},
        {"poly", {"x^3 - 2x^2 + 3x - 1", [](double x) { return x*x*x - 2*x*x + 3*x - 1; }}},
        {"exp", {"exp(-x^2)", [](double x) { return std::exp(-x*x); }}},
        {"gaussian", {"N(0,1) PDF", [](double x) {
            return (1.0 / std::sqrt(2*PI)) * std::exp(-0.5*x*x);
        }}},
        {"oscillating", {"sin(10x)*exp(-x)", [](double x) {
            return std::sin(10*x) * std::exp(-x);
        }}},
        {"rational", {"1/(1+x^2)", [](double x) { return 1.0 / (1.0 + x*x); }}},
    };
}

// ---------------------------------------------------------------------------
// Custom trapezoidal rule
// ---------------------------------------------------------------------------

double trapezoidal_rule(const Function& 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_rule(const Function& 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++) {
        double x = a + i * h;
        sum += (i % 2 == 0) ? 2 * f(x) : 4 * f(x);
    }
    return sum * h / 3.0;
}

// ---------------------------------------------------------------------------
// Romberg integration
// ---------------------------------------------------------------------------

double romberg_integration(const Function& f, double a, double b, int max_order) {
    std::vector<std::vector<double>> R(max_order + 1, std::vector<double>(max_order + 1, 0.0));
    R[0][0] = trapezoidal_rule(f, a, b, 1);

    for (int i = 1; i <= max_order; i++) {
        int n = static_cast<int>(std::pow(2, i));
        R[i][0] = trapezoidal_rule(f, a, b, n);

        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);
        }
    }
    return R[max_order][max_order];
}

// ---------------------------------------------------------------------------
// Convergence analysis
// ---------------------------------------------------------------------------

void convergence_analysis(const Function& f, double a, double b) {
    std::cout << "\n--- Convergence Analysis ---\n";
    std::cout << std::left << std::setw(12) << "Intervals"
              << std::setw(22) << "Trapezoidal"
              << std::setw(22) << "Simpson"
              << std::setw(15) << "Trap Error"
              << std::setw(15) << "Simp Error" << "\n";
    std::cout << std::string(86, '-') << "\n";

    double prev_trap = 0, prev_simp = 0;
    for (int n = 4; n <= 4096; n *= 2) {
        double trap = trapezoidal_rule(f, a, b, n);
        double simp = simpsons_rule(f, a, b, n);

        std::cout << std::left << std::setw(12) << n
                  << std::fixed << std::setprecision(15)
                  << std::setw(22) << trap
                  << std::setw(22) << simp;

        if (n > 4) {
            std::cout << std::scientific << std::setprecision(2)
                      << std::setw(15) << std::abs(trap - prev_trap)
                      << std::setw(15) << std::abs(simp - prev_simp);
        } else {
            std::cout << std::setw(15) << "---" << std::setw(15) << "---";
        }
        std::cout << "\n";

        prev_trap = trap;
        prev_simp = simp;
    }
}

// ---------------------------------------------------------------------------
// Method comparison using Boost.Math
// ---------------------------------------------------------------------------

void method_comparison(const Function& f, double a, double b) {
    std::cout << "\n--- Method Comparison ---\n";
    std::cout << std::left << std::setw(32) << "Method"
              << std::setw(22) << "Result" << "\n";
    std::cout << std::string(54, '-') << "\n";

    // Custom trapezoidal
    double trap = trapezoidal_rule(f, a, b, 2048);
    std::cout << std::left << std::setw(32) << "Trapezoidal (n=2048)"
              << std::fixed << std::setprecision(15) << trap << "\n";

    // Custom Simpson
    double simp = simpsons_rule(f, a, b, 2048);
    std::cout << std::left << std::setw(32) << "Simpson (n=2048)"
              << std::fixed << std::setprecision(15) << simp << "\n";

    // Romberg
    double romb = romberg_integration(f, a, b, 10);
    std::cout << std::left << std::setw(32) << "Romberg (order=10)"
              << std::fixed << std::setprecision(15) << romb << "\n";

    // Boost trapezoidal
    double boost_trap = boost::math::quadrature::trapezoidal(f, a, b);
    std::cout << std::left << std::setw(32) << "Boost Trapezoidal"
              << std::fixed << std::setprecision(15) << boost_trap << "\n";

    // Boost Gauss-Legendre (various orders)
    double g7 = boost::math::quadrature::gauss<double, 7>::integrate(f, a, b);
    std::cout << std::left << std::setw(32) << "Boost Gauss-Legendre (7)"
              << std::fixed << std::setprecision(15) << g7 << "\n";

    double g15 = boost::math::quadrature::gauss<double, 15>::integrate(f, a, b);
    std::cout << std::left << std::setw(32) << "Boost Gauss-Legendre (15)"
              << std::fixed << std::setprecision(15) << g15 << "\n";

    double g30 = boost::math::quadrature::gauss<double, 30>::integrate(f, a, b);
    std::cout << std::left << std::setw(32) << "Boost Gauss-Legendre (30)"
              << std::fixed << std::setprecision(15) << g30 << "\n";
}

// ---------------------------------------------------------------------------
// Multi-function comparison
// ---------------------------------------------------------------------------

void multi_function_comparison() {
    std::cout << "\n--- Multi-Function Comparison (Romberg order=10) ---\n";
    std::cout << std::left << std::setw(25) << "Function"
              << std::setw(15) << "Interval"
              << std::setw(22) << "Romberg Result" << "\n";
    std::cout << std::string(62, '-') << "\n";

    auto funcs = get_test_functions();
    struct Config { std::string key; double a; double b; };
    std::vector<Config> configs = {
        {"sin", 0, PI}, {"poly", 0, 3}, {"exp", 0, 2},
        {"gaussian", -3, 3}, {"oscillating", 0, 5}, {"rational", 0, 1},
    };

    for (const auto& cfg : configs) {
        double result = romberg_integration(funcs[cfg.key].fn, cfg.a, cfg.b, 10);
        std::ostringstream interval;
        interval << std::fixed << std::setprecision(1) << "[" << cfg.a << ", " << cfg.b << "]";
        std::cout << std::left << std::setw(25) << funcs[cfg.key].name
                  << std::setw(15) << interval.str()
                  << std::fixed << std::setprecision(15) << result << "\n";
    }
}

// ---------------------------------------------------------------------------
// Main
// ---------------------------------------------------------------------------

int main(int argc, char** argv) {
    CLI::App app{"Numerical Integration Calculator (boost-math + CLI11)"};

    std::string func_name = "sin";
    double lower = 0.0;
    double upper = PI;

    app.add_option("-f,--function", func_name, "Function to integrate (sin, poly, exp, gaussian, oscillating, rational)")
       ->default_val("sin");
    app.add_option("-a,--lower", lower, "Lower bound")->default_val(0.0);
    app.add_option("-b,--upper", upper, "Upper bound")->default_val(PI);

    CLI11_PARSE(app, argc, argv);

    auto funcs = get_test_functions();
    if (funcs.find(func_name) == funcs.end()) {
        std::cerr << "Unknown function: " << func_name << "\n";
        return 1;
    }

    std::cout << std::string(60, '=') << "\n";
    std::cout << "Numerical Integration Calculator\n";
    std::cout << "  (boost-math + CLI11)\n";
    std::cout << std::string(60, '=') << "\n";

    const auto& tf = funcs[func_name];
    std::cout << "\nFunction: " << tf.name << "\n";
    std::cout << std::fixed << std::setprecision(4);
    std::cout << "Interval: [" << lower << ", " << upper << "]\n";

    convergence_analysis(tf.fn, lower, upper);
    method_comparison(tf.fn, lower, upper);
    multi_function_comparison();

    std::cout << "\n" << std::string(60, '=') << "\n";
    std::cout << "Integration complete.\n";
    std::cout << std::string(60, '=') << "\n";

    return 0;
}

/**
 * Numerical Integration - Definite integrals via trapezoidal/Simpson/Gauss/Romberg.
 *
 * Uses GSL (GNU Scientific Library) for numerical integration and CLI11 for CLI parsing.
 */

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

#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <CLI/CLI.hpp>

const double PI = M_PI;

// === Custom Integration Methods ===

double trapezoidalRule(std::function<double(double)> 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;
}

double simpsonsRule(std::function<double(double)> 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++) {
        sum += ((i % 2 == 1) ? 4.0 : 2.0) * f(a + i * h);
    }
    return sum * h / 3.0;
}

double midpointRule(std::function<double(double)> f, double a, double b, int n) {
    double h = (b - a) / n;
    double sum = 0;
    for (int i = 0; i < n; i++) {
        sum += f(a + (i + 0.5) * h);
    }
    return sum * h;
}

double rombergIntegration(std::function<double(double)> f, double a, double b,
                          int maxOrder = 12, double tol = 1e-14) {
    std::vector<std::vector<double>> R(maxOrder + 1, std::vector<double>(maxOrder + 1, 0));
    R[0][0] = (b - a) / 2.0 * (f(a) + f(b));

    for (int i = 1; i <= maxOrder; i++) {
        int nInt = 1 << i;
        double h = (b - a) / nInt;
        double sum = 0;
        for (int k = 1; k <= nInt / 2; k++) {
            sum += f(a + (2 * k - 1) * h);
        }
        R[i][0] = 0.5 * R[i - 1][0] + h * sum;

        for (int j = 1; j <= i; j++) {
            double fac = std::pow(4.0, j);
            R[i][j] = (fac * R[i][j - 1] - R[i - 1][j - 1]) / (fac - 1);
        }
        if (i > 1 && std::abs(R[i][i] - R[i - 1][i - 1]) < tol) {
            return R[i][i];
        }
    }
    return R[maxOrder][maxOrder];
}

double gaussLegendre5(std::function<double(double)> f, double a, double b) {
    static const double nodes[] = {-0.906179845938664, -0.538469310105683, 0.0,
                                    0.538469310105683, 0.906179845938664};
    static const double weights[] = {0.236926885056189, 0.478628670499366, 0.568888888888889,
                                      0.478628670499366, 0.236926885056189};
    double sum = 0;
    for (int i = 0; i < 5; i++) {
        double x = (b - a) / 2.0 * nodes[i] + (a + b) / 2.0;
        sum += weights[i] * f(x);
    }
    return sum * (b - a) / 2.0;
}

double compositeGauss(std::function<double(double)> f, double a, double b, int panels) {
    double h = (b - a) / panels;
    double total = 0;
    for (int i = 0; i < panels; i++) {
        total += gaussLegendre5(f, a + i * h, a + (i + 1) * h);
    }
    return total;
}

double adaptiveSimpson(std::function<double(double)> f, double a, double b,
                       double tol = 1e-12, int maxDepth = 50) {
    std::function<double(double, double, double, double, int)> rec =
        [&](double a, double b, double tol, double whole, int depth) -> double {
        double mid = (a + b) / 2;
        double left = simpsonsRule(f, a, mid, 4);
        double right = simpsonsRule(f, mid, b, 4);
        double total = left + right;
        if (depth >= maxDepth || std::abs(total - whole) < 15 * tol) {
            return total + (total - whole) / 15.0;
        }
        return rec(a, mid, tol / 2, left, depth + 1) + rec(mid, b, tol / 2, right, depth + 1);
    };
    double whole = simpsonsRule(f, a, b, 4);
    return rec(a, b, tol, whole, 0);
}

// === GSL Integration Wrapper ===

struct GslFuncWrapper {
    std::function<double(double)> func;
};

double gslCallback(double x, void* params) {
    auto* wrapper = static_cast<GslFuncWrapper*>(params);
    return wrapper->func(x);
}

double gslQAGS(std::function<double(double)> f, double a, double b) {
    GslFuncWrapper wrapper{f};
    gsl_function gslF;
    gslF.function = &gslCallback;
    gslF.params = &wrapper;

    gsl_integration_workspace* ws = gsl_integration_workspace_alloc(1000);
    double result, error;
    gsl_integration_qags(&gslF, a, b, 1e-12, 1e-12, 1000, ws, &result, &error);
    gsl_integration_workspace_free(ws);
    return result;
}

double gslQAG(std::function<double(double)> f, double a, double b) {
    GslFuncWrapper wrapper{f};
    gsl_function gslF;
    gslF.function = &gslCallback;
    gslF.params = &wrapper;

    gsl_integration_workspace* ws = gsl_integration_workspace_alloc(1000);
    double result, error;
    gsl_integration_qag(&gslF, a, b, 1e-12, 1e-12, 1000, GSL_INTEG_GAUSS61, ws, &result, &error);
    gsl_integration_workspace_free(ws);
    return result;
}

// === Test Functions ===

struct TestFunction {
    std::string name;
    std::function<double(double)> f;
    double a, b, exact;
};

std::vector<TestFunction> getTestFunctions() {
    return {
        {"sin(x)", [](double x) { return std::sin(x); }, 0, PI, 2.0},
        {"1/(1+x^2)", [](double x) { return 1.0 / (1 + x * x); }, 0, 1, PI / 4.0},
        {"exp(-x^2)", [](double x) { return std::exp(-x * x); }, 0, 3, 0.8862073482595214},
        {"x*sin(x^2)", [](double x) { return x * std::sin(x * x); }, 0, std::sqrt(PI),
            0.5 * (1 - std::cos(PI))},
        {"cos(x)*exp(-x)", [](double x) { return std::cos(x) * std::exp(-x); }, 0, PI,
            0.5 * (1.0 - std::exp(-PI))},
    };
}

// === Display ===

void printRow(const std::string& method, int n, double result, double exact, double timeUs) {
    std::string nStr = n > 0 ? std::to_string(n) : "auto";
    double error = std::abs(result - exact);
    std::cout << "  " << std::left << std::setw(24) << method
              << std::right << std::setw(8) << nStr
              << std::setw(21) << std::fixed << std::setprecision(15) << result
              << "  " << std::scientific << std::setprecision(2) << std::setw(12) << error
              << std::fixed << std::setprecision(1) << std::setw(10) << timeUs << "\n";
}

void evaluateFunction(const TestFunction& tf, int baseIntervals) {
    std::cout << "\n--- " << tf.name << " on [" << tf.a << ", " << tf.b << "] ---\n";
    std::cout << std::fixed << std::setprecision(15) << "  Exact: " << tf.exact << "\n";
    std::cout << "  " << std::left << std::setw(24) << "Method"
              << std::right << std::setw(8) << "N"
              << std::setw(21) << "Result"
              << std::setw(14) << "Error"
              << std::setw(10) << "Time(us)" << "\n";
    std::cout << "  " << std::string(77, '-') << "\n";

    for (int n : {baseIntervals / 10, baseIntervals, baseIntervals * 10}) {
        if (n < 4) n = 4;
        auto t0 = std::chrono::high_resolution_clock::now();
        double r = trapezoidalRule(tf.f, tf.a, tf.b, n);
        auto t1 = std::chrono::high_resolution_clock::now();
        printRow("Trapezoidal", n, r, tf.exact,
                 std::chrono::duration<double, std::micro>(t1 - t0).count());

        t0 = std::chrono::high_resolution_clock::now();
        r = simpsonsRule(tf.f, tf.a, tf.b, n);
        t1 = std::chrono::high_resolution_clock::now();
        printRow("Simpson 1/3", n, r, tf.exact,
                 std::chrono::duration<double, std::micro>(t1 - t0).count());
    }

    auto t0 = std::chrono::high_resolution_clock::now();
    double r = rombergIntegration(tf.f, tf.a, tf.b);
    auto t1 = std::chrono::high_resolution_clock::now();
    printRow("Romberg", -1, r, tf.exact,
             std::chrono::duration<double, std::micro>(t1 - t0).count());

    t0 = std::chrono::high_resolution_clock::now();
    r = compositeGauss(tf.f, tf.a, tf.b, 20);
    t1 = std::chrono::high_resolution_clock::now();
    printRow("Gauss-Legendre(5x20)", -1, r, tf.exact,
             std::chrono::duration<double, std::micro>(t1 - t0).count());

    t0 = std::chrono::high_resolution_clock::now();
    r = adaptiveSimpson(tf.f, tf.a, tf.b);
    t1 = std::chrono::high_resolution_clock::now();
    printRow("Adaptive Simpson", -1, r, tf.exact,
             std::chrono::duration<double, std::micro>(t1 - t0).count());

    t0 = std::chrono::high_resolution_clock::now();
    r = gslQAGS(tf.f, tf.a, tf.b);
    t1 = std::chrono::high_resolution_clock::now();
    printRow("GSL QAGS", -1, r, tf.exact,
             std::chrono::duration<double, std::micro>(t1 - t0).count());

    t0 = std::chrono::high_resolution_clock::now();
    r = gslQAG(tf.f, tf.a, tf.b);
    t1 = std::chrono::high_resolution_clock::now();
    printRow("GSL QAG (Gauss61)", -1, r, tf.exact,
             std::chrono::duration<double, std::micro>(t1 - t0).count());
}

int main(int argc, char** argv) {
    CLI::App app{"Numerical Integration - GSL + CLI11"};

    int intervals = 100;
    double tolerance = 1e-12;
    bool verbose = false;

    app.add_option("-n,--intervals", intervals, "Base number of intervals")->default_val(100);
    app.add_option("-t,--tolerance", tolerance, "Convergence tolerance")->default_val(1e-12);
    app.add_flag("-v,--verbose", verbose, "Verbose output");

    CLI11_PARSE(app, argc, argv);

    std::cout << "=== Numerical Integration (gsl + CLI11) ===\n";
    std::cout << "Intervals: " << intervals << ", Tolerance: " << std::scientific << tolerance << "\n";

    auto functions = getTestFunctions();
    for (const auto& tf : functions) {
        evaluateFunction(tf, intervals);
    }

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

# Numerical Integration Calculator (C++ - Trial 2)

Computes definite integrals via trapezoidal, Simpson, Gauss-Legendre, and Romberg methods with convergence comparison.

## Dependencies

- **boost-math** - Boost.Math quadrature (trapezoidal, Gauss-Legendre) and mathematical constants
- **CLI11** (2.4.1) - Command-line argument parsing library

## Build

```bash
mkdir build && cd build
cmake ..
make
```

### Prerequisites

- Ubuntu/Debian: `sudo apt install libboost-dev`
- macOS: `brew install boost`
- CLI11 is fetched automatically via CMake FetchContent

## Usage

```bash
./numerical_integrator
./numerical_integrator --function sin --lower 0 --upper 3.14159
./numerical_integrator -f poly -a 0 -b 3
```

## Features

- Trapezoidal rule with convergence tracking
- Simpson's 1/3 rule implementation
- Romberg integration with Richardson extrapolation
- Boost.Math trapezoidal quadrature
- Boost.Math Gauss-Legendre quadrature (7, 15, 30 nodes)
- Multi-function comparison (sin, polynomial, exponential, Gaussian, oscillating, rational)
- CLI with CLI11 for configurable parameters