数理解析研究所講究録920 巻 1995 年 62-73

627.
Fraction-free による行列式の計算効率
工学院大電気系 木下 孝
(Takashi Kinoshita)
工学院大数学 牧野 野夫(Isao Makino)
工学院大電子 三好 和憲(Kazunori Miyoshi)

We consider to calculate determinants of matrices whose coefficients are integers or polynomials
with integral coefficients. Speed and precision are important factors for mathematical computing
system. On the calculation of determinants, Gauss elimination method is well known. But themethod takes much time and causes significant errors, because the calculation deals with rationalor decimal fraction numbers.

Our paper presents Sylvester’s Identity and Chinese remainder theorem as the methods tocalculate determinants which avoid computation of rational numbers, and the close analysis of
the methods.

7.1 はじめに整数に関する線形計算を行う際数式処理システムにおいて, 任意精度でかつ高速な計算は不可欠である. また科学技術計算では, 線形計算の解を求めるのに殆ど実行時間を費やしている現状である.

そこで本報告では, 線形計算でよく使われている計算の中から, 行列式と連立 1 次方程式を扱った.
またその行列式計算を応用し, 終結式を用いた, 互いに素である 1 変数多項式の extended
1 $f,$ $g$ が互いに素な多項式のとき,$f\cdot X+g\cdot Y=1$$\mathrm{g}\mathrm{c}\mathrm{d}1$も取となる多項式 $X,$ $Y$ をこのように呼ぶこととする.

63り扱った. 特に我々は, これらの基礎である整数要素の行列式計算を主に取り上げ, 効率よく求める方法を調査し, 評価したので報告する.

7.2 行列式の計算法
一般的な行列式の求め方に掃き出し法がある. しかしこの方法では, 整数要素の演算において結果が整数であることは自明であるにもかかわらず, 計算時間のかかる有理数演算を行う. さらに, 有理数を実数として処理するシステム上では, 誤差を生じる問題も発生する. 従って, 有理数や実数を用いた掃き出し法をそのまま利用するのは, 数式処理システム上において適した方法でない.
そこで有理数演算を用いずに, 整数要素や整数係数の多項式要素の行列式を整数演算によって求めるアルゴリズムとして, 次のような方法がある.
(J)

掃き出し法において分母を別に計算する (Sylvester’s Identity).

(II)

互いに素である複数の法を掃き出し法で計算し, 最後にまとめて中国剰余定理と Hadamardの不等式により行列式を求める.

7.2.1 計算法 (I)$\mathrm{S}\mathrm{y}1_{\mathrm{V}\mathrm{e}\mathrm{s}}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s})$Identity による方法

ここでは, 以下に示す Sylvester の恒等式 (Sylvester’s Identity) を用いて, 行列式の値を求める.

$\blacksquare$

Sylvester’s Identity$A=[a:,j1$ を$n$次正方行列とする.

各 $k(=0,1, \star. . , n-1)$ に対し, $A^{(k)}=\mathrm{r}_{a}!_{j}^{k)}.,1(k+1\leq i\leq n, k+1\leq i\leq n)$ を次のように定める(1)$|a_{i,j}^{(0\rangle}|=A$$a_{0.0^{1)}}^{\mathrm{t}-}=1$とし,$a_{i,j}^{\langle k)}=\det\{$$\mathrm{L}$$a_{1,1}$$a_{1,2}$$a_{1,k}$$a_{1,j}$$a_{2,1}$$a_{2,2}$$a_{2,k}$$a_{2,j}$.

.

.

$a_{k.1}$$a_{k,2}$$a_{1}.,1$$a:,2$:.
$a_{k},j$$a_{k.k}$$\mathrm{u}_{l,k}$と定義したとき,$[$$\det(A)[a_{p,p}^{()}p-\iota]n-l-1=\det\{$$a_{p+1,n}^{(p)}$$a_{\ell+1,\ell+}^{\mathrm{t}^{p_{)}}}1$.
$(\ell)$$a_{n,\ell+1}$

..

:

(3)
$\rfloor$$(\ell)$$a_{\mathfrak{n},n}$64となる $(1\leq P\leq n-1)$ . 特に式 (2) の右辺にあらわれる行列を改めて $A$ とおいたとき, $A$ は $(k+1)$次正方行列であるので, ここで式 (3) を適用すると以下のように応用できる.
$.\ell=k-1$ としたとき single-step fraction-freeと呼ばれ, 次のように表される.

$a_{i,j}^{(k)}= \frac{1}{a_{k-}^{(k2}\overline{\iota},k-1)}$

(4)

つまり
$a_{i,j}^{(k)}=(a_{k,k}^{(k)}-1a!_{j}^{k1)}.,--a^{(k-}ak,j\cdot.,k1)(k-1))/a(k-1k-2,)k-1$’$(k+1\leq i\leq n, k+1\leq j\leq n;k=1,2, \ldots, n-1 )$となる. またこのとき $A^{(k\iota)}-arrow A^{(k)}$ の計算量は以下の通りである.

乗算:$2(n-k)^{2}$ ,$.\ell=k-2$ としたとき

加減算 :

除算

:

(5)

$(n-k)^{2}$

two-step fraction-free と呼ばれ, 次のように表される.

$a_{k,k-1}^{1^{k}}a_{k,k-}^{(k}a_{i}^{()}-,-k\mathit{2}1k^{-}-1-2)2)1$$a^{(k)}.

