from numpy import *
import sympy as sym
from scipy.linalg import hilbert

n = 10

# Numerically
Pnum = eye(n+1)
altsign = matrix((-1)**(arange(0, n+1)))
Nnum = hilbert(n+1) + array(transpose(altsign)*altsign)*hilbert(n+1)
for k in range(1,n+1):
    Pnum[:k,k] = linalg.solve( -Nnum[:k,:k], Nnum[:k,k] )

# Symbolically
hilbertmatr = sym.zeros(n+1)
for i in range(n+1):
    for j in range(n+1):
        hilbertmatr[i,j] = sym.Rational(1, i+j+1)
P = sym.eye(n+1)
altsign = sym.zeros(1,n+1)
for k in range(n+1):
    altsign[0,k] = (-1)**k

newmatr = altsign.T*altsign
toadd = newmatr.multiply_elementwise(hilbertmatr)
N = hilbertmatr + toadd
for k in range(1,n+1):
    P[:k,k] = -N[:k,:k].inv()*N[:k,k]