繰り返し二乗法で行列のべき乗の高速化

【例題】 以下の行列の$10^8$乗を剰余 100007で求めよ

\large 
\begin{pmatrix} 
2 & 2 & -2 & 1\\
1 & 0 & 0 & 0\\ 
0 & 1 & 0 & 0\\ 
0 & 0 & 1 & 0
\end{pmatrix}

まずは素直にnumpyのmatrixを用いて書いてみますが、これでは$10^6$までが精いっぱい。

import numpy as np
N, M = 10**6, 100007
m = np.matrix([[2,2,-2,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]], dtype='object')
print(m**N % M)
#[[87688 87212 96356 67590]
# [67590 52515 52039 31522]
# [31522 4546 89478 15076]
# [15076 1370 74401 19623]]

行列の高速化としては対角化等の技法があるようですが、今回はNが大きいのと剰余を取るので繰り返し二乗法を使います。アルゴリズムは以下のような簡単なものです。

1. 与えられた行列を$2^k$がnを超える直前まで$k$回2乗を行う
2. 残った$n-2^k$に対して1.を繰り返す、これを残りが１になるまで続ける。

関数powmtxは引数に、行列mtx、べき乗の数n、剰余Mを与えると答えの行列を返します。

import numpy as np
def powmtx(mtx,n,M):   # mtx: matrix, n: power, M: modulo
  if n == 0: return np.eye(len(mtx), dtype='object')
  if n == 1: return mtx
  m2, p2 = mtx, 2
  while p2 < n:
    m2 = (m2*m2)%M
    p2 *= 2
  return powmtx(mtx,(n-p2//2),M)*m2 % M

このpowmtxを使って例題を解くと$10^8$も問題なく解くことが出来ました。どこまで出来るか試してみましたが$10^{100}$位までは楽勝でした。

N, M = 10**8, 100007
m = np.matrix([[2,2,-2,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]], dtype='object')

print(powmtx(m, N, M))
#[[81734 3148 15339 71658]
# [71658 38425 59846 58648]
# [58648 54369 21136 77135]
# [77135 4385 106 75399]]

（開発環境：Google Colab）