$(n-k)^{2}$ ,$a_{k.j}^{(.)}a_{k}^{\langle)}a^{(k^{-1}2}|,k,-\mathit{2}kjj-2-)$$a_{k}^{1,)}a_{k,k}^{(k}a_{\mathfrak{i}}^{(}k-2k-2)-1,kk^{-2})$

$|$

(6)

\cdot,j=\frac{1}{[a_{k-2}^{(k3)}-,k-\mathit{2}]^{2}}$つまり$=(a_{k-1}^{(-},-a_{k}^{(}k\mathit{2})k1k,-k2)-a_{k-1k}^{(\mathrm{z})(\mathit{2})}k-,a_{k},-1)k/a_{k}k-(k-\mathrm{s})-2,k-2$$c_{0^{k-}}^{(\mathit{2})}$$c_{i1}^{(k-2)}$$=(a_{k-1}^{(-}k\mathit{2},)k\cdot.,k^{-}-1O(k2)-a-a(k2)(.k-2))/a(k-k-1,k-1\cdot.kk-2,ks)-2$$c_{i\mathit{2}}^{(k-\mathit{2})}$$=(a_{k,k-1}^{(-}ak2)!_{k}.,k-2)-a_{k,k}^{(-}k2)a-)\mathrm{t}k\mathit{2}i,k^{-}1)/a_{k-}^{(}k-3\mathit{2},k-2)$

’

’

,

$\{$

$=(c_{0^{k2}}^{(-}a!^{k}).,-2)(.k-2)(k-2)++\mathrm{C}_{1}acj\cdot k,j!_{2^{-}}.k2)a_{k-}^{(2)}-,)1j/ka_{k}(k-s_{k-2}-\mathit{2}|)$$a^{(k)}.\cdot,\mathrm{j}$,,$a_{k,k}^{(k-1)}$$=c_{\mathrm{o}^{k-}}^{(2)}$$a_{k,m}^{(k-1)}$$=a_{k,m}^{(k)}=(a_{k-1}^{(-},Ok2)k-1k(k,-m2)-a-1,ma-1)(k-\mathit{2})(k-2kk,k)/a_{k}(k-s)-2,k-\mathit{2}$’$(k+1 \leq i\leq n, k+1\leq j\leq n, k+1\leq m\leq n;k=2,4, \ldots, 2\lfloor\frac{n-1}{2}\rfloor)$となる. またこのとき A(k-2)\rightarrow A(りの計算量は以下の通りである.

:$3(n-k)2+4(n-k)+2(n-k)+2$

加減算 :$2(n-k)^{2}+2(n-k)+(n-k)+1$

岡づ n$(n-k)^{\mathit{2}}+2(n-k)+(n-k)+1$

乗算
$\{$

:

$.P=k-3$ としたとき three-step fraction-free と呼ばれ, 次のように表される.

(7)

65

$(k-3)$

a $(k-3)$a$a_{k-1}^{(}k-\mathit{3},)k-2$$a_{k-1}^{\langle k3)}-,$$k-\mathit{2},$ $k-\mathit{2}$$k-\mathit{2},$$k-$ ]

