%typeset_mode False
from scipy import special as sp
from scipy import integrate
from sage.symbolic.integration.integral import indefinite_integral
from sage.symbolic.integration.integral import definite_integral
import numpy as np
import math

# The error for numerical integration is 0.0034

numerical_integral(9.096e12 * x ** 58, -1, 1)

# Definite integral and then applying limits

sage_integral = definite_integral(9.096e12 * x**58, x, -1, 1)
print(sage_integral)

poly_x = np.poly1d([1, 0])
poly_y = np.poly1d([1, 0])

for i in range(0,57):
    poly_x = poly_x * poly_y

print(poly_x)

poly1d_integral = np.poly1d.integ(9.096e12 * poly_x)(1) - np.poly1d.integ(9.096e12 * poly_x)(-1)
print(poly1d_integral - numerical_integral(9.096e12 * x ** 58, -1, 1)[0])

print((sage_integral- poly1d_integral) * 1e14)