- N桁のパスワードを数える問題を、再帰・メモ化・2配列DP・行列累乗の順に実装し、それぞれの速度を比較しました。
- 隣り合う桁の更新規則が線形変換であることに着目すると、遷移を9×9の行列で表せます。
- 行列の累乗を二分累乗で計算すると、N=1,000,000でも約20回の行列積で済み、2配列DPの約900万回から大幅に削減できます。
- 状態数が小さくNが大きい問題では、行列累乗DPが数百倍の速度差を生み出します。
1. 問題
AtCoder の「1111gal password」を Common Lisp で解きながら、アルゴリズムによって計算時間がどう変わると見てみました。
N桁のパスワードで、各桁が1〜9、かつ隣り合う桁の差が1以下になるものを数えます。
答えは 998244353 で割った余りです1。
制約は 2 ≤ N ≤ 1,000,000 です。
たとえば N=2 のとき、11, 12, 21, 22, 23 … のように隣の差が0か1であればよく、13 や 19 はNGです。
2. まず再帰で書く(TLE 7)
この問題を、1〜9の数直線上をN回移動する経路の数を求める問題と、読み替えて考えました。
各ステップでは、「左に1」「そのまま」「右に1」の3種類の移動ができますが、1と9は端でそれ以上外には出られません。
このように見ると、再帰の構造が自然に出てきます。
「末尾桁が d で長さが N のパスワードの数」を rec(N, d) とすると、長さ N のパスワードは、長さ N-1 のパスワードの末尾に1桁追加したものです。
つまり、末尾桁 d に来られるのは、直前が d-1、d、d+1 のいずれかだったときなので、
rec(N, d) = rec(N-1, d-1) + rec(N-1, d) + rec(N-1, d+1)
基底ケースとなる N =1桁の場合では、どの桁も1通りです。
ただし、端の処理は d=1 なら d-1 への移動がなく、d=9 なら d+1 への移動がありません。
コードにするとこうなります。
(defparameter *prime* 998244353)
(defun satisfied-patterns/rec (N)
(labels ((rec (N d)
(mod (cond
((= N 0) 0)
((= N 1) 1)
((= d 1) (+ (rec (1- N) 1)
(rec (1- N) 2)))
((= d 9) (+ (rec (1- N) 8)
(rec (1- N) 9)))
(t (+ (rec (1- N) (1- d))
(rec (1- N) d)
(rec (1- N) (1+ d)))))
*prime*)))
(loop for digit from 1 to 9
sum (rec N digit))))Code language: Lisp (lisp)
2.1. メイン関数
実際に動かしてみます。
(defun main/rec ()
(let ((N (read)))
(princ (satisfied-patterns/rec N))
(terpri)))
#-swank (main/rec)Code language: Lisp (lisp)
コードは簡潔なのですが、同じ (N, d) の組み合わせを何度も計算するため、N が少し大きくなるだけで指数的に遅くなり2、N=1,000,000 ではスタックオーバーフローになります。