a $(k-3)$$k-2,$

a$a_{k-\overline{\iota},k}^{(k\mathit{3}})$$a_{k-1}^{(3)}k-,J$

$k$

$k-1$

$a_{i,j}^{(k)}= \frac{1}{[a_{k-3}^{1,)}k-4k-3]^{3}}$$(k-\mathit{3})$$k-\mathit{2},j$$a_{k,k^{-}-2}^{()}k\mathit{3}$$a_{k,k^{-3}1}^{(k}-)$$a_{kk}^{1,-}k3)$$a_{k,j}^{1^{k3})}-$$k^{-}-2)$$a^{(k3}.\cdot,$$a^{(k-3)}|.,k-1$$\mathrm{a}_{\mathfrak{i},k}^{(3)}k-$$a^{(k-3)}$.
$(8^{\cdot}\dagger$

つまり$Det3(a, b, C, d, e, f, g, h, i, j)=\{a(ei-fh)+b(fg-di)+c(dh-eg)\}/j$としたとき$c_{0^{k-}}^{(3)}$$a_{k,k-}^{(k-s)(}2’ a_{k}k,-s)k-1’ ka^{1},-\mathit{3})akk’(k-k-3,k-4)3)/a_{k-}(k-4)3,k-3$$c_{i\mathit{3}}^{(k-3)}$,$Det3(a_{k\mathit{2},k}^{1^{k}}--2’ k-2,k-1’ a-3)a(k-3)(k-3)k-2,k’ k-a-1(k\mathit{3}.)k-2’ a^{(}k-k-13,)k-1’ a_{k,k-1}^{(-}k\mathit{3})$$=$’$a^{\mathrm{t}k3)}k\mathit{3})k-\mathit{2},k-1a-ka_{k}k2,-1.k-2a_{k}-\iota.k-1k,k--1$$Det3(a_{k}^{(}-2,k-2’
a^{\mathrm{t})(k3)(k}-k-3,-,3),(k-3),$$=$$a_{i,k^{-}}-1’ a(k)\mathrm{t}k\mathit{3})(.\cdot,kk^{-}’ k3)a-3,k-\mathit{3})(k-4)/a\mathrm{t}k-4)k-3,k-3$,$a_{\mathfrak{i}.k-\mathit{2}}-3,$$c_{i2}^{(k-\mathit{3})}$$a_{k}k\mathit{3})kk3)(k-\mathit{3}2,k-1’
k,k)(ak-3\rangle$$Det3(a-2,k-2’ a\mathrm{t}k-3)\mathrm{t}k-\mathit{3})kk-2,k-1’ k-2a^{(}-,’
a_{k,k-}^{(-},$$=$$a^{(k-3}.\cdot.-,$$a^{(k}-1’$\{$$c_{i1}^{1-}k3)$’ak2)i.k^{-}3)!^{k-}..k’ k-3.k-3)/3)a^{1-}k4)(k-ak-3,k-4)\mathit{3}$’$Det3(a-1,’ a1^{k}-kk-2k-1,k3)(k-\mathit{3})-1’ a_{k1k}^{(3)1^{k})(k-s)}k--,’ a_{k},-\mathit{3}a-1’ a_{k}k-2’ k,k(k,-3)k$$=$$a^{(k\mathit{3})}k2)|k-1\cdot.,k^{-}’
a\mathrm{t}k-k-3,k-4)\mathit{3})/a_{k-\mathit{3}}(k-4,)k-3’$$a^{(k-3}.\cdot,-,$$a^{(k3)}.,-,$)$(c_{0^{k-3}..j}^{\mathrm{t})}a^{(k}.-\mathit{3})-Ca-\mathit{3}+i3k,jc(k-3)\langle k$$1.\cdot 2k-3$

)

,

’

.

)

