#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 10 14:11:57 2024

@author: sam
"""

def qtrapn(f,a,b,N):
    h = (b-a)/N
    I = f(a) + f(b)
    for i in range(1,int(N)):
        xi = a + i*h
        I += 2*f(xi)
    I *= h/2
    return I

def qsimpn(f,a,b,N):
    h = (b-a)/N
    I = f(a) + f(b)
    for i in range(1,int(N)):
        if i % 2 == 0:
            xi = a + i*h
            I += 2*f(xi)
        elif i % 2 != 0:
            xi = a + i*h
            I += 4*f(xi)
    I *= h/3
    return I



def I(f,a,b,n):
    I = 0
    x,w = p_roots(n)
    for i in range(len(x)):
        I += w[i]*f((b-a)/2*x[i]+(a+b)/2)
    I *= (b-a)/2
    return I

def quadz(f,a,b,tol):
    Iold = 0
    for i in range(1,100):
        x,w = p_roots(i)
        I = (b-a)/2*sum(w*f((b-a)/2*x+(a+b)/2))
        if i > 3 and abs(I - Iold) < tol:
            return I
        Iold = I
    return I

def squad2d(f,a,b,c,d):
    def yint(f,x,c,d):
        def fy(y):
            return f(x,y)
        return I(fy,c,d,13)
    def fx(x):
        return yint(f,x,c,d)
    return I(fx,a,b,13)