mata
function pow(a, n) {
	x = a
	for(p=1; n>0; n=floor(n/2)) {
		if(mod(n,2)==1) p = p*x
		x = x*x
	}
	return(p)
}
end