3. ハッシュテーブルでメモ化(1514 ms)
再帰の再計算を避けるため、dp[N][d] をハッシュテーブルに記録してみました。
これは「トップダウンDP」の典型的な手法で、関数の返り値をメモしておくだけでかんたんに実装でき、再帰の形をそのまま残しながら重複計算を避けられます3。
キーは n * 10 + d の整数にすれば重複しません。n が最大 10^6、d が 1〜9 なので、この組み合わせは一意になります。
ハッシュテーブル memo は、関数呼び出しのたびにローカルに定義されますが、結果を返すまでは保持され、再帰関数 rec からクロージャで直接参照できます。
(defun satisfied-patterns/memo (N)
(let ((memo (make-hash-table)))
(labels ((rec (n d)
(cond
((= n 1) 1)
(t (let ((key (+ (* n 10) d)))
(or (gethash key memo)
(setf (gethash key memo)
(mod (+ (if (> d 1) (rec (1- n) (1- d)) 0)
(rec (1- n) d)
(if (< d 9) (rec (1- n) (1+ d)) 0))
*prime*))))))))
(mod (loop for d from 1 to 9
sum (rec N d))
*prime*))))Code language: Lisp (lisp)
再計算はなくなりますが、N=1,000,000 の場合、エントリ数が約 9,000,000 になります。
(defun main/memo ()
(let ((N (read)))
(princ (satisfied-patterns/memo N))))
#-swank (main/memo)Code language: Lisp (lisp)
メモリは、約720MB。
メモリ制限の 1024 MiB までギリギリです。
4. 2配列DPに切り替える(197 ms)
単純なメモ化では、全エントリを保持するので、空間計算量を圧迫してしまいます。
そこで、2配列DPに変えてみました。
遷移をよく観察すると、長さ i の結果は、長さ i-1 だけで計算できます。
なので、順番に i を上がっていくなら、直前の状態を保持しておけば、それ以前のメモは不要です。
dp_next[d] = dp_prev[d-1] + dp_prev[d] + dp_prev[d+1]
これは、ボトムアップDP4の典型的な実装です。
prev と now の2本の配列を用意して、長さを1ずつ伸ばしていきます。
配列は :element-type 'fixnum で特化配列として確保します5。
(defparameter *prime* 998244353)
(defun satisfied-highest (N d prev now)
(setf (aref now d)
(mod
(cond
((= N 0) 0)
((= N 1) 1)
((= d 1) (+ (aref prev 1)
(aref prev 2)))
((= d 9) (+ (aref prev 8)
(aref prev 9)))
(t (+ (aref prev (1- d))
(aref prev d)
(aref prev (1+ d)))))
*prime*)))
(defun satisfied-patterns (N)
(let ((now (make-array 10 :element-type 'fixnum :initial-element 0))
(prev (make-array 10 :element-type 'fixnum :initial-element 0)))
(loop for i from 1 to (1- N)
do (loop for d from 1 to 9
do (satisfied-highest i d prev now))
do (loop for d from 1 to 9
do (setf (aref prev d) (aref now d))))
(mod
(loop for d from 1 to 9
sum (satisfied-highest N d prev now))
*prime*)))Code language: Lisp (lisp)
計算量は O(9N)、約900万回の更新ですが、積み上げ式にしたことで、計算途中が減り 197 ms になりました。
メモリ消費量も 20MBまで削減できました。

(defun main ()
(let ((N (read)))
(princ (satisfied-patterns N))))
#-swank (main)Code language: Lisp (lisp)
4.1. 2配列をフリップで少し効率化(161 ms)
毎回、prev から now にコピーするのが無駄な気がしたので、odd-now と even-now の2配列を N の偶奇で切り替えることにしました。
初期値は odd-now[1]〜[9] = 1 にセットします。
また、prev[0] = 0、prev[10] = 0 に固定しておけば、d=1 のとき d=9 のときの比較演算も不要になります。
(defun satisfied-highest/effi (N d odd-now even-now)
(let* ((now (if (evenp N) even-now odd-now))
(prev (if (evenp N) odd-now even-now)))
(setf (aref now d)
(mod (+ (aref prev (1- d))
(aref prev d)
(aref prev (1+ d)))
*prime*))))
(defun satisfied-patterns/effi (N)
(let ((even-now (make-array 11 :element-type 'fixnum :initial-element 0))
(odd-now (make-array 11 :element-type 'fixnum :initial-element 0)))
(loop for d from 1 to 9
do (setf (aref odd-now d) 1))
(loop for i from 2 to (1- N)
do (loop for d from 1 to 9
do (satisfied-highest/effi i d odd-now even-now)))
(mod
(loop for d from 1 to 9
sum (satisfied-highest/effi N d odd-now even-now))
*prime*)))Code language: Lisp (lisp)
実行時間は最大 161 ms でした。
コピーのコストを省くことで、18.3% 速くなりました。
ただし、コードは読みにくくなり、N の偶奇で odd-now と even-now を切り替える設計は境界条件で混乱しやすいです。
いわば「最後の手段」の最適化で、これぐらいならもとのボトムアップDPで十分だったかもしれません。
5. 線形変換をベクトルと行列の積で表す(1047 ms)
ある程度、解答ができたところで、ほかの参加者の解答を見ていると、1ms台の提出があることに気づきました。
- 解説 – AtCoder Beginner Contest 242
- 提出 #65946762 – AtCoder Beginner Contest 242
- 提出 #47604007 – AtCoder Beginner Contest 242
この差は、アルゴリズムの細かな無駄を削る「定数倍最適化」ではなく、解法の根本を変えることで計算量のオーダーを改善しているはずです。
こういうときに役立つのが「行列(matrix)」です。
5.1. 線形変換は行列で表せる
「ある値の集まりを、別の値の集まりに変換する操作のうち、各出力が入力のスカラー倍の和で書けるもの」を「線形変換」といいます。
2配列DPの更新式はまさに線形変換で、各出力 dp_next[d] は 入力 dp_prev の隣接する成分の和になっています。
dp_next[d] = dp_prev[d-1] + dp_prev[d] + dp_prev[d+1]
線形変換は、「基底」を選ぶことで必ず行列の掛け算として表せます6。
行列 Mat を使うと、更新式は成分ごとに掛け算で定義できます。
dp_next[j] = Σ_i dp_prev[i] × Mat[i][j]
つまり Mat[i][j] は「桁 i から桁 j への寄与の重み」を表しています。
今回の問題では、隣り合う桁への遷移だけを許すので、|i - j| ≤ 1 のとき Mat[i][j] = 1、それ以外は 0 です。
Mat を具体的に書くとこうなります。
1 1 0 0 0 0 0 0 0
1 1 1 0 0 0 0 0 0
0 1 1 1 0 0 0 0 0
0 0 1 1 1 0 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 0 1 1 1 0 0
0 0 0 0 0 1 1 1 0
0 0 0 0 0 0 1 1 1
0 0 0 0 0 0 0 1 1
この行列 Mat を使うと、dp_prev から dp_next の更新は、このように書けるわけです。
dp_next = dp_prev × Mat
この遷移行列は、以下のコードで作ります。
(defun make-matrix (size)
(make-array (list size size)
:element-type 'fixnum
:initial-element 0))
(defun make-transition-matrix (size)
(let ((Mat (make-matrix size)))
(loop for i below size do
(loop for j below size do
(setf (aref Mat i j)
(if (<= (abs (- i j)) 1) 1 0))))
Mat))Code language: Lisp (lisp)
(<= (abs (- i j)) 1)の部分が、問題の条件に相当しています。
5.2. ベクトルと行列の掛け算
状態の遷移は、ベクトルと行列の掛け算で計算できます。
(defun multiply-vector-matrix (size dp Mat)
(let ((tmp (load-time-value
(make-array 9 :element-type 'fixnum
:initial-element 0))))
(loop for j below size do
(setf (aref tmp j)
(mod (loop for i below size
sum (* (aref dp i) (aref Mat i j)))
*prime*)))
(loop for j below size do
(setf (aref dp j) (aref tmp j)))))Code language: Lisp (lisp)
いったん、tmp[j] = Σ dp[i] * Mat[i][j] の計算して、それを最後に dp[j]に書き込みます。
これが、更新式に対応しているわけです。
5.3. まず素直に実装する
行列で更新を表せたので、N-1 回 Mat を掛けていく実装を書いてみます。
(defun satisfied-patterns/mat-naive (N)
(let* ((size 9)
(dp (make-array size
:element-type 'fixnum
:initial-element 1))
(Mat (make-transition-matrix size)))
(loop repeat (1- N) do
(multiply-vector-matrix size dp Mat))
(mod (loop for d below size
sum (aref dp d)) *prime*)))
(defun main/mat-naive ()
(let ((N (read)))
(princ (satisfied-patterns/mat-naive N))))
#-swank (main/mat-naive)Code language: Lisp (lisp)
初期ベクトル dp を [1,1,...,1] から始めて、Mat を1回ずつ掛けて更新します。
これを計測すると、実行時間は 1,047msになりました。

ここまでなら、本質的にはさっきのDPと変わりません。
むしろ、ループ回数は同じ N-1 回、各ステップで 9×9 の積和を計算しているので、2配列DPより重くなります。
行列式が効果を発揮するのは、「累乗」の性質を利用できることです。
6. 行列の累乗を先に計算する
長さ1の初期状態は、各桁がそれぞれ1通りなので dp_1 = [1,1,1,1,1,1,1,1,1] です。
dp_N は、これに Mat を書けていくことで求めました。
そのため、dp_N = dp_1 × Mat^(N-1) と表すことができます。
dp_2 = dp_1 × Mat
dp_3 = dp_2 × Mat = dp_1 × Mat × Mat = dp_1 × Mat²
dp_N = dp_1 × Mat^(N-1)
ここで重要なのは、行列の掛け算は結合則を満たすことです7。
(dp_1 × Mat) × Mat = dp_1 × (Mat × Mat)
つまり、dp_1 に Mat を1回ずつ掛けていく代わりに、Mat どうしを先に掛け合わせて Mat^(N-1) を作ってから、最後に dp_1 と合わせることができます。
というのも、Mat は各桁の更新規則に対応したものなので、N の値に関わらず変わりません。
だから Mat^(N-1) だけを先に計算しておくことが可能です。
このような形にすると、「二分累乗」というアルゴリズムで効率化できます。
6.1. 二分累乗で O(log N) に落とす
Mat^(N-1) をどう計算するか、という問題ですが、素直に Mat を N-1 回掛けると結局 O(N) で変わりません。
ここで使うのが「二分累乗」です。
指数を2進数で分解すると、必要な累乗だけを掛け合わせて目的の累乗が作れます。
二分累乗は、通常の数値の累乗で使われるアルゴリズムですが、行列積でも利用できます。
たとえば、Matの13乗は、以下の掛け算に分解できます。
Mat^13 = Mat^8 × Mat^4 × Mat^1 (13 = 1101₂)
Mat^1 から始めて Mat^2、Mat^4、Mat^8 と順に二乗していけば、各ステップは前の結果を使い回せます。
N=13 なら、行列積は4回で済みます。
必要な行列積の回数は log₂(N) 程度になります。
つまり、N=1,000,000 でも約20回です。
行列積1回の計算量は 9³ = 729 回の積和なので、全体では:
9³ × log₂(N) ≒ 729 × 20 ≒ 15,000 回
2配列DPでは約900万回の計算がありましたが、約1.5万回と数百倍の差が期待できます。
そこで、行列の累乗の準備をしていきます。
6.2. 一時的な行列から結果をコピーする
また、行列積の結果を安全に書き込むためには、一時的な行列が必要です。
そこで、行列をコピーする copy-matrix を用意します。
(defun copy-matrix (size S D)
(loop for i below size do
(loop for j below size do
(setf (aref D i j) (aref S i j))))
D)Code language: Lisp (lisp)
setfと同じように、最後に D 自身を返しています。
6.3. 行列積
まず、行列 X と 行列 Y の積を計算する multiply-matrix 関数は、数学の定義をそのまま実装しています。
(defun multiply-matrix (size X Y D)
(let ((out (if (or (eq X D) (eq Y D))
(load-time-value (make-matrix 9))
D)))
(loop for i below size do
(loop for j below size do
(setf (aref out i j) 0)
(loop for k below size do
(setf (aref out i j)
(mod (+ (aref out i j)
(* (aref X i k)
(aref Y k j)))
*prime*)))))
(when (not (eq out D))
(copy-matrix size out D))
D))Code language: Lisp (lisp)
入力と出力に同じ配列を渡してそのまま計算すると、値が壊れてします。
そこで、出力する配列が入力にあるときには、いったん out に一時的に保存して copy-matrix で 書き戻しています。
multiply-matrix では積和の途中で毎回 mod を取っています8。
6.4. N回の行列の累乗を計算する
さて、ここまでの準備で、行列の累乗を実装できます。power-matrix/natural では、Mat に対して Matを power - 1回 掛けます。
(defun power-matrix/natural (size Mat power)
(let ((result (copy-matrix size Mat (make-matrix size)))
(times (1- power)))
(loop while (plusp times) do
(when (oddp times)
(multiply-matrix size result Mat result))
(multiply-matrix size Mat Mat Mat)
(setf times (ash times -1)))
result))Code language: Lisp (lisp)
(oddp times) で最下位ビットが 1 かどうかを判定し9、立っているときだけ result に掛け込みます。(ash power -1) は右に1ビットシフト10、つまり2で割る操作で、ビットを1つずつ処理していきます。
6.5. 結果を返す
最後に sum-matrix で全要素を足して答えにします。
初期ベクトルが全桁 1 なので、dp_1 × Mat^(N-1) の総和は Mat^(N-1) の全要素の和と等しくなります。
(defun satisfied-patterns/mat (N)
(let* ((size 9)
(Mat (make-transition-matrix size))
(mat-pow (power-matrix/natural size Mat (1- N)))
(dp (make-array size :element-type 'fixnum
:initial-element 1)))
(multiply-vector-matrix size dp mat-pow)
(mod (loop for d below size
sum (aref dp d)) *prime*)))Code language: Lisp (lisp)
6.6. 完成したコード(5 ms)
(defparameter *prime* 998244353)
(defun make-matrix (size)
(make-array (list size size)
:element-type 'fixnum
:initial-element 0))
(defun copy-matrix (size S D)
(loop for i below size do
(loop for j below size do
(setf (aref D i j) (aref S i j))))
D)
(defun multiply-matrix (size X Y D)
(let ((out (if (or (eq X D) (eq Y D))
(load-time-value (make-matrix 9))
D)))
(loop for i below size do
(loop for j below size do
(setf (aref out i j) 0)
(loop for k below size do
(setf (aref out i j)
(mod (+ (aref out i j)
(* (aref X i k)
(aref Y k j)))
*prime*)))))
(when (not (eq out D))
(copy-matrix size out D))
D))
(defun power-matrix/natural (size Mat power)
(let ((result (copy-matrix size Mat (make-matrix size)))
(tmp (load-time-value (make-matrix 9)))
(times (1- power)))
(loop while (plusp times) do
(when (oddp times)
(multiply-matrix size result Mat tmp)
(copy-matrix size tmp result))
(multiply-matrix size Mat Mat tmp)
(copy-matrix size tmp Mat)
(setf times (ash times -1)))
result))
(defun make-transition-matrix (size)
(let ((Mat (make-matrix size)))
(loop for i below size do
(loop for j below size do
(setf (aref Mat i j)
(if (<= (abs (- i j)) 1) 1 0))))
Mat))
(defun multiply-vector-matrix (size dp Mat)
(let ((tmp (load-time-value
(make-array 9 :element-type 'fixnum
:initial-element 0))))
(loop for j below size do
(setf (aref tmp j)
(mod (loop for i below size
sum (* (aref dp i) (aref Mat i j)))
*prime*)))
(loop for j below size do
(setf (aref dp j) (aref tmp j)))))
(defun satisfied-patterns/mat (N)
(let* ((size 9)
(Mat (make-transition-matrix size))
(mat-pow (power-matrix/natural size Mat (1- N)))
(dp (make-array size :element-type 'fixnum
:initial-element 1)))
(multiply-vector-matrix size dp mat-pow)
(mod (loop for d below size
sum (aref dp d)) *prime*)))
(defun main/mat ()
(let ((N (read)))
(princ (satisfied-patterns/mat N))))
#-swank (main/mat)Code language: Lisp (lisp)
7. 実装の比較
| 解法 | 時間計算量 | 空間計算量 | 最大実行時間 |
|---|---|---|---|
| 再帰 | O(3^N) | O(N) | TLE・スタック溢れ |
| メモ化再帰 | O(9N) | O(9N) | 1514 ms |
| 2配列DP | O(9N) | O(9) | 197 ms |
| 2配列DPフリップ | O(9N) | O(9) | 161 ms |
| 行列ナイーブ | O(9²N) | O(9²) | 1047 ms |
| 行列累乗DP | O(9³ log N) | O(9²) | 5 ms台 |
メモ化再帰は時間計算量こそ 2配列DPと同じ O(9N) ですが、空間計算量が O(9N) になります11。
N=1,000,000 では約900万エントリをハッシュテーブルに保持するため、アクセスのオーバーヘッドが積み重なって TLE になりました。
2配列DPは常に2本の配列しか持たないので、空間計算量は O(9) で一定です。
行列ナイーブは、行列という道具を導入したにもかかわらず 2配列DPより遅くなっています。
それは、1回の更新で 9×9 の積和を計算するためで、行列が速さを発揮するのはあくまで累乗との組み合わせによります。
速度面では行列累乗が圧倒しています。
ただ、Common Lispの標準関数には、行列演算がないので実装量は増えます。
7.1. 行列累乗が有効な条件
行列累乗が速いのは、今回のように「状態数が小さく、N が大きい」場合です。
状態数を S とすると計算量は O(S³ log N) なので、S が大きくなると逆に遅くなります。
今回は S=9 と非常に小さかったため、効果が出ました。
「1桁ずつ進めるDP」を「2桁分、4桁分、8桁分とまとめて進むDP」に変換する、というのがこの解法の核心です。
更新規則が毎回同じ線形変換であることに気づければ、行列累乗という道具が使えます12。
7.2. 【補足】行列の累乗で 0 を許容する場合
power-matrix/naturalでは、power>=1 を仮定していて、power-1回かけています。
これを、power>=0にするなら、「単位行列」からresult を計算し始める必要があります。
「単位行列」とは、何かに掛けても相手を変えない行列で、数の 1 に相当します13。
(defun make-identity-matrix (size)
(let ((Mat (make-matrix size)))
(loop for i below size do
(setf (aref Mat i i) 1))
Mat))Code language: Lisp (lisp)
その場合、power-matrixは、このように実装します。
(defun power-matrix (size Mat power)
(let ((result (make-identity-matrix size)))
(loop while (plusp power) do
(when (oddp power)
(multiply-matrix size result Mat result))
(multiply-matrix size Mat Mat Mat)
(setf power (ash power -1)))
result))Code language: Lisp (lisp)
このようにすると、powerが 0 の場合でも動作します。
- 998244353 は素数で、かつ 998244353 = 2²³ × 119 + 1 という形を持つため、数論変換(NTT)に使いやすい性質があります。競プロでは 10⁹+7 とともによく登場します。 – NTTのやさしい解説
- 各桁から最大3方向に分岐し、深さ N まで展開するので時間計算量は O(3^N) です。メモ化なしではN=20程度でも数秒かかります。
- 動的計画法には、大きな問題から再帰的に小さな問題へ分解するトップダウン方式(メモ化再帰)と、小さな問題から順に積み上げるボトムアップ方式の2種類があります。 – 動的計画法 – Wikipedia
- 小さな部分問題(長さ1)から順に大きな問題(長さN)へ積み上げるため、ボトムアップ方式と呼ばれます。再帰呼び出しのスタックコストがなく、必要な状態だけを保持できます。 – 動的計画法入門
:element-type 'fixnumを指定すると、SBCLは配列をタグなし整数の特化配列として確保します。SBCLの fixnum は63ビット符号付き整数で、most-positive-fixnumは 2⁶²-1 = 4611686018427387903 です。型宣言により aref のボックス化解除やバウンドチェック省略などの最適化が入ります。 – 競技プログラミングでCommon Lispを使っている人のために- 有限次元ベクトル空間上の線形変換は、基底を選ぶことで行列と一対一に対応します。行列の積は線形変換の合成に対応し、これが行列累乗を「繰り返し適用」に使える理由です。
- 行列積は一般に可換ではありません(AB ≠ BA)が、結合則 (AB)C = A(BC) は常に成立します。これは行列積が線形写像の合成に対応するためで、写像の合成が結合則を満たすことから従います。 – Associativity of Matrix Multiplication
- 998244353² ≈ 10¹⁸ は fixnum の上限(約 4.6×10¹⁸)以内ですが、9項を蓄積すると超過するリスクがあります。毎ステップで mod を取ることでオーバーフローを防いでいます。SBCLは fixnum を超えると自動的に bignum に昇格しますが、bignum 演算はヒープ確保を伴うため速度が落ちます。
- 整数が奇数であることと、2進数表現の最下位ビットが1であることは等価です。
(oddp n)は(= (mod n 2) 1)と同じ結果ですが、SBCLでは単純なビット検査にコンパイルされます。 ashは Arithmetic Shift の略で、符号ビットを保持したままビットシフトします。(ash n -1)はnを2で割った商(切り捨て)と等価です。C言語のn >> 1に相当します。- 空間計算量(space complexity)はアルゴリズムが使用するメモリ量を入力サイズの関数として表したものです。メモ化再帰はすべての (N,d) の組み合わせを保持するため O(9N) になります。再帰呼び出しのスタックも O(N) 分消費します。
- 行列累乗の最もよく知られた例はフィボナッチ数列です。F(n+1) = F(n) + F(n-1) という漸化式を2×2行列で表すと、F(N) を O(log N) で求められます。競プロでは「K回操作後の状態を求める」問題に広く使えます。 – 行列累乗についてのメモと出題例
- n×n 単位行列は対角成分がすべて1、それ以外が0の行列です。任意の行列 A に対して AI = IA = A が成立します。
make-identity-matrixでは(setf (aref Mat i i) 1)で対角成分だけを1にセットし、残りはmake-matrixのゼロ初期化に任せています。