$|1^{-}k/(k\mathit{3})(k-\mathit{3})a(k-4)k-3.k-3$$a_{kj}^{(3)}k--1,-C\cdot
a-2,j$$a^{(k)}.\cdot,j$$=$$a_{k.k}^{(k)}-\mathrm{l}$$=$$c_{0^{k-}}^{(3)}$$a_{k1h}^{(k2)}--$$=$$a_{k-1,h}^{(k)}=(a_{k-2}^{(-}.ak3)k-2(k-3)k-1.h-a_{k-\overline{\iota}}^{(k}3,)k-\mathit{2}(a_{k-2}-,)h/a^{(k-}k3)k-3,k-4)3$$a_{k,m}^{(k-1})$$=$$a_{k,m}^{(k)}=Det3(a_{k}^{(}-2_{1}k-\mathit{2}’a-\mathit{3})(k-3),$$a-1k-\mathit{2},m’a-3)(k-\mathit{3}),$$a-,k-\mathit{2}k-2,k-1’
akk-\mathit{2},k\mathrm{t}kk1(k-3)(k-3)k-1,m$

.

’

,
’

’

$a_{k,k-}^{(-}a_{k}^{(},-\mathit{3})a^{()}k\mathit{3})\mathit{2}’ kk-\iota’ km’ k-|3,k-a^{\langle}\mathit{3}k-\mathit{3}k-4))/a_{k-}(k-4)3,k-3$’$(k+1 \leq i\leq n, k+1\leq j\leq n, k+1\leq m\leq n, k-1\leq h\leq n;k=3,6, \ldots, 3\lfloor\frac{\mathfrak{n}-1}{3}\rfloor)$となる. ただし,$n\equiv 2$(mod 3)のときは最後に以下の手順を行う.

$a_{\mathfrak{n},n}^{(\mathfrak{n}-1)}=(a_{n-}^{(\mathfrak{n}-\mathit{3})}2,n-2n-a^{()}-1,nn\mathit{3}-1-a_{\mathfrak{n}-}^{(}n-3)\mathrm{t}n2,\mathfrak{n}-1n-a-1\mathit{3},)\mathfrak{n}-2)/a^{(}\mathfrak{n}n-\overline{s}_{n-}4,)\mathit{3}$

またこのとき A(k-3)\rightarrow A(りの計算量は以下の通りである.

:$4(n-k)\mathit{2}7(+2n-k)+9(n-k)+2(n-k+2)+9$

加減算 :$3(n-k)^{\mathit{2}}+15(n-k)+5(n-k)+(n-k+2)+5$

乗算
$\{$

隔づ

:

となる. ただし, いずれの

(9)$(n,-k)^{2}+6(n-k)+2(n-k)+(n-k+2)+2$$\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}$の fraction-free$a_{k_{J}}^{(1)}k,-=0$

においても

,

$J<k$

,

(10)

である. 式 (1), (10) の条件に基づき, 式 (4), (6), (8) のいずれかの掃き出し法による漸化式を利用して計算していくと, 最終的に元の行列は, 式 (11) に示すような上三角行列に変形される (ただし, 空

白は

$0$

). なお本プログラムでは掃き出し法において, 整数要素を扱うことを前提にしているので, 数

66

要素のときのみ掃き出しの軸成分を絶対値最小になるように選んで掃き出しを行っている.

. ..

$a_{1,1}^{(0)}$)$a_{2,2}^{(1)}a_{1}^{()}0,2$$n\mathfrak{n}$$a_{k,k}^{(k-1)}\ldots$$A^{(n-1)}=[$$a_{\mathfrak{n}’}^{(\mathfrak{n}-1}a_{k}^{(k}a^{(}a_{2}^{\langle)}10,’$$nn1-1$))$’(11)$

このとき$a_{n,n}^{(n-}$りが元の行列の行列式になるアルゴリズムである. すなわち
$a_{n,\mathfrak{n}}^{(\mathfrak{n}-1)}=\det(A)$.

つまり本アルゴリズムは, 式 (1), (10) に基づき, 式 (4), (6), (8) のいずれかの掃き出し法による漸化式を利用し, 元の行列$A$を順に帰納的に変形させて, 果 n,n-]) である行列式を計算する方法である.

