HammingP:=function(n,q)
 R<x>:=PolynomialRing(Rationals());
 v:=[x,1/2*((x-(q-2))*x-n*(q-1))];
 if n gt 2 then 
  for i in [2..n-1] do
   Append(~v,1/(i+1)*((x-i*(q-2))*v[i]-(n-i+1)*(q-1)*v[i-1]));
  end for;
 end if;
 v:=[1] cat v;
 return Matrix(Rationals(),n+1,n+1,
  [[Evaluate(v[j],n*(q-1)-q*i):j in [1..n+1]]:i in [0..n]]);
end function;