#%% imports
import numpy as np
from numpy.polynomial import polynomial as P
#%% integrand & grenser

a = 0       # start
b = 1       # stopp

mu = 0      # middelverdi
sg = 1      # standardavvik

def f(x):      #    normal distribution
    return np.exp(-(((x-mu)/sg)**2)/2) / (sg*((2*np.pi)**(1/2)))


#%% riemannsum

def RS(a,b,f,n):
    #    startverdi for oppsamlet areal = 0
    A_hoyre = 0
    A_venstre = 0
    A_midt = 0
    A_trapes = 0
    
    h = (b - a)/n

    for i in range(n):
        A_hoyre += h * f(a + i*h)       #    hoyre riemannsum
        A_venstre += h * f(a + (i+1)*h)      #    venstre riemannsum
        A_midt += h * f(a + (2*i+1)*h/2)      #    midtpunktsummer
        A_trapes += h/2 * (f(a+i*h) + f(a+(i+1)*h)) #    trapesmetoden
    
    return [A_hoyre, A_venstre, A_midt, A_trapes]


#%% kvadratur

def kvadratur(xn):
    def l_coeff(xi,xn):
        r = []
        s = 1
        for i in range(len(xn)):
            if xn[i] != xi:
                r.append(xn[i])
                s *= (xi-xn[i])
        p = P.polyfromroots(r)
        return [x/s for x in p]

    def integrate_polynomial(p_coeff,a,b):
        hv = [0]
        hv2 = 0
        for i in range(len(p_coeff)):
            hv.append(p_coeff[i]/(i+1))
        for i in range(len(hv)):
            hv2 += (b**i - a**i)*hv[i]
        return hv2
    
    hv = 0
    for i in range(len(xn)):
        hv += f(xn[i])*integrate_polynomial(l_coeff(xn[i],xn),a,b)
    return hv



## Newton-Cotes-kvadratur
def NC(a,b,f,n=2):
    # Punkter jevnt fordelt
    xn = []
    for i in range(n):          
        xn.append(a + i*(b-a)/(n-1))
    
    return kvadratur(xn)
    
    
## Clenshaw-Curtis-kvadratur
def CS(a,b,f,n=2):    
    # https://en.wikipedia.org/wiki/Chebyshev_nodes
    xn = []
    for i in range(n):          # lager punktene
        xn.append((a+b)/2 + (b-a)*np.cos((2*(i+1)-1)*np.pi/(2*n))/2)
    
    return kvadratur(xn)


## Gauss-Legendre-kvadratur
def GL(a,b,f,n=2):
    # Approksimasjon. Se https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial
    xn = []
    for i in range(n):
        xn.append((1 - 1/(8*n**2) + 1/(8*n**3)) * np.cos((4*(i+1)-1)*np.pi/(4*n+2)))
    
    return kvadratur(xn)


#%% kall paa funksjoner

print('\nP(0<Z<1) har løsning:  0.34134474606854294')

RSn = 100000        # Beste verdi: hoyest mulig
print("\n====== Riemannsummer ======")
print(f'Valgte parametre: start={a}, stopp={b}, n={RSn}')
print(f'Hoyre riemannsum: {RS(a,b,f,RSn)[0]}')
print(f'Venstre riemannsum: {RS(a,b,f,RSn)[1]}')
print(f'Midtpunktmetoden: {RS(a,b,f,RSn)[2]}')
print(f'Trapesmetoden: {RS(a,b,f,RSn)[3]}')
print('Feil (trapes):',abs(0.34134474606854294-RS(a,b,f,RSn)[3]))


NCn = 10            # Beste verdi: 10
print('\n====== Newton-Cotes ======')
print(f'Valgte parametre: start={a}, stopp={b}, n={NCn}')
print(f'Areal N-C: {NC(a,b,f,NCn)}')
print('Feil:',abs(0.34134474606854294-NC(a,b,f,NCn)))

CSn = 9             # Beste verdi: 9
print('\n====== Clenshaw-Curtis ======')
print(f'Valgte parametre: start={a}, stopp={b}, n={CSn}')
print(f'Areal C-S: {CS(a,b,f,CSn)}')
print('Feil:',abs(0.34134474606854294-CS(a,b,f,CSn)))

GLn = 15            # Beste verdi: 15
print('\n====== Gauss-Legendre ======')
print(f'Valgte parametre: start={a}, stopp={b}, n={GLn}')
print(f'Areal N-C: {GL(a,b,f,GLn)}')
print('Feil:',abs(0.34134474606854294-GL(a,b,f,GLn)))