⏰️

クロックパズルを単因子論で解く

に公開

これは「Haskell Advent Calendar 2025」12日目の記事です。


タイトルにクロックパズルと銘打ちましたが、このパズルの本当の名前を私は知りません(知ってたら教えて下さい)。パズルというのは例えば、以下のように

左、中央、右にそれぞれ3時間と4時間と5時間を計る時計があり、時計の針は 1,2,3 を指しているとします。それぞれの押しボタンを押すと時計の針は以下のように動きます。

  • 押しボタン a を押すと、左の時計の針が2時間進み、中央の時計の針が1時間進む
  • 押しボタン b を押すと、中央の時計の針が2時間進み、右の時計の針が3時間進む
  • 押しボタン c を押すと、左の時計の針が1時間進み、右の時計の針が2時間進む

この時、初期状態から始めて押しボタンを押し、時計の針が全て0を指すように揃えてください、というようなものです。

一般に、それぞれの時計が何時間計れるかは揃っててもいいしバラバラでも構いません、また時計と押しボタンの数が一致していなくても構いません。どの押しボタンにどの時計が対応しているのかも基本的には自由です(もちろん解けるかどうかは設定によります)。

このようなパズルは実際のゲームでも

  • ソニックフロンティア - トゥームストーンの謎(ネタバレ動画
  • 杖と剣の伝説 - 石碑の遺跡(ネタバレ動画
  • (FFXでもあったような気がする)

と見かけたのですが、他にも色々あると思います(知ってたら教えて下さい)。古くはライツアウトと呼ばれるパズルがあり、これも2時間の時計が升目上に並んだ同様のパズルと考えることができます。またルービッククロックも同様のパズルだと教えてもらいました。

牛刀割鶏な方法にはなりますが本稿ではこのパズルを単因子論を使って解いていきます。

解き方

上記例のパズルに焦点を当てて考えますが、以下の議論は一般の場合にもそのまま拡張可能です。

まず時計の状態は、0時を基準に時計の針が左右に何時間分回転したかを数えることで(時計回りなら+1、反時計回りなら-1)整数の三つ組

{\mathbb Z}\times{\mathbb Z}\times{\mathbb Z} = {\mathbb Z}^3

の元として表現できます。また文字盤の数字に注目すれば巡回群の直積

{\mathbb Z}/3{\mathbb Z}\times{\mathbb Z}/4{\mathbb Z}\times{\mathbb Z}/5{\mathbb Z}

の元だと考えることもできます。これらの間にはそれぞれの時計で計れる時間で割った余りを考えるという自然な全射

\pi: {\mathbb Z}^3\to{\mathbb Z}/3{\mathbb Z}\times{\mathbb Z}/4{\mathbb Z}\times{\mathbb Z}/5{\mathbb Z}

が存在します。例えば時計の針が全て0を指す状態は

(\overline{0},\overline{0},\overline{0})\in{\mathbb Z}/3{\mathbb Z}\times{\mathbb Z}/4{\mathbb Z}\times{\mathbb Z}/5{\mathbb Z}

を表します。

初期状態が

(\overline{y_1},\overline{y_2},\overline{y_3})\in{\mathbb Z}/3{\mathbb Z}\times{\mathbb Z}/4{\mathbb Z}\times{\mathbb Z}/5{\mathbb Z}

だったとしましょう。パズルを解くために必要なのは押しボタンを押すことで初期状態の逆元 (-\overline{y_1},-\overline{y_2},-\overline{y_3}) を作り出すことです。いま初期状態は (\overline{1},\overline{2},\overline{3}) なので逆元は (-\overline{1},-\overline{2},-\overline{3})、すなわち (\overline{2},\overline{2},\overline{2}) となります。これを {\mathbb Z}^3 の世界に引き戻すと集合 \pi^{-1}((\overline{2},\overline{2},\overline{2})) のいずれかの元を作ればよいということになります。

押しボタンaを押すことは {\mathbb Z}^3 の世界で({\mathbb Z}^3 の元を縦ベクトルで表すことにする)

\begin{pmatrix} 2 \\ 1 \\ 0 \end{pmatrix}

を、押しボタンbを押すことは

\begin{pmatrix} 0 \\ 2 \\ 3 \end{pmatrix}

を、押しボタンcを押すことは

\begin{pmatrix} 1 \\ 0 \\ 2 \end{pmatrix}

を時計の状態に足すことを表します。これらの縦ベクトルを並べた行列を

A = \begin{pmatrix} 2 & 0 & 1 \\ 1 & 2 & 0 \\ 0 & 3 & 2 \end{pmatrix}

と置きましょう。

押しボタンa,b,cを押した回数を並べたベクトル x = (a, b, c)^{\sf T} を考えると、押しボタンを押した結果は

Ax

となります。

パズルを解くためには

Ax\in\pi^{-1}((\overline{2},\overline{2},\overline{2}))

つまり

\begin{pmatrix} 2 & 0 & 1 \\ 1 & 2 & 0 \\ 0 & 3 & 2 \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} 2 \\ 2 \\ 2 \end{pmatrix} + \begin{pmatrix} 3 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 5 \end{pmatrix} \begin{pmatrix} n \\ m \\ l \end{pmatrix}

を満たすような a, b, c, n, m, l を求める必要があります。

\~{A} = \begin{pmatrix} 2 & 0 & 1 & 3 & 0 & 0 \\ 1 & 2 & 0 & 0 & 4 & 0 \\ 0 & 3 & 2 & 0 & 0 & 5 \end{pmatrix},\ \ w = \begin{pmatrix} a \\ b \\ c \\ -n \\ -m \\ -l \end{pmatrix},\ \ b = \begin{pmatrix} 2 \\ 2 \\ 2 \end{pmatrix}

と置くと、元の式は

\~{A}w=b

となり、このような線形方程式系を解けば良いことがわかりました。

線形方程式系を解くだけなら係数行列を掃き出し法で階段形に整理できればよいのでエルミート標準形に帰着するのが素直な方法です。しかし本稿ではパズルの構造をもう少し詳しく調べるために単因子論、具体的にはスミス標準形を用いた解き方を考えることにします。

単因子論

単因子論は単因子環 R 上の行列のスミス標準形と単因子によって R 上の行列や有限生成加群の構造を調べる理論です。主イデアル整域(PID)ならば単因子環なので具体例として整数 {\mathbb Z} が含まれます。今回必要なのは整数行列のスミス標準形だけなので、整数 {\mathbb Z} を前提に話を進めます。

まずは整数行列の基本変形について見ていきましょう。

【整数行列の基本変形】
<左基本変形>
(左1)2つの行を入れ換える
(左2)ある行に -1 を掛ける
(左3)ある行に整数 z を掛けて、他の行に加える
<右基本変形>
(右1)2つの列を入れ換える
(右2)ある列に -1 を掛ける
(右3)ある列に整数 z を掛けて、他の列に加える

これらの基本変形は行列を掛けることで実現できます。例えば

\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}

