# Simpson's 1/3 Rule implementation in Python
import numpy as np
 
def simpsons_one_third_rule(func, a, b, n):
    """
    Approximate the integral of a function using Simpson's 1/3 Rule.
 
    Parameters:
        func (function): The function to integrate.
        a (float): The lower limit of integration.
        b (float): The upper limit of integration.
        n (int): The number of sub-intervals (must be even).
 
    Returns:
        float: Approximation of the integral.
    """
    # Check if n is even
    if n % 2 != 0:
        raise ValueError("Number of sub-intervals (n) must be even.")
 
    h = (b - a) / n  # Step size
    x = np.linspace(a, b, n + 1)  # Generate x values
    y = func(x)  # Evaluate the function at x values
 
    # Apply Simpson's 1/3 Rule
    integral = y[0] + y[-1]  # First and last terms
    integral += 4 * sum(y[1:n:2])  # Odd terms
    integral += 2 * sum(y[2:n-1:2])  # Even terms
 
    return (h / 3) * integral
 
# Example usage
if __name__ == "__main__":
    # Define the function to integrate
    def f(x):
        return np.sin(x)  # Example: Integrating sin(x)
 
    # Define limits and number of intervals
    a = 0  # Lower limit
    b = np.pi  # Upper limit
    n = 6  # Number of sub-intervals (must be even)
 
    # Perform integration
    try:
        result = simpsons_one_third_rule(f, a, b, n)
        print(f"The approximate value of the integral is: {result}")
    except ValueError as e:
        print(e)
 

IyBTaW1wc29uJ3MgMS8zIFJ1bGUgaW1wbGVtZW50YXRpb24gaW4gUHl0aG9uCmltcG9ydCBudW1weSBhcyBucAoKZGVmIHNpbXBzb25zX29uZV90aGlyZF9ydWxlKGZ1bmMsIGEsIGIsIG4pOgogICAgIiIiCiAgICBBcHByb3hpbWF0ZSB0aGUgaW50ZWdyYWwgb2YgYSBmdW5jdGlvbiB1c2luZyBTaW1wc29uJ3MgMS8zIFJ1bGUuCgogICAgUGFyYW1ldGVyczoKICAgICAgICBmdW5jIChmdW5jdGlvbik6IFRoZSBmdW5jdGlvbiB0byBpbnRlZ3JhdGUuCiAgICAgICAgYSAoZmxvYXQpOiBUaGUgbG93ZXIgbGltaXQgb2YgaW50ZWdyYXRpb24uCiAgICAgICAgYiAoZmxvYXQpOiBUaGUgdXBwZXIgbGltaXQgb2YgaW50ZWdyYXRpb24uCiAgICAgICAgbiAoaW50KTogVGhlIG51bWJlciBvZiBzdWItaW50ZXJ2YWxzIChtdXN0IGJlIGV2ZW4pLgoKICAgIFJldHVybnM6CiAgICAgICAgZmxvYXQ6IEFwcHJveGltYXRpb24gb2YgdGhlIGludGVncmFsLgogICAgIiIiCiAgICAjIENoZWNrIGlmIG4gaXMgZXZlbgogICAgaWYgbiAlIDIgIT0gMDoKICAgICAgICByYWlzZSBWYWx1ZUVycm9yKCJOdW1iZXIgb2Ygc3ViLWludGVydmFscyAobikgbXVzdCBiZSBldmVuLiIpCgogICAgaCA9IChiIC0gYSkgLyBuICAjIFN0ZXAgc2l6ZQogICAgeCA9IG5wLmxpbnNwYWNlKGEsIGIsIG4gKyAxKSAgIyBHZW5lcmF0ZSB4IHZhbHVlcwogICAgeSA9IGZ1bmMoeCkgICMgRXZhbHVhdGUgdGhlIGZ1bmN0aW9uIGF0IHggdmFsdWVzCgogICAgIyBBcHBseSBTaW1wc29uJ3MgMS8zIFJ1bGUKICAgIGludGVncmFsID0geVswXSArIHlbLTFdICAjIEZpcnN0IGFuZCBsYXN0IHRlcm1zCiAgICBpbnRlZ3JhbCArPSA0ICogc3VtKHlbMTpuOjJdKSAgIyBPZGQgdGVybXMKICAgIGludGVncmFsICs9IDIgKiBzdW0oeVsyOm4tMToyXSkgICMgRXZlbiB0ZXJtcwoKICAgIHJldHVybiAoaCAvIDMpICogaW50ZWdyYWwKCiMgRXhhbXBsZSB1c2FnZQppZiBfX25hbWVfXyA9PSAiX19tYWluX18iOgogICAgIyBEZWZpbmUgdGhlIGZ1bmN0aW9uIHRvIGludGVncmF0ZQogICAgZGVmIGYoeCk6CiAgICAgICAgcmV0dXJuIG5wLnNpbih4KSAgIyBFeGFtcGxlOiBJbnRlZ3JhdGluZyBzaW4oeCkKCiAgICAjIERlZmluZSBsaW1pdHMgYW5kIG51bWJlciBvZiBpbnRlcnZhbHMKICAgIGEgPSAwICAjIExvd2VyIGxpbWl0CiAgICBiID0gbnAucGkgICMgVXBwZXIgbGltaXQKICAgIG4gPSA2ICAjIE51bWJlciBvZiBzdWItaW50ZXJ2YWxzIChtdXN0IGJlIGV2ZW4pCgogICAgIyBQZXJmb3JtIGludGVncmF0aW9uCiAgICB0cnk6CiAgICAgICAgcmVzdWx0ID0gc2ltcHNvbnNfb25lX3RoaXJkX3J1bGUoZiwgYSwgYiwgbikKICAgICAgICBwcmludChmIlRoZSBhcHByb3hpbWF0ZSB2YWx1ZSBvZiB0aGUgaW50ZWdyYWwgaXM6IHtyZXN1bHR9IikKICAgIGV4Y2VwdCBWYWx1ZUVycm9yIGFzIGU6CiAgICAgICAgcHJpbnQoZSkK