f(x)=1/(x^2+1)
xs=0:0.5:4
N=length(xs)
M=zeros(N,N)
for i=1:N 
    M[i,1]=f(xs[i])
end
for j=2:N
    for i=1:N-j+1
        M[i,j]=M[i+1,j-1]-M[i,j-1]
    end
end

function p(alpha,d)
    b=1.0
    s=0
    for j=1:d+1
        s+=b*M[1,j]
        b*=(alpha-j+1)/j
    end
    return s
end