euler 101

euler101.py

def solve(polynomial):
    polymatrix = []
    for i in range(len(polynomial)):
        polymatrix.append([])
        for j in range(len(polynomial)):
            polymatrix[i].append(polynomial[i][0] ** j)
        polymatrix[i].append(polynomial[i][-1])
    for i in range(len(polymatrix)):
        divisor = polymatrix[i][i]
        for j in range(len(polymatrix[0])):
            polymatrix[i][j] /= 1.0 * divisor
        for j in range(1 + i, len(polymatrix)):
            divisor = polymatrix[j][i]
            for k in range(len(polymatrix[0])):
                polymatrix[j][k] -= divisor * polymatrix[i][k]
    for i in range(len(polymatrix) - 1, -1, -1):
        divisor = polymatrix[i][i]
        for j in range(i):
            polymatrix[i][j] /= 1.0 * divisor
        for j in range(i):
            divisor = polymatrix[j][i]
            for k in range(len(polymatrix[0])):
                polymatrix[j][k] -= divisor * polymatrix[i][k]
    answer = []
    for i in range(len(polymatrix)):
        answer.append(polymatrix[i][-1])
    return answer

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

boptotal = 0
polynomial = []
for i in range(1, 20):
    polynomial.append([i, u(i)])
    subtotal = 0
    for j in range(1, i + 2):
        subtotal = 0
        for k in range(len(polynomial)):
            subtotal += solve(polynomial)[k] * j ** k
    if int(subtotal) != int(u(i + 1)):
        boptotal += subtotal
print boptotal