QUAD

=QUAD(function_table, a, b, [epsabs], [epsrel])

=QUAD({0,0;1,1;2,4;3,9;4,16}, 0, 4)

=QUAD({0,0;1.5708,1;3.1416,0}, 0, 3.1416, 1E-9)

=QUAD({0,1;2.5,0.0821;5,0.0067}, 0, 5, 1E-9, 1E-9)

=QUAD({1,3;2,5;3,7}, 1, 3, , 1E-9)

from bisect import bisect_left
import math
from scipy.integrate import quad as scipy_quad
 
 
def quad(function_table, a, b, epsabs=1.49e-8, epsrel=1.49e-8):
    """
    Numerically integrate a function defined by a table of x, y values over [a, b] using adaptive quadrature.
 
    Args:
        function_table (2D list): Each row is [x, y]; must contain at least two rows with strictly increasing x values.
        a (float): Lower limit of integration.
        b (float): Upper limit of integration (must be greater than or equal to a).
        epsabs (float, optional): Absolute error tolerance (default: 1.49e-8).
        epsrel (float, optional): Relative error tolerance (default: 1.49e-8).
 
    Returns:
        The estimated value of the definite integral (float), or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
 
    def _convert_float_any(value, name):
        try:
            numeric = float(value)
        except Exception:
            return f"Invalid input: {name} must be a number."
        if math.isnan(numeric) or math.isinf(numeric):
            return f"Invalid input: {name} must be finite."
        return numeric
 
    def _convert_tolerance(value, name, *, allow_zero: bool):
        converted = _convert_float_any(value, name)
        if isinstance(converted, str):
            return converted
        if allow_zero:
            if converted < 0.0:
                return f"Invalid input: {name} must be non-negative."
        else:
            if converted <= 0.0:
                return f"Invalid input: {name} must be greater than 0."
        return converted
 
    # Normalize 2D-list args coming from Excel: single-element 2D lists may be passed as scalars.
    def to2d(x):
        return [[x]] if not isinstance(x, list) else x
 
    function_table = to2d(function_table)
 
    if not isinstance(function_table, list) or len(function_table) < 2:
        return "Invalid input: function_table must be a 2D list with at least two rows."
 
    processed_rows = []
    for row in function_table:
        if not isinstance(row, list) or len(row) < 2:
            return "Invalid input: each row in function_table must contain at least two values."
        try:
            processed_rows.append([float(row[0]), float(row[1])])
        except Exception:
            return "Invalid input: function_table must contain numeric values."
 
    xs = [row[0] for row in processed_rows]
    ys = [row[1] for row in processed_rows]
 
    if any(x2 <= x1 for x1, x2 in zip(xs, xs[1:])):
        return "Invalid input: x values in function_table must be strictly increasing."
 
    converted_a = _convert_float_any(a, "a")
    if isinstance(converted_a, str):
        return converted_a
    converted_b = _convert_float_any(b, "b")
    if isinstance(converted_b, str):
        return converted_b
 
    if converted_b < converted_a:
        return "Invalid input: a must be less than or equal to b."
 
    x_min, x_max = xs[0], xs[-1]
    if converted_a < x_min or converted_b > x_max:
        return "Invalid input: integration limits must fall within the range of function_table x values."
 
    converted_epsabs = _convert_tolerance(epsabs, "epsabs", allow_zero=True)
    if isinstance(converted_epsabs, str):
        return converted_epsabs
    converted_epsrel = _convert_tolerance(epsrel, "epsrel", allow_zero=False)
    if isinstance(converted_epsrel, str):
        return converted_epsrel
 
    def _interpolated_function(x_value):
        if math.isnan(x_value) or math.isinf(x_value):
            raise ValueError("Integration points must be finite.")
        index = bisect_left(xs, x_value)
        if index == 0:
            return ys[0]
        if index >= len(xs):
            return ys[-1]
        x0, x1 = xs[index - 1], xs[index]
        y0, y1 = ys[index - 1], ys[index]
        if x1 == x0:
            raise ValueError("Interpolation failed due to duplicate x values.")
        weight = (x_value - x0) / (x1 - x0)
        return y0 + weight * (y1 - y0)
 
    try:
        result, _ = scipy_quad(
            _interpolated_function,
            converted_a,
            converted_b,
            epsabs=converted_epsabs,
            epsrel=converted_epsrel,
        )
    except Exception as exc:  # noqa: BLE001 - repackage SciPy errors
        return f"scipy.integrate.quad error: {exc}"
 
    if not isinstance(result, (int, float)):
        return "scipy.integrate.quad error: non-numeric result."
    numeric_result = float(result)
    if math.isnan(numeric_result) or math.isinf(numeric_result):
        return "scipy.integrate.quad error: result is not finite."
    return numeric_result