最後に式 (5), (7), (9) により, $A^{(0\rangle}arrow A^{(n-1)}$ の総計算量をまとめると表 1 となる.

7.2.2 計算法 (II) 中国剰余定理による方法$\blacksquare$

中国剰余定理$m_{1},$ $m_{2},$ $\ldots,$$m_{1}$のどの 2 つも互いに素である連立合同式$r_{1}$$r_{2}$$x\equiv\{$

(mod

(mod

$m_{1}$$m_{\mathit{2}}$)).
$r_{t}$$(\mathrm{m}\mathrm{o}\mathrm{d} m_{\iota)}$

(12)

$]$

67

に解は存在し, その 1 つの解を

$x_{t}$とすると -般解は(mod$x\equiv x_{l}$である.

$m_{1}m_{\mathit{2}}\ldots m_{1}$

(13)

)

$\blacksquare$

ここでは上記の定理を利用して, 以下のように行列式を計算する.

加減乗算はある整数を法として計算することによって, 小さな桁数の数で計算ができ計算量を減らせる.

さらにそのような法のもとで, 逆元が存在する場合は, 除算も有理数演算を利用せずに乗算によって同様に計算できる. すなわち行列式は, ある整数を法とした掃き出し法で計算できる場合がある. しかし法の値が小さい時そのような法のもとで条件を満たす行列式は, 法を周期とした比較的に短周期で複数個存在するので, -意に定まらない.
そこで, 互いに素なる法で行列式$r_{1},$$r_{2},$$\ldots,$$t$$\Gamma_{l}$個の比較的小さな整数$m_{1},$ $m_{2},$ $\ldots,$$m_{1}$

を法として, 各法について掃き出しを計算する. それからこれらを式 (12) に当てはめて, 連立合同式を解けば,

式 (13) のように解の周期を長くすることができる. さらに与えられた行列式の限界値が既知のとき,その範囲を超えるように式 (13) の法 $m_{1}m_{2}\ldots m_{t}$ を定めれば, 行列式の値は–意に求まる.

従って本アルゴリズム (は, まず行列式の限界値を計算し, この値を超えるように$t$個の法を定める.

ここで定めた各法における行列式を掃き出し法によって計算し, これらの連立合同式を解いて行列式を求める方法である.
その行列式の値の限界を求める方法として, 次に示す Hadamard の不等式がある. また, 連立合同式の解法は, 次に示す Gauss の方法, Garner の方法の 2 通りある.
さらに我々は, 大きな次数や桁数要素の行列を扱うことを前提にしている. 従ってそのような場合においても, 行列式を計算できるような法を抽出せねばならない. そこでそのような法の抽出方法として, 次の 3 通りの方法を利用する.

1)

システム内の素数テーブルを最大の数から利用する方法.

(本開発システム上では 10 進 8桁の素数 999 個:99, $981,793\sim 99,999,989$ )2)

円分数 $M(35, X),$ $M(39, X)$ ;$2\leq x\leq$

1,000 により 50 桁以上の素数 508 個を予め (手作業で) 抽出した素数テーブルを読み込んで計算する.
3)$n$番目の–般フィボナッチ数$\mathrm{g}\mathrm{c}\mathrm{d}(P, q)=1$,$fo=0$ ,$f_{n}$の値の相関性$f_{1}=1$とし,$f_{n}=p\cdot f_{n-1}+q\cdot f_{\mathfrak{n}-2}$と定めたとき$\mathrm{g}\mathrm{c}\mathrm{d}(m, n)=1\Rightarrow \mathrm{g}\mathrm{c}\mathrm{d}(f_{m}, f_{\mathfrak{n}})=1$

を利用して法$f_{n}$を抽出する. ただし計算の簡略化のため $m,$$n$は素数, $p=q=1$ とした.

68$\blacksquare$

Hadamard の不等式
次正方行列 $A=1^{a}\cdot,j1$ に対し$n$$(1 \leq i\leq n, 1\leq j\leq n)$$| \det(A)|\leq.\cdot\prod_{=1}^{n}..\sqrt{a^{2}+1a^{2}2^{++a^{2}}i\mathfrak{n}}$が成立する. この式により, $|\det(A)|$ の上限値がわかる.

$\blacksquare$

Gauss の方法,

$M=m_{1}m_{2}\ldots m_{t}$とする. このとき$\mathrm{g}\mathrm{c}\mathrm{d}(M_{k}, m_{k})=1$$M_{k}=M/m_{k}$である. そして(mod$\ell_{k}M_{k}\equiv 1$なる$p_{k}$$m_{k}$)をとったとき, 式 (12) の解は以下の通りである.
$x\equiv r_{1}P_{1}M1+r_{2}l_{2}M_{2}+\ldots+r_{t}\ell \mathrm{c}M_{t}$$\blacksquare$$(1 \leq k\leq t)$(mod $M$ )Garner の方法

