連立1次ディオファントス方程式

2019/05/21
2019/05/31 更新
2019/06/17 更新(bug fix and shape up in PDF)
2020/11/19 追加
2020/11/22 追加
2020/11/26 プログラム更新と追加 (shape up)
2020/12/01 プログラム更新と追加 (shape up)
2020/12/07 リニューアル (shape up)

Abstract: Another solution to a system of linear Diophantine equations is presented.
The idea stands on the work of Gilbert and Pathria eliminating some defects
in their solution.

連立1次ディオファントス方程式

連立1次ディオファントス方程式とは

連立1次ディオファントス方程式、すなわち下記の式(1)
\[ \begin{equation} %\label{eq:bi:2} \left. \begin{array}{lcl} a_{1,1}x_1 + a_{1,2}x_2+\cdots + a_{1,n}x_n &=& b_1\\ a_{2,1}x_1 + a_{2,2}x_2+\cdots + a_{2,n}x_n &=& b_2\\ \cdots &=& \cdots\\ a_{m,1}x_1 + a_{m,2}x_2+\cdots + a_{m,n}x_n &=& b_m \end{array} \right\} \end{equation} \] において整数 $a_{i,j} \;(i=1,...,m;\; j=1,...,n)$ と整数 $b_i \;(i=1,...,m)$ を与え、整数解 $x_1,...,x_n$ を求める問題である。
例えば
\[ \begin{equation} %\label{eq:bi:2} \left. \begin{array}{rcr} 5x_1 + 6x_2 + 8x_3 &=& 1\\ 6x_1 - 11x_2 + 7x_3 &=& 9 \end{array} \right\} \end{equation} \] の解は
\[ \begin{equation} x_1 = 5 + 10k,\quad x_2=k,\quad x_3=-3 - 7k \end{equation} \] である。ここに $k$ は任意の整数である。

式(1)を連立一次ディオファントス方程式と言う。これの、効率の良い、しかも易しい解法を紹介する。
以下、これを文献[1]とする。

ダウンロード: http:lde.pdf

主なアイデアは Gilbert and Pathria の記事[2]から来ているが、改良が加えられている。
この解法は、筆者の知る限り、どの本にも載っていないし、またネット検索でも見つからない。
素朴に計算すると、嫌になるような計算が、巧みに避けられている。
1次ディオファントス方程式にチャレンジした経験のある方は是非とも読んでみてください。
新しい解法の威力が分かると思います。

読者としては、行列の初等的な知識を持っていることが想定されています。

文献

[1] 有澤健治:『連立1次ディオファントス方程式』
http://ar.nyx.link/lde/lde.pdf
[2] Willliam J. Gilbert, Anu Pathria: “Linear Diophantine Equations”
https://www.math.uwaterloo.ca/~wgilbert/Research/GilbertPathria.pdf 1990

文献[1]の正誤表

p.7
誤: また $\hat e_{i,j}=(\hat M)_{i,j}$ である。
正: また $\hat e_{i,j}=(\hat E)_{i,j}$ である。

アルゴリズムの要点

式(1)は次のように表現できる。
\[ \begin{equation} %\label{eq:bi:3} A= \begin{pmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n}\\ a_{2,1} & a_{2,2} & \cdots & a_{2,n}\\ \cdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{pmatrix},\quad {\boldsymbol x}=\begin{pmatrix} x_1\\ x_2\\ \dots\\ x_n \end{pmatrix},\quad {\boldsymbol b}=\begin{pmatrix} b_1\\ b_2\\ \dots\\ b_m \end{pmatrix},\quad A {\boldsymbol x} = {\boldsymbol b} \end{equation} \] これから次の表を作る。
\[ \begin{equation} \begin{array}{rrrr|r} a_{1,1} & a_{1,2} & \cdots & a_{1,n} & b_1\\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} & b_2\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} & b_m\\ \hline 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & 1 \end{array} \end{equation} \]

この表に対して次の変形規則を導入する。(これを基本変形と言う)
行の変形規則: 1行目から$m$行目について
(a) 任意の行$i$に、$i$ と異なる任意の行の整数倍を加えてよい
(b) 任意の行は符号を反転してもよい
(c) 任意の$2$つの行は交換してもよい
列の変形規則: 1列目から$n$列目について行と同様な操作が許される。

変形の目標は $A$ の部分を対角化することである。対角化すると次のようになるであろう。

