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