式 (12) の解は以下のように定められる.
$v_{\mathit{0}}\equiv r_{1}$とする. このとき(mod$m_{1}$$\mathrm{g}\mathrm{c}\mathrm{d}(L_{k-}1, mk)=1$),$L_{k}=.

\cdot\prod_{=1}^{k}m$.

$(1 \leq k\leq t)$である. そして$v_{k-1}L_{k-1}\equiv r_{k}-(v_{\mathit{0}}+v_{1}L_{1}+\ldots+v_{k-2}L_{k2}-)$によって$m_{k}$$v_{1},$(mod$m_{k}$)$v(0\leq v_{k}<m_{k+1})$ を順に帰納的に定めたとき (このとき途中の計算は法$v_{2\cdot\cdot \mathrm{t}-1},.,$の最大値内で行われる), 解は$x\equiv v_{\mathit{0}}+v_{1}L_{1}+\ldots+v_{1-1}L_{\iota-1}$(mod$L_{1}$)

である. この 1 回の計算でのみ大きな値を扱うため, 計算が高速化される利点がある.

7.3 応用プログフ- ム
72 節で述べた行列式計算アルゴリズムを応用して, 整数または整数多項式の係数の連立 1 次方程式,
1 変数多項式の extended$\mathrm{g}\mathrm{c}\mathrm{d}$解法プログラムも開発した.

連立 1 次方程式は, アルゴリズムとして Sylvester の fraction-free を応用した Gauss の消去法を用

69

いた. この方法で不定解, 不能解も求められ, また文字係数の場合も解けるようにした.

1 変数多項式の extended$\mathrm{g}\mathrm{c}\mathrm{d}$は, 整数係数の多項式演算を扱わねばならないので, Sylvester の

fraction-free を応用した. 入力した 2 つの多項式が互いに素 (Resolutant$\mathrm{g}\mathrm{c}\mathrm{d}$$\neq 0$) であるときに extendedを求めるようにした.

なおいずれのプログラムにおいても, 入力に対して正しい出力が得られることを確認した.

7.4 実験内容ならびに結果
はじめに開発環境は以下の通りである.
ソフトウェア

富士通情報社会科学研で開発中の数式処理処理システム/ ライブラリ Risa の計算エンジンの言語インタ一フェース asir
ハードウエア$\mathrm{H}\mathrm{P}9000/735$($\mathrm{O}\mathrm{S}$:HP-UX, メモリ$80\mathrm{M}\mathrm{B}$)

この環境下で 72 節で述べた 3 種のアルゴリズムによる行列式の演算効率の比較を下記の要領で行った.
比較に用いた行列式の演算プログラムは (I) の single, two, three-step fraction-free, (II) の Gaussの方法, Garner の方法, 有理数演算を用いる掃き出し法, asir 内の行列式を求める組み込み関数
して asir のライブラリ内の行列式を求める関数 deter の 8 種である. なお, 関数$\mathrm{d}\mathrm{e}\mathrm{t}$$\mathrm{d}\mathrm{e}\mathrm{t}$, そ, deter は single-

step fraction-free に基づいて作成されている. これらの関数を用いて以下の実験, 計測を行った.
A)

次数は 9\sim 15 次で, 要素は 100\sim 1,000 桁までの整数要素の正方行列を用いて行列式を求め, 計算時間の計測を行う.

B)

Sylvester の fraction-free アルゴリズムについて, 次数は 10\sim 300 次で, 要素は 5 または

10 桁までの整数要素の正方行列を用いて行列式を求め, 計算時間の計測をする.
C)

Sylvester の fraction-free アルゴリズムについて, 多項式演算ができることを確認し, 1 次多項式要素 (係数は 5 桁までの整数) の正方行列を用いて行列式を求め, 計算時間の計測を行う.
D)

中国剰余定理のアルゴリズムついて, 722 節で述べた法により, 8\sim 11 次の正方行列の行列式を求め, 計算時間の計測を行う.
上記の方法により, 以下のような行列式の演算効率がわかる. なおこれらの実験結果の–部を 77 節にまとめておいた.

