RosettaCodeData/Task/Matrix-arithmetic/PARI-GP/matrix-arithmetic-3.pari

matperm(M)=
{
	my(n=matsize(M)[1],innerSums=vectorv(n));
	if(n==0, return(1));
	sum(x=1,2^n-1,
		my(k=valuation(x,2),s=M[,k+1],gray=bitxor(x, x>>1));
		if(bittest(gray,k),
			innerSums += s;
		,
			innerSums -= s;
		);
		(-1)^hammingweight(gray)*factorback(innerSums)
	)*(-1)^n;
}