\[ \begin{equation} \begin{array}{rrrr|r} \hat a_{1,1} & 0 & \cdots & 0 & \hat b_1\\ 0 & \hat a_{2,2} & \cdots & 0 & \hat b_2 \\ \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & 0 & \cdots & 0 & \hat b_m \\ \hline \hat e_{1,1} & \hat e_{1,2} & \cdots & \hat e_{1,n} \\ \hat e_{2,1} & \hat e_{2,2} & \cdots & \hat e_{2,n} \\ \cdots & \cdots & \cdots & \cdots \\ \hat e_{n,1} & \hat e_{n,2} & \cdots & \hat e_{n,n} \end{array} \end{equation} \] ここに $\hat a_{i,i},\hat b_{i}, \hat e_{i,j}$ などは整数である。
特に $\hat a_{i,i}$ は $i=1,...,k$ で正整数、$i=k+1,...,m$ で $0$ になるようにできる。
なお後に述べる単因子定理によれば、$\hat a_{i+1,i+1}$ は $\hat a_{i,i}$ で割り切れるようになるまで変形できるが、今の問題ではそこまで変形する必要はない。

対角化するために $A$ を次の形にもっていく。これを $A^{(1)}$としよう。
\[ \begin{equation} \begin{array}{ccccc} \hat a_{1,1} & 0 & 0 & \cdots & 0\\ 0 & * & * & \cdots & *\\ 0 & * & * & \cdots & *\\ 0 & * & * & \cdots & *\\ \end{array} \end{equation} \] ここに $*$ は何かの整数である。これができれば残った $*$ の部分を同様に処理すればよい。
文献[1]にはそのための手順が書かれているが、ここでは Maxima に適した方法を紹介する。
STEP1: $A$ の中の $0$ でない絶対値の最小の要素を $(1,1)$ の位置に移動する。
行と列の交換で達成できる。
STEP2: $1$ 列目の $(1,k) \;(k=2,..,m)$ を $0$ にすべく基本変形を繰り返す。
STEP3: $1$ 行目の $(k,1) \;(k=2,..,m)$ を $0$ にすべく基本変形を繰り返す。

以上のステップを繰り返すと式(7)の $*$ の部分が減少して行き、最終的には式(6)の形にまで行き着く。
文献[1]との違いは、文献[1]では行の処理を先に行っているが、ここでは列の処理を先に行っている。
また STEP1 では $A$ の全体を対象にしているところが違う。(この方が実行が早いし、プログラムも簡単になる)
正確には Maxima のプログラムコードを参照すべし。

\[ \begin{equation} \hat A= \begin{pmatrix} \hat a_{1,1} & 0 & \cdots & 0 \\ 0 & \hat a_{2,2} & \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & 0 \\ \end{pmatrix} ,\quad \hat{\boldsymbol b}= \begin{pmatrix} \hat b_1 \\ \hat b_2 \\ \cdots \\ \hat b_m \\ \end{pmatrix} ,\quad \hat E=\begin{pmatrix} \hat e_{1,1} & \hat e_{1,2} & \cdots & \hat e_{1,n} \\ \hat e_{2,1} & \hat e_{2,2} & \cdots & \hat e_{2,n} \\ \cdots & \cdots & \cdots & \cdots \\ \hat e_{n,1} & \hat e_{n,2} & \cdots & \hat e_{n,n} \end{pmatrix} \end{equation} \] と置くと次の方程式
\[ \begin{equation} \hat A \hat{\boldsymbol x}=\hat{\boldsymbol b} \end{equation} \] を解く問題に帰着する。$\hat{\boldsymbol x}$ についてのこの方程式は$\hat A$ が対角化されているので容易に求められる。${\boldsymbol x}$ は
\[ \begin{equation} {\boldsymbol x}=\hat E \hat{\boldsymbol x} \end{equation} \] で求められる。

