import numpy as np


def rev(xs):
    return list(reversed(xs))


def bop(degree, f):

    if degree == 1:
        return 1
    A = [rev([n ** x for x in range(degree)])
         for n in range(1, degree + 1)]
    y = np.array([f(n) for n in range(1, degree + 1)])
    x = np.linalg.solve(A, y)
    ns = rev([(degree + 1) ** x for x in range(degree)])
    return np.dot(x, ns)


def euler_101():
    # max_degree = 4
    # def u(x):
        # return x * x * x

    max_degree = 10
    def u(n):
        return 1 - n + n**2 - n**3 + n**4 - n**5 + \
               n**6 - n**7 + n**8 - n**9 + n**10

    return int(sum([bop(d, u) for d in range(1, max_degree + 1)]))


if __name__ == "__main__":
    solution = euler_101()
    print("e101.py: " + str(solution))
    assert(solution == 37076114526)