という3次行列を行列の左から掛けると2行目と3行目を入れ換える基本変形(左1)を実現します。反対にこの行列を右から掛けると2列目と3列目を入れ換える基本変形(右1)を実現します。この行列は交換を表すことからも分かるように自分自身が逆行列です。

次に

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

を行列の左から掛けると2行目に -1 を掛ける基本変形(左2)を実現します。反対にこの行列を右から掛けると2列目に -1 を掛ける基本変形(右2)を実現します。この行列もまた自分自身が逆行列となります。

最後に

\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & z \\ 0 & 0 & 1 \end{pmatrix}

を行列の左から掛けると2行目に3行目を z 倍して加える基本変形(左3)を実現します。反対にこの行列を右から掛けると3列目に2列目を z 倍して加える基本変形(右3)を実現します。この行列の逆行列は

\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & -z \\ 0 & 0 & 1 \end{pmatrix}

となります。

基本変形を表現する行列を基本行列と呼びます。可逆な整数行列、すなわち行列式が \pm 1 となる整数行列をユニモジュラ行列といいます。基本行列は可逆なのでユニモジュラ行列であり、反対に任意のユニモジュラ行列は基本行列の積で表すことができます(後で証明します)。

行列 A に基本変形を繰り返して行列 B が得られる時、行列 AB対等であると言います。定義より行列 AB が対等である時、あるユニモジュラ行列 U, V を用いて

B = UAV

と書くことができます。この対等という関係は同値関係です。

さて、準備ができたのでスミス標準形の存在性を証明しましょう。

【スミス標準形の存在性】
整数を要素に持つ任意のm\times n行列 A は次のような対角成分にのみ値を持つ m\times n行列(標準形)と対等である