\[ \begin{array}{rrr|r} 5 & 6 & 8 & 1 \\ 6 & -11 & 7 & 9 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|r} 5 & 6 & 8 & 1 \\ 1 & -17 & -1 & 8 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|r} 1 & -17 & -1 & 8 \\ 5 & 6 & 8 & 1 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|r} 1 & -17 & -1 & 8 \\ 0 & 91 & 13 & -39 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \\ \Longrightarrow \begin{array}{rrr|r} 1 & 0 & 0 & 8 \\ 0 & 91 & 13 & -39 \\ \hline 1 & 17 & 1 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|r} 1 & 0 & 0 & 8 \\ 0 & 13 & 91 & -39 \\ \hline 1 & 1 & 17 & \\ 0 & 0 & 1 & \\ 0 & 1 & 0 & \end{array} \Longrightarrow \begin{array}{rrr|r} 1 & 0 & 0 & 8 \\ 0 & 13 & 0 & -39 \\ \hline 1 & 1 & 10 & \\ 0 & 0 & 1 & \\ 0 & 1 & -7 & \end{array} \] 従って
\[ \begin{equation} \hat A=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 13 & 0 \end{pmatrix} ,\quad \hat{\boldsymbol b}=\begin{pmatrix} 8 \\ -39 \end{pmatrix} ,\quad \hat E=\begin{pmatrix} 1 & 1 & 10 \\ 0 & 0 & 1 \\ 0 & 1 & -7 \end{pmatrix} \end{equation} \] であり、式(9)と式(10)により
\[ \begin{equation} \hat{\boldsymbol x}=\begin{pmatrix} 8 \\ -3 \\ \hat x_3 \end{pmatrix} ,\quad {\boldsymbol b}=\hat E \hat{\boldsymbol b}=\begin{pmatrix} 5+10\hat x_3 \\ \hat x_3 \\ -3 - 7\hat x_3 \end{pmatrix} \end{equation} \] として式(3)の結果が得られる。ここに $\hat x_3$ は任意の整数であり、式(3)では $k$ になっている。

Maxima のプログラム

Maxima のプログラム
http:lde.mac
を紹介する。

実行例

Maxima の会話的実行は既知とする。
Maxima を立ち上げ

directory("");

で示されるディレクトリに "lde.mac" を置き

batchload("lde.mac");

を実行する。その下で、式(2)を例題として、式(3)を得る実行例を示す。ただし $k$ は $x_3$ に置き換わっている。

