from math import sqrt

# Increase the allowed recursion depth:
from sys import setrecursionlimit
setrecursionlimit(1000000)

# The adaptive integration function. g is the integrand, c/d the lower/upper
# lomtit and 未/蔚 are the absolute and relative accuracy.
def adapt(g, c, d, 未, 蔚):
    # The recursive subfunction:
    def adapt24(f, a, b, 未, 蔚, f2, f3, nrec=0):
        f1 = f(a+(b-a)/6)
        f4 = f(a+5*(b-a)/6)
        Q = (2*f1 + f2 + f3 + 2*f4)/6*(b-a)
        q = (f1 + f2 + f3 + f4)/4*(b-a)
        tol = 未 + 蔚*abs(Q)
        err = abs(Q-q)
        if err < tol:
            return (Q, nrec)
        else:
            (Q1, nrec) = adapt24(f, a, (a+b)/2, 未/sqrt(2), 蔚, f1, f2, nrec+1)
            (Q2, nrec) = adapt24(f, (a+b)/2, b, 未/sqrt(2), 蔚, f3, f4, nrec+1)
            return (Q1 + Q2, nrec)
    # Performing variable transformation in case of infinity limits:
    inf = float('inf')
    if c == -inf and d == inf:
        (a, b) = (-1, 1)
        f = lambda t: g(t/(1-t**2)) * (1+t**2)/(1-t**2)**2
    elif c == -inf:
        (a,b) = (0,1)
        f = lambda t: g(d - (1-t)/t) / t**2
        #(a, b) = (-1, 0)
        #f = lambda t: g(d - t/(1+t)) / -(1+t)**2
    elif d == inf:
        (a, b) = (0, 1)
        f = lambda t: g(c + t/(1-t)) / (1-t)**2
    else:
        (a, b) = (c, d)
        f = g
    # The integration begins:
    f2 = f(a+2*(b-a)/6)
    f3 = f(a+4*(b-a)/6)
    return adapt24(f, a, b, 未, 蔚, f2, f3)