7.4.1 行列式の演算効志
まず, 上記 8 種の関数の計算時間の優劣を調査するため, A) の実験を行った. 実験結果より次のことがわかった. なお 12 次の測定結果のみ表 2 に示しておいた.

. 少ない桁数要素の時, Sylvester の fraction-free は中国剰余より計算時間を要さない. しかし,
◇ 9 次正方行列で 1,000 桁以上の要素において,

70$\text{◇}10$

次正方行列で 700 桁以上の要素において,$\text{◇}11$

次正方行列で 500 桁以上の要素において,$\text{◇}12$

次正方行列で 400 桁以上の要素において,$\text{◇}13$次正方行列で 350 桁以上の要素において,◇ 14 次正方行列で 300 桁以上の要素において,$\text{◇}15$

次正方行列で 250 桁以上の要素においていかなる step の Sylvester の fraction-free より中国剰余の方が, 逆に速くなる.

. 有理数演算は, Sylvester や中国剰余定理による整数演算アルゴリズムより約 10 倍以上計算時間を要する.

.Sylvester の fraction-free に着目すると, $9\sim 15$ 次といった低次において, two-step が最速であり, 次いで single-step, three-step の順である. このとき two-step は, single-step より 2 割程度
計算時間が短い. ただし次数が 3 の倍数の時, single-step と three-step の形成が逆転する.

. 組み込み関数やライブラリ関数は, 我々の single-step より若干計算時間を要する. とくに, 要
素の桁数が大きくなるにつれて, これらの計算時間の差も大きくなる.
実験 A) では, three-step の効果が発揮されなかった. しかし, 次数を増やすにつれて, three-stepが効果を発揮することが考察できる. そこで, 効果を発揮する次数を調べるため, 実験 B) を行った.

10 桁の整数要素における 100 次までの結果を表 3 に示し, この結果から次のようなことがいえる.
. 先に述べたように, 10$\text{次程度の低次}\mathrm{t}\text{こおいて}$

two-step が最速であり, 次いで single-step, three-step の順の計算スピードである. ところが, 20 次前後を超えると, three-step と single-step の計算時間が逆転する. さらに次数を増やし, 70 次前後以上に至ると two-step より three-step が速くな.

り, three-step が最速となる.
このことから, 次数を増すにつれて, 高 step の fraction-free が効果を発揮し, 有用となることが理解できる.

これまで整数要素について調査してきたが, C) の実験により, 多項式要素についての調査も行った.
結果は表 4 であり, これより以下のようなことが理解できる.
$.\mathrm{S}\mathrm{y}1_{\mathrm{V}\mathrm{e}\mathrm{S}}\mathrm{t}\mathrm{e}\mathrm{r}$の fraction-freeを用いたプログラムでは, 多項式の演算もできる.

. 整数要素と同様に, single-step より two-step は, 低次においても有効なアルゴリズムである.
さらに, single-step と three-step の計算時間の逆転する次数も, 整数要素の時とほぼ同じ 20 次前後である.

. 従って, 次数が大きくなるにつれて, 整数要素のとき同様に高 step の fraction-free アルゴリ.

ズムが有効であることが推測できる.
またこのとき, single-step は組み込み関数とほぼ同じ計算時間であり, GC time の分だけ僅かながら我々の single-step の方が速い. 従って整数要素の時には, 掃き出しの軸成分を絶対値最小になるように選んで掃き出しを行うと, 計算時間を短縮できることが考察できる.
これまで Sylvester の fraction-free アルゴリズムについてのみ着目してきたが, ここで中国剰余定

71

理によるアルゴリズムに着目する. 実験 A) では, 法を 2) によって抽出したが, 他の抽出方法について調査比較するため D) の実験を行った. 8 次の結果のみ表 5 に示しておいた. 実験結果より次のよ
うなことがいえる.
.1) の方法, つまり組み込み関数 lprime$()$

によって抽出した素数を法として利用する方法の実

験結果について考察する.
$\text{◇}$$\text{◇}$抽出された法は 1word 整数であるので, 1 つの法についての行列式計算の時間は短い.

しかし利用する法は, 2), 3) の方法に比べて多く必要とする. 従って, 他の 2 方法に比べて計算時間がかかる.
◇ lprime$()$の参照できる素数の数も 999 年分少ないため, 約 8,000 桁未満の行列式しか計算できない.

