Project Euler 101

http://projecteuler.net/index.php?section=problems&id=101

これは、線型方程式を解くだけなので簡単です。
しかし、これって必ず整数になるのでしょうか。重ね合わせの原理が効くので単項式だけ考えればよいのです。Vandermonde行列がでてきます。10乗までの計算では禁じた公式の係数は全て整数になるようです。

from fractions import Fraction

def inverse(A):
    def swap(m):
        for m1 in range(m + 1, n):
            if A[m1][m] != 0:
                A[m], A[m1] = A[m1], A[m]
                return
    
    def div(row, d):
        for k in range(n):
            row[k] /= d
    
    def sub(row, row0, m):
        for k in range(n):
            row[k] -= row0[k] * m
    
    n = len(A)
    B = [ [ 1 if i == j else 0 for j in range(n) ] for i in range(n) ]
    for m in range(n):
        if A[m][m] == 0:
            swap(m)
        a = A[m][m]
        div(A[m], a)
        div(B[m], a)
        for k in range(m + 1, n):
            a = A[k][m]
            sub(A[k], A[m], a)
            sub(B[k], B[m], a)
    
    for m in range(n - 2, -1, -1):
        for k in range(m + 1, n):
            sub(B[m], B[k], A[m][k])
    
    return B

def make_matrix(n):
    return [ [ Fraction(m ** e) for e in range(n) ]
                            for m in range(1, n + 1) ]

def mul(A, v):
    return [ sum(a[k] * v[k] for k in range(len(v))) for a in A ]

def BOP(f, n):
    A = [ [ Fraction(m ** e) for e in range(n) ]
                            for m in range(1, n + 1) ]
    B = inverse(A)
    v = mul(B, [ f(m) for m in range(1, n + 1) ])
    return sum(v[m] * (n + 1) ** m for m in range(n))

def f(n):
    return reduce(lambda x, y: 1 - n * x, xrange(N), 1)

N = 10
print sum(BOP(f, n) for n in range(1, N + 1))