real matrix matpow(real matrix a, real scalar n) {
	real matrix p, x
	real scalar i, s
	s = n<0
	n = abs(n)
	x = a
	p = I(rows(a))
	for (i=n; i>0; i=floor(i/2)) {
		if (mod(i,2)==1) p = p*x
		x = x*x
	}
	return(s?luinv(p):p)
}