. 次に 2) の方法の実験結果について考察する.
$\text{◇}$法は 50 桁以上の整数を用いているため, 1 つの法についての計算時間は, 1word 整数による計算よりかなり要する.

◇使用する法は 1) よりかなり少ないため, 2 割程度計算時間が短くなっている.
$\text{◇}$しかしながら, 素数テーブルを手作業で作らねばならないといった大変な手間を要するため,

数式処理システム上において適した方法でない.

◇さらにテーブルの大きさに応じた行列式計算の限界 (本素数テーブルでは約 30,000 桁程度
までしか計算できない) がある.

. 最後に 3) の方法の実験結果について考察する. この方法は上記 2 方法の弱点を克服している.
$\text{◇}$$\text{◇}$

計算時間は, 2) の方法より多少要するが, 比較的短時間で計算できる.

法の計算方法も簡単であり, 数式処理システムを有効に利用できる.

◇計算できる行列式の桁数も, 500 番目までのフィボナッチ数の法を用いた場合で 172,176 桁であり, 相当な桁数まで計算できる. 本プログラムにおいて, 1,229 番目の素数 $(9,997)$ 番目のフィボナッチ数まで, つまり$f_{\mathit{2}}\sim f_{9997}$の 1,229 個の法をもてる.

. 従って, 中国剰余定理による行列式計算に用いる法は, フィボナッチ数を用いるのがよいと考えられる.

ことが確認できた.

7.5 評価本報告では, 整数要素による行列式演算の効率について述べてきた. 実験結果より次のようなことが確認できた.

oSylvester の fraction-free アルゴリズムは多項式要素の演算もでき, 次数を増やすにつれて高step のアルゴリズムが有用となる. また two-step は, 低次でも single-step より計算時間を要さな72

いため, かなり有効なアルゴリズムといえる. しかし, 要素の桁数が大きくなると不適である. また

整数要素の行列式計算は, 絶対値最小になるように軸成分を抽出すれば, 計算時間を短縮できる.

. 逆に中国剰余定理の方法は, 要素の桁数が大きいときに有効な手段である. また法としてフィボナッチ数を用いると, 計算時間も比較的短時間で計算でき, さらにかなり大きい桁となる行列式の計算ができる. しかしながら, 本プログラムでは, 多項式の演算ができない.

. 今回は整数整数係数の多項式要素について扱ったが, 要素の lcm を乗ずることによって有理数のみならず, Sylvester の fraction-free では, 多項式も扱えるので整数以外の要素のときにも応用できる.

. 京都の発表以後, 富士通情報社会研の野呂正行氏に HP の asir へ karatsuba による多倍長乗
算を impliment していただきました.

この karatsuba 法を用いると, 今回検討した中でかなり有用な fraction-free アルゴリズムである Sylvester の two-step において, 約 2\sim 3 割の計算時間の改善がはかれました. 他の step の計算時間も改善されています (表 2 参照).

7.6 参考文献[1]

Erwin H. Bareiss, ”Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elim-ination”, Math. Comp., 22(103)[2]$\mathrm{p}\mathrm{p}.565-578(1968)$.

Erwin H. Bareiss, ”Computational Solutions of Matrix Problems Over an Integral Do-

main”, J.Inst.Maths Applics, 10, $\mathrm{p}\mathrm{p}68-104(1972)$ .

[3]

K.Geddes, S.Czapor and G.Labahn: ” Algorithms for Computer Algebra”, Kluwer Aca-

demic Publishers, 1992.
[4]

木田祐司, 牧野潔夫: ” UBASIC によるコンピュータ整数論”, 日本評論社.

[5]

森本光生, 木田祐司, 小林美千代: ” 円分数の素因数分解 (その 3)”, 上智大学数学講究録.

[6]

長嶋秀世: ” 数値計算法”, 槙書店.

7.7 参考以下に実験結果のデータを –部示す. なお各データとも骨内の数字は, 特に断りのないかぎり計算時

time $(\sec)$ ]” である. また表 5 の間であり, ” [CPU time $(\sec)$ ]$+$目のフィボナッチ数, つまりのことである.

$f_{2}$[$\mathrm{G}\mathrm{C}$$f_{p1}$とは, 1 番目の素数番

73