import numpy as np
import quadrature
import polynomial_new as poly

""" Settings: """
# Create an instance of the quadrature rule:
x = np.array([-1.0, 0.0, 1.0])
w = np.array([1.0/3, 4.0/3, 1.0/3])
S = quadrature.QuadratureRule(x, w)

# Boundary values of interval:
a = -1.0
b = 1.0
# Number of test:
ntest = 10

""" Function to integrate polynomial: """
def int_poly(P, a, b):
    # Coefficients of the primitive function:
    pc = P.coeff / np.arange(1, P.deg+2, dtype=float)
    # Evaluate powers at both boundary points:
    bx = np.array([b**i for i in range(1, P.deg+2)])
    ax = np.array([a**i for i in range(1, P.deg+2)])

    return (pc @ bx - pc @ ax)

""" Main Test: """
for ii in range(ntest):
    # Generate random third order polynomial:
    #c = np.array([1.0, 0.0, 0.0, 1.0])
    c = np.random.randn(4)
    print("Test %d: coeffs: "%ii, c)
    P = poly.polynomial(c)
    # Compute integral:
    I = int_poly(P, a, b)
    # Approximation with Simpson rule:
    ISimp = S.integrate(P)
    print("Test %d: diff = "%ii, np.abs(I - ISimp))