(%i3) M:matrix([5,6,8],[6,-11,7]); b:[1,9];
(%o3) matrix([5,6,8],[6,-11,7])
(%o4) [1,9]
(%i5) lde(M,b);
step0 0
     [matrix([1],[9]),matrix([5,6,8],[6,-11,7]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step1 1
     [matrix([1],[9]),matrix([5,6,8],[6,-11,7]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step2 1
     [matrix([1],[8]),matrix([5,6,8],[1,-17,-1]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step1 1
     [matrix([8],[1]),matrix([1,-17,-1],[5,6,8]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step2 1
     [matrix([8],[-39]),matrix([1,-17,-1],[0,91,13]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step1 1
     [matrix([8],[-39]),matrix([1,-17,-1],[0,91,13]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step3 1
     [matrix([8],[-39]),matrix([1,0,0],[0,91,13]),
      matrix([1,17,1],[0,1,0],[0,0,1])]
step1 1
     [matrix([8],[-39]),matrix([1,0,0],[0,91,13]),
      matrix([1,17,1],[0,1,0],[0,0,1])]
step1 2
     [matrix([8],[-39]),matrix([1,0,0],[0,13,91]),
      matrix([1,1,17],[0,0,1],[0,1,0])]
step3 2
     [matrix([8],[-39]),matrix([1,0,0],[0,13,0]),
      matrix([1,1,10],[0,0,1],[0,1,-7])]
step1 2
     [matrix([8],[-39]),matrix([1,0,0],[0,13,0]),
      matrix([1,1,10],[0,0,1],[0,1,-7])]
[matrix([8],[-39]),matrix([1,0,0],[0,13,0]),matrix([1,1,10],[0,0,1],[0,1,-7])]

(%o5) [10*x3+5,x3,(-7*x3)-3]
(%i6)

最後の出力

[10*x3+5,x3,(-7*x3)-3]

が式(3)に対応する。

プログラムは式(6)に至るまでの中間結果も出力する。
ステップごとの表は
\[
\begin{pmatrix}
A^{(l)}&{\boldsymbol b}^{(l)} \\
E^{(l)}&O
\end{pmatrix}
\]
の形をしている。$O$ は式(5)から式(6)に至る表では空欄であるが、それを形式的に $O$ で表している。
各ステップでは ${\boldsymbol b}^{(l)},A^{(l)},E^{(l)}$ がこの順でリスト形式で表示されている。
例えば

[matrix([8],[-39]),matrix([1,0,0],[0,13,0]),
matrix([1,1,10],[0,0,1],[0,1,-7])]

は式(11)を表している。

プログラムコード

lde.mac の主関数の一つ lde(M,b) の定義を次に示す。
列と行変形の基本的なアルゴリズムが正確に解るはずである。詳しくは lde.mac を直接参照されたい。

lde(M,b):=block([B,L,n,k],
  B:tab_gen2(M,b),
  if B = false then return("# ERROR,inconsistent b"),
  n:min(length(M),length(M[1])),
  dprint("step0",0,B),
  /* don't replace B[2] below with A, which behaves differently! */
  for k:1 while k <= n and chk_zero(B[2],k)=false do (
    X:step1(B,k),
    if X=false then (k:n+1, return(false)),
    B:X,
    dprint("step1",k,B),
    loop,
    unless chk_col(B[2],k) do (
      B:step2(B,k),
      dprint("step2",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    unless chk_row(B[2],k) do (
      B:step3(B,k),
      dprint("step3",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    if chk_col(B[2],k)=false then go(loop)
   ),
  print(B),
  L:lde_solution(B),
  if lde_final(L,M,b) then return(expand(L))
  else return("no solution")
);

B は3つの行列 S,A,T を要素とするリストである。
tab_gen2(M,b) によってそのようなリストが生成される。
その場合、 S はリスト b から生成される1列から成る行列、A は M のコピー、TはMを基にした単位行列(Mの行数と列数の小さい方を採用する)である。
S,A,T は基本変形によって変形を受ける。
その途中経過が

dprint("step2",k,B)

などによって示されているのである。
k は処理中の行番号(同じことであるが処理中の列番号)である。

lde(M,b) の内部コードは行列 $M$ の単因子標準形を求めるコードと多くを共有している。

SymPy のプログラム

2020-12-01

SymPy のプログラム
http:mlde.py
を紹介する。
名称が lde.py ではなく、mlde.py としたのは、モジュール化した時に中に含まれている関数 lde がモジュール名と一致していることが何かと不便なので、ファイル名に "m" を添えたのである。

実行例

Python の会話的実行は既知とする。
Python3 を立ち上げ

import os
os.getcwd()

で示されるディレクトリに "mlde.py" を置き

import mlde
from mlde import *

を実行する。その下で、式(2)を例題として、式(3)を得る実行例を示す。ただし $k$ は $x
_3$ に置き換わっている。

>>> M=Mat([5,6,8],[6,-11,7]); b=[1,9]
>>> lde(M,b)
step0 0 [Matrix([
[1],
[9]]), Matrix([
[5,   6, 8],
[6, -11, 7]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
step1 1 [Matrix([
[1],
[9]]), Matrix([
[5,   6, 8],
[6, -11, 7]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
中略
step1 2 [Matrix([
[  8],
[-39]]), Matrix([
[1,  0, 0],
[0, 13, 0]]), Matrix([
[1, 1, 10],
[0, 0,  1],
[0, 1, -7]])]
[10*x3 + 5, x3, -7*x3 - 3]
>>>

SymPy における行列は正式には Matrix であり、Maxima の matrix に比べて引数の与え方が少し複雑である。
つまり matrix で行列を表すには

matrix(行1,行2,...)

の形式で済むが Matrix では

Matrix([行1,行2,...])

の形式を要する。つまり "[ ]" が余分に必要である。
このように設計したのは、Matrix の使い方を柔軟にするためである。
しかし、他方ではコーディングが面倒で間違えやすくなるので、mlde では新たに関数 Mat を定義して、用法を matrix に合わせている。

Maxima と Python の違いの一つに、リストや行列のインデックスの開始番号がある。Maxima は 1 から始めるが Python は 0 から始める。 Maxima から Python への翻訳において、最も間違えやすい部分である。

プログラムコード

mlde.py の主関数の一つ lde(M,b) の定義を次に示す。
列と行変形の基本的なアルゴリズムが正確に解るはずである。詳しくは mlde.py を直接参照されたい。

def lde(M,b):
  B = tab_gen2(M,b)
  if B == False: error("# ERROR,inconsistent b")
  n = min(M.shape)
  dprint("step0",-1,B)
  # don't replace B[1] below with A, which behaves differently!
  for k in range(0,n):
    if chk_zero(B[1],k)==True: break
    X = step1(B,k)
    if X == False: break
    B = X
    dprint("step1",k,B)
    while True:
      while chk_col(B[1],k)==False:
        B = step2(B,k)
        dprint("step2",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      while chk_row(B[1],k)==False:
        B = step3(B,k)
        dprint("step3",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      if chk_col(B[1],k)==True: break
  L = lde_solution(B)
  if lde_final(L,M,b) : return list(expand(L))
  else: return "no solution"

dprint("step2",k,B)

などによって示されているのである。
k は処理中の行番号(同じことであるが処理中の列番号)である。

出力は Maxima 版と完全に同じになるように調整されている。

lde(M,b) の内部コードは行列 $M$ の単因子標準形を求めるコードと多くを共有している。

整数行列の単因子標準形

lde.mac には整数行列 $M$ の単因子標準形を求めるコードが含まれている。
多項式行列の単因子標準形も同じ考え方で求められるが、実際に必要になるコードは遥かに複雑であり、lde.mac ではサポートされていない。elmdiv.mac でサポートされている。

整数行列

整数を要素とする行列を整数行列と言う。正方行列でなくても構わない。
整数行列 $M$ が逆行列を整数行列の中に持つには、$M$ は正方行列で、行列式 $|M|$ が $\pm 1$ であることが必要十分条件である。逆行列を持つ整数行列は正則であると言う。

整数行列の単因子定理

全ての整数行列は、或る正則な整数行列 $S$ と $T$ によって次の形の整数行列 $\hat M=SMT$ に変換できる。
(a) $\hat M$ は対角行列で
(b) 対角要素は非負整数、
(c) 対角要素を $d_1,d_2,...$ とすれば、$d_{k+1}$ は $d_k$ の倍数である。
$0$ は $0$ の倍数であることに留意する。

従って対角要素は $d_1,d_2,...,0,0,...$ の形をとる。$0,0,...$ は発生しない場合もある。

$\hat M$ を $M$ の単因子標準形と言う。
例えば $ M=\begin{pmatrix}
5&6&8 \\
6&-11&7
\end{pmatrix} $ は $ S=\begin{pmatrix}-1&1\\6&-5\end{pmatrix},
\quad T=\begin{pmatrix}1&1&10\\0&0&1\\0&1&-7\end{pmatrix} $ によって単因子標準形 $ SMT=\hat M= \begin{pmatrix}1&0&0 \\ 0&13&0 \end{pmatrix} $ に変換される。

整数行列の基本変形

行の変形規則:
(a) 任意の行$i$に、$i$ と異なる任意の行の整数倍を加えてよい
(b) 任意の行は符号を反転してもよい
(c) 任意の$2$つの行は交換してもよい
列の変形規則: 行と同様な操作が許される。

行の変形規則と列の変形規則を合わせて整数行列の基本変形という。これらは可逆な変形規則であり、正則な整数行列で表現できる。

変形の目標は $M$ を単因子標準形になるよう対角化することである。最初の目標は\[
\begin{equation}
\left(
\begin{array}{ccccc}
d_1&0&0&\cdots&0\\
0&*&*&\cdots&*\\
0&*&*&\cdots&*\\
0&*&*&\cdots&*\\
\end{array}
\right)
\end{equation}
\]の形に持っていくことである。ここに $*$ の部分は全て $d_1$ の倍数とする。これに成功すれば $d_2,d_3,...$ と進んで行けばよい。ディオファントス方程式を解く場合には単に対角化すればよかったが単因子標準形の場合にはそれでは不十分である。

対角化アルゴリズム

以下では単因子標準形への変換行列も同時に求める方法を紹介する。すなわち $SMT$ が単因子標準形 $\hat M$ となる正則行列 $S$ と $T$ をも求める。この場合には次のようにやればよい。

$M$ を $n\times m$ の行列とする。その下で $(n+m)\times(n+m)$ の正方行列 \[
\begin{pmatrix} M&E_n \\ E_m&O \end{pmatrix} \] を作る。ここに $E_n$ は $n\times n$ の単位行列、$E_m$ は $m\times m$ の単位行列である。$O$ は $m\times n$ の零行列である。

基本変形は $O$ の部分に作用しないように範囲を限定する。従って行変形は $1$ から $n$ に、列変形は $1$ から $m$ に限定する。その下で変形を繰り返して対角化を進めるのである。最終的な結果は $$ \begin{pmatrix} \hat M&S \\ T&O \end{pmatrix} $$ になっている。

例えば $M=\begin{pmatrix}
5&6&8 \\
6&-11&7
\end{pmatrix}$ であれば \[
\begin{pmatrix}
5&6&8&1&0\\
6&-11&7&0&1\\
1&0&0&0&0\\
0&1&0&0&0\\
0&0&1&0&0
\end{pmatrix} \]を作るのであるが、実際の手計算では、これを
\[
\begin{array}{rrr|rr}
5&6&8&1&0 \\
6&-11&7&0&1 \\
\hline
1&0&0&\\
0&1&0&\\
0&0&1 &
\end{array}
\]と書き換えた方がわかりやすいであろう。

\[ \begin{array}{rrr|rr} 5 & 6 & 8 & 1 & 0 \\ 6 & -11 & 7 & 0 & 1 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|rr} 5 & 6 & 8 & 1 & 0 \\ 1 & -17 & -1 & -1 & 1 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|rr} 1 & -17 & -1 & -1 & 1 \\ 5 & 6 & 8 & 1 & 0 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|rr} 1 & -17 & -1 & -1 & 1 \\ 0 & 91 & 13 & 6 & -5 \\ \hline 1 & 0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \\ \Longrightarrow \begin{array}{rrr|rr} 1 & 0 & 0 & -1 & 1 \\ 0 & 91 & 13 & 6 & -5 \\ \hline 1 & 17 & 1 & \\ 0 & 1 & 0 & \\ 0 & 0 & 1 & \end{array} \Longrightarrow \begin{array}{rrr|rr} 1 & 0 & 0 & -1 & 1 \\ 0 & 13 & 91 & 6 & -5 \\ \hline 1 & 1 & 17 & \\ 0 & 0 & 1 & \\ 0 & 1 & 0 & \end{array} \Longrightarrow \begin{array}{rrr|rr} 1 & 0 & 0 & -1 & 1 \\ 0 & 13 & 0 & 6 & -5 \\ \hline 1 & 1 & 10 & \\ 0 & 0 & 1 & \\ 0 & 1 & -7 & \end{array} \] これから
\[ S=\begin{pmatrix} -1 & 1 \\ 6 & -5 \end{pmatrix} ,\quad \hat M=\begin{pmatrix} 1 & 0 & 0\\ 0 & 13 & 0 \end{pmatrix} ,\quad T=\begin{pmatrix} 1 & 1 & 10\\ 0 & 0 & 1 \\ 0 & 1 & -7 \end{pmatrix} \] を得る。

連立1次ディオファントス方程式の解法との類似は明瞭である。すなわち ${\boldsymbol b}$ が $S$ に、$E$ が $T$ に置き換わったに過ぎない。なお文献[1]での $S^{(1)}S^{(2)}\cdots$ が、ここでは $S$ で、$T^{(1)}T^{(2)}\cdots$ が $T$ とされている。

対角化のアルゴリズムを正確に自然言語で説明することは難しい。正確なアルゴリズムはプログラムコードを見るのが一番良い。

Maxima のプログラム

実行例1

行列 $M=\begin{pmatrix}5&6&8\\ 6&-11&7\end{pmatrix}$ の単因子標準形 $\hat M$ を求める。

(%i49) M:matrix([5,6,8],[6,-11,7]);

(%o49) matrix([5,6,8],[6,-11,7])
(%i50) elmdiv(M);

step0 0
              [matrix([1,0],[0,1]),matrix([5,6,8],[6,-11,7]),
               matrix([1,0,0],[0,1,0],[0,0,1])]
step1 1
     [matrix([1,0],[0,1]),matrix([5,6,8],[6,-11,7]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step2 1
     [matrix([1,0],[-1,1]),matrix([5,6,8],[1,-17,-1]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
中略
step1 2
     [matrix([-1,1],[6,-5]),matrix([1,0,0],[0,13,0]),
      matrix([1,1,10],[0,0,1],[0,1,-7])]
(%o50) [matrix([-1,1],[6,-5]),matrix([1,0,0],[0,13,0]),
        matrix([1,1,10],[0,0,1],[0,1,-7])]

ここで得られた

matrix([-1,1],[6,-5])
matrix([1,0,0],[0,13,0])
matrix([1,1,10],[0,0,1],[0,1,-7])

が $S,\hat M, T$ を与えている。すなわち
\[
S=\begin{pmatrix}
-1&1 \\
6&-5
\end{pmatrix}
,\quad
\hat M=\begin{pmatrix}
1&0&0\\
0&13&0
\end{pmatrix}
,\quad
T=\begin{pmatrix}
1&1&10\\
0&0&1 \\
0&1&-7
\end{pmatrix}
\]
である。

実行例2

行列 $M=\begin{pmatrix}4&3&6&2\\ 2&3&6&4 \\6&6&13&5\end{pmatrix}$ の単因子標準形 $\hat M$ を求める。

(%i51) M:matrix([4,3,6,2],[2,3,6,4],[6,6,13,5]);
(%o51) matrix([4,3,6,2],[2,3,6,4],[6,6,13,5])
(%i52) elmdiv(M);
step0 0
              [matrix([1,0,0],[0,1,0],[0,0,1]),
               matrix([4,3,6,2],[2,3,6,4],[6,6,13,5]),
               matrix([1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1])]
step1 1
     [matrix([1,0,0],[0,1,0],[0,0,1]),matrix([2,3,6,4],[4,3,6,2],[5,6,13,6]),
      matrix([0,0,0,1],[0,1,0,0],[0,0,1,0],[1,0,0,0])]
step2 1
     [matrix([1,0,0],[-2,1,0],[-2,0,1]),
      matrix([2,3,6,4],[0,-3,-6,-6],[1,0,1,-2]),
      matrix([0,0,0,1],[0,1,0,0],[0,0,1,0],[1,0,0,0])]
中略
step1 3
     [matrix([-2,0,1],[8,1,-4],[13,1,-6]),
      matrix([1,0,0,0],[0,1,0,0],[0,0,6,0]),
      matrix([0,0,0,1],[0,1,-2,-4],[0,-1,3,1],[1,1,-3,1])]
(%o52) [matrix([-2,0,1],[8,1,-4],[13,1,-6]),
        matrix([1,0,0,0],[0,1,0,0],[0,0,6,0]),
        matrix([0,0,0,1],[0,1,-2,-4],[0,-1,3,1],[1,1,-3,1])]

ここで得られた

matrix([-2,0,1],[8,1,-4],[13,1,-6])
matrix([1,0,0,0],[0,1,0,0],[0,0,6,0])
matrix([0,0,0,1],[0,1,-2,-4],[0,-1,3,1],[1,1,-3,1])

が $S,\hat M, T$ を与えている。すなわち
\[
S=\begin{pmatrix}
-2&0&1 \\
8&1&-4 \\
13&1&-6
\end{pmatrix}
,\quad
\hat M=\begin{pmatrix}
1&0&0&0\\
0&1&0&0 \\
0&0&6&0
\end{pmatrix}
,\quad
T=\begin{pmatrix}
0&0&0&1\\
0&1&-2&-4 \\
0&-1&3&1 \\
1&1&-3&1
\end{pmatrix}
\]
である。

プログラムコード

lde.mac の主関数の一つ elmdiv(M) の定義を次に示す。
列と行変形の基本的なアルゴリズムが正確に解るはずである。詳しくは lde.mac を直接参照されたい。

elmdiv(M):=block([B,n,k,r],
  n:min(length(M),length(M[1])),
  B:tab_gen(M),
  dprint("step0",0,B),
  /* don't use 'for' loop */
  k:1,
  /* don't replace B[2] with A, which behaves differently! */
  while k <= n and chk_zero(B[2],k)=false do (
    X:step1(B,k),
    if X=false then (k:n+1, return(false)),
    B:X,
    dprint("step1",k,B),
  	loop, /* we should avoid go(loop) if possible. this is merely an example for go(loop) */
    unless chk_col(B[2],k) do (
      B:step2(B,k),
      dprint("step2",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    unless chk_row(B[2],k) do (
      B:step3(B,k),
      dprint("step3",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    if chk_col(B[2],k)=false then go(loop),
    /* we may assume here both chk_col(B[2],k) and chk_row(B[2],k) are true */
    if chk_mod(B[2],k) then k:k+1
    else (
      X:step4(B,k),
      if X=false then k:n+1
      else (
        B:X,
        dprint("step4",k,B)
      )
    )
  ),
  r:chk_final(B,M),
  if r # 0 then error("# ERROR: bug!",r),
  return(B)
);

step4(B,k) が追加されたことになる。

SymPy のプログラム

実行例1

行列 $M=\begin{pmatrix}5&6&8\\ 6&-11&7\end{pmatrix}$ の単因子標準形 $\hat M$ を求める。

>>> M=Mat([5,6,8],[6,-11,7])
>>> elmdiv(M)
step0 0 [Matrix([
[1, 0],
[0, 1]]), Matrix([
[5,   6, 8],
[6, -11, 7]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
step1 1 [Matrix([
[1, 0],
[0, 1]]), Matrix([
[5,   6, 8],
[6, -11, 7]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
中略
step1 2 [Matrix([
[-1,  1],
[ 6, -5]]), Matrix([
[1,  0, 0],
[0, 13, 0]]), Matrix([
[1, 1, 10],
[0, 0,  1],
[0, 1, -7]])]
[Matrix([
[-1,  1],
[ 6, -5]]), Matrix([
[1,  0, 0],
[0, 13, 0]]), Matrix([
[1, 1, 10],
[0, 0,  1],
[0, 1, -7]])]
>>>

ここで得られた

Matrix([[-1,1],[6,-5]])
Matrix([[1,0,0],[0,13,0]])
Matrix([[1,1,10],[0,0,1],[0,1,-7]])

実行例2

行列 $M=\begin{pmatrix}4&3&6&2\\ 2&3&6&4 \\6&6&13&5\end{pmatrix}$ の単因子標準形 $\hat M$ を求める。

>>> M=Mat([4,3,6,2],[2,3,6,4],[6,6,13,5])
>>> elmdiv(M)
step0 0 [Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]), Matrix([
[4, 3,  6, 2],
[2, 3,  6, 4],
[6, 6, 13, 5]]), Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])]
step1 1 [Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]), Matrix([
[2, 3,  6, 4],
[4, 3,  6, 2],
[5, 6, 13, 6]]), Matrix([
[0, 0, 0, 1],
[0, 1, 0, 0],
[0, 0, 1, 0],
[1, 0, 0, 0]])]
中略
step1 3 [Matrix([
[-2, 0,  1],
[ 8, 1, -4],
[13, 1, -6]]), Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 6, 0]]), Matrix([
[0,  0,  0,  1],
[0,  1, -2, -4],
[0, -1,  3,  1],
[1,  1, -3,  1]])]
[Matrix([
[-2, 0,  1],
[ 8, 1, -4],
[13, 1, -6]]), Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 6, 0]]), Matrix([
[0,  0,  0,  1],
[0,  1, -2, -4],
[0, -1,  3,  1],
[1,  1, -3,  1]])]
>>>

ここで得られた

Matrix([[-2,0,1],[8,1,-4],[13,1,-6]])
Matrix([[1,0,0,0],[0,1,0,0],[0,0,6,0]])
Matrix([[0,0,0,1],[0,1,-2,-4],[0,-1,3,1],[1,1,-3,1]])

プログラムコード

mlde.py の主関数の一つ elmdiv(M) の定義を次に示す。
列と行変形の基本的なアルゴリズムが正確に解るはずである。詳しくは mlde.py を直接参照されたい。

def elmdiv(M):
  n = min(M.shape)
  B = tab_gen(M)
  dprint("step0",-1,B)
  k = 0
  # don't replace B[1] with A, which behaves differently!
  while k < n and chk_zero(B[1],k)==False:
    X = step1(B,k)
    if X == False: break
    B = X
    dprint("step1",k,B)

    while True:
      while chk_col(B[1],k) == False:
        B = step2(B,k)
        dprint("step2",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      while chk_row(B[1],k) == False:
        B = step3(B,k)
        dprint("step3",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      if chk_col(B[1],k) == True: break

    # we may assume here both chk_col(B[1],k) and chk_row(B[1],k) are true
    if chk_mod(B[1],k): k += 1
    else:
      X = step4(B,k)
      if X==False: break
      B = X
      dprint("step4",k,B)
  r = chk_final(B,M)
  if r != 0 : error("# ERROR: bug!",r)
  return B

step4(B,k) が追加されたことになる。

正則な整数行列の生成

2020-11-26

教材を作るのに、行列式が $\pm 1$ の整数行列 $P$ をランダムに生成できると何かと良い。
例えば与えられた Jordan 行列を $J$ とすると $P^{-1}JP$ によって、$J$ と相似な整数行列を作ることができる。
lde.mac および mlde.py には

mat_rand(n,m,nr)

が付属している。n は正方行列のサイズ、m と nr には適当な正整数を与えればよい。
m が大きいと、生成される行列要素も大きくなる傾向がある。
nr が大きいとランダム性が増す。