\begin{pmatrix} d_1 & \\ & d_2 \\ & & \ddots \\ & & & d_r \\ & & & & 0 \\ & & & & & \ddots \\ & & & & & & 0 \\ \end{pmatrix}

d_1,\dots,d_r は0でない整数であり、d_{i-1}d_i を割り切る(2\leq i\leq r

証明
もし A が零行列 O ならA 自身が標準形である。

A\neq O と対等な行列 D
(1) d_{11}\neq 0
(2) d_{i1}=d_{1j}=0\ (i>1,j>1)
(3) d_{11}d_{ij}\ (i>1,j>1) を割り切る
という以下のような行列

D = \begin{pmatrix} d_{11} & 0 & \cdots & 0\\ 0 & d_{22} & \cdots & d_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ 0 & d_{m2} & \cdots & d_{mn} \end{pmatrix}

の存在を証明できれば、同様の議論を右下のブロック行列を対象として再帰的に繰り返すことにより有限回の操作で標準形にたどり着き証明が完了する。

A の0でない絶対値の最も小さな要素が a_{ij} ならば第1行と第i行、第1列と第j列を入れ替えた行列 B を考えることで、A と対等で b_{11}=a_{ij} を満たす行列が得られる。

② もし b_{11} が 第1行の全ての要素を割り切るならば、第1列を -b_{1j}/b_{11} 倍して第j(j>1) に加えることで (1,j) 成分が全て0となる対等な行列が得られる。同様にもし b_{11} が第1列の全ての要素を割り切るならば、第1行を -b_{i1}/b_{11} 倍して第i(i>1) に加えることで (i,1) 成分が全て0となる対等な行列が得られる。

③ 反対に第1行の中に b_{11} で割り切れない要素があるとき、すなわち (1,j) 成分において

b_{1j} = qb_{11}+r,\ \ 0<r<|b_{11}|

となるなら、第1列を -q 倍して第j列に加えることで (1,j) 成分が r となる対等な行列が得られる。そしてこの行列を A と見做して①から議論をやり直す。第1列の中に b_{11} で割り切れない要素があるときも同様。0<r<|b_{11}| であるためこの議論は必ず有限回で停止し、②と合わせると最終的に (1,1) 成分以外の第1行及び第1列の要素が0となるような行列 C が得られる。

c_{ij}\ (i>1,j>1)c_{11} で割り切れない要素が存在した場合、第i行を第1行にそのまま(1倍して)加えた行列を B と見做して③から議論をやり直す。この議論が繰り返されるたびに |c_{11}| は真に小さくなるので、この議論は必ず有限回で停止し、(i,j) 成分 (i>1,j>1) が全て (1,1) 成分で割り切れるような行列 D が得られる。

この行列 D は性質(1)(2)(3)を満たすので、あとは同様の議論を再帰的に繰り返せば良い。


実は標準形は対角成分の単元倍(-1倍)の違いを除いて一意に定まることも証明できます[1]。ある行列 A に対等で一意に存在するこの標準形を スミス標準形 と呼び、0でない対角成分 d_1,\dots,d_r を行列 A単因子 と呼びます。

ユニモジュラ行列 A のスミス標準形 S を考えると AS は対等なので基本行列の積(ユニモジュラ行列)で表せる行列 U,V が存在して S=UAV と書けます。ユニモジュラ行列の行列式は \pm 1 なので

\begin{matrix} \operatorname{det}(S) &=& \operatorname{det}(U)\operatorname{det}(A)\operatorname{det}(V) \\ &=& \pm 1\cdot\pm 1\cdot\pm 1 \\ &=& \pm 1 \end{matrix}

となります。S は対角行列なので対角成分は全て \pm 1 でありそのような行列は基本行列の積で表せるので、A=U^{-1}SV^{-1} を考えればユニモジュラ行列 A 自身も基本行列の積で表せることがわかりました。

実装

パズルの話に戻りましょう。今考えている3次行列ならスミス標準形を手で求めることも簡単ですが、一般の場合になると大変なのでプログラムを用いて計算できるようにしましょう。

実装には Haskell を用います。型レベル自然数を用いれば行と列の操作を混同した処理を書いてしまったとしてもコンパイル時に検出されるので安心です。

プラグマの宣言とライブラリのインポート
{-# LANGUAGE GHC2024 #-}
{-# LANGUAGE BlockArguments #-}

import Control.Monad (forM_, when)
import Data.Function (fix)
import Data.Maybe (fromMaybe, listToMaybe)
import Data.Proxy (Proxy(..))
import GHC.TypeLits (KnownNat, natVal)

import Control.Monad.State
import Data.Finite (Finite, packFinite)
import Data.Vector.Sized (Vector)
import qualified Data.Vector.Sized as V

まずは純粋な線形代数に関する関数の実装です。ベクトルや行列の型にはサイズを型レベル自然数で扱えるよう vector-sizedVector を使います。

-- | 内積
dot :: Num a => Vector n a -> Vector n a -> a
dot = (V.sum .) . V.zipWith (*)

-- | 行列
type Matrix n m a = Vector n (Vector m a)

-- | 単位行列
eye :: KnownNat n => Matrix n n Int
eye = V.generate \i -> V.generate \j -> if i == j then 1 else 0

-- | 行列の転置
transpose :: KnownNat n => Matrix m n a -> Matrix n m a
transpose = sequenceA

-- | 行列積
(!*!) :: (KnownNat r, Num a) => Matrix m n a -> Matrix n r a -> Matrix m r a
a !*! b = fmap (flip fmap (transpose b) . dot) a

もっと色々な線形代数に関する関数の実装に興味がある方は以下の記事もご参照ください。

https://zenn.dev/lotz/articles/6b0d8081ad2f8a

次にスミス標準形の計算に欠かせない整数行列の基本変形を実装しましょう。基本行列の積としてではなく基本変形に対応した操作として実装します。

-- | 整数行列の基本変形(左1)
swapRows :: (KnownNat n, KnownNat m) => Finite n -> Finite n -> Matrix n m Int -> Matrix n m Int
swapRows i j mat
  | i == j = mat
  | otherwise = V.generate \r ->
      if r == i then mat `V.index` j
      else if r == j then mat `V.index` i
      else mat `V.index` r

-- | 整数行列の基本変形(右1)
swapCols :: (KnownNat n, KnownNat m) => Finite m -> Finite m -> Matrix n m Int -> Matrix n m Int
swapCols i j mat
  | i == j = mat
  | otherwise = flip V.map mat \row ->
      V.generate \r ->
        if r == i then row `V.index` j
        else if r == j then row `V.index` i
        else row `V.index` r

-- | 整数行列の基本変形(左2)
negateRow :: (KnownNat n, KnownNat m) => Finite n -> Matrix n m Int -> Matrix n m Int
negateRow i = V.imap (\r row -> if r == i then V.map negate row else row)

-- | 整数行列の基本変形(右2)
negateCol :: (KnownNat n, KnownNat m) => Finite m -> Matrix n m Int -> Matrix n m Int
negateCol j mat = flip V.map mat \row -> flip V.imap row (\c z -> if c == j then negate z else z) 

-- | 整数行列の基本変形(左3)
addRow :: (KnownNat n, KnownNat m) => Finite n -> Finite n -> Int -> Matrix n m Int -> Matrix n m Int
addRow i j z mat = V.generate \r ->
  if r == i
    then V.zipWith (\x y -> x + z * y) (mat `V.index` i) (mat `V.index` j)
    else mat `V.index` r

-- | 整数行列の基本変形(右3)
addCol :: (KnownNat n, KnownNat m) => Finite m -> Finite m -> Int -> Matrix n m Int -> Matrix n m Int
addCol i j z mat = flip V.map mat \row ->
    V.generate \c ->
      if c == i
        then (row `V.index` i) + z * (row `V.index` j)
        else row `V.index` c

行列 A のスミス標準形 S と対等を与えるユニモジュラ行列 U,V

S = UAV

の3つ組 (S,U,V) を表すデータ構造を定義しましょう。

-- | 整数行列のスミス標準形
-- S = UAV となるS,U,Vを表す
data SNFState n m = SNFState
  { snfS :: Matrix n m Int
  , snfU :: Matrix n n Int
  , snfV :: Matrix m m Int
  } deriving Show

スミス標準形を計算するアルゴリズムの中では SNFState を頻繁に参照したり更新したりを行います。なので State モナドを用いた管理を採用し、必要なヘルパー関数も用意しておきます。

-- | スミス標準形の計算に用いるモナド
type SNFMonad n m a = State (SNFState n m) a

getS :: SNFMonad n m (Matrix n m Int)
getS = gets snfS

getVal :: (KnownNat n, KnownNat m) => Finite n -> Finite m -> SNFMonad n m Int
getVal i j = do
  s <- getS
  pure $ (s `V.index` i) `V.index` j

modifyS :: (Matrix n m Int -> Matrix n m Int) -> SNFMonad n m ()
modifyS f = modify \s -> s { snfS = f (snfS s) }

modifyU :: (Matrix n n Int -> Matrix n n Int) -> SNFMonad n m ()
modifyU f = modify \s -> s { snfU = f (snfU s) }

modifyV :: (Matrix m m Int -> Matrix m m Int) -> SNFMonad n m ()
modifyV f = modify \s -> s { snfV = f (snfV s) }

SNFState に対する整数行列の基本変形を定義します。例えば行列 S の左側から基本行列 P を掛けると

PS = PUAV

となるので行列 S に対する基本変形は行列 U あるいは V に対する基本変形も同時に引き起こします。

-- | `SNFState` に対する整数行列の基本変形
swapRowsM :: (KnownNat n, KnownNat m) => Finite n -> Finite n -> SNFMonad n m ()
swapRowsM i j = do
  modifyS $ swapRows i j
  modifyU $ swapRows i j

swapColsM :: (KnownNat n, KnownNat m) => Finite m -> Finite m -> SNFMonad n m ()
swapColsM i j = do
  modifyS $ swapCols i j
  modifyV $ swapCols i j

negateRowM :: (KnownNat n, KnownNat m) => Finite n -> SNFMonad n m ()
negateRowM i = do
  modifyS $ negateRow i
  modifyU $ negateRow i

negateColM :: (KnownNat n, KnownNat m) => Finite m -> SNFMonad n m ()
negateColM j = do
  modifyS $ negateCol j
  modifyV $ negateCol j

addRowM :: (KnownNat n, KnownNat m) => Finite n -> Finite n -> Int -> SNFMonad n m ()
addRowM i j m = do
  modifyS $ addRow i j m
  modifyU $ addRow i j m

addColM :: (KnownNat n, KnownNat m) => Finite m -> Finite m -> Int -> SNFMonad n m ()
addColM i j m = do
  modifyS $ addCol i j m
  modifyV $ addCol i j m

いよいよスミス標準形を求めるアルゴリズムを実装しましょう。少し長いですが、そもそもスミス標準形の存在性の証明が構成的であり、基本的にはそのまま実装されているので証明と対応付けながらだと読みやすいと思います。単因子は単元倍の自由度があるので最後に -1 を掛けて0以上の値に揃えているところだけ追加しています。

-- | スミス標準形を求めるアルゴリズム
smithNormalForm :: forall n m. (KnownNat n, KnownNat m) => Matrix n m Int -> SNFState n m
smithNormalForm mat0 = execState algorithm initialState
  where
    initialState = SNFState mat0 eye eye
    rows = fromIntegral $ natVal (Proxy :: Proxy n)
    cols = fromIntegral $ natVal (Proxy :: Proxy m)
    limit = min rows cols

    algorithm :: SNFMonad n m ()
    algorithm = do
      forM_ [0 .. limit - 1] \k -> do
        processPivot k
        ensurePositiveDiagonal k

    toFinite :: KnownNat k => Int -> Finite k
    toFinite x = fromMaybe (error "Index out of bounds") $ packFinite (fromIntegral x)

    processPivot :: Int -> SNFMonad n m ()
    processPivot k = do
      fix \restart -> do
        -- ① 最小の非ゼロ要素を探して (k,k) に移動
        maybeMin <- findMinNonZero k
        case maybeMin of
          Nothing -> pure () -- 残りが零行列ならこの処理は完了
          Just (r, c) -> do
            let k_row = toFinite @n k
                k_col = toFinite @m k
            swapRowsM k_row r
            swapColsM k_col c

            -- ② 行方向
            forM_ [k + 1 .. rows - 1] \i -> do
              let i_row = toFinite @n i
              val <- getVal i_row k_col
              pivot <- getVal k_row k_col
              when (val /= 0) do
                let q = val `div` pivot
                addRowM i_row k_row (-q)
                rem' <- getVal i_row k_col
                when (rem'/= 0) do
                  swapRowsM k_row i_row
                  restart

            -- ② 列方向
            forM_ [k + 1 .. cols - 1] \j -> do
              let j_col = toFinite @m j
              val <- getVal k_row j_col
              pivot <- getVal k_row k_col
              when (val /= 0) do
                let q = val `div` pivot
                addColM j_col k_col (-q)
                rem' <- getVal k_row j_col
                when (rem' /= 0) do
                  swapColsM k_col j_col
                  restart

            -- ③
            pivot <- getVal k_row k_col
            badIndices <- findBadIndices k pivot
            case badIndices of
              Just (_, c') -> do
                addColM k_col c' 1
                restart
              Nothing -> pure ()

    -- (k, k)成分から右下のブロック行列における0でない絶対値が最小の要素の位置を探す関数
    findMinNonZero :: Int -> SNFMonad n m (Maybe (Finite n, Finite m))
    findMinNonZero k = do
      s <- getS
      let candidates = [ (abs val, (r, c))
                       | r <- map (toFinite @n) [k .. rows - 1]
                       , c <- map (toFinite @m) [k .. cols - 1]
                       , let val = (s `V.index` r) `V.index` c
                       , val /= 0
                       ]
      if null candidates
        then pure Nothing
        else pure $ Just (snd $ minimum candidates)

    -- (k, k)成分から右下のブロック行列において与えられた整数で割り切れない要素の位置を探す関数
    findBadIndices :: Int -> Int -> SNFMonad n m (Maybe (Finite n, Finite m))
    findBadIndices k pivot = do
      s <- getS
      let indices = [ (r, c)
                    | r <- map (toFinite @n) [k + 1 .. rows - 1]
                    , c <- map (toFinite @m) [k + 1 .. cols - 1]
                    , ((s `V.index` r) `V.index` c) `mod` pivot /= 0
                    ]
      pure $ listToMaybe indices

    -- 対角成分を正にする
    ensurePositiveDiagonal :: Int -> SNFMonad n m ()
    ensurePositiveDiagonal k = do
      let k_row = toFinite @n k
          k_col = toFinite @m k
      val <- getVal k_row k_col
      when (val < 0) $ negateRowM k_row

fix を使ったループの実装が気になる方は以下の記事もご参照ください

https://qiita.com/lotz/items/409112751a23689d7d85

https://qiita.com/lotz/items/0894079a44e87dc8b73e

解き方(続き)

スミス標準形を計算する関数が実装できたので早速パズルを解くために使いましょう。我々は

\~{A} = \begin{pmatrix} 2 & 0 & 1 & 3 & 0 & 0 \\ 1 & 2 & 0 & 0 & 4 & 0 \\ 0 & 3 & 2 & 0 & 0 & 5 \end{pmatrix},\ \ w = \begin{pmatrix} a \\ b \\ c \\ -n \\ -m \\ -l \end{pmatrix},\ \ b = \begin{pmatrix} 2 \\ 2 \\ 2 \end{pmatrix}

として

\~{A}w=b

という線形方程式系を考えていたのでした。この係数行列 \~{A} に対してスミス標準形を計算してみましょう。

main :: IO ()
main = do
  let matList :: [[Int]]
      matList = [ [2, 0, 1, 3, 0, 0]
                , [1, 2, 0, 0, 4, 0]
                , [0, 3, 2, 0, 0, 5]]
      maybeMat :: Maybe (Matrix 3 6 Int)
      maybeMat = V.fromListN =<< (sequence $ map V.fromListN matList)

  case maybeMat of
    Nothing -> putStrLn "行列生成エラー"
    Just matA -> do
      putStrLn "Original Matrix A:"
      V.mapM_ print matA
      putStrLn "--------------------"

      let res = smithNormalForm matA

      putStrLn "Smith Normal Form S:"
      V.mapM_ print $ snfS res
      putStrLn "--------------------"

      putStrLn "Left Transform U:"
      V.mapM_ print $ snfU res
      putStrLn "--------------------"

      putStrLn "Right Transform V:"
      V.mapM_ print $ snfV res
      putStrLn "--------------------"

      putStrLn "Verification (U * A * V):"
      let s_calc = snfU res !*! matA !*! snfV res
      V.mapM_ print s_calc
      putStrLn $ "Match: " ++ show (snfS res == s_calc)
      putStrLn "--------------------"

出力

$ cabal run
Original Matrix A:
Vector [2,0,1,3,0,0]
Vector [1,2,0,0,4,0]
Vector [0,3,2,0,0,5]
--------------------
Smith Normal Form S:
Vector [1,0,0,0,0,0]
Vector [0,1,0,0,0,0]
Vector [0,0,1,0,0,0]
--------------------
Left Transform U:
Vector [1,0,0]
Vector [0,1,0]
Vector [-2,4,1]
--------------------
Right Transform V:
Vector [0,1,0,0,-4,-2]
Vector [0,0,0,0,0,1]
Vector [1,-2,3,-15,-40,-29]
Vector [0,0,-1,5,16,11]
Vector [0,0,0,0,1,0]
Vector [0,0,-1,6,16,11]
--------------------
Verification (U * A * V):
Vector [1,0,0,0,0,0]
Vector [0,1,0,0,0,0]
Vector [0,0,1,0,0,0]
Match: True
--------------------

上記より

S = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ \end{pmatrix}

\~{A} のスミス標準形であり、

U = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -2 & 4 & 1 \\ \end{pmatrix},\ V = \begin{pmatrix} 0 & 1 & 0 & 0 & -4 & -2 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & -2 & 3 & -15 & -40 & -29 \\ 0 & 0 & -1 & 5 & 16 & 11 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & -1 & 6 & 16 & 11 \\ \end{pmatrix}

として

S = U\~{A}V

という関係が成り立つことがわかりました。

元の線形方程式系を変形すると

\begin{matrix} && \~{A}w &=& b \\ &\iff& U\~{A}VV^{-1}w &=& Ub \\ &\iff& SV^{-1}w &=& Ub \end{matrix}

となり、 Ub

Ub= \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ -2 & 4 & 1 \end{pmatrix} \begin{pmatrix} 2\\ 2\\ 2 \end{pmatrix} = \begin{pmatrix} 2 \\ 2 \\ 6 \end{pmatrix}

となります。

S は対角行列なので S の対角成分がそれぞれ対応する Ub の要素を割り切らなければそもそも整数解は存在しません。つまり \~{A} の単因子が Ub の成分を割り切るかどうかでパズルが解けるかどうかがわかるのです。今回のように単因子が全て1となるように設計すればどのような初期値を与えても解けるパズルになります。

計算を続けましょう。

SV^{-1}w = Ub

より

V^{-1}w = \begin{pmatrix} 2 \\ 2 \\ 6 \\ p \\ q \\ r \end{pmatrix}

となります。但し p, q, r は任意の整数です。

VV^{-1}w=w を計算すると

\begin{pmatrix} 0 & 1 & 0 & 0 & -4 & 2 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & -2 & 3 & -15 & -40 & -29 \\ 0 & 0 & -1 & 5 & 16 & 11 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & -1 & 6 & 16 & 11 \\ \end{pmatrix} \begin{pmatrix} 2 \\ 2 \\ 6 \\ p \\ q \\ r \\ \end{pmatrix} = \begin{pmatrix} 2-4q+2r \\ r \\ 16-15p-40q-29r \\ -6+5p+16q+11r \\ q \\ -6+6p+16q+11r \\ \end{pmatrix}

となります。

w = \begin{pmatrix} a \\ b \\ c \\ -n \\ -m \\ -l \end{pmatrix}

だったので p,q,r に好きな数字を入れて a,b,c が0以上の整数になるように調整すると p=1,q=0,r=0 のとき a=2,b=0,c=1 となります。よってパズルの答えは

  • 押しボタンaを2回押す
  • 押しボタンbを0回押す
  • 押しボタンcを1回押す

であることがわかりました。

クロックパズルが解けましたね👏

おわりに

本稿では特定のパズルの解き方を具体的に考察しましたが、一般的な問題設定の場合でも問題なく同様の方法を適用することが可能です。
本稿で実装したスミス標準形を計算するアルゴリズムは証明に即した単純な方法ではありますが、計算の中間段階で現れる数値の桁数が入力行列の成分の桁数に対して指数関数的に増大する係数爆発という現象を引き起こすことも知られています。こういった計算の困難さを回避するために Kannan-Bachemアルゴリズム を始めとする様々な改良アルゴリズムが研究されています。
みなさんもクロックパズルとみなせるような連動して動く形のパズルを見かけたらぜひ行列のスミス標準形を使って解くことを試してみてください!


\読んでいただきありがとうございました!/
この記事が面白かったら いいね♡ をいただけると嬉しいです☺️
バッジを贈っていただければ次の記事を書くため励みになります🙌

脚注
  1. 例えば「代数学基礎 B 講義ノート (環論)」を参照 ↩︎

Discussion