有理数計算による対称行列の固有値問題における特性多項式の
因子探し
寒川 光, 早稲田大学理工学術院
平成 27 年 8 月 7 日
目次
1
対称行列の固有値の多重度解析アルゴリズム
4
1.1
4
5
1.2
1.3
2
特性多項式を求めるアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1
ヘッセンベルグ変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2
1.1.3
重複固有値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
フロベニウス変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
整数行列の数値例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1
内部 2 × 2 メッシュ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2
1.2.3
1.2.4
分割数と固有値の多重度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
6
8
8
有理数行列の整数行列変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
11
12
1.2.5 包含（inclusion） . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
プログラムの概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
13
根と係数の関係公式の応用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
テストプログラム sehfv1.cpp
2.1
2.2
2.3
15
処理の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
根と係数の関係公式と組合せのプログラミング . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 根と係数の関係公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 組合せのプログラミング（十進 BASIC による準備） . . . . . . . . . . . . . . . . . . . .
17
17
18
2.2.3
2.2.4
整数性判定と有理数での整数の抽出 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
22
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
多項式の格納形式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
22
2.3.2 多項式の除算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
包含チェックと 1 次の因子探し . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
24
2.4.1
2.4.2
2.4.3
1 次の因子探し . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
25
27
2 次以上の因子探し . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
多項式の格納と除算
2.3.1
2.4
2.5
根と係数の公式と組合せのプログラミング（C++） . . . . . . . . . . . . . . . . . . . . .
包含チェック . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ホーナー法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
3
4
有理行列の解析プログラム sehfv2.cpp
31
3.1
31
32
3.1.1
固有値の区間での取り扱い . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2
3.1.3
関数多重定義による askinty 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.4
3.1.5
精度要求の判定値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
整数性の判定値 epsint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
33
包含の判定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
34
3.2
3.3
3.1.6 精度改良 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
計算例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
零固有値のある問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
38
39
3.4
チェックポイントとリスタート . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
倍精度浮動小数点行列の解析プログラム sehfv3int.cpp
42
4.1
4.2
区間算術演算の応用
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
倍精度浮動小数点数行列の読込み . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
45
4.3
4.2.1 浮動小数点行列のファイルへの出力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 浮動小数点行列のファイルからの入力 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
チェックポイントとリスタート . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
47
48
4.4
近接根の分離 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
49
4.5
5
有理算術演算による精度改良の sehfv1.cpp に対する変更点 . . . . . . . . . . . . . . . . . . . . .
4.4.1
ループ構成の変更 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2
4.4.3
近接する固有値の分離（2 根の場合） . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.4
4.4.5
零を含むグループに対する包含 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
近接する固有値の分離（3 根の場合） . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
51
精度改良の sehfv2.cpp に対する変更点 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
52
4.4.6 その他の sehfv2.cpp に対する変更点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
計算例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 最適化オプションなしで CT2D を作成した場合 . . . . . . . . . . . . . . . . . . . . . . .
53
53
54
4.5.2
4.5.3
A = K を選択した場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A = M −1 K を選択した場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
57
4.5.4
4.5.5
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = S −1 KS −1 とした場合 . . .
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = K とした場合 . . . . . . .
59
61
4.5.6
4.5.7
4.5.8
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = M −1 K とした場合 . . . .
多重精度算術演算のサポートに関して . . . . . . . . . . . . . . . . . . . . . . . . . . . .
デバッグ例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
61
61
高速化チューニング
5.1
因子探しループでの反復のスキップ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
64
5.2
5.3
整数性判定の高速化
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
66
インライン展開 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
根と係数の関係公式
5.3.1
5.3.2
5.4
一時配列の利用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
69
5.3.3 再帰呼出し . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.4 計測結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 分法の桁数削減 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
76
77
5.4.1
78
最良近似分数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
5.4.2
プロファイルの取得 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
5.4.3
5.4.4
素因数分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
80
中点の有理数の分子の移動による桁数削減 . . . . . . . . . . . . . . . . . . . . . . . . . .
A ヤコビ法
83
B シェルプログラム
86
B.1 検証用データの作成用シェルプログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 実行して検証するためのシェルプログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
概 要
有理算術演算は，
「紙と鉛筆」による計算の延長で，扱う計算対象の数値が有理数の範囲であれば誤差なしで
計算できる．このため，数値計算で誤差理論を学ぶ前の学生に，計算機による正確な計算を体験させることがで
きる．十進 BASIC の有理数モードはこれにあたる．しかし，扱う数値の範囲が無理数に広がると，区間演算な
どの工夫が必要になったり，あるいは計算を諦めることになる．現在開発中の「有理数計算プログラミング環
境」は，浮動小数点計算の誤差解析や数式処理との融合，数学教育などを目的として，
（C++に有理数の変数型
rational を追加することで）有理算術演算と浮動小数点算術演算をひとつのプログラムで混在して使用できる
（十進 BASIC では混在使用ができない）．現在は便利な解析システムやツール類が広く行きわたり，計算機を
使用して種々の解析を行うにも，汎用プログラミング言語を使用する必要がほとんどなくなった．そこで大学の
初年度の，力学をカリキュラムに持たない学科でも理解できる HPC プログラミングの例題「対称行列の固有値
の多重度解析」を作成した．対称行列の特性多項式を，有理算術演算によって，無平方多項式の形で正確に求め
る．これに浮動小数点計算によって得られた固有値の近似値を組み合わせて当てはめることで，特性多項式の
因数分解された形を作り上げる．前半は，ヘッセンベルグ変換とフロベニウス変換を使用するが，これらは基本
変換を使用するだけなので，線形代数の授業で習う内容と考える．後半の当てはめは，根と係数の関係公式なの
で，高校時代から馴染み深い．プログラミング技術に長けていれば，実質的に数式処理システムの因数分解を代
替できる．この演習を通して，数学知識をプログラミング技法で生かせれば何ができるかを知り，プログラムの
高速化の楽しみを体験できる内容を目指した．この例題はたんに教育用の題材としてだけでなく，有理数計算プ
ログラミング環境に，多重精度算術演算の必要性や，有理数を形成する分母・分子と gcd アルゴリズムを調査
する上で有意義な示唆を与えた．
3
86
86
1
対称行列の固有値の多重度解析アルゴリズム
固有値問題は，次の線形同次方程式の n 個の未知数 λ を決定する．
Ax = λx.
(1)
本稿で扱う行列 A は対称で，その要素は有理数とする1 ．要素が有理数の行列は，全要素の分母の最小公倍数
（Least Common Multiplyer，LCM）を掛けることによって整数行列に変換できる．固有方程式 (1) の係数を整
数に変換すれば，整数を係数行列の要素とする対称行列の固有値問題となる．方程式 (1) は次式が成立したとき，
非自明解をもつ．
det(A − λI) = 0.
(2)
式 (2) の左辺の行列式は展開されて次式を得る．
(−1)n λn + pn−1 λn−1 + pn−2 λn−2 + · · · + p0 = 0
(3)
これを係数行列 A の特性方程式（characteristic equation）と呼ぶ2 ．
本資料で扱う例題のポイントを述べる．
• 式 (3) の左辺の多項式の全係数は，固有方程式 (1) の行列要素 aij に対する乗算と加減算だけで得られるの
で，行列要素が整数なら，特性方程式 (3) の係数 pi も整数である．
• 多重度が高い（重根がある）場合は，ヘッセンベルグ（Hessenberg）変換の副対角項が零になる．
• 多重度が高い場合は，ヘッセンベルグ行列は複数のヘッセンベルグ小行列に分離され，各小行列をフロベニ
ウス（Frobenius）変換して得られる特性方程式は重根を持たない（特性多項式が無平方（squarefree）で
ある）．
• 特性方程式の最高次の係数は 1 か −1 なので，係数の数と近似固有値の数が等しい．整数行列なら「根と
係数の関係公式」を用いて，近似固有値から「整数になるはずの係数」を求められる3 ．計算誤差を含む近
似固有値から根と係数の関係公式で求める係数が，許容範囲内に整数があるかどうかの判定を「整数性」の
判定と呼ぶことにする．
• 「整数性」の判定は絶対誤差評価で行うので，整数性判定には近似固有値の精度が十分になくてはならな
い．本稿で用いるアルゴリズムは，近似固有値を「包含」（inclusion）する区間を見つけ，2 分法で精度改
良する．これらは有理算術演算による正確な計算で可能である．
これらの項目のうち，未説明のものの数学的説明を本章で行い，2 章以降で例題について解説する．
1.1
特性多項式を求めるアルゴリズム
A をヘッセンベルグ行列 H に変換すると，多重度が高い固有値が存在する場合は，副対角項に零要素が現れ
る．その後の計算は各ヘッセンベルグ小行列に対して独立に行える．
1 計算機では有限のビット数で実数を近似するので，実数といっても，数学的には有理数である．
˛
˛
˛
˛ a11 − λ
˛
a12
a13
˛
˛
˛
a12
˛=
˛ = λ2 − (a11 + a22 )λ + a11 a22 − a12 a21 ，n = 3 の場合 ˛
a21
a22 − λ
a23
˛
˛
˛
a22 − λ
˛
a31
a32
a33 − λ ˛
−λ3 + (a11 + a22 + a33 )λ2 − (a12 a21 + a23 a32 + a13 a31 )λ + det(A) になる．n 次と n − 1 次の項は対角項の積の項から出てくる．
pn = (−1)n ，pn−1 = (−1)n−1 (a11 + a22 + · · · + ann )，p0 = det(A) である．特異行列 det(A) = 0 は特性多項式の定数項が零で，零
固有値をもつ．
3 有理数係数の特性多項式を通分して，整数係数の多項式を作ると，近似固有値の数よりも決めなくてはならない係数の数が 1 つ多くな
るので，根と係数の関係公式が使えない．
2n
˛
˛ a −λ
= 2 の場合 ˛˛ 11
a21
4
ヘッセンベルグ変換
1.1.1
次数 n の対称行列 A に基本変換行列 R を (n − 2) 回，相似変換（similarity transformation）の形で掛ける4 ．
−1
H = Rn−2 · · · R2 R1 AR1−1 R2−1 · · · Rn−2
.
(4)
浮動小数点計算では，このアルゴリズムは非対称行列に対して使用されることはあるが，ヘッセンベルグ行列が
非対称行列であるため，対称行列に対して用いられることはめったにない5 ．対称行列でも，計算が誤差なしで
行われる場合は，対称性保持の必要はない．これは有理数計算でアルゴリズム選択の重要な特長である．
次数 5 のヘッセンベルグ行列を示す．
⎛
h11
⎜
⎜ k2
⎜
H=⎜
⎜
⎜
⎝
h12
h22
h13
h23
h14
h24
k3
h33
k4
h34
h44
k5
⎞
h15
⎟
h25 ⎟
⎟
h35 ⎟
⎟.
⎟
h45 ⎠
h55
(5)
対角項の下の副対角項 hs+1,s を ks+1 で表す．ヘッセンベルグ行列への変換アルゴリズムは次のようになる [1,
p. 478];
第 s 段階までに，元の行列 A ≡ A1 は As になっており，その最初の (s − 1) 行および列は上ヘッセンベルグ
形である．第 s 段階は次の一連の操作からなる．
• 第 s 列の対角要素より下側で，絶対値最大の要素を捜す．それが零なら，次の 2 項目を飛ばし，この段階
は終わる．そうでないなら，その最大要素のあった行の番号を s と置く．
• 第 s 行と第 (s + 1) 行を交換する．これは軸選択の手順である．この交換を相似変換にするため，列 s と
列 (s + 1) も交換する．
• i = s + 2, s + 3, . . . , n に対して，乗数
ri,s+1 =
ais
as+1,s
(6)
を計算する．行 (s + 1) に ri,s+1 を掛けて行 i から引く．この消去を相似変換にするため，列 i に ri,s+1
を掛けて列 (s + 1) に加える6 ．第 s 列の要素 ais を消去するために R の s + 1 列目に非零項を使うとこ
ろが，ガウス消去法で用いる基本変換とは異なる．
この段階が全体で (n − 2) 回必要である．
関数 elmhes は，数値計算分野では，非対称行列用に固有値計算として一般的に使用されるものに類似してお
り，Numerical Recipes には ELMHES サブルーチンとして紹介されている．
4 Ax
= λx の両辺に H −1 を掛けると H −1 Ax = λH −1 x になるが，H −1 H = I = HH −1 なので
H −1 A HH −1 x = λH −1 x
| {z }
=I
(H −1 AH)H −1 x
λH −1 x
となる．H −1 AH
すなわち
=
= B ，H −1 x = y とおけば，By = λy の固有方程式になる．λ は同じなので，
−1
が掛けられる．
相似変換により固有値は保存され，固有ベクトルは H
5 実数の固有値を持つはずの問題でも，非対称行列の固有値問題を，浮動小数点計算で解くと（固有値計算の反復計算で使用するアルゴ
リズムにもよるが）
，丸め誤差のために虚数部分に計算誤差が現れるので，
complex
0
10
1
10
1 型で解く必要がある．
0
1
0
0 0
1
0
0 0
a11 a12 a13 a14
h11 h12 a13 a14
B
CB
B
B
1
0 0 C
1
0 0 C
h22 a23 a24 C
6B 0
C
C B a21 a22 a23 a24 C B 0
C B k2
@ 0 −r32 1 0 A @ a31 a32 a33 a34 A @ 0 r32 1 0 A = @ 0
a32 a33 a34 A
0 −r42 0 1
0 r42 0 1
a41 a42 a43 a44
0
a42 a43 a44
5
重複固有値
1.1.2
行列が縮重（derogatory）していて，複数の重複する固有値が存在する場合は，ヘッセンベルグ行列の副対角
項に零要素が現れる．この現象は，対称行列を対称 3 重対角行列に変換すると，非対角項に零要素が現れる現象
に対応している．Wilkinson の古典的教科書から引用する．
行列 A が 3 重対角化された状態を考える．
⎛
SAS −1
α1
⎜
⎜ β
⎜ 1
=⎜
⎜
⎝
⎞
β1
α2
..
.
..
.
..
.
βn−1
βn−1
⎟
⎟
⎟
⎟.
⎟
⎠
αn
対称行列が多重度 k を持つ場合，ギブンス法かハウスホルダー法によって 3 重対角化された行列は，
少なくとも (k − 1) 個の副対角項を持つ．一方，ギブンス法かハウスホルダー法によって 3 重対角化
された行列が零副対角項を持っても，それは多重度の高い固有値の存在を意味するとは限らない．
この引用に続く証明は，スツルム列の性質を使用し，非線形方程式を解く漸化式の計算に沿った方法でなされて
いる [2, p. 300]．ヘッセンベルグ行列のスツルム列の計算式は 3 重対角行列の場合よりも複雑なので，ここでは
ヘッセンベルグ行列の階数から証明する．
Proof. 対称行列 A が縮重していると仮定する．この場合，多重度の高い固有値を λj とすると，この固有値
に対応するジョルダン細胞が複数存在する．A を変換して得られるヘッセンベルグ行列 H は， λj に対応す
る固有ベクトルを複数もつ．これは (H − λj I) の階数が (n − 2) 以下であることを意味する．もし，副対角項
ki = 0, i = 2, · · · , n とするなら，行列 (H − λj I) の第 1 列から 第 (n − 1) 列は線形独立なので [2, p. 406]，
(H − λj I) の階数が (n − 1) 以上であることを意味する．これは仮定に反する．したがって，βj = 0 となる行列
A に対しては，必ず kj+1 = 0 となる．
対称行列の 3 重対角化の場合と同様，ヘッセンベルグ行列が零副対角項を持っても，それは多重度の高い固有
値の存在を意味するとは限らない．
なお，浮動小数点計算で 3 重対角化を行う場合は，βj = 0 の判定は丸め誤差のために，閾値の設定に苦慮す
るが，計算誤差を含まない有理数計算では「イコール零」でチェックできる．浮動小数点演算による「丸め誤差
を考慮した零判定」のようなプログラミングの難しさはない．
kr+1 = 0 の場合は次のように書ける．
H=
Hr
0
Br
Hn−r
,
Hr は次数 r，Hn−r は次数 (n − r) で，どちらも上ヘッセンベルグ形である．行列式 (2) は次のように計算で
きる．
det(H − λI) = det(Hr − λI) det(Hn−r − λI),
つまり，kr+1 = 0 が検出されれば，後続の計算は 2 つの小行列に対して独立に進められる．
1.1.3
フロベニウス変換
特性多項式を求めるために，上ヘッセンベルグ行列 Hr または Hn−r をフロベニウス標準形 F に変換する．
この標準形は，有理数行列の要素間の四則演算（有理演算）のみで計算可能なので，有理標準系とも呼ばれる7 ．
7 ジョルダン標準系は，一般に，体の拡大を必要とする
[2, p. 16]．
6
この段落で，行列の添字の表記を変更する．Hr に付けられた “r” は行列のサイズを表していた．ここではサ
イズは r か n − r で，省略し，替わりに消去の段階を表記する．消去は (n − 1) 段の主たるステップからなる．
元の行列は H ≡ H1 とする．n = 5, s = 3, の場合には R3 H3 R3−1 を計算する．
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎞⎛
−r14
1
1
1
0
⎟⎜
⎟ ⎜ k2
⎟⎜
⎟⎜
⎟⎜
⎟⎜
⎠⎝
−r24
−r34
1
0
h13
h14
h15
0
k3
h23
h33
h24
h34
h25
h35
k4
h44
k5
h45
h55
1
⎞⎛
⎟⎜
⎟⎜
⎟⎜
⎟⎜
⎟⎜
⎟⎜
⎠⎝
1
⎞
r14
1
⎟
⎟
⎟
⎟
⎟
⎟
⎠
r24
1 r34
1
1
i が 1 から s について，次の操作を行う [2, p. 406];
• ri,s+1 =
his
を求める．
ks+1
• 行 i から，ri,s+1 × 行 (s + 1) を引く．
• i が 1 から s について，ki+1 ri,s+1 を hi,s+1 に加える．
第 3 段の終了時点では次のようになる．
⎛
0 0 0
⎜
⎜ k2 0 0
⎜
⎜
k3 0
⎜
⎜
k4
⎝
h14
h24
h34
h44
k5
⎞
h15
⎟
h25 ⎟
⎟
h35 ⎟
⎟,
⎟
h45 ⎠
h55
hij = hij + ki+1 rij .
この変換操作を (n − 1) 回繰り返し，相似変換を完結させる．行列 X が得られる．
⎛
0
⎜
⎜ k2
⎜
X=⎜
⎜
⎜
⎝
0
0
0
0
0
0
k3
0
k4
0
0
k5
⎞
x0
⎟
x1 ⎟
⎟
x2 ⎟
⎟.
⎟
x3 ⎠
x4
フロベニウス行列は対角行列を相似変換 F = D−1 XD の形で掛けて得られる．ここに対角行列 D は d1 =
1, d2 = k2 , d3 = k2 k3 , d4 = k2 k3 k4 ,
⎛
0
⎜
⎜ 1
⎜
F =⎜
⎜
⎜
⎝
d5 = k2 k3 k4 k5 , である8 ．
⎞
0 0 0 p0
⎟
⎞
⎛
0 0 0 p1 ⎟
n−2−i
⎟
pi = ⎝
kn−j ⎠ xi .
1 0 0 p2 ⎟
⎟,
⎟
j=0
1 0 p3 ⎠
1 p4
(7)
次数 5 の特性多項式が −λ5 + p4 λ4 + p3 λ3 + p2 λ2 + p1 λ + p0 として得られる．最高次の係数は “−1” であるが
n が偶数の場合は “+1” として，全係数 pi の符号を反転させる．次数 n の特性方程式は式 (3) である．
8 D を右から掛ける操作は，1 列目はそのまま，2 列目は k 倍，3 列目は k k 倍，. . .，n 列目は k k · · · k 倍する．D −1 を左から掛け
n
2
2 3
2 3
る操作は，1 行目はそのまま，2 行目は k2 で割り，3 列目は k2 k3 で割り，. . .，n 列目は k2 k3 · · · kn で割る．したがって p0 = k2 k3 k4 k5 x0 ，
p1 = k3 k4 k5 x1 ，p2 = k4 k5 x2 ，p3 = k5 x3 ，p4 = x4 になる．なお，式 (7) の pi は教科書の (52.6) 式とは異なる [2, p. 407]．
7
100◦ C
13
9
100◦ C
？
14
10
15
16
11
43 44 45 46 47 48 49
12
100◦ C
5
6
7
8
36
42
29
35
22
28
15
21
8
40◦ C
(a) boundary condition
1
2
3
4
14
1
(b) 3 × 3 mesh, 2 × 2 matrix
2
3
4
5
6
7
(b) 6 × 6 mesh, 5 × 5 matrix
図 1: 四角柱の内部の温度分布
ヘッセンベルグ行列を経由してフロベニウス行列を求める相似変換は，1 つの行列 T にまとめられ，これを用
いて相似変換の形で F = T AT −1 と記述できる．多重度が高い場合は F はブロック対角行列となる．
⎛
T AT
−1
⎜
⎜
=⎜
⎜
⎝
⎞
F1
⎟
⎟
⎟.
⎟
⎠
F2
..
.
(8)
Fm
特性多項式は各フロベニウス小行列 Fi から得られる特性多項式の積になる．
p(λ) =
m
pi (λ)
(9)
i=1
各特性多項式 pi (λ) は，重根をもたないので無平方（squarefree）である．
1.2
整数行列の数値例
正方形の断面をもつ四角柱で，表面の温度を与えられた場合の熱伝導問題を，断面領域で考える [3, p. 6]．図
1 に，内側のメッシュが 2 × 2 節点と 5 × 5 節点の場合を示す．1 つの差分要素を中図と右図に示した．熱量は上
下左右の要素から，温度勾配と熱伝導係数に比例して流入する，5 点差分スキームを用いる．解析領域を m × m
に分割すると，内部は (m − 1) × (m − 1) にメッシュ切りされ，行列のサイズは n = (m − 1) × (m − 1) になる．
未知数は内側の節点にあるので，行列の行・列番号は，境界節点の番号を詰めて振られる．3 × 3 のメッシュで
は内側は 2 × 2 になるので，行列の行・列番号では節点 6 が最初の行になり，節点 7 が 2 番目，節点 10 が 3 番
目，節点 11 が 4 番目になる．
1.2.1
内部 2 × 2 メッシュ
熱伝導係数を 1 にすると係数行列は次のようになる．
⎛
4 −1
⎜
⎜ −1 4
A=⎜
⎜ −1 0
⎝
0 −1
8
−1
0
4
−1
⎞
0
⎟
−1 ⎟
⎟
−1 ⎟
⎠
4
(10)
⎛
⎞
⎛
⎞
1 0 0 0
1 0 0 0
⎜
⎟
⎜
⎟
⎜ 0 1 0 0 ⎟
⎜ 0 1 0 0 ⎟
−1
⎜
⎟
⎜
⎟
第 1 段の消去は a3,1 = −1 を消去する．R1 = ⎜
⎟ と R1 = ⎜ 0 1 1 0 ⎟ を用いて
⎝ 0 −1 1 0 ⎠
⎝
⎠
0 0 0 1
0 0 0 1
⎞
⎛
4 −2 −1 0
⎟
⎜
⎜
−1
4
0 −1 ⎟
−1
⎟ となる．第 2 段の消去は， a4,2 と a4,3 を交換すると a4,3 = 0
⎜
R1 AR1 により A = ⎜
0
4
0 ⎟
⎠
⎝ 0
0 −2 −1 4
であるため行われない9 ．
行列は物理的な形状の対称性のために縮退しており，ヘッセンベルグ行列は次のようになる．
⎛
⎞
4 −2 0
0
⎜
⎟
⎜ −1 4 −1 0 ⎟
⎜
⎟
(11)
H=⎜
−2 4 −1 ⎟
⎝
⎠
0
4
k4 = 0 なので，H は H3 と H1 に分かれる．大きいほうの行列から
⎛
⎞
⎛
4 −2 0
0
⎜
⎟
⎜
H3 = ⎝ −1 4 −1 ⎠ , X3 = ⎝ −1
−2
4
0
24
⎞
⎟
0 22 ⎠
−2 12
(12)
が得られ，特性多項式が得られて，それは整数の範囲で因数分解できる．
p3 (λ) = 48 − 44λ + 12λ2 − λ3 = (2 − λ)(4 − λ)(6 − λ).
(13)
小さいほうの小行列 H1 から特性多項式 p1 (λ) = −λ + 4 が得られる．
計算は H1 から始める．p3 (λ) は p1 (λ) で割り切れ，p3 (λ) ÷ p1 (λ) = λ2 − 8λ + 12 の因子を探す．ここで特
性多項式がその定義式 (2) から λ4 − 16λ3 + 92λ2 − 224λ + 192 と得られても，それが (2 − λ)(4 − λ)2 (6 − λ)
と 2 乗の項を含む形に因数分解できるまでは，重根の存在を言えない．重根の難問を解決するのはヘッセンベル
グ変換である．これは有理算術演算の正確な計算によってもたらされるのであって，浮動小数点演算では実現で
きない．
1.2.2
分割数と固有値の多重度
メッシュ分割を細かくすると行列は大きくなってゆく．図 1 の右のように，内部を 5 × 5 のメッシュに切ると，
行列のサイズは n = 25 になる．ヘッセンベルグ行列は 5 つの小行列 H13 , H9 , H1 , H1 , H1 に分かれる．固有値
の分布を表 1 に示す．多重度 5（λ11 = · · · = λ15 ）が存在する．11 個の固有値が有理数で，14 個の固有値が無
理数である10 ．
5 × 5 分割の場合は，小行列の特性多項式が次の小行列の特性多項式で割り切れる（H9 の特性多項式が H13
の特性多項式を割り切る）が，これは一般的には成立しない．4 × 4 の場合を表 2 に示す．2 つの H3 の特性多
項式は異なっているので，割り切る関係は成立しない11 ．
ヘッセンベルグ行列が Ha ，Hb などのヘッセンベルグ小行列の積に分かれる状態を表 3 に示した．表では 4 − λ
の 1 次式の因子が複数現れる場合を “3*1” のように表した．すべてのケースで λ = 4 の固有値が分割数に一致す
9 HesFrb.BAS を実行すると分かりやすい．2 進モードでは x87 FPU の命令セットの倍精度浮動小数点演算で実行される．5 × 5 分割
では特性多項式の係数は精度が不足して正しく計算されない．
10 固有値が有理数か無理数であるかは，対応する固有ベクトルを求める場合は重要になる．有理数に対応する固有値は正確に計算できる
ので，無理数に対応する固有ベクトルよりも先に計算するほうがよい場合がある．
11 このため，プログラムを作成する場合には注意が必要になる．
9
表 1: 5 × 5 の場合の固有値
λi
λ3
.535898384862
1.26794919243
exact λ
√
2(2 − 3)
√
3− 3
λ4
λ5
λ6
2.0
2.26794919243
2
√
4− 3
λ7
λ9
λ8
λ10
3.0
3.26794919243
3
√
5− 3
λ11
λ16
λ18
λ12
λ17
λ19
4.0
4.73205080757
5.0
λ20
λ22
λ21
5.73205080767
6.0
λ23
λ25
λ24
6.73205080757
7.46410161514
4
√
3+ 3
5
√
4+ 3
6
√
5+ 3
√
2(2 + 3)
H13
H9
λ1
λ2
approx. ei
λ13
λ14
H1
H1
λ15
H1
表 2: 4 × 4 の場合の固有値
λi
λ1
λ2
λ4
λ5
λ7
approx. ei
λ3
λ8
λ11
λ13
λ6
λ9
λ15
H9
H3
exact λ
√
3− 5
√
4− 5
√
5− 5
3.0
4.0
3
4
5.0
5.2360679775
5
√
3+ 5
√
4+ 5
√
5+ 5
λ10
λ12
λ14
λ16
0.7639320225
1.7639320225
2.7639320225
6.2360679775
7.2360679775
H3
H1
表 3: 分割数 m と小行列サイズと因子多項式の最高次数 d
m
2
3
4
5
6
7
8
9
10
11
12
13
n
4
9
16
25
36
49
64
81
100
121
144
169
Ha
3
5
9
13
19
25
33
41
51
55
73
85
Hb
Hc
Hd
1
3
1
3
3
1
9
3*1
13
4*1
19
5*1
25
6*1
33
7*1
41
8*1
37
15
7
61
10*1
73
11*1
6
?
He
d
7*1
2
2
2
2
3
4
3
10
4
5
4
る多重度を持つ．また分割数 m の最大の小行列 Ha のサイズと m + 1 の 2 番目の小行列 Hb のサイズが，m = 4
と m = 11 を例外として一致している．なお，本資料で使用するプログラムでは m = 13 の場合は組合せ数が大
きすぎて，long long 型の変数の最大値を超えるため計算できない（後述）．
1.2.3
根と係数の関係公式の応用
「有理数計算プログラミング環境」
倍精度計算で得られた固有値の近似値 ei を小数で表 1 の第 2 列に示した．
は，有理数と倍精度浮動小数点数を，1 つのプログラムで使用できる点が十進 BASIC とは異なる．行列 A が整
数行列のときは特性方程式の係数も整数になる．特性多項式が次数の低い多項式の積に分解される場合，浮動小
数点計算で得られた固有値を用いて，因数分解ルーチンを簡単なフィルタープログラムで置き換えられそうに見
える．
特性多項式が 2 次式に因数分解されるとする．モニック 2 次方程式 “x2 + rx + s = 0” に対する根と係数の関
係公式 “α + β = −r” と “αβ = s” を利用する．2 つの候補 ei と ej を近似固有値から選び，r = −(ei + ej ) と
s = ei × ej を作成し，
Floor(|r|) − |r| or Ceil(|r|) − |r|
and
Floor(|s|) − |s| or Ceil(|s|) − |s|
が許容範囲に入るかどうかを調べる．これを整数性判定ということにする12 ．範囲内であれば多項式の除算 pr (λ)÷
(x2 + rx + s) を試み，剰余が零ならこの 2 次多項式は特性多項式の因子（divisor）であることが分かる．すべ
ての i と j の組み合わせをテストする．
たとえば表 1 の 2 つめと 8 つめを選んで α = 1.2679，β = 4.7321 と 5 桁で計算すると，r = 6.0000，s =
5.9998 · · · となるので，許容範囲が 0.0002 で整数が係数の 2 次多項式 x2 −6x+6 が得られる．つまり “x2 −6x+6 =
√
√
(x − (3 − 3))(x − (3 + 3))” のようなペアを見つければよい13 ．
ここでの例題は，表 3 に示したように，内部分割が 5 × 5 以下は，整数の範囲での因数分解が 2 次以下の因子
に分解される．6 × 6 の場合は 3 次式を含む．この場合はモニック 3 次方程式 “x3 + rx2 + sx + t” に対する公
式 “α + β + γ = −r,
αβ + βγ + γα = s,
αβγ = −t” によって 3 つ子組を探す．36 個の近似固有値から 3 つ
選んで有効数字 5 桁，α = 0.95108，β = 5.3569，γ = 5.6920 で計算すると，r = −11.99998，s = 40.9998 · · · ，
t = −28.9998 · · · が得られる14 ．
7 × 7 の場合は 4 次式を含む．pr (λ) から次数の低い多項式を分離することのメリットは，正確な計算を続行
できる点にある．一度丸め誤差の混入を許すと，解の精度を保証するために，真の解と近似解の距離に気をくば
らなければならなくなる．7 × 7 の場合，陰的ではあるが，解はフェラーリの公式から
5
5
√
8 − 22
22 + 8 √
− 2 + 4, ±
+ 2 + 4.
λ=±
2
2
を正確な解と考えることができる．この形のほうが，解を挟む区間の上限と下限を示すよりも好ましい場合があ
る．また，どの近似固有値がどの因子に繋がるかが分かるので，特性多項式の構成を考えることができる（固有
値を数値として何 10 桁求めてみても，これは分からない）．表 3 の最下段から，9 分割以下，および 11 分割で
は 4 次多項式以下に分解されるので，べき根による公式を正確な解と考えれば，正確な固有値を得ることができ
る．10 分割は 5 次多項式以下に，12 分割は 6 次多項式以下に分解される．
12 自然数が素数であることを英語では primality というが，数がある範囲内に整数をもつかどうかを，本稿では「整数性」という．この
語に対する英語が見つからなかったので，primality に類似した造語 integer に ty を付して短縮した inty をコメントや関数名に使用する．
13 {x − (p + √q)}{x − (p − √q)} = x2 − 2px + p2 − q
14 x3 − 12x + 41x − 29 = 0 の 3 根を十進 BASIC のプログラム CCubicEqn.BAS で計算するともとの近似固有値に近い値が得られる．
このプログラムは十進 BASIC の複素数モードを使用している．3 次方程式の根がすべて実数の場合でも，途中に虚数を用いないと実根は
得られない．
11
1.2.4
有理数行列の整数行列変換
代数方程式 a0 xn + · · · an−1 x + an = 0 の解を α1 , · · · , αn とするとき，その定数 k 倍 kα1 , · · · , kαn を解とす
る方程式は，
a0 xn + a1 kxn−1 + a2 k 2 xn−2 + · · · an−1 k n−1 x + an k n = 0
(14)
となる．行列要素が有理数の行列の特性多項式は有理数係数である．この多項式を通分して，整数係数の多項式
に変換すると，特性多項式の係数は，pi が pi k i に変わる．しかし最高次の係数が 1 でなくなるので，未知数が
ひとつ増えて n + 1 個になり，係数を決めることができなくなる．そこで特性多項式を通分するのではなく，も
との行列のほうを通分（分母の最小公倍数倍）して，最高次の係数が 1 の整数係数の特性多項式を求めることに
する．
熱伝導係数を 0.1 にすると，式 (10) の係数行列は次のように変わる．
⎞
⎛
0.4 −0.1 −0.1
0
⎟
⎜
⎜ −0.1 0.4
0
−0.1 ⎟
⎟
⎜
A=⎜
0
0.4 −0.1 ⎟
⎠
⎝ −0.1
0
−0.1 −0.1 0.4
(15)
全要素の分母の LCM は 10 であるが，これを掛ければ整数行列に変換できる．固有方程式 (1) の係数行列 A を
10 倍した行列 A = 10A の固有値は，元の行列の固有値の 10 倍である．この行列 A は式 (10) の A に一致
する．
0.1 は 10 進数だから循環小数にならなかったが，2 進数では循環小数になる．数値計算するには，どこかで丸
めなくてはならない．この丸め誤差の影響を受けて，倍精度浮動小数点数を有理数に変換すると，たいていの場
合は分母と分子のそれぞれが 16 桁程度の数になる．例えば 0.4 と 0.1 は次の数になる．
0.4 3, 602, 879, 701, 896, 397
,
9, 007, 199, 254, 740, 992
0.1 3, 602, 879, 701, 896, 397
36, 028, 797, 018, 963, 968
(16)
0.4 は「分母の 4 倍に 2 を加えると分子の 10 倍に一致する」有理数になる．0.1 は「分子の 10 倍から 2 を引く
と分母に一致する」有理数になる．15 または 16 桁の数に対して 2 だけずれている．これが倍精度浮動小数点数
の丸め誤差である．倍精度浮動小数点数の有意数（significand）は 53 ビットあり，丸め誤差は最後の桁に入る
ので， 253 1016 より，分母と分子の桁数は 16 桁程度の数になる．
熱伝導係数を 10 進で 0.1 にした，式 (15) の係数行列は，2 進では 0.4 と 0.1 に式 (16) の有理数を入れた行列
に変わる．この行列の分母の最小公倍数 LCM = 36, 028, 797, 018, 963, 968 倍した整数行列は次のようになる．
⎞
⎛
14, 411, 518, 807, 585, 588
⎜ −3, 602, 879, 701, 896, 397
A =⎝
−3, 602, 879, 701, 896, 397
0
−3, 602, 879, 701, 896, 397
14, 411, 518, 807, 585, 588
0
−3, 602, 879, 701, 896, 397
−3, 602, 879, 701, 896, 397
0
14, 411, 518, 807, 585, 588
−3, 602, 879, 701, 896, 397
0
−3, 602, 879, 701, 896, 397 ⎟
−3, 602, 879, 701, 896, 397 ⎠
14, 411, 518, 807, 585, 588
(17)
式 (16) の 0.4 の分母の 4 倍が LCM に一致するので，対角項は 0.4 の分子の 4 倍で済み，17 桁で収まっている．
固有値も LCM 倍されるので，根と係数の関係公式から作成する多項式係数も大きな値になるため，近似固有
値もそれに合わせて高精度化しないと整数性を判定できない．この行列の多重度はやはり 2 である．分割数を増
やしていっても，多重度は整数行列の場合と同じに m である．このことは，対角項が非対角項のマイナス 4 倍
であることは，多重度 m の必要条件ではなさそうだ．また，特性多項式の因数分解後の次数も，表 3 の最下段
の形になるだろうか．これらの疑問は，実際に計算することによって確かめられるだろうか．
1.2.5
包含（inclusion）
有理算術演算を用いると，対象の数値が有理数の範囲では，素朴なプログラミングによって正確に行うこと
ができるが，無理数は区間演算の考え方を取り入れる必要がある．区間演算では，存在（existence），唯一性
12
（uniqueness），包含（inclusion）の 3 つを保証することで行う15 ．特性多項式が重根を 2 個と数えれば n 個の実
根をもつことは，係数行列が対称であることから（正確な計算では）保証される．残る 2 つは次のようにする．
まず，解を 1 つだけ挟む区間 (a, b) を見つける．固有値 λ1 から λn を昇順にとると，
· · · < λi−1 < a < λi < b < λi+1 < · · ·
(18)
を満たせば，λi は区間 (a, b) に存在し，そしてこの区間には他の根は存在しない．この状態を「包含」と呼ぶ．
包含は inclusion の訳で，区間 (a, b) 内に p(λ) = 0 の根が 1 つだけ存在する状態を指す．包含が成立したら，こ
の区間 (a, b) を縮小するアルゴリズムを用いて区間幅を，精度要求を満たすまで「縮小反復」する．因子探しの
前半は，n 次多項式に対して n 個の区間を見つけることで包含を達成する．後半で，近似値の桁数と見つけるべ
き因子の次数から精度要求を決定して，固有値の存在区間を，2 分法で縮小反復してから，候補を選び，整数性
を判定し，多項式の除算を行う．これによって式 (17) の A の特性多項式の因数分解された形を得られる．
1.3
プログラムの概要
本資料では sehfv1,2,3 の 3 つのプログラムを解説する．これらのプログラムは，生成された行列から有理算術
演算を用いて特性多項式を得るが，同時に多重度の高い固有値がある場合はその存在を知る．これは浮動小数点
演算ではできない．
テストプログラム sehfv1.cpp は，条件のよいデータに対して，正確な計算が最低限必要となる部分（行列変
換と多項式の除算）だけを有理算術演算を用い，整数性判定は倍精度浮動小数点演算で行う．次の版 sehfv2.cpp
は行列を有理数とするので，特性多項式の桁数が増え，倍精度では不足する固有値の精度を有理算術演算で補う．
3 番目の版 sehfv3.cpp は有理数で区間算術演算を導入して，倍精度浮動小数点演算で作成された物理的・工学的
に意味の付随した問題の行列の解析を行う．
テストプログラム 熱伝導行列を，m × m の内部分割で解析する．この行列は多重度 m の固有値をもつ．有理算
術演算で正確な計算を行うのは，対称行列のヘッセンベルグ変換，フロベニウス変換，多項式除算で，整
数性判定は倍精度浮動小数点演算で行う．これにより，整数行列に対してアルゴリズムが稼働するかどう
かを確かめる．固有値は 0.1 から 10 の間に分布するので，整数性は倍精度浮動小数点演算で間に合う．特
性多項式は表 3 の下段の d で示したように，低い次数の因子多項式に分解されるので，計算時間も短い．
有理行列の解析 行列の要素を有理数にして，全要素の分母の最小公倍数 LCM を求めて掛けることで整数行列
にする．行列の桁数が増えるので，特性多項式の係数の桁数も大きくなる．また固有値も LCM 倍される
ので大きい．したがって整数性判定は倍精度浮動小数点演算ではできない．倍精度演算で得られた近似固
有値を有理数型の変数に変換して格納し，これを有理算術演算により必要な精度まで精度改良する．改良
された近似固有値を用いて根と係数の関係公式を使用して，数 10 桁に及ぶ多項式係数を作成する．桁数は
増えるがヘッセンベルグ変換には時間はかからない．これは対角項と非対角項の 2 種類の数値があるだけ
で，行列の構造は同じだからである．特性多項式は表 3 のように，低い次数の因子多項式に分解されるの
で，計算時間も短い．
熱伝導行列の境界条件を外して，零固有値をもつグラフ・ラプラシアン行列（整数行列）も解析する．こ
の行列は多重度 4 の固有値をもつ．零固有値は特性多項式の定数項がゼロとして現れるので，簡単に見つ
けられる．因子多項式は高次になるので，組合わせの数が膨大になり，計算時間は長くなる．
倍精度浮動小数点行列を有理行列に変換して解析 数値シミュレーションで使用されるような構造解析プログラ
ム CT2D の行列の最も小さなものをファイル経由で読込み，有理行列に変換する．理論的には 3 つの零固
有値をもつ問題であるが，浮動小数点演算で得られた行列は，厳密な零固有値は持たない．このため必要
15 ドイツ語で
Existenz，Eindeutigkeit，Einschließung の 3 語の頭文字をとって 3E 法と呼ばれた．
13
な精度を決定する部分を改訂する必要がある．根と係数の関係公式で作成する多項式係数を，有理数を用
いた区間算術演算を導入することで，閾値の設定に曖昧性をもつ整数性判定値を不要にする（区間内部に
整数が存在するかどうかで判定できるので，閾値なしで整数性を厳密に判定できる）．解析結果は理論的に
零となるはずの固有値は零近傍の小さな値をもち，構造物の形状対称性に由来する重根は重根として存在
する（重根か近似固有値かは，倍精度浮動小数点演算では分からない）．係数行列の桁数は sehfv2.cpp の
行列の桁数と同程度であるが，ヘッセンベルグ行列の桁数は増え，ヘッセンベルグ変換とフロベニウス変
換に計算時間がかかる．これは行列要素が種々の数値で構成され，それらの要素の演算結果が，ヘッセン
ベルグ行列要素になるからである．構造解析プログラム CT2D を 32 ビットのアーキテクチャである x87
FPU 命令セットを使用する環境でコンパイルするとき，最適化オプションを “-O3” にすると，多重度の分
布が変わり，重根が消滅する．この解析を行う．
次章でテストプログラムの作成を解説する．3 章で有理数行列をを分析するプログラムを解説する．4 章で倍精
度浮動小数点行列を分析するプログラムを解説する．テストに使用する行列の特性多項式が，次数の低い因子多
項式の積に因数分解できる場合は，比較的短時間で因子探しが完了するが，グラフ・ラプラシアン行列では次数
が高くなるため，因子探しに時間がかかるようになる．また，倍精度浮動小数点行列ではヘッセンベルグ変換の
時間は長くなる（ヘッセンベルグ行列要素の桁数の増加のため）．プログラムが採用している組合せを網羅的に
探すアルゴリズムは指数時間アルゴリズムであるため，問題を少し大きくすると組合せ数が爆発して計算時間の
観点から計算不能になる．したがってチューニングの題材としては長く使用できる．5 章では，簡単なプログラ
ム変更で可能なチューニングを施す．有理算術演算は浮動小数点演算よりもチューニングのマージンが大きく，
簡単なプログラム変更で計算時間は 1 桁近く短縮できる．このチューニング作業を通して，有理算術演算に必要
なユーティリティ・ルーチンの洗い出しと，有理算術演算の浮動小数点演算とは異なるチューニングの着目点を
知る．
14
2
テストプログラム sehfv1.cpp
有理算術演算（rational arithmetic）は，有理数を対象に正確な計算を行うが，対象の数値が無理数になると，
要求される精度で妥協して誤差を含む近似解を求めるか，解の存在する区間を求めるのが普通の流れである．こ
こで行う「特性多項式の因子探し」はこの流れとは反対に，誤差を含む近似解から出発して正確な因子多項式を
探す．近似解を測定されたデータ，根と係数の関係公式を理論式と見なせば，測定データに合致する理論式を見
つけるデータ解析，あるいは逆問題の形になっている．プログラム sehfv1.cpp では素朴に倍精度浮動小数点演
算で得られた近似解から候補の組合せを選んで，根と係数の関係公式によって多項式係数を作る．その係数の整
数性判定を行い，小数点以下の小さな値を切り捨てるか切り上げるかして整数の係数を抽出し，整数係数の多項
式で特性多項式を割ってみる．抽出された係数と特性多項式は有理数型で扱い，除算を有理算術演算で正確に行
うので，ここに計算誤差は入らない．割り切れれば因子であるから，特性多項式を商多項式で置き換えて，特性
多項式の次数を減らしてゆく．プログラム sehfv1.cpp で，有理算術演算で正確な計算を行うのは，対称行列の
ヘッセンベルグ変換，フロベニウス変換，多項式除算で，整数性判定は倍精度浮動小数点演算で行うことでプロ
グラミングを簡単にし，整数行列に対してアルゴリズムが稼働することを確かめる．
2.1
処理の概要
メッシュ分割数 m を引数としてとる 250 行ほどの関数 eigvals が処理の流れを制御する．処理の概要は以下
のように進む．
• 係数行列 A を倍精度浮動小数点演算で 2 次元配列 matrix<double> a(n,n) に生成する（SetFtherm）
16
．有理数 2 次元配列 rational matrix ar(n+1,n+1) に cnvmat 関数によって有理数変換し，ヘッセン
ベルグ変換 A → H する（elmhes）17 ．
• 近似固有値を倍精度計算で求め（jacobi2）ソートする（bsort）18 ．近似固有値は 1 次元配列 est に格納
され，グループ化し（selectev），各グループの下限 eigvl と上限 eigvh，そのグループに属する近似固有
値の多重度 mult を抽出する．近似固有値に関する数値は倍精度浮動小数点数の配列に置かれる．グルー
プ化は，ファイル全体をスコープ（有効範囲）とする変数 eps によって ε = 10−9 で行う19 ．
• ヘッセンベルグ行列の副対角項 ki のゼロ要素を探してヘッセンベルグ小行列の数を得る（nmhes）．
5 分割の場合は，H13 ，H9 ，H1 ，H1 ，H1 の 5 つの小行列が得られるので，以下，各小行列に対して逆順に H1 ，
H1 ，H1 ，H9 ，H13 と進める．
• 小行列 Hk をフロベニウス変換して特性多項式 p(λ) を生成する（hesfrb）20 ．次数は変数 nsize に格納
する．
• 既知の因子多項式がある場合はこれを getpolynomial 関数で取り出し，特性多項式 p(λ) を割り（polydivchk），割り切れれば商多項式で p(λ) を置き換え，次数 nsize を小さくする21 ．
• 近似固有値の p(λ) に対する包含を調べる（chkinc）．p(λ) が 1 次式なら，p(a) · p(b) < 0 が条件になる．
n 次多項式の場合，p(a) · p(b) < 0 で，かつこのような区間が n 個見つかればよい．これらのケースで
16 プログラミング言語 C は，他の言語で「配列の配列」と呼んでいるものを，2 次元配列と呼ぶ場合があるが，Fortran や Pascal のよ
「有理数計算プログラミング環境」では 2 次元配列として使用するために matrix<double> 型を用意した．
うな 2 次元配列は扱えない．
17 cnvmat 関数は rational クラス，elmhes 関数は rblas クラスにある．
18 jacobi2 のアルゴリズムであるヤコビ法の説明を付録に記す．
19 10 ページの表 1 の 5 分割の例の場合は 13 グループに分かれる．
20 hesfrb 関数は rblas クラスにある．
21 表 1 の 5 分割の例で H を扱う場合には，特性多項式は 9 次であるが，既知の (4 − λ) で割り切れるため 8 次式 (λ8 − 32λ7 + 436λ6 −
9
3296λ5 + 15079λ4 − 42608λ3 + 72312λ2 − 67008λ + 25740) になる．
15
は，それぞれの区間に 1 つづつ根が存在するので，各区間で「インクルージョンが成立している」という．
sehfv1.cpp では，扱う行列の特質から，包含チェックを通過しなければ計算を停止する22 ．
• 近似固有値 ei の整数性を askinty 関数で調べ，整数らしければ ei を有理数型の整数 e として抽出し，1
次式 (−λ + e) で p(λ) を割り，割り切れれば putpolynomial 関数によって因子として保存し，商多項式
で p(λ) を置き換え，nsize をデクリメントする23 ．
• 次数 d = 2, 3, . . . について以下の 2 次以上の因子探しを行う．この処理は while ループで書かれているが，
内側に for ループによる「組合せループ」があり，因子が見つかると，内側のループを終える場合と，継
続する候補の組合わせを探す場合に分岐する（ goto 文で戻る場合がある）．組合せループの前処理として，
包含が成立している近似固有値だけを setecand 関数で evcan に詰めてセットし，その属するグループ添
字を配列 usev にセットする24 ．この時点で因子多項式の候補となる近似固有値の数が ncandi 個と決ま
り，作業用配列 idx を初期化し，組合せループの反復回数を得る（nn=comb(ncandi,d)）25 ．この後，組
合せループを最大 nn 回繰り返すが，その前に Top ラベルを置く（ここへ下から goto 文で飛ぶ）．
– nxtcmb 関数によって「組合せ」を作業用配列 idx に得る26 ．この配列に従って候補固有値を eigvl
d
から figvl に詰めて格納する．根と係数の関係公式の 1 次項
ej を求める．
j
– 1 次項の整数性だけを先行して調べ（先行判定），整数と見なせれば，vietaterm 関数によって他の
係数を cof[1] から cof[d] に求め，それらの整数性を調べる27 ．
– 全係数が整数と見なせれば，係数を有理数型の整数として抽出し，それらによって整多項式 v(λ) を
作り，これで p(λ) を割る．p(λ) = v(λ)q(λ) なら割り切れて，v(λ) を因子多項式とし保存し，商多
項式 q(λ) で p(λ) を置き換え，nsize を減らす28 ．特性多項式も取り出された因子多項式の候補も有
理数型なので，桁数が多くなっても除算は正確に行える．
見つかった近似固有値を候補から除外して，配列 evcan を詰め直し，ncandi と作業用配列 idx もリ
セットする．配列 idx の状態から，d 個の組合わせの因子探しが完了したか途中かを判定し，組合わ
せループを抜けるか，あるいは，comb 関数で nn を再設定して Top に帰る29 ．
倍精度浮動小数点演算で得られた近似固有値のグループ化を行う関数を示す．
22 5 分割の例で H を扱う場合には，グループ 2,4,5,6,8,9,10,12 で包含が成立するので，これらのグループの近似固有値を組み合わせて，
9
因子多項式を探す．
23 5 分割で H を扱う場合には，グループ 5 の e = 2.999 · · · と，グループ 9 の e = 5 で割り切れて，特性多項式は 6 次式
9
5
9
(λ6 − 24λ5 + 229λ4 − 1104λ3 + 2812λ2 − 3552λ + 1716) になる．
24 usev は組合わせの番号からグループ番号へのポインタで， usev[idx[j]] が j 番目のグループを指す．
25 初期化は，idx[1] から idx[d] に 1 から d を格納し，最後の idx[d] からは 1 を引く（後述）．5 分割で H を扱う場合には，6 つ
9
の候補があり nn = 6 C2 = 15 である．
26 組合せは d=2 なら 12,13,· · · , ncandi,23,24,· · · のように idx[1] と idx[2] に入る．usev は組合わせの番号からグループ番号へのポ
インタで， usev[idx[j]] が j 番目のグループを指す．
27 全係数を作成してから整数性のチェックを行うと，プログラムは単純になるが計算時間がかかる．そこで総和で計算できる 1 次項だけ
を先に判定し，これを通過した場合だけ，他の係数を求める．この先行判定は有理算術演算では大きなチューニング効果を与える．なお第 d
次係数を総積で求めるより，1 次の係数を求めるほうが，桁数が少ないので速い場合が多い．
28 5 分割で H を扱う場合には，6 つの候補から idx[1]=1, idx[2]=4 のとき，usev[1]=2，usev[4]=8 で，e = 1.2679 · · · と
9
2
e8 = 4.7320 · · · が選択されると，係数 −6 と 6 が作られ，(λ6 − 24λ5 + 229λ4 − 1104λ3 + 2812λ2 − 3552λ + 1716) ÷ (λ2 − 6λ + 6) =
λ4 − 18 ∗ λ3 + 115 ∗ λ2 − 306 ∗ λ1 + 286 と割り切れる．
29 5 分割で H を扱う場合には，2 次式 (λ2 − 6λ + 6) で割り切れると，6 つの候補が 4 つに減り，配列 evcan から 2 つを取り除いて，
9
グループ 4,6,10,12 を evcan にリセットする．
16
selectev 関数
void selectev(const vector<double>& a, int n, int& ngrp, vector<double>& eigvl,
vector<double>& eigvh, vector<int>& mult) {
int i,ii,j=0;
double x,am;
am=fabs((a[n-1]-a[0])/((double)n)); /* abs(mean) */
std::cout << "selectev Amean=" << am << std::endl;
eigvl[j+1]=a[0]; eigvh[j+1]=0; x=a[0]/am; ii=0; mult[j+1]=0;
for(i=1; i < n; i++){
if(fabs(a[i]/am-x) > eps){ j++; x = a[i]/am; eigvl[j+1]=a[i]; ii=i; mult[j+1]=0; }
if(ii != i){ eigvh[j+1]=a[i]; mult[j+1]++; }
}
ngrp = j;
}
複数ある固有値の真ん中のものを代表として選び，eps を相対誤差判定で使用し，4 分割の表 2 では 16 個の近
似固有値が 9，5 分割の表 1 では 25 個の近似固有値が 13 のグループに分類される．この関数は後述の章でも使
用する．
2 次以上の因子を探す部分は，2 重のループ構造（次数 2 から特性多項式の次数の while ループと，組合わせ
のループ）の中に 3 つの if ブロックがあり，goto 文で制御が戻るところもあるのでプログラムは複雑である．
この理由は，候補数 ncandi が減ってゆき，これがループ構成で組合せループの終端を決めるからである．候補
が減っても，常に候補の組合せをはじめから調べれば，プログラムは単純になるが，すでにチェック済みの組合
せを再びチェックすることになる30 ．これを避けるために組合せループの前処理をして goto 文で制御を Top に
戻している．
2.2
根と係数の関係公式と組合せのプログラミング
誤差を含む近似固有値から，誤差を除外した整数の多項式係数を抽出するところが鍵である．本節では，根と
係数の関係公式と組合せのアルゴリズムとそのプログラミングを説明する．
2.2.1
根と係数の関係公式
公式（Vieta’s formula）を，次数 n = 5 の多項式 p0 x5 + p1 x4 + p2 x3 + p3 x2 + p4 x + p5 を例に説明する．根
を α1 から α5 とすると次のようになる（最高次の係数は 1 か −1 である．最高次の係数 p0 = 1 の多項式をモ
ニック多項式という）．
p0 = −1
p1 = α1 + α2 + α3 + α4 + α5
p2 = −(α1 α2 + α1 α3 + α1 α4 + α1 α5 + α2 α3 + α2 α4 + α2 α5 + α3 α4 + α3 α5 + α4 α5 )
p3 = α1 α2 α3 + α1 α2 α4 + α1 α2 α5 + α1 α3 α4 + α1 α3 α5 + α1 α4 α5 + α2 α3 α4 + α2 α3 α5 + α2 α4 α5 + α3 α4 α5
p4 = −(α1 α2 α3 α4 + α1 α2 α3 α5 + α1 α2 α4 α5 + α1 α3 α4 α5 + α2 α3 α4 α5 )
lack α5
lack α4
lack α3
lack α2
lack α1
p5 = α1 α2 α3 α4 α5
符号は，次数が奇数なら負，偶数なら正である．
30 候補が 5 つあって，組合わせ (2,3) で因子が見つかると，候補は 3 つに減る．このとき 1 番目の候補に対するチェックは終わっている
ので，新しい組合わせ (2,3) から因子探しを始める．
17
プログラミングを行うにあたり，添字の並びのルールを調べる．各係数を形成する項 αk の添字の並び方に着
n
αk である．p1 と p4 は互いに補（compliment）の関係がある．つまり，
目する．定数項は総積 pn = p5 =
p1 は k = 1, 2, · · · , n の総和
n
k=1
αk で，総和の第 1 項は α1 である．p4 は 4 つの α の積の総和だが，その第 1
k=1
項は 1 から n までの α の項から α1 が抜ける．p1 の総和の第 2 項は α2 で，p4 の第 2 項は 1 から n までの α
の項から α2 が抜ける．· · · ．どちらも項数は n である．
p2 と p3 にも互いに補の関係があり，p3 の項を後ろから逆に読むと，抜けた添字が p2 の項を前から読む添字
10
10
αi1 αi2 であり，p3 =
αi1 αi2 αi3 である．総和記号の上と下につく式を一般的に書く
に一致する．p2 =
i1,i2
i1,i2,i3
と次のようになる．
n Cr
αi1 αi2 · · · αir
(19)
i1,i2,··· ,ir
2.2.2
組合せのプログラミング（十進 BASIC による準備）
もっとも基本的な「1 から 100 まで数える」ループプログラムを，C++ と 十進 BASIC で示す31 ．BASIC は
1960 年代に生まれた．当時の文字は 6 ビットコードだったので，プログラミングに小文字は使えなかった．7 ビッ
トコードの ASCII コードが一般化した後に生まれた C や C++ では大文字と小文字は別の文字として区別する
が，Fortran や BASIC では区別しない．十進 BASIC ではカーソルを移動すると，制御文字は大文字に自動的
に変更され，ユーザの書いた変数名はそのままになる（字下げも自動的にブランクを挿入してくれる）．本稿で
は変数名は主に小文字を使用する．
C などのコンパイラ型プログラミング言語では変数に型を指定するが，BASIC はインタープリータ型のプログ
ラミング言語なので，数値は電卓のように一通りに扱われる（型を指定する必要はない）．C や C++ では for
と初期化式 int i=1; を区切る「セパレータ」に括弧を使用するが，BASIC はセパレータにブランクを用いる
ので，FOR の次にブランクが必要である．1 から「1 つづつ増やして」を指定する部分は「再初期化式」という
が，BASIC は STEP 1 と指定する．増分が 1 の場合は省略可能である．C や C++ ではステートメントの終わ
りはセミコロンで示すが，60 年代のパンチカードにプログラムのステートメントを記述した時代の言語ではセミ
コロンは使わない（改行コードでステートメントが終わる）．十進 BASIC は，セミコロンで区切って 1 行に複
数のステートメントを書くことはできない．
int n=100;
for(int i=1; i <= n; i++){
.....; .....;
}
正規化ループ
n=100
FOR i=1 TO n STEP 1
.....
.....
NEXI i
C や C++ では変数 i や n の変数型を int にすると 32 ビット整数型の最大値，2,147,483,647= 2G − 1 まで
数えられ，それより大きな数まで数えたいときは変数型を，例えば long long 型などに変更する．十進 BASIC
では p/q のアイコンをクリックして「有理数モード」にすれば，1000 桁の数まで数えられる．
順列，組合せ，重複順列，重複組合せの 4 つは，数え上げプログラミングのループ構成の基本である．5 個の
ものから 2 個を選んで作る組合せは n = 5 として 2 重のループ構成で，5 個のものから 3 個を選んで作る組合
せは 3 重のループ構成で作ることができる．
31 最も基本的なループで「正規化ループ」
（normalized
loop）と呼ぶ．
18
Combination.BAS
LET n=5
FOR i1=1 TO n
FOR i2=i1+1 TO n
PRINT i1;i2
NEXT i2
NEXT i1
LET n=5
FOR i1=1 TO n
FOR i2=i1+1 TO n
FOR i3=i2+1 TO n
PRINT i1;i2;i3
NEXT i3
NEXT i2
NEXT i1
左のプログラムは次の出力を与える．
i1
1
1
1
1
2
2
2
3
3
4
i2
2
3
4
5
3
4
5
4
5
5
i1,i2 を 2 桁の数値として読むと，昇順になっている．右のプログラムは次の出力を与える．
i1
i2
1
2
1
2
1
2
1
3
1
3
1
4
2
3
2
3
2
4
3
4
i3
3
4
5
4
5
5
4
5
5
5
根と係数の関係公式では，与えられた n に対して，r = 1 から r = n までを変数にして動かす必要がある．
プログラムがコード生成して，生成されたコードに制御を渡してそのコードを実行する方式をコード生成（code
gen）方式という32 ．この方式を使用すれば上記のループ構成でもよいが，コード生成方式を採用するまでもな
く，ここでは Combination.BAS の r 重にネストしたループ構成を 1 重にループ変換すればよい33 ．そこでまず，
呼び出されるたびに表の i1,i2 や i1,i2,i3 の数列を引数の配列に返すサブルーチン nxtcmb （「次の組合せ」
の意）を用意する．n = 5 で r = 2 なら，ith(1)=1，ith(2)=1 に初期化して呼ぶと配列 ith には 1,2 が返る．
次は 1,3 が返る．· · · ．
CombSet.BAS の nxtcmb サブルーチン
EXTERNAL SUB nxtcmb (n,r,ith())
IF ith(r) = n THEN
FOR i=r-1 TO 1 STEP -1
IF ith(i) <= n-(r-i)-1 THEN
LET ith(i)=ith(i)+1
FOR j=i+1 TO r
LET ith(j)=ith(i)+j-i
NEXT j
EXIT FOR
END IF
NEXT i
ELSE
LET ith(r)=ith(r)+1
END IF
END SUB
!
!
!
!
!
!
最後の桁が n だったら
前の第 i 桁を探す．
第 i 桁が大きくできれば
第 i 桁をインクリメントして
それ以降の桁を
第 i 桁に従ってセットする．
! 最後の桁が n でなかったら
! 最後の桁をインクリメントする．
十進 BASIC には組合せの総数を与える comb 関数が用意されている．これを使って n と r を変数にして，組
合せの添字を順番に得ることができる．プログラム CombSet.BAS は，n を与えると，n C1 ，n C2 ，n C3 ，n C4 ，
n C5
それぞれについて組合せを表示するが，各組合わせは次の SetInd サブルーチンで作成して表示している．
32 現代のプロセッサは，命令とデータを分けて記憶域からフェッチする．コード生成はデータ領域に生成し，これを実行するときは命令と
してフェッチすることでコード生成方式は実行される．
33 「ループ変換技術（loop transformation technique）
」は HPC プログラミングの根幹をなすプログラミング技術である．
19
CombSet.BAS の SetInd サブルーチン
EXTERNAL SUB SetInd (n,r)
DIM ith(r)
FOR i=1 TO r
LET ith(i)=i
NEXT i
LET ith(r)=ith(r)-1
PRINT "comb(n,r)=";comb(n,r)
FOR i=1 TO comb(n,r)
CALL nxtcmb(n,r,ith)
FOR j=1 TO r
PRINT ith(j);
NEXT j
PRINT
NEXT i
END SUB
!
!
!
!
配列の初期化
配列の初期化
配列の初期化
配列の初期化
! _n C_r 回反復
! 次の組合せを ith[1] から [r] に得る
CombSet.BAS によって 1 から n までの組合せループの 1 重化の方法を確認できる．
2.2.3
根と係数の公式と組合せのプログラミング（C++）
十進 BASIC の comb 関数を作成する．戻り値は 64 ビット整数（long long 型）とした．この変数型で扱え
る最大数は 263 − 1 で，n > 61 では桁溢れするので停止するようにした．
n Cr
=
n(n − 1) · · · (n − r + 1)
r(r − 1) · · · 2
(20)
を計算するので，r が 2 なら分子のどちらかは偶数である．r が 3 なら分子の内の１つは 3 の倍数である．した
がって分子の 1 項と分母の 1 項を交互に計算すれば，中間桁溢れは起こらない34 ．
rational.cpp の comb 関数（long long 型）
long long comb(const int n, const int r){
long long t=(long long)n;
for(int i=1; i<r-1;i++){
t=t*((long long)(n-i));
t=t/((long long)(i+1));
}
return t;
}
ただし，64 ビット整数の除算は時間がかかるので，式 (20) の分子と分母を longint 型の変数で独立して乗算のみ
で計算してから，除算を 1 回としたほうが速い．ここでは longint 型の紹介を兼ねて，型変換の方法と longint
クラスの除算関数 ldivide の使用例を示す．
34 最終結果が桁溢れしなくても，それを算出する過程で桁溢れすることを「中間桁溢れ」という．
20
rational.cpp の comb 関数（longint 型）
int i;
longint t("1"), s("1"),x,y; // how to set a number into longint variable
luint32_t ni;
long long result;
for(i=0; i<r;i++){
ni=n-i;
x=ni;
// cast int to longint
t=t*x;
}
for(i=1; i<=r;i++){
ni=i;
if(ni != 0){ x=ni; s=s*x; }
}
x=longint::ldivide(t,s,&y);
result=x;
return result;
}
nxtcmb 関数を示す．
rational.cpp の nxtcmb 関数
void nxtcmb(const int n, const int r, int ith[]){
int i,j;
if(ith[r]==n){
for(i=r-1; i>0;i--){
// search position to add
if(ith[i]<=n-(r-i)-1){
ith[i]++;
for(j=i+1; j<=r;j++) ith[j]=ith[i]+j-i;
break;
}
}
}else{ ith[r]++; }
}
倍精度配列 e に近似固有値が n 個格納されており，これから r 個の組合せを n Cr 選び出し，それらの積の和 pr
を倍精度浮動小数点数で返す関数 vietaterm を示す．
rational.cpp の vietaterm 関数
double vietaterm(const double e[], int n, int r){ // <== rational_vector
int ith[r+1];
for(int i=1; i<=r;i++) ith[i]=i;
ith[r]=ith[r]-1;
double t=0;
// <== rational
long long nn=comb(n,r);
for(long long i=1; i<=nn; i++){
nxtcmb(n,r,ith);
double s=1;
// <== rational
for(int j=1; j<=r; j++) s=s*e[ith[j]];
t=t+s;
}
if(r%2==1) t=-t;
! 次数が奇数は符号反転
return t;
}
同名の vietaterm で，コメントでマークした 2 行の変数型を rational とした rational 型で計算する関数も多
重定義した（後述）．
“for(初期化式; 継続条件式; 再初期化式)” 文は，セミコロンで区切られた 3 つの式は独立に評価される．した
21
がって，初期化式を long long i を用いているので，反復回数は 32 ビットを超えて 263 − 1 まで実行できる．
さらに多くの反復を実行する場合は，変数を rational とする．
2.2.4
整数性判定と有理数での整数の抽出
倍精度浮動小数点数 t が，絶対誤差で epsint の範囲内に整数が存在するか否かを判定する．t が絶対値の小
さな数なら問題ないが，次章で扱う大きな数ではこの方法に限界が出る．
rational クラスに次の倍精度数の場合の整数性チェックのための askinty 関数を用意した．第 2 引数の epsint
以内にある整数を有理数として抽出して第 3 引数に返す．
rational.cpp の askinty 関数（倍精度数）
int askinty(const double t, const double epsint, rational& vrat) {
double u;
int rtnval=0;
if(fabs(floor(t)-t) < epsint){
u=floor(t); vrat=Rdset(&u); rtnval=1;
}
if(fabs(ceil(t)-t) < epsint){
u=ceil(t); vrat=Rdset(&u); rtnval=1;
}
return rtnval;
}
同名の関数で引数 t を rational 型で計算する関数も多重定義した（後述）．
2.3
多項式の格納と除算
多項式の表記は， p(x) = a3 x3 + a2 x2 + a1 x + a0 のように冪と係数の添字を一致させる（降順の）書き方と，
逆に p(x) = a0 x3 + a1 x2 + a2 x + a3 のように最高次の係数添字をゼロに，定数項の添字を次数に一致させる（昇
順の）書き方とがある．高校の数学では昇順を用いるが，多項式の乗算や除算を行う場合には降順のほうが自然
にプログラミングできる．7 ページのフロベニウス変換式 (7) でも冪と行番号を一致させる降順が用いられる．本
資料では，17 ページの根と係数の関係公式や，12 ページの方程式の解の変換公式 (14) では係数の添字とべきを
一致させるために昇順を使用した．一般に両方の使用法が用いられる．
n 次多項式の n + 1 個の係数を配列 a[0] から a[n] に与えてこの多項式係数を putpolynomial 関数で格納
し，また getpolynomial 関数で配列に取り出す．これらの操作では，多項式の項が昇順か降順かを決める必要
はなく，使用者が昇順か降順かを理解していればよい．本例題プログラムでは降順で使用している．
2.3.1
多項式の格納形式
項数可変の複数の多項式を格納するのに，2 つの配列を用いて，一方をポインタ配列として，他方を多項式の
項を詰めて格納する配列として使用する．各多項式の項数は，連続するポインタの差として得られる35 ．係数配
列の先頭を指すポインタの差 ppolynomial[k]-ppolynomial[k-1]-1 が多項式の次数になる．したがって，次
の putpolynomial 関数で格納でき
35 このデータ形式は，行列計算では疎行列をスカイライン形式で扱う場合のスカイライン添字配列と行列要素の扱い方と同じである．
22
rational.cpp の putpolynomial 関数
int putpolynomial(const int nsize, const rational_vector& cof, int npolynomial,
vector<int>& ppolynomial, rational_vector& cpolynomial){
int i,i1,i2;
i1=ppolynomial[npolynomial]; i2=i1+nsize+1;
for(i=i1; i<i2; i++) cpolynomial[i]=cof[i-i1];
npolynomial++;
ppolynomial[npolynomial]=i2;
return npolynomial;
}
逆の getpolynomial 関数で取り出せる．引数 npolynomial は値渡しである．何番目の多項式かを与えて取り出
し，戻り値は取り出した関数の次数である．
rational.cpp の putpolynomial 関数
int getpolynomial(const int k, rational_vector& cof, const int npolynomial,
const vector<int>& ppolynomial, const rational_vector& cpolynomial){
int i,i1,i2;
if(k>npolynomial){
cout << "getpolynomial error npolynomial=" << npolynomial <<
" smaller than the first argument=" << k << endl;
exit(1);
}
i1=ppolynomial[k-1]; i2=ppolynomial[k];
for(i=i1; i<i2; i++) cof[i-i1]=cpolynomial[i];
return i2-i1-1;
}
使用例を示す．
put/getpolynomial 関数の使用
/********* get, put polynomials
*********************************/
int npolynomial=0;
/* number of polynomials found */
vector<int> ppolynomial(n+2);
/* pointer to polynomials
*/
rational_vector cpolynomial(n*2); /* coefficients of polynomials */
/******************************************************************/
ppolynomial[0]=0;
// 初期化
int ksz=getpolynomial(iod,vrat,npolynomial,ppolynomial,cpolynomial);
polytexformr(vrat,ksz);
取り出した関数を polytexformr で表示するが，表示する関数は昇順の場合は polytexform を使用する．
polytexformr 関数も polytexform 関数も rational クラスに含まれる．また，有理数の分母が 1 の場合は “/1”
を取り除く polytexformint 関数も用意した．特性多項式 (13) の polytexformint 関数による表示は次のよう
に得られる．
polytexformint 関数の出力
-1*x^3 +12*x^2 -44*x^1 +48
これをコピー&ペーストで数式処理システムに入力できる．
2.3.2
多項式の除算
多項式の除算は次の rational クラスの関数で行う（高校の数学で「整式の除法」）．
23
rational.cpp の polydivchk
int polydivchk(rational_vector& u, const int m, const rational_vector& v,
const int n, rational_vector& q, rational_vector& r){
for (int k=m-n; k>=0; k--){
q[k] = u[n+k] / v[n];
for (int j=n+k-1; j>=k; j--) u[j] = u[j] - q[k] * v[j-k];
}
int rtnval=1;
! 剰余のチェック
for (int j=0; j<=n-1; j++){
r[j] = u[j];
if(r[j] != rational::ZERO) rtnval=0;
}
return rtnval;
}
u は変更されるので注意のこと．割り切れた場合は戻り値は 1 で，剰余が非零の場合は 0 を返す．
2.4
包含チェックと 1 次の因子探し
包含は，区間両端での関数値を求めて，符号を調べる．関数値はホーナー法で計算するが，導関数の値もホー
ナー法で調べられる．十進 BASIC と C++で説明する．
2.4.1
包含チェック
selectev 関数でグループ化された近似固有値に対し，各グループに属する近似固有値の下限値が引数 e に，上
限値が引数 h に，相対誤差 ε = 10−9 を考慮して，特性多項式 p(λ) に対して包含を調べる．近似固有値は倍精
度，多項式係数は有理数で得られているので，近似固有値を rational 型の変数 a と b に有理数変換して格納し，
有理算術演算によるホーナー法 Horner で関数値 p(a) と p(b) を正確に計算し，符号変化を調べることで包含を
判定する．
sehfv1.cpp の chkinc 関数
int chkinc(const int n, const rational_vector& p, const double e,
const double h, const int mult){
double ta,tb;
int rtnval=0;
rational a,b,fa,fb;
if(fabs(e) > eps){
ta=e-fabs(e)*eps; tb=h+fabs(h)*eps;
if(mult==0) tb=e+fabs(e)*eps;
}else{
ta=e-eps; tb=h+eps;
// consider eigenvalues close to zero
if(mult==0) tb=e+eps;
}
a = Rdset(&ta); fa=Horner(n,p,a);
b = Rdset(&tb); fb=Horner(n,p,b);
if(fa*fb<rational::ZERO) rtnval=1;
return rtnval;
}
包含が成立していれば戻り値は 1 で，これを配列 included に格納する．sehfv1.cpp では，この時点での多項
式次数 nsize と，包含の成立数を比較して，一致していた場合のみ計算を続行する．
24
2.4.2
ホーナー法
n 次の多項式
pn (x) = a0 xn + a1 xn−1 + a2 xn−2 + · · · + an−2 x2 + an−1 x + an
(21)
を = 0 と置いた方程式を代数方程式（algebraic equation）と呼ぶ．係数 ak は数学の教科書では複素数か実数
として議論されるが，本資料では計算機で扱える数を扱うので有理数とする．多項式 (21) の x = ξ における値
pn (ξ) は，ホーナー法によって計算する．
pn (ξ) = (· · · ((a0 ξ + a1 )ξ + a2 )ξ + · · · + an−1 )ξ + an
(22)
多項式の各項を，式 (21) の左から順に a0 × x × x × x · · · ，次に a1 × x × x × x · · · のように計算すると，n(n + 1)/2
回の乗算が必要になる．ホーナー法ではこれを n 回で行える．十進 BASIC で示す．
EXTERNAL FUNCTION Horner(n,a(),x)
OPTION BASE 0
LET r=a(0)
FOR i=1 TO n
LET r=r*x+a(i)
NEXT i
LET Horner=r
END FUNCTION
pn (x) を x − ξ で割ったとき，その商と剰余を Q(x) と R で表せば
pn (x) = (x − ξ)Q(x) + R
(23)
pn (ξ) = R
(24)
となる．これに x = ξ を代入すれば
が得られる．つまり，pn (ξ) の値は pn (x) を x − ξ で割った剰余としても求められる．これは剰余定理（remainder
theorem）として知られている．例えば p3 (x) = x3 − 15x − 4 を x − 3 で割った場合，組立除法（synthetic
division）によって商と剰余を求めてみる．
1
0
3
–15
9
–4
–18
1
3
–6
–22
+)
3
上段には多項式の係数を並べ，右端には x − 3 の 3 を記す．次に最高次の係数 1 をそのまま下段に下ろし，それ
に上段右端の 3 を掛けた値 3 を（矢印に従って）中段右隣に記す．上段の値と中段の値の和を下段に記す．こ
の操作を繰り返すと，商多項式の係数と剰余 p3 (3) = (x2 + 3x − 6)(x − 3) − 22 が得られる．
組立除法の一般の場合を示す．
a0
a1
a2
···
an−2
an−1
an
a0
a0 ξ
b1
b1 ξ
b2
···
···
bn−3 ξ
bn−2
bn−2 ξ
bn−1
bn−1 ξ
R
となる．ここで，最下段の値 bi ,
b1
ξ
i = 0, 1, . . . , n − 1 および R は，
= a0 ξ + a1
b2 = b1 ξ + a2 = (a0 ξ + a1 )ξ + a2
..
.
bn−1
R
= bn−2 ξ + an−1 = (· · · ((a0 ξ + a1 )ξ + a2 )ξ + · · · + an−2 )ξ + an−1
=
(· · · ((a0 ξ + a1 )ξ + a2 )ξ + · · · + an−1 )ξ + an
25
となる．すなわち組立除法による剰余 R の計算は，ホーナー法による pn (ξ) の計算順序に一致している．pn (ξ)
は pn (x) を x − ξ で割ったとき，その商を Q(x)，剰余を R で表わせば，
pn (x) = (x − ξ)Q(x) + R,
pn (ξ) = R
(25)
である．両辺を x で微分すると
pn (x) = (x − ξ)Q (x) + Q(x)
なので pn (ξ) = Q(ξ) を得る．つまり，多項式 pn (x) の導関数 pn (x) の x = ξ での値 pn (ξ) は， Q(x) を x − ξ
で割ったとき，その剰余として求められる．
Q(x) = (x − ξ)T (x) + S,
pn (ξ) = Q(ξ) = S
(26)
1 階導関数 pn (x) の x = ξ での値（接線の傾き）を求めるサブルーチンを作っておく．
EXTERNAL SUB Horner2(n,a(),x,r,s)
OPTION BASE 0
LET r=a(0)
LET s=a(0)
FOR i=1 TO n-1
! i=1
LET r=r*x+a(i) ! b1=a(0)*x+a(1)
LET s=s*x+r
! s =a(0)*x+b1
NEXT i
LET r=r*x+a(n)
END SUB
i=2
i=3
b2=b1*x+a(2)
b3=b2*x+a(3)
s=(a(0)*x+b1)*x+b2 s=((a(0)*x+b1)*x+b2)*x+b3
2 階導関数まで考慮したサブルーチン Horner3 を用意しておこう．
EXTERNAL SUB Horner3(n,a(),x,r,s,t)
OPTION BASE 0
LET r=a(0)
LET s=a(0)
LET t=a(0)
FOR i=1 TO n-2
! i=1
LET r=r*x+a(i) ! b11=a(0)*x+a(1)
LET s=s*x+r
! b21=a(0)*x+b11
LET t=t*x+s
! t =a(0)*x+b21
NEXT i
LET r=r*x+a(n-1)
LET s=s*x+r
LET r=r*x+a(n)
LET t=t*2
END SUB
i=2
i=3
b12=b11*x+a(2)
b13=b12*x+a(3)
b22=b21*x+b12
b23=b22*x+b13
t=(a(0)*x+b21)*x+b22 t=t*x+b23
Horner2 サブルーチンを用いる例 NewtonSDFig.BAS を例題として加えた．このプログラムは，初期値（例
えば 3）を与えると，x3 − 15x − 4 = 0 の近似解をニュートン法反復で求め，収束までの過程を作画する．
Horner と Horner3 の C++ 版を，変数型を rational で，rational クラスに含めた．
rational Horner(const int n, const rational_vector& coef, const rational& x) {
rational r = coef[0];
for (int i=1; i <= n; ++i) r = r*x + coef[i];
return r;
}
void Horner3(const int n, const rational_vector& coef, const rational& x,
rational& r, rational& s, rational& t) {
if(n >= 2){
r = coef[0]; s = coef[0]; t = coef[0];
for (int i=1; i < n-1; ++i) { r = r*x + coef[i]; s = s*x + r; t = t*x + s; }
26
r = r*x + coef[n-1]; s = s*x + r; r = r*x + coef[n]; t = t*rational::TWO;
}else if(n ==1){
r = coef[0]*x + coef[1]; s = coef[1]; t = rational::ZERO;
}else if(n ==0){
r = coef[0]; s = rational::ZERO; t = rational::ZERO;
}
}
2.4.3
1 次の因子探し
包含が成立しているグループに対して36 ，近似固有値を askinty 関数で評価して，整数と見なせればこの整
数を rational 型の変数に抽出し，1 次式 v(λ) の定数項として格納し， v(λ) で p(λ) を割ってみる（polydivchk
関数）．
1 次の因子探し
for(i=1; i<=ngrp1;i++){
包含が成立していれば
if(included[i]){
if(askinty(eigvl[i],epsint,vrat[0])){
整数性が満たされていれば
for (j=0; j<=nsize; j++) urat[nsize-j]=p[j]; 特性多項式を urat に
vrat[1]=-rational::ONE;
割るほうの vrat を作り
if(polydivchk(urat,nsize,vrat,1,qrat,rrat)){ 割り切れれば
included[i]=0; nsize --;
for (i=0; i<=nsize; i++) p[nsize-i]=qrat[i]; 特性多項式を商多項式で置換え
npolynomial=putpolynomial(1,vrat,npolynomial,ppolynomial,cpolynomial);
} } } }
割り切れれば戻り値が 1 で，このとき商多項式で特性多項式を置き換え，次数 nsize をデクリメントし，1 次
式を因子多項式として putpolynomial 関数で保存する．
2.5
2 次以上の因子探し
ループ構成は，探索する多項式の次数 d を 2 次から nsize 次まで探索するが，因子多項式が見つかるたびに，
特性多項式が因子で割られて，商多項式によって置き換えられるので，特性多項式の次数 nsize が小さくなる．
このようなループ端が変化する場合のループ構成 は，注意してプログラミングする必要がある．Fortran の DO
ループや BASIC の FOR ループでは，ループの反復終端が反復中に変化しても，反復は反復開始前に設定した回
数を実行する．これに対し，C 言語の for ループでは，ループの反復終端が反復中に変化すると，変化に応じて
反復は終了する．このようなプログラミング言語の差異を考慮して，sehfv1.cpp では while(d<=nsize) ループ
を用いた37 ．
探索は，まず setecand 関数で包含が成立している近似固有値だけを選び，詰めて配列 evcan に格納し，ま
た，次数 d 以下の番号でグループ番号を間接参照するための配列 usev を用意する．戻り値は選択された固有値
の数である．
36 for(i=1;
i<= ngrp1; i++) ループを if(included[i]) で条件選択する．
では do while(...) 文 と end do 文で，十進 BASIC では DO WHILE ... 文と LOOP 文で行う．BASIC ではブランクがセ
パレータなので，条件を括弧で囲わない．
37 Fortran
27
setecand
int setecand(const vector<double>& eigvl, const int ncan,
const int included[], double evcan[], int usev[]){
int j=0;
for(int i=1; i<=ncan;i++){
if(included[i]){ j++; usev[j]=i; evcan[j]=eigvl[i]; };
}
return j;
}
戻り値を ncandi に保存すると，探索は，ncandi 個の候補から，次数 d 個を選ぶ組合せ
ncandi Cd に対して行う．
nn=comb(ncandi,d) によって組合せ数が決められ，配列 idx を 1 から d までの最初に返すために idx[d]=d-1
に初期化する．
2 重のループ構造をもち，内側のループの内部に 3 つの if ブロックがある（1 つは後述する高速化の先行判
定のため）．探索は内側の for(ii=1; ii<=nn;ii++) のループで行うが，このループ反復中に多項式因子が見
つかると，候補は除外されて，候補数 ncandi が減る38 ．この時点で配列 idx の状態によって，同じ次数 d の候
補探しがまだ残っていて継続するか，あるいは終了するかで探索方法が変わる．継続の場合は，idx 配列の初期
化を行って，goto Top で飛ぶが，終了の場合はループを抜けて新たな次数 d の探索に移る．飛ぶか終了するか
の判定は，リセットされた idx[d] が ncandi より大きいかどうかで判定する39 ．
38 変数
ii は long long 型
ページの 2.2.1 項の脚注の例で述べる．5 つの候補に対し，2,3 で割り切れると，ncandi は 3 になり，noff = idx[1] − 1 = 0 で
idx は 1,2 にリセットされる．d = 2 なので ncandi よりも小さく，nn をリセットして goto 文で飛ぶ．
39 17
28
2 次以上の因子探し
d=2;
while(d<=nsize){
// search for polynomial of degree 2 to nsize
ncandi=setecand(eigvl,ngrp1,included,evcan,usev); // set number of candidate eigenvalues
for(i=1; i<=d;i++) idx[i]=i;
// initialize for nxtcmb
idx[d]=idx[d]-1;
// initialize for nxtcmb
nn=comb(ncandi,d);
// since nn is type "long long", ngrp1-nuse should be "< 62"
Top:
for(ii=1; ii<=nn;ii++){ // loop over nn combinations of candidates
nxtcmb(ncandi,d,idx);
// select one of _(ncan-nused) C_d candidates
u=0;
for(j=1; j<=d; j++){ figvl[j]=evcan[idx[j]]; u=u+figvl[j];}
if(askinty(u,epsint,vrat[0])){ // pass the first cofficient check, then ...
int ok=1;
for(int r=1; r<=d; r++){
cof[r]=vietaterm(figvl,d,r);
if(! askinty(cof[r],epsint,vrat[d-r])){ ok=0; break; }
}
if(ok){
// OK, try trial division
vrat[d]=rational::ONE;
for (j=0; j<=nsize; j++) urat[nsize-j]=p[j];
if(polydivchk(urat,nsize,vrat,d,qrat,rrat){; // trial division
for (j=0; j<=nsize-d; j++) p[nsize-d-j]=qrat[j]; // set quatient for next p(\lambda)
nsize = nsize-d;
// and decrease nsize
npolynomial=putpolynomial(d,vrat,npolynomial,ppolynomial,cpolynomial);
for (j=1; j<=d; j++) included[usev[idx[j]]]=0;
ncandi=setecand(eigvl,ngrp1,included,evcan,usev); // reset candidate eigenvalues
noff=idx[1]-1;
// number of OFFed indecies
for(j=1; j<=d;j++) idx[j]=j+noff; // initialize for nxtcmb
if(idx[d]>ncandi) break;
// there are no degree d left, then break for loop
idx[d]=idx[d]-1;
// initialize for nxtcmb
nn=comb(ncandi-noff,d);
// reset nn (number of combinations)
goto Top;
/////////////// to loop Top //////////////////
} } } }
d++;
}
5 分割で H9 を扱う場合には，グループ 2,4,6,8,10,12 の 6 つの候補から idx[1]=1, idx[2]=4 のとき，usev[1]=2，
√
√
usev[4]=8 で，e2 = 1.2679 · · · 3 − 3 と e8 = 4.7320 · · · 3 + 3 が選択されると，係数 −6 と 6 が作られ
て特性多項式が 2 次式 (λ2 − 6λ + 6) で割り切れる．この組合せは 6 C2 = 15 のうちの 3 番目である．2 つのグ
ループが除外され，setecand で 4 グループがリセットされるが，割り切れたときの idx[1]=1 なので，noff = 0
より，リセットされる idx は 1, 1 である（4 候補を始めから探索）．
次にグループ 4,6,10,12 の 4 つの候補から idx[1]=1, idx[2]=3 のとき，usev[1]=4，usev[3]=10 で，e4 =
√
√
2.2679 · · · 4− 3 と e10 = 5.7320 · · · 4+ 3 が選択されると，係数 −8 と 13 が作られて 2 次式 (λ2 −8λ+13)
で割り切れる．この組合せは 4 C2 = 15 のうちの 2 番目である．2 つのグループが除外され，setecand で 2 グ
ループがリセットされるが，割り切れたときの idx[1]=1 なので，noff = 0 より，リセットされる idx は 1 で
ある（2 候補を始めから探索）．
最後にグループ 6,12 の 2 つの候補から idx[1]=1, idx[2]=2 のとき，usev[1]=6，usev[3]=12 で，e6 =
√
√
3.2679 · · · 5− 3 と e12 = 6.7320 · · · 5+ 3 が選択されると，係数 −10 と 22 が作られて 2 次式 (λ2 −10λ+22)
で割り切れる．この組合せは 4 C2 = 15 のうちの 1 番目である．ncandi がなくなるので探索は終わる．この場
合，直観的に H13 から
13 Cr
による多くの反復を想像しがちだが，このように次数の低い因子が見つかると，組
合せループの反復は 3 + 2 + 1 で 6 回で済む．
H13 は H9 の特性多項式で割り切れるため，4 次式となり，これが 1 次式で割り切れるためすぐに計算が終
わる．
29
熱伝導行列は構造が単純で高い多重度をもち，低い次数の因子多項式に分解されるため，組合せループを探索
する回数も少なく，このアルゴリズムでも分解は短時間で終了する．
30
3
有理行列の解析プログラム sehfv2.cpp
前章のテストプログラムで扱った行列の固有値は 0.1 と 10 の間に収まっていたので，浮動小数点計算で得ら
れた固有値を，そのまま近似値として使用し，倍精度浮動小数点演算で整数性も判定できた．倍精度浮動小数点
計算で得られた固有値の精度を 10 進で 15 桁とすると，近似固有値の小数点以下に 14 桁の有効数字があり，ま
た因子多項式の係数も少ない桁数で収まっていたからである．しかし目標の倍精度浮動小数点演算で求められた
行列の解析には，近似固有値の精度の問題を解決しなくてはならない．倍精度浮動小数点数は数学的には有理数
である．行列要素を有理数変換し，さらに整数係数の特性多項式を得るには，有理行列の全要素の分母の LCM
を掛けるので，特性多項式の係数は大きくなる．また固有値も LCM 倍されるため，倍精度計算で得られた近似
値を使って整数性から因子を探すには，近似固有値の精度が不足する．sehfv1 は，データの好条件に支えられて
偶然稼働していた，と言っても過言ではない．
数学で扱う「整数」は数の大小に関係なく整数としての性質を備えているが，計算機で行う「計算」で扱える
整数は 231 − 1 以下とか 263 − 1 以下などと，大きさに制限がある．既存の多くの汎用プログラミング言語で扱
える整数はごく小さな数だけである40 ．また，浮動小数点演算は，限られた精度に丸めて計算するため，正確な
計算はできないことが多い（ごく小さな数の計算は正確にできることがある）．浮動小数点演算では，精度の解
決策として，倍精度を 4 倍精度に拡張することも考えられるが，このような方法は中途半端である．4 倍精度で
も不足すると 8 倍精度に，さらに 16 倍精度にと拡張しなくてはならない．どこまで拡張すれば正しい計算がで
きるのかは，問題に依存する．問題に応じてコンパイラやライブラリといったシステムソフトを改訂しつつ問題
を解く方法は誰もとりたくないだろう．
要求される有効数字の桁数は，整数性判定が小数点以下の数桁を見るので（絶対誤差評価的になり），固有値
の小数点よりも下に何桁の有効な数字が残るかと，因子多項式の次数に依存する41 ．特性多項式の次数が与えら
れても，因数分解後の因子多項式の次数は分からない．そこで近似固有値を 2 つの有理数型の変数で挟み，区間
で捉える．探索すべき因子多項式の次数に応じて近似固有値の精度を，有理算術演算によって必要なところまで
伸ばす方法を試す．このプログラムは，
「有理数計算プログラミング環境」が，多倍精度算術演算をサポートする
必要性を判断する上で有効な示唆を与える．また，零固有値を持つ問題（グラフ・ラプラシアン行列）の多重度
解析も行う．この行列はヘッセンベルグ変換によって複数の小行列に分割できるところまでは熱伝導行列に似て
いるが，特性多項式が低次の因子多項式に分解されにくいので，計算時間の配分が熱伝導行列の場合とかなり異
なる．
3.1
有理算術演算による精度改良の sehfv1.cpp に対する変更点
包含が成立した固有値の近似値は，区間の上限と下限が正確に抑えられる．ここで 2 分法を用いると，固有値
の存在する区間を，反復ごとに半分に狭めてゆける．テストプログラム sehfv1.cpp では倍精度型の配列 eigvl
に浮動小数点数で格納された近似値で計算を進めてきたが，これを有理数型の配列 eigvr を用意して精度改良
された固有値を格納することにする．2 分法で挟み込まれた区間の上下限をそれぞれ，i 番目のグループに属す
る固有値の下限（lower）を eiglow[i] に，上限（upper）を eigup[i] に格納する42 ．これが有理算術演算で
無理数を必要な精度まで求めて精度保証する第 1 歩である．sehfv1.cpp に対する改訂箇所を列記する．
• 包含の判定で用いた chkinc 関数のローカル変数 rational a,b を仮引数（パラメータ）に出して，挟み
込んだ区間の上限と下限として，呼出し側で保存できるようにする．つまり，chkinc 関数は，包含チェッ
40 十進
BASIC の有理数モードや Java の BigInteger クラスを用いると，大きな整数も計算することができる．
次の因子多項式を d 個の近似固有値か
ら推定するとき，
「整数かどうか」を使うし，これを使えるために整数行列にしたところがアイデアなので，絶対誤差評価から逃れられない．
42 浮動小数点数は，データ形式から小数点の位置が分かり，それが表す数値の小数点以下の桁数（精度）を判定できるが，有理数はデー
タから精度を判定することができない．有理数としての数値は（分母・分子の桁数に関係なく）つねに正確なので，その有理数が何かを近
似している場合には，何かを挟み込む区間の上限と下限が与えられなければ精度は分からない．
41 浮動小数点演算では相対誤差評価を用いることが多い．しかしここで用いる整数性の判定は，d
31
クに使用した区間 (a, b) を ratoinal 型の引数とする．ここでは参照渡しを使用する43 ．
• 組合せで使用する固有値を格納する一時配列も，倍精度型の evcan，figvl から有理数型の ercan，figvr
に変更する（変数名に rational の “r” を入れた）．
• 精度が改良された近似固有値を用いて，根と係数の関係公式で作成される係数も，倍精度型の cof から有
理数型の cor に変更する．
これらの変数型の変更に伴い，rational クラスの askinty 関数，vietaterm 関数は，多重定義によって rational
型の引数を扱う版を追加する．なお，区間は包含成立時に作成する．
3.1.1
固有値の区間での取り扱い
sehfv2.cpp で扱うデータは，浮動小数点演算で得られた近似固有値で作成したグループに，2 つ以上の根が存
在することがない（多重度の高い固有値は，別の小行列に分かれて現れる）ので，グループの下限と上限を使っ
て包含を調べる．
chkinc 関数
int chkinc(const int n, const rational_vector& p, const double e, const double h, const int mult,
引数に a と b を回す
const double eps, rational& a, rational& b){
double ta,tb;
int rtnval=0;
rational fa,fb;
a と b は引数に
以下，sehfv2.cpp の場合と同じなので省略する
根を挟み込む区間 (a, b) が引数で返るので，次のように呼び出し，eiglow[i] と eigup[i] として抑える．
区間
included[i] = chkinc(nsize,p,eigvl[i],eigvh[i],mult[i],eps,eiglow[i],eigup[i]);
−k
2 分法の反復で，区間端を引数に置けば，反復によって区間幅は反復 k 回で 2
に縮小される．次の bisect
関数呼び出しは，1 次の因子を探すために，収束判定値に，askinty 関数で使用する整数性判定値 epsint の
0.00001 倍を直接指定して精度を改良する（bisect 関数は後述する）．
区間縮小
int irtn=bisect(p,nsize,epsint*0.00001,100,eiglow[i],faeig,eigup[i],fbeig,ceig);
3.1.2
関数多重定義による askinty 関数
整数性を判定する askinty 関数は，近似固有値を倍精度浮動小数点数から有理数に変更になるので，多重定
義により有理数版を作成する．次のプログラムは 22 ページの 2.2.4 項に示した askinty 関数の有理数版である．
floor(t) と書かないで t.floor() と書くのは，floor 関数が rational クラスのメンバー関数だからである44 ．
この書き方はオブオブジェクト指向言語の書き方である．abs や numerator，denominator も同様．
43 実引数と仮引数の関係は，Fortran ではアドレス渡し（call by address）を標準として使用した．C では，誤って呼び出し側の定数な
どを変えてしまうことのないように，スカラー変数は値渡し（call by value）を取り入れ，値渡しとポインタ渡しの両者を使い分けるよう
にした．参照渡し（call by reference）は，実際にはアドレスを渡すが，更新したくない仮引数には const を付けることで，アドレス渡し
とポインタ渡しの良さを兼ね備えたものとなっている．
44 rational.h でメンバー関数として定義し，public としている．
32
rational.cpp の askinty 関数
int askinty(const rational& t, const double epsint, rational& vrat) {
rational x;
longint xd,xn,xq,xr,xh;
int rtnval=0;
rational t1=t.floor()-t;
rational t2=t1.abs();
if((double)t2 < epsint) rtnval=1;
t1=t.ceil()-t;
t2=t1.abs();
if((double)t2 < epsint) rtnval=1;
if(rtnval==1){
if(t < rational::ZERO){
x=-t;
}else{
x=t;
}
xd = x.denominator();
xh = longint::ldivide(xd,longint::TWO,&xr);
xn = x.numerator();
xq = longint::ldivide(xn,xd,&xr);
if( xr > xh) xq=xq+longint::ONE;
vrat=RRset(xq,longint::ONE);
if(t < rational::ZERO) vrat=-vrat;
}
return rtnval;
}
Fortran や C 言語などのコンパイラ型の汎用プログラミング言語では，浮動小数点数の小数部分を切り捨てた整
数を取り出す場合は，floor 関数を使うか，浮動小数点数から整数への「型変換」を用いる．これらの型変換処
理が，コンパイラによって生成されて機械語になるが，その機械語は計算機アーキテクチャ（命令セットとデー
タ形式）によってさまざまで，計算機の下位の階層である計算機アーキテクチャに強く依存する．
「有理数」からその数にもっとも近い整数を抽出するための方法も，下位の階層である longint クラスを使用
する．ここでは，longint クラスの除算を用いる．longint 型の除算は，計算機の整数除算命令に似て，被除数を
除数で割り，整数の商と剰余を返す．そこで有理数の分母と分子を .denominator() と .numerator() で取り
出し，longint 型の除算関数 ldivide で割って，剰余 xr が分母の半分 xh よりも大きければ商 xq プラス 1 を，
半分以下ならば商を有理数として返す．ldivide 関数は longint クラスの関数なので，上位の階層からこのクラ
スの関数を呼び出す場合には，スコープ解決演算子 “::” を使用してクラス名を指定する．
vietaterm 関数は，20 ページの 2.2.3 項に示したプログラムのコメントに記した変数型に変更する．変更は変
数の型だけなので省略する．
3.1.3
整数性の判定値 epsint
近似固有値をグループ化したとき，隣接するグループの固有値が近いと，2 通りのセットで同じ多項式を作ら
れてしまう危険性がある．e1 e2 のとき，e1 + e3 e2 + e3 であるから，本来は e2 + e3 で多項式係数を作ると
ころを，同じ係数を e1 + e3 で作ってしまうと，グループ 1 が包含から外れてしまうエラーが生ずる．sehfv2.cpp
では隣接グループの近似固有値の最小値を求めて，初期設定した epsint よりもこの値が小さい場合は，これに
よって置き換える45 ．これはグループ化の直後に入れる．
45 グラフ・ラプラシアン行列で “G 9” の場合に，グループ 7,12,22,28,35 で 5 次多項式 (λ5 − 14λ4 + 69λ3 − 142λ2 + 112λ − 29) が生
成されるが，e11 = .7560 · · · で e12 = .7561 · · · と 0.01 以内にあるため，epsint=0.01 ではこの処理なしではエラーとなる．
33
epsint の固有値の接近による置換え
int ismal=0;
for (i=1; i<=ngrp1; i++){
if(i==1){
eveps=epsint;
}else{
if(fabs(eigvl[i]-eigvl[i-1]) < eveps){
eveps = fabs(eigvl[i]-eigvl[i-1]); ismal=i;
} } }
if(epsint > eveps) epsint=eveps;
3.1.4
精度要求の判定値
2 次以上の因子多項式については，近似固有値の候補を選び，根と係数の関係公式を使用して係数を作成して，
その整数性を判定する．複数の e の演算結果で係数ができるので，e に対する直接的な精度要求値は見えなくな
る．このときに必要となる近似固有値の精度は，因子多項式の項の次数に比例して増える．これに誤差の伝搬の
影響を考慮する．
√
3 の例で，有効数字を 3 桁として，e2 1.27，e16 4.73 として，近似固有値の精
度の問題を考えてみよう．この例では，抽出されるべき 2 次式の係数は −6 と 6 である．近似固有値の和は
前述した固有値 3 ±
−r = e2 + e16 = 6.00 で積は s = e2 · e16 = 6.0071 になる．係数行列を 10 倍すると固有値も 10 倍されて 12.7
と 47.3 になり，和は 60.0，積は 600.71 になる．次の桁に誤差が混入していて，12.74 と 47.34 だとしたら，和
は 60.08 だが積は 603.1116 で，積の整数性は正しく得られない．つまり，整数性判定は，近似固有値の精度が，
固有値の小数点よりも上の桁数と，根と係数の関係公式で使用する根の総積の項数に依存する．倍精度浮動小
数点数を有理数化した場合は LCM の桁数が大きくなり，固有値自身も桁数の大きな数になる．行列 (17) の場
合は e1 = 7, 205, 759, 403, 792, 791 と 16 桁の数になるので，少なくとも 16 桁以上の精度が要求される．この
ときの因子多項式の定数項がもっとも厳しい条件を与える．因子多項式の次数を d とした場合，定数項は総積
d
pd =
αk である．
k=1
プログラムでは，固有値の近似値を ei から，小数点よりも上のビット数（2 進数表現したときの桁数）を
log2 ei + 1 によって得る．計算誤差の伝搬を考慮してこれに 2 を加えてから d 倍した値を用いる．
d · (log2 ei + 3)
(27)
因子多項式の次数は特性多項式からは分からないので，既知の因子で割った後の多項式に対して，まず 1 次式用
に設定する．固有値が有理数なら 1 次式で割りきれる．
次に候補の因子多項式の次数を，2 次から 3 次，4 次と上げてゆくループの中で，次数に応じて精度を高める．
次数 d に関するループ中に精度改良を入れるので，不必要に桁数を増加させることを食い止められる．
判定式 (27) が機能しているかどうかを確かめる目的で，プリプロセッサの変数 ASKINTY を定義（#define）し
てコンパイルすると，判定値の余裕を見ることができる46 ．
3.1.5
包含の判定
4 × 4 分割を例とする（表 2）．はじめに 16 個の近似固有値をグループ化して 9 つの近似値が得られる．ヘッ
センベルグ小行列をフロベニウス変換して特性多項式を得た段階で，近似固有値の各グループの上限に ε を加
え，下限から ε を引いて，区間 (ai , bi ) を定める．特性多項式は既知の因子で割ってから，包含を調べられる．
2 番目に処理する H3 の特性多項式は，4 − λ で割りきれて，商多項式で置き換えられ p(λ) = λ2 − 8λ + 15 に
46 判定式
(27) を d · (log2 ei + 1) にすると，“s 12” で正しく 12 次多項式を 2 つの 6 次多項式の積に分解できない．
34
なっている．p(ai ) · p(bi ) < 0 なら区間 i に 1 つ，あるいは奇数個の根がある．グループをすべて調べて，この
次数が 2 の特性多項式 p(λ) に対して，このような区間が 2 つ見つかれば，包含が成立する．4 × 4 分割の 2 番
目に計算す小行列 H3 の場合は， e4 = 3 と e6 = 5 だけが包含条件を満たす．
もとの行列が整数の場合は，浮動小数点計算で得られた精度で因子探しは間に合う場合があるが，倍精度浮
動小数点計算では精度が不足する．例えば “s 4” の場合は，2 番目に処理する H3 の特性多項式は，既知の因子
(14, 411, 518, 807, 585, 588 − λ) で割られて
p(λ) = λ2 − 28, 823, 037, 615, 171, 176λ + 194, 711, 132, 195, 056, 057, 687, 171, 823, 724, 135
になる47 ．
包含は桁数が増えても同様に成立する．有理算術演算の威力は，数学的に成立する事柄を，数の大小に関係な
くプログラムで確かめられる点にある．
3.1.6
精度改良
必要な精度が与えられ，包含が成立すれば，近似固有値の精度は 2 分法反復で改良できる．2 分法は高校の数
a+b
学 B の教科書にも載っている素朴なアルゴリズムで，根の存在する区間 (a, b) の中点 c =
で関数値を計
2
算して，p(a) · p(c) < 0 なら区間 (a, c) を，p(c) · p(b) < 0 なら区間 (c, b) を選択する．したがって，反復ごとに
区間幅を半分にしてゆける．
収束判定には，関数値を使うか，区間幅を使うかの選択肢がある．正確な計算により多項式の零点を挟む区間
を縮小するには，区間幅を収束判定値に用いるほうが使いやすい．この版の bisect 関数を rational クラスに含
めた．これを用いると，区間 (a, b) に根が包含されたら，区間幅を指定した値まで確実に縮小できる．戻り値は
反復回数とした．指定された反復回数でも区間幅が tol 以下にならなかったら，反復回数の (−1) 倍を返し，根
が有理数で根に到達したら 0 を返す．なお，Horner も rational クラスの関数である（前述）．
47 28, 823, 037, 615, 171, 176
÷ 8 = 3, 602, 879, 701, 896, 397 で，式 (16) の 0.1 の分子である．
35
rational.cpp の bisect 関数
int bisect(const rational_vector& p, const int n, const double tol, const int maxit, rational& a,
rational& fa, rational& b, rational& fb, rational& c){
int it=0, rtnval;
fa=Horner(n,p,a); fb=Horner(n,p,b);
while( it < maxit){
if(fa*fb > rational::ZERO){
std::cout << " bisect error a=" << setprecision(16) << (double)a << " b=" << (double)b << " fa="
exit(3);
}
it++;
if(it >= maxit){
std::cout << " bisect MAXIT reached it=" << it << std::endl;
rtnval=-it; break;
}
c = (a + b) / rational::TWO; rational fc=Horner(n,p,c);
if(fc == rational::ZERO){
cout << " bisect converged it=" << it << " fc==rational::ZERO" << endl;
a=c; fa=fc; b=c; fb=fc;
rtnval=0; break;
}
if(fa*fc > rational::ZERO){
a=c; fa=fc;
}else{
b=c; fb=fc;
}
if(it%50==0) cout << " bisect it=" << it << " b-a=" << (double)(b-a) << endl;
if((double)(b-a) < tol){
cout << " bisect converged it=" << it << " b-a=" << setprecision(16) << (double)(b-a) << endl;
rtnval=it; break;
}
}
return rtnval;
}
bisect 関数は，1 次の因子を見つけるための精度改良として，区間幅が epsint よりも十分に小さくなるように
指定して使用する．
sehfv2.cpp の bisect の使用
faeig=Horner(nsize,p,eiglow[i]); fbeig=Horner(nsize,p,eigup[i]);
if(faeig*fbeig > rational::ZERO){
cout << " ERROR check Horner .....faeig=" << faeig << " fbeig=" << fbeig << endl;
polytexform(p,nsize); exit(3);
}
//********* bisect to refine eigenvalue *********************
int irtn = bisect(p,nsize,epsint*0.00001,100,eiglow[i],faeig,eigup[i],fbeig,ceig);
//***********************************************************
cout << " Return from bisect irtn=" << irtn << endl;
included[i] = abs(irtn);
eigvr[i]=ceig; mxitr[i]=mxitr[i]+abs(irtn);
因子多項式が 2 次以上になると，while ループの中で，多項式の次数 d に応じて精度改良を行う．収束判定値
（区間幅）を零として使用し，常に指定した反復回数で最大反復回数となって，区間を 2 のマイナス反復回数乗
に狭める．浮動小数点演算と異なり，有理算術演算では誤差なしで計算でき，分母と分子の桁数も拡張されるの
で，区間幅は厳密に反復ごとに半分になる．したがって，精度要求を 2 進数の桁数で求めているのであるから，
不足するビット数と同じ回数だけ反復すれば，要求精度を満たす．
36
sehfv2.cpp の精度改良
for(i=1; i<=ngrp1; i++){
if(included[i]>=1){
log2ev=log2(eigvl[i])+1;
int reqndeig=((int)log2ev+2)*d+accdouble; // <============
int niteration=reqndeig-mxitr[i];
if(niteration<2) niteration=2;
cout << endl << "required accuracy i=" << i << " for eigval=" << eigvl[i] << " log2ev=" << log2ev << "
int iwas = chkonce(nsize,p,eiglow[i],eigup[i],faeig,fbeig);
if(iwas==0){
cout << "chkonce ERROR then STOP" << endl;
cout << " REfine accuracy of eigval" << " Width of interval=" << (double)(eigup[i]-eiglow[i]) << e
exit(3);
}
double dzero=0;
bisect(p,nsize,dzero,niteration,eiglow[i],faeig,eigup[i],fbeig,ceig);
cout << " Refined accuracy of eigval eigvr=" << eigvr[i] << " ceig=" << ceig << " aeig=" << (double)ei
eigvr[i]=ceig; mxitr[i]=reqndeig;
cout << " Refined accuracy of eigval eigvr=" << eigvr[i] << " " << (double)eigvr[i] << " Width of inter
cout << "included[" << i << "]=" << included[i] << " mxitr[i]=" << mxitr[i] << " maxit=" << maxit << e
}
}
この精度改良された近似固有値を使用して 2 次以上の因子探しを行う48 ．
48 chkonce 関数は，2 つの変数に対する関数値を Horner 関数で得て，2 つの関数値の積が正のときは戻り値を 0 にする．これはエラー
検出が目的である．プログラムの開発段階では，エラーの早期発見は，デバッグの効率を高める効果がある．
37
sehfv2.cpp の 2 次以上の因子探し
while(d<=nsize){
// search for polynomial of degree 2 to nsize
上記の精度改良
ncandi=setecand(eigvr,ngrp1,included,ercan,usev); //sehfv2:
for(i=1; i<=d;i++) idx[i]=i;
// initialize for nxtcmb
idx[d]=idx[d]-1;
// initialize for nxtcmb
nn=comb(ncandi,d);
// since nn is type "long long", ngrp1-nuse should be "< 62"
Top:
for(ii=1; ii<=nn;ii++){
// loop over nn combinations of candidates
nxtcmb(ncandi,d,idx);
// select one of _(ncan-nused) C_d candidates
ur=rational::ZERO;
//sehfv2: u ==> ur
for(j=1; j<=d; j++){ figvr[j]=ercan[idx[j]]; ur=ur+figvr[j];}//sehfv2:
if(askinty(ur,epsint,cor[1])){ //sehfv2: pass the first cofficient check, then ...
int ok=1;
for(int r=1; r<=d; r++){
cor[r]=vietaterm(figvr,d,r);
//sehfv2: cof,figvl ==> cor,figvr
if(! askinty(cor[r],epsint,vrat[d-r])){ ok=0; break; }
}
if(ok){
// OK, try trial division
vrat[d]=rational::ONE;
for (j=0; j<=nsize; j++) urat[nsize-j]=p[j];
if(polydivchk(urat,nsize,vrat,d,qrat,rrat){;
// trial division divided up
for (j=0; j<=nsize-d; j++) p[nsize-d-j]=qrat[j]; // set quatient for next p(\lambda)
nsize = nsize-d;
// and decrease nsize
cout << " Divisor = " ; polytexformintr(vrat,d);
npolynomial=putpolynomial(d,vrat,npolynomial,ppolynomial,cpolynomial);
for (j=1; j<=d; j++) included[usev[idx[j]]]=0;
ncandi=setecand(eigvr,ngrp1,included,ercan,usev); //sehfv2: reset number of candidate eigenvalues
noff=idx[1]-1; // number of OFFed indecies
for(j=1; j<=d;j++) idx[j]=j+noff; // initialize for nxtcmb
if(idx[d]>ncandi) break;
// there are no degree d left, then break for loop
idx[d]=idx[d]-1;
// initialize for nxtcmb
nn=comb(ncandi-noff,d);
// reset nn (number of combinations)
goto Top;
/////////////// to loop Top //////////////////
}
// if(divup)
}
// if(ok){
}
//
if(areyouint(u,epsint))
}
//
for(i=1; i<=nn;i++){
d++;
}
//
while(d<nsize){
3.2
計算例
実行時の入力データは，コマンド行引数に “T 4” のように英文字で行列の種類を，数字で分割数を指定する．
T が sehfv1 と同じ熱伝導行列，小文字の t が 10 進数の 10 で割った熱伝導率，s が 2 進数の 10 で割った熱伝導
率を使用する熱伝導行列を指定する．G が後述するグラフ・ラプラシアン行列を指定する．
“s 11” の場合を調べてみる．H55 , H37 , H15 , H7 と 7 つの H1 に分割される．この小行列のサイズは 10 ペー
ジの表 3 の m = 11 の場合に一致する．H7 の特性多項式は H1 の特性多項式で割り切れて 6 次式になり，2 次
式と 4 次式で割り切れる．H15 の特性多項式は H7 の特性多項式で割り切れて 8 次式になり，2 つの 1 次式と 3
つの 2 次式で割り切れる．H37 の特性多項式は既知の 1 次，2 次，4 次多項式因子で割り切れて 30 次式になり，
3 つの 2 次式と 6 つの 4 次式で割り切れる．30 C2 = 435 で，割り切れる組合せは，131 番目，12 番目に現れる．
24 C3 = 2, 024，24 C4 = 10, 626 で，割り切れる組割り切れる組合せは，1587 番目，733 番目，383 番目，63 番目
に現れる．H55 の特性多項式は既知の 3 つの 1 次，7 つの 2 次，6 つの 4 次多項式因子で割り切れて 10 次式に
なり，2 つの 1 次式と 2 つの 2 次式と 4 次式で割り切れる．この場合は，特性多項式の次数は高いが，実際の除
算は 4 次式までで済むので，計算時間も短い．
“s 10” の場合は，H51 , H41 と 8 つの H1 に分割される．H41 の特性多項式は H1 の特性多項式で割り切れて
40 次式になり，8 つの 5 次式で割り切れる．40 C5 = 658, 008 で，割り切れる組合せは，66,903 番目，36,980 番
目，18,032 番目，5,680 番目，624 番目に現れる．H51 の特性多項式は H41 の特性多項式で割り切れて 10 次式
38
11
7
3
12
8
9
4
41 42 43 44 45
10
5
6
34
40
27
33
20
26
13
19
6
1
2
12
1
(a) 4 × 4 分割
2
3
4
5
(b) 7 × 7 分割
図 2: グラフ・ラプラシアン行列
になり，2 つの 5 次式で割り切れる．次数 40 を相手に，5 次多項式で割ってゆくので，“s 11” の場合よりも組
合せの数が多く，また 5 次式なので近似固有値の精度も高くなり，桁数の多い計算をしなけらばならない．計算
時間を比較すると，11 の場合よりも 5 倍以上を要する．プログラムでは cof[r] に pi が得られるが，これが許
容範囲内で整数に収まるかどうかを調べる．次数が大きくなるとこの計算量は膨大になるので， p1 に関してだ
け先行して行い，これに通らないものは調べない．p1 の先行判定を行わないと，m = 10 のケースでは，パソコ
ンで約 1.6 秒の処理が 10.4 秒に伸びた．
3.3
零固有値のある問題
“G m” の入力で，熱伝導行列の境界条件を指定しない行列を扱う．m は分割数で図 2 のように節点番号を数
える．節点数は (m − 2)2 + 4(m − 2) になる．m = 4 の場合は 12 節点で，図左のように接続したグラフ・ラプ
ラシアン行列が生成される．
⎛
−1
1
⎜
⎜
⎜
⎜ −1
⎜
⎜
⎜
⎜
A=⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
1
−1
1
−1
−1
4
−1
⎞
−1
−1
4
−1
−1
−1
−1
1
−1
1
−1
−1
4
−1
−1
−1
−1
4
−1
−1
−1
1
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
−1
⎟
−1 ⎟
⎟
⎟
⎠
1
(28)
1
グラフ・ラプラシアン行列に，境界条件として境界固定を与えると熱伝導率が 1 の熱伝導行列になる49 ．
固有値の分布は表 4 のようになる．
H2 から得られる特性多項式は (λ2 − 5λ + 2) である．H7 から得られる特性多項式は (−24λ + 134λ2 − 267λ3 +
234λ4 − 92λ5 + 16λ6 − λ7 ) と，定数項がなく，零固有値を含むことが分る．これは H1 と H2 の特性多項式で
割り切れて (λ4 − 10λ3 + 25λ2 − 12λ) が得られ，λ(λ − 3)(λ2 − 7λ + 4) と分解される．境界条件を設定すると固
有値は 4 になり，多重度は分割数に増えるところが興味深い．
表 5 に m を 4 から 12 までの結果を示すが，m = 7 を例外として，多重度 4 は最後の 3 つの小行列と Ha に
現れる．表では 1 − λ の 1 次式の因子が複数現れる場合を “3*1” のように表した．熱伝導行列と比較すると，次
数の高い因子があるため，計算時間が長くなる．
49 行列
(28) の対角項が 4 の行と列だけ残して，それ以外の行と列を取り除くと，式 (10) の行列になる．
39
表 4: 4 × 4 分割のグラフ・ラプラシアン行列の固有値
λi
λ1
λ2
λ4
λ5
λ9
λ10
λ3
exact λ
2.37e-16
0
0.438447187
0.627718673
λ6
λ7
λ8
λ11
λ12
H7
approx. ei
0.99999999
2.99999999
4.561552812
6.392281323
H2
H1
H1
√
0.5(5 − 17)
√
0.5(7 − 33)
1
3
√
0.5(5 + 17)
√
0.5(7 + 33)
H1
表 5: グラフ・ラプラシアン行列の分割数 m と小行列サイズと因子多項式の最高次数 d
m
4
5
6
7
8
9
10
11
12
n
12
21
32
43
60
77
96
117
140
Ha
7
12
18
25
33
42
52
63
75
Hb
Hc
2
3*1
6
3*1
7
4
15
1
14
10
32
3*1
23
18
42
9
34
28
3*1
1
2
3*1
3*1
3*1
3*1
7
10
23
42
Hd
He
Hf
d
1
2
4
14
40
18
m = 9 を見ると，H42 ，H32 と 3 つの H1 に分かれ，m = 2 の場合と同様に H32 は既知の因子で割りき
れない．32 次の特性多項式は 5 次，9 次，18 次の多項式に分解される．32 C5 = 201, 376 で，因子は組合せが
111,182 番目である．27 C6 = 296, 010，27 C7 = 888, 030，27 C8 = 2, 220, 075，27 C9 = 4, 686, 825 で，組み合わせ
が 2,264,208 番目である．H42 は 1 次の因子 (1 − λ) と H32 の因子で割り切れて 9 次となり，これが 1 次（λ），
4 次，4 次と分解される．このように，熱伝導行列の場合に比較すると，因子多項式の見つかるまでの反復回数
が増える．
m = 10 を見ると，H52 ，H23 ，H18 と 3 つの H1 に分かれる．H18 は 2 つの 9 次式に分解されるが，H23 は分
解されない．H52 は 1 次の因子 (1 − λ) と H18 の 2 つの 9 次式と H23 の因子で割り切れて 10 次となり，これ
が 1 次（λ），4 次，5 次と分解される．そのため，計算時間の大部分は 23 次多項式の係数の作成にかかる．こ
の計算で，近似固有値の精度を高めるところを e2 = 0.08284674618487071 に着目して見てみる．はじめは区間
幅は 10−11 程度であるが，
eiglow=
eigup =
0.0828467461848711311800386
0.0828467462262947110240852
|
23 次多項式の係数の整数性を求めるためには 10−28 程度にまで精度改良されている．
eiglow=
0.082846746184871488843615432956758127913444
eigup =
0.082846746184871488843615433531625779229182
|
この精度は 4 倍精度演算の限界に近い．
3.4
チェックポイントとリスタート
この問題では，ヘッセンベルグ変換やフロベニウス変換にはあまり計算時間がかからないので，putpolynomial
関数で格納した因子多項式だけをチェックポイントをとった．“#define CHKPTRST” の指定によって，チェック
ポイントとリスタート機能の使用例を見られる．“G 9” の場合を例にする．ヘッセンベルグ小行列が H42 , H32 ,
H1 , H1 , H1 の 5 つに分かれた場合に，H42 だけを再計算したい場合に，H32 , H1 , H1 , H1 に関する計算を行わ
ずに，それらの特性多項式の因子を，チェックポイント・ファイルから読み込み，計算時間の大半をスキップす
ることができる．変更箇所を示す．
• チェックポイント・ファイル名と変数 krst の指定（コマンド行引数）
• ファイルのオープン
• 既知の因子の読み込み
• 因子多項式のチェックポイント・ファイルへの書き込み
コンパイルは，CHKPTRST を指定し，リンクに archive.o を加える．実行は，“G 9” を例にすると，“G 9 filename
0” でチェックポイントがとられる．最後の引数が変数 krst に格納され，この値が零で，小行列が最後のもので
なければ，指定されたファイルに見つかった多項式因子を書く．“G 9 filename 1” でリスタートすると，最後以
外の小行列の特性多項式の因子探しをスキップして，チェックポイント・ファイルから H32 , H1 , H1 , H1 の 4 つ
の小行列の因子探しの結果を読み込んで，H42 についてだけ計算する．
oarchive
template<class T> oarchive& operator<<(oarchive& oa, const vector<T>& v) {
oa.getOs() << " " << v.len();
for (int i=0; i < v.len(); ++i) oa << v[i];
return oa;
}
41
4
倍精度浮動小数点行列の解析プログラム sehfv3int.cpp
2 次元トラス構造解析プログラム CT2D から抽出した，倍精度浮動小数点演算で作成された対称行列の固有値
の多重度を解析する．前章の行列との違いは，ヘッセンベルグ小行列から得られる特性方程式が病的に近接した
根を持つ点にある．前章の行列の多重度の問題は，ヘッセンベルグ変換で解決され，小行列の特性方程式は近接
根はなかった．本章で扱う CT2D の行列は，ヘッセンベルグ変換後の小行列の特性方程式が，1 つのグループに
複数の根をもつ難しさがある．sehfv2.cpp では近似固有値を 2 つの rational 型の変数で挟み込み，この区間を縮
小することで，根と係数の関係を使用できるところまで精度改良すれば，因子多項式を見つけられることを確か
めた．これは有理算術演算で無理数を扱う第 1 歩である．sehfv3int.cpp ではさらに次の項目を試みる．
• 区間算術演算の応用範囲を広げ，vietaterm 関数が生成する多項式の係数も rational 型の変数で扱う「点」
ではなく，2 つの rational 型の変数で表される「区間」として扱う．区間で表された 2 つの数の加算や乗
算といった「区間算術演算」（interval arithmetic）も使用して，因子多項式の係数も区間演算で扱う．
• 精度改良のための包含では，関数値を使うだけでは分離が難しく，導関数も使用して近接根を分離する．
• 2 分法を反復回数で制御するのではなく，区間幅そのものを 2 分法の収束判定値に使用することで，2 分法
以外の精度改良法を使用できるようにする．
• 「有理数計算プログラミング環境」が，倍精度浮動小数点演算以外に，4 倍精度，8 倍精度といった多倍精
度演算をサポートすべきかどうかを評価する．
プログラムの開発では，信頼できるプログラムの結果と比較して正しい稼働を確かめながら進められると効率
がよい50 ．前章の因子探しプログラムの作成では，テストプログラム sehfv1.cpp と sehfv2.cpp の実行結果を比
較しつつ進め，また数式処理システム Maxima も利用した．sehfv3int.cpp では，プリプロセッサに対する変数
GENFTHERM で切り替えて，sehfv2.cpp で稼働した有理行列に対する計算が，プログラムの改変で影響を受けな
いことを確認しつつ，CT2D が生成した倍精度浮動小数点行列の解析機能を組み込む．
本章でははじめに区間算術演算を説明する．次節で浮動小数点演算による倍精度の行列をテキストファイル経
由で読み込む方法を説明する．この行列のヘッセンベルグ変換とフロベニウス変換は計算時間がかかるので，3
節で，チェックポイントとリスタートについて触れる．4 節で，1 グループに最大で 3 つの根をもつ特性多項式
の包含方法を説明する．5 節で計算例を示すが，浮動小数点行列が，コンパイルオプションの違いによって固有
値の多重度に影響を与える現象を解明する．なお次章で，このプログラムに対する簡単なチューニングを行って
計算時間の短縮を試みる．
4.1
区間算術演算の応用
sehfv2 で近似固有値を挟む区間の下限と上限を得たが，sehfv3int ではこれをコンパイラに対して，区間とし
てまとめて interval 型の数として認識させる．有理数型の 2 つの配列 eiglow(n+1) と eigup(n+1) を区間型
の配列 eigint(n+1) に置き換える．
//
rational_vector eiglow(n+1), eigup(n+1);
// rational 型の下限と上限を
vector<interval> eigint(n+1);
// interval 型に
中略
rational tempa,tempb;
included[i] = chkinc(nsize,p,eigvl[i],eigvh[i],mult[i],tempa,tempb); 包含が成立したら
nexist = nexist + included[i];
eigint[i]=interval(tempa,tempb);
//
interval 型の配列要素に格納
50 信頼できる参照プログラムがない場合は，数学的に計算可能な理論解や，筆算や電卓で行う計算結果と比較する．
42
これ以降，下限と上限を 2 つの rational 型の変数で扱っていた区間を，1 つ interval 型の変数で扱う．したがっ
て rational クラスの Horner と bisect 関数は，区間を引数とする interval クラスの版を呼ぶ．
Horner(nsize,p,eigint[i],faeig,fbeig);
int irtn = bisect(p,nsize,eps,100,eigint[i],faeig,fbeig,ceig);
formatprn(eigint[i]);
// 引数は interval
// 引数は interval
formatprn(eigint[i]) は interval 型の数（区間）を表示する（後述）．
1 次の因子探しの整数性判定は次のように書ける．interval 型変数の上限値を interval クラスのメンバー関数
upper() で抽出し，rational クラスのメンバー関数 floor() で小数点以下を捨てる．
if(eigint[i].upper().floor() == eigint[i].lower().ceil()){
精度改良された近似固有値も区間を引数として eab に渡し，rational 型の近似固有値の下限 ercana と上限
ercanb に入れる（sehfv2 では setecand 関数で rational 型の配列 ercan に詰め替えただけ）．
int setecanint(const rational_vector& eigvr, const vector<interval>& eab, const int ncan,
const int included[], rational_vector& ercan, rational_vector& ercana, rational_vector& ercanb,
rational_vector& efcana, rational_vector& efcanb, int usev[]){
int j=0;
for(int i=1; i<=ncan;i++){
if(included[i]){
// exclude eigvl[i] according to included[i]
j++; usev[j]=i; ercan[j]=eigvr[i]; ercana[j]=eab[i].lower(); ercanb[j]=eab[i].upper();
efcana[j]=ercana[j]-ercana[j].floor(); efcanb[j]=ercanb[j]-ercanb[j].floor();
};
}
return j;
// return number of candidates set in array evcan[1] to [j]
}
以上の変更は，2 つの rational 型を 1 つの interval 型に詰めるだけの改訂作業であるが，区間をコンパイラが数
として扱うので，演算子多重定義により，vietaterm の “s=s*e[ith[j]] のステートメントは変数型が interval
になっているのでそのまま使用できる．nxtcmb で d 個の候補が決まった時点で，改めて区間配列 figv に格納
し，1 次の係数の整数性判定を行い，パスすれば根と係数の関係から係数を作成する．
for(ii=1; ii<=nn;ii++){
// loop over nn combinations of candidates
nxtcmb(ncandi,d,idx);
// select one of _(ncan-nused) C_d candidates
for(j=1; j<=d; j++){
figvr[j]=ercan[idx[j]];
figv[j]=interval(ercana[idx[j]],ercanb[idx[j]]);
区間に代入
figvaf[j]=efcana[idx[j]]; figvbf[j]=efcanb[idx[j]];
}
ura=figvaf[1]; urb=figvbf[1];
for(j=2; j<=d; j++){
ura=ura+figvaf[j]; urb=urb+figvbf[j];
先行判定のための 1 次の係数
}
if(urb.floor() == ura.ceil()){
for(int r=1; r<=d; r++){
rational ca,cb;
//**** make coefficieants terms of Vieta’s formula from figvr ***********
vietaterm(figv,d,r,ca,cb);
// 引数は interval
根と係数の関係公式で，近似固有値が区間で与えられた場合は，それらを根とする方程式の係数を作成する．
n 個の近似固有値の候補から r 個を選ぶ組合せを作成して，r 次の多項式係数を作成する vietaterm の区間算
術演算版を示すが，近似固有値の変数型が interval なので，
（実行文は同じだが）係数も interval 型で得られる．
43
区間算術演算（interval クラス）の vietaterm 関数
int vietaterm(const vector<interval>& e, const int n, const int r, rational& ta, rational& tb) {
interval s;
int ith[r+1];
for(int i=1; i<=r;i++){ith[i]=i;}
ith[r]=ith[r]-1;
interval t=interval(rational::ZERO,rational::ZERO); // (0,0) の区間
long long nn=comb(n,r);
for(long long i=1; i<=nn; i++){
nxtcmb(n,r,ith);
s=interval(rational::ONE,rational::ONE);
// (1,1) の区間
for(int j=1; j<=r; j++) {
if(e[ith[j]].lower().s * e[ith[j]].upper().s < 0){
std::cout << "rational.vietaterm ERROR STOP" << std::endl;
exit(3);
}
区間乗算
s=s*e[ith[j]];
}
区間加算
t=t+s;
}
if(r%2==1) { t=t.neg(); }
ta=t.lower(); tb=t.upper();
return 1;
}
結果は ta と tb の有理数で返し，これを次のように整数性判定で使用する．
区間演算での整数性判定
if(cb.floor() == ca.ceil()){
cor[r]=cb.floor();
}else{
std::cout << "cof[" << r << "] != its a and b
ok=0; break;
}
exit loop" << std::endl;
係数を区間に変更することで，整数性判定に，epsint のように閾値の設定が難しい値を不要にできた．
ただし，有理算術演算では floor() 関数は，longint 型の除算を行うため短時間で処理できない．このため，
d
αk
実行回数が多いカーネルではチューニングが必要になる．例えば，根と係数の関係公式の先行判定 p1 =
k=1
では，小数点以下だけを加えればよいが，近似固有値の小数点以下をあらかじめ用意しておくと，組合わせの数
が多い場合に効果が大きい．sehfv2.cpp で setecand 関数で準備した包含された近似固有値の詰め替えを，次の
setecanint 関数に改訂して，小数点以下のみの近似固有値を floor() 関数を用いて作成して efcana と efcanb
に返す．
44
setecanint 関数
int setecanint(const rational_vector& eigvr, const rational_vector& ea, const rational_vector& eb,
const int ncan, const int included[], rational_vector& ercan, rational_vector& ercana,
rational_vector& ercanb, rational_vector& efcana, rational_vector& efcanb, int usev[]){
int j=0;
for(int i=1; i<=ncan;i++){
if(included[i]){
j++; usev[j]=i; ercan[j]=eigvr[i]; ercana[j]=ea[i]; ercanb[j]=eb[i];
efcana[j]=ercana[j]-ercana[j].floor(); efcanb[j]=ercanb[j]-ercanb[j].floor();
};
}
return j;
// return number of candidates set in array evcan[1] to [j]
}
この関数を呼び出し，小数点以下に桁数を詰めた候補近似固有値を使用して整数性の先行判定を行う．
sehfv3.cpp の区間演算による整数性判定
ncandi=setecanint(eigvr,eiglow,eigup,ngrp1,included,ercan,ercana,ercanb,efcana,efcanb,usev);
for(i=1; i<=d;i++) idx[i]=i;
// initialize for nxtcmb
idx[d]=idx[d]-1;
// initialize for nxtcmb
nn=comb(ncandi,d); nok=0;
// since nn is type "long long", ngrp1-nuse should be "< 62"
Top:
if(nok != 0) cout << endl << "%reset comb(" << ncandi << "," << d << ")=" << nn << endl << endl;
cout << "Vieta’s formula check for d=" << d << " ngrp1=" << ngrp1 << " ncandi=" << ncandi << endl;
for(ii=1; ii<=nn;ii++){
// loop over nn combinations of candidates
nxtcmb(ncandi,d,idx);
// select one of _(ncan-nused) C_d candidates
for(j=1; j<=d; j++){
figvr[j]=ercan[idx[j]]; figva[j]=ercana[idx[j]]; figvb[j]=ercanb[idx[j]];
figvaf[j]=efcana[idx[j]]; figvbf[j]=efcanb[idx[j]];
}
ura=figvaf[1]; urb=figvbf[1];
for(j=2; j<=d; j++){ ura=ura+figvaf[j]; urb=urb+figvbf[j]; }
if(urb.floor() == ura.ceil()){
4.2
倍精度浮動小数点数行列の読込み
浮動小数点計算の途中の行列を抽出して有理数変換して正確な計算を実行することが「有理数計算プログラミ
ング環境」の目的である．抽出された行列要素の下位の桁は，浮動小数点計算のさまざまな条件の組合せで変化
する．これらの条件による摂動を分析することで，誤差の分布の原因を探す．浮動小数点計算プログラムとして，
構造解析プログラム CT2D（C program for Truss 2 Dimensional analysis）を使用する．CT2D は C 言語で書
かれた単純はプログラムで，入手も可能（ダウンロード可能）である [3, p. 68]．
剛性行列 K と質量行列 M を CT2D で作成する．ただし質量行列は，CT2D プログラムに機能追加して，固
有振動モードの解析を可能にして作成した．solveSTATIC.c を変更する形で NormalMode.c プログラムを作成
し，このプログラムの NormalMode 関数が massmat 関数（NormalMode.c 中）を呼び出し，そこで質量行列
を生成する．質量行列 M は，トラス要素の質量（断面積，密度，部材長の積）を両節点に等配分する集中質量
行列（対角行列）とした．したがって 1 次元配列に格納される．さらに，NormalMode は，eigv.c プログラム
の genamat 関数を呼び出す．genamat は n と 行列 K と M と A = S −1 KS −1 を，非零要素だけを，行番号，
列番号，行列要素の組でファイルに出力する．その後，jacobev 関数で固有値と固有ベクトルを求め，固有値は
ソートする51 ．
51 jacobev
は sehfv1.cpp などに入れた jacob2 に固有ベクトルの計算を加えた．ヤコビ法の説明を付録に記した．
45
調和振動（定常解）を表す 2 階の微分方程式 M ü(t) + Ku(t) = 0 は一般化固有値問題
ω 2 M u = Ku,
(29)
となる52 ．M は対角行列なので，M の対角項の平方根を対角項とする対角行列 S によって
A = S −1 KS −1 ,
x = Su,
(30)
と対称性を保存したまま固有方程式 (1) に変形できる53 ．
なお sehfv3int.cpp では行列は次の 3 通りを解析できる（後述するコマンド行入力データの mattype に 0，1，
2 のいずれかを指定する）．
• A = S −1 KS −1 を倍精度で計算してから有理数変換する（mattype = 0）．
• A = K として有理数変換する（mattype = 1）．
• K と M を有理数変換してから，A = M −1 K を有理数計算する（mattype = 2）．
4.2.1
浮動小数点行列のファイルへの出力
例題集の CT2Dnorm ディレクトリに，トラスデータ作成プログラム trusgenorm.f があり54 ，gfortran でコ
ンパイルして nx と ny に 2 を与えると，最小の 13 節点のデータが作成される．これを S13Free.dat とする．2
つのコマンド “gcc -c *.c” と “gcc -lm *.o” で作成した a.exe を実行すれば，S13free.txt が作成される．ファイ
ル名は CT2D のメインプログラム Clientmain.c に埋め込んだので，コンパイルオプションを変更する場合
は，その都度，Clientmain.c を変更して使用する．
IEEE 標準の倍精度浮動小数点数は，書式 “%25.17le” を指定したテキストファイル経由で渡す．10 進数で 17
桁の有効桁数は，IEEE 倍精度浮動小数点数を正確に復元できる [4]．CT2D の eigv.c の genamat 関数を示す．
void genamat(double * mvec, double *b, double *a, int n, int nb, FILE *fp_mat){
int i, j, j1;
double scal;
for(j = 0; j < n ; j++){
for(i = 0; i < n; i++) a[j*n+i]=0;
}
for(j = 1; j <= n ; j++){
j1=MAX(j-nb+1,1);
for(i=j1; i<=j; i++) a[(j-1)*n+i-1]=b[j*nb-j+i-1];
}
/*** fill lower portion by transpose
***/
for(j = 0; j < n ; j++){
for(i=j+1; i<n; i++) a[j*n+i] = a[i*n+j];
}
/*** file output n and K matrix ***/
fprintf(fp_mat,"%d\n",n);
for(i=0; i<n; i++){
for (j=0; j<n; j++){
if(a[n*j+i] != 0){;
fprintf(fp_mat," %d",i);
fprintf(fp_mat," %d",j);
fprintf(fp_mat,"%25.17le\n",a[n*j+i]);
} } }
52 調和振動では
u(t) = u cos ωt のようにおくことができる（u は固有振動モードの振幅）．
= S 2 とおき，Ku = λSSu の両辺に S −1 を掛けて S −1 KS −1 ( Su ) = λ( Su ) と変形し x = Su とおくことで固有方程式
|
{z
} |{z}
|{z}
x
x
A
Ax = λx が得られる．
54 このデータ生成プログラムは，外周にない垂直部材（節点 2 と 7，7 と 12）と水平部材（節点 6 と 7，7 と 8）は 2 重になっている．
拘束条件がないので，3 つの剛体モードと，上方向と横方向の対称性のために，複数組の重根をもつ．
53 M
46
/*** file output M matrix ***/
for (j=0; j<n; j++) fprintf(fp_mat,"%25.17le",mvec[j]);
fprintf(fp_mat,"\n");
for(j = 0; j < n ; j++){
scal = 1.0 / sqrt(mvec[j]);
for(i=0; i<n; i++){
a[j*n+i] = scal * a[j*n+i]; a[i*n+j] = scal * a[i*n+j];
} }
/*** file output A matrix ***/
for(i=0; i<n; i++){
for (j=0; j<n; j++){
if(a[n*j+i] != 0){;
fprintf(fp_mat," %d",i);
fprintf(fp_mat," %d",j);
fprintf(fp_mat,"%25.17le\n",a[n*j+i]);
} } } }
ファイル名は，例題集の中では S13Free.txt とした．
4.2.2
浮動小数点行列のファイルからの入力
対応する C++ ルーチンでの読込みは，まず自由度 n を読み込む．
sehfv3int.cpp の CT2D ファイルの n の読込み
strcpy(matfile,matname); // required #include <string.h>
strcat(matfile, ".txt");
cout << "Reading CT2D output matfile =" << matfile << endl;
fin.open(matfile);
if(!fin){
cout << "cannot open input file\n";
exit(EXIT_FAILURE);
}
fin >> n;
std::cout << "sehfv3 n=" << n << std::endl;
n × n の 2 次元配列を定義する．
sehfv3int.cpp の読込み領域の定義
matrix<double> a(n,n);
1 次元配列 ka，ma，aa に C の 1 次元配列の倍精度データを読み込む．非零要素のみを読み込むので，配列を零
クリアして，行と列の番号に従って行列要素を読み込む．
/*************** read matrix from file made by CT2D program *************/
matrix<double> ak(n,n);
rational_matrix kar(n+1,n+1);
ka = new double[n*n];
ma = new double[n];
aa = new double[n*n];
for (j=0; j<n; j++){
for (i=0; i<n; i++){ aa[j*n+i] = 0; ka[j*n+i] = 0; }
ma[j] = 0;
}
std::cout << " Now read file=" << matfile << std::endl;
for (;;){
fin >> i; fin >> j; fin >> tmp;
ka[n*j+i]=tmp;
if(i==n-1 && j==n-1) break;
}
for (j=0; j<n; j++){
fin >> tmp;
ma[j]=tmp;
47
}
for (;;){
fin >> i; fin >> j; fin >> tmp;
aa[n*j+i]=tmp;
if(i==n-1 && j==n-1) break;
}
fin.close();
std::cout << "sehfv3 read file end" << " n*n=" << n*n << std::endl;
1 次元配列から 2 次元配列に移す．
sehfv3int.cpp の行列の形式の変換
/******** format conversion from C type one dimensional array to matrix format *********/
for (j=0; j < n; ++j) {
for (i=0; i < n; ++i){ a[i][j] = aa[j*n+i]; }
}
倍精度から有理数に変換する．sehfv3int の実行時の指定で，A か K か M
sehfv3int.cpp の行列の型変換
−1
K かを選択する．
/*********** convert matrix from double to rational ****************/
cout << "matk=" << matk << " mattype=" << mattype << endl;
printf("mattype=%x\n",mattype);
printf("matk=%s\n",matk);
if(mattype == 2){
//MinvK
cout << " MinvK selected for A matrix" << endl;
for (j=0; j < n; ++j) {
for (i=0; i < n; ++i){ ak[i][j] = ka[j*n+i]; }
}
cnvmat(ak,kar,n,n);
// kar(n,n) <-- K matrix rational
for (i=0; i < n; ++i) {
rational mforinv = Rdset(&ma[i]);
for (j=0; j < n; ++j){ h[i][j] = kar[i][j] / mforinv; }
}
// h(n,n) <-- M^{-1}K matrix rational
}else if(mattype == 1){
cout << " K selected for A matrix" << endl;
for (j=0; j < n; ++j) {
for (i=0; i < n; ++i){ a[i][j] = ka[j*n+i]; }
}
cnvmat(a,h,n,n);
// h(n,n) <-- K matrix rational
}else{
cout << " SinvKSinv selected for A matrix" << endl;
cnvmat(a,h,n,n);
// h(n,n) <-- S^{-1}KS^{-1} matrix rational
}
cout << n << "x" << n << " A matrix" << endl;
有理数の行列になれば，固有値の計算は熱伝導行列やグラフ・ラプラシアン行列のように計算できるはずである．
4.3
チェックポイントとリスタート
CT2D の構造解析行列は，sehfv2.cpp の行列よりも，ヘッセンベルグ変換とフロベニウス変換に計算時間がか
かる．これはヘッセンベルグ行列の MatDgtCount から得られる平均の桁数が多い（234 に及ぶ，グラフ・ラプラ
シアン “G 9” は 7.3 “G 10” は 9.2 ）ことからも分かる．そこで変換された行列のチェックポイントをとり，リス
タート可能とした．チェックポイント・ファイル名は，入力行列のファイル名が S13free.txt の場合は，S −1 KS −1
の場合は S13freeA.chk に，K の場合は S13freeK.chk に，M −1 K の場合は S13freeMinvK.chk になる．
if(chkres==0){
rblas::elmhes(ar,n); /***
// ----- checkpoint H
std::stringstream ofn;
if(mattype==0){
ELMHES
*******************/
48
ofn << matname << "A" << ".chk";
}else if(mattype==1){
ofn << matname << "K" << ".chk";
}else{
ofn << matname << "MinvK" << ".chk";
}
ofs.open(ofn.str().c_str(),std::fstream::trunc);
if (!ofs) {
throw std::runtime_error("Write open failed for " + ofn.str());
}
oa << ar;
}else{
//
read H from checkpoint file
std::stringstream ifn;
if(mattype==0){
ifn << matname << "A" << ".chk";
}else if(mattype==1){
ifn << matname << "K" << ".chk";
}else{
ifn << matname << "MinvK" << ".chk";
}
ifs.open(ifn.str().c_str(),std::ofstream::in);
if (!ifs) {
throw std::runtime_error("Read open failed for " + ifn.str());
}
ia >> ar;
std::cout << "Read Hessenberg matrix from checkpoint file" << std::endl;
}
フロベニウス行列についての説明は，ヘッセンベルグ行列と同様なので省略する．
4.4
近接根の分離
1 つのグループに複数の近似固有値が近接して存在し，ヘッセンベルグ変換後の小行列でも近接固有値の状態
で存在する場合は，分離が困難になる．グループ化は相対誤差 ε = 10−9 で行っている（selectev）が，区間 [a, b]
に 2 根が存在すると， p(a)p(b) > 0 になる．これを分離するには接線の傾き p (a) と p (b) を使用して，導関数
p (x) の零点を含む区間を求める．近接する 2 根の間に存在する点 c を求めれば，区間 (a, b) を p(a)p(c) < 0 と
p(c)p(b) < 0 とに分解できる．3 根が存在する場合は曲率 p (a) と p (b) を判定に使用して，2 階導関数 p (x)
の零点を含む区間を求めて区間を分解して，2 根を求める問題に帰着する．
4.4.1
ループ構成の変更
sehfv2.cpp では，グループ数 ngrp1 を制御変数として使用して，グループ番号を包含フラッグ included[i]
でマスクしたループ構成を用いたが，sehfv3int.cpp ではグループに複数の根があるので，包含した根の数 nexist
でループ構成する．プログラムの開発では，追加や変更の機能を，熱伝導行列のケースと比較して確かめながら
行うために，プリプロセッサの変数 GENFTHERM を #define し，指定された場合は熱伝導行列やグラフ・ラプ
ラシアン行列，指定されていない場合は CT2D の行列を扱うように，ソースコードを 2 通りに切り替える．こ
のため，ループや if ブロックの字下げが異なるが，字下げはプログラミング言語の文法で決められているもので
はないので無視する．プリプロセッサのマクロプログラミングに不慣れで，プリプロセッサ変数を処理した後の
コードを想起しずらい人は，–E オプションを付けてコンパイルすることで，コメントやマクロプログラミング
を処理した後のコードを参照することができるので利用されたい．
49
p (a1 )p (b1 ) < 0
p(a1 )p(b1 ) < 0
p(a)p(b) > 0
p (a)p (b) < 0
a
p (a)p (b) > 0
p (a2 )p (b2 ) < 0
p (a2 )p (b2 ) < 0
p (a)p (b) < 0
a1
b1 b
λ
a2
a
b2
b
図 3: bisect2 による [a, b] から [a1 , b1 ] の絞り込み（左）と bisect3 による [a, b] から [a2 , b2 ] の絞り込み（右）
sehfv3int.cpp のループ構成，熱伝導行列の場合と CT2D の場合
#ifdef GENFTHERM
for(i=1; i<=ngrp1;i++){ // グループ数
if(included[i]){
#else
for(i=1; i<=nexist;i++){ // inclusion が成立した数
#endif
略
#ifdef GENFTHERM
}
//
if(includea[i])
#endif
}
4.4.2
近接する固有値の分離（2 根の場合）
1 つのグループに最大で 3 根が存在する場合を想定した．引数 mult[i] は，15 ページに記した「処理の概要」
に示したグループ化を行う selectev 関数で得た，各グループに属する根の数マイナス 1 である．chkinc 関数
で区間を探す（引数の am は固有値の平均の値で，selectev 関数の引数に加えた）．countinc 関数は，示唆さ
れた数の根を包含して，それらの区間を，それぞれ上限用の btemp と下限用の atemp 配列に返す．countinc 関
数は，包含が成立した数を戻り値に返す．呼び側は制御変数 nexist をインクリメントして，包含の成立した数
を得る．
sehfv3int.cpp の包含カウント
rational tempa,tempb;
int idummy, nmult, maxit = 53;
if(eigvl[i]*eigvh[i] >= 0){ //////////////////////////////////////////////////////////////////
idummy = chkinc(nsize,p,eigvl[i],eigvh[i],mult[i],am,tempa,tempb); // sehfv2 と異なる （am 追加）
nmult = countinc(nsize,p,tempa,tempb,mult[i],epsint,maxit,abtemp);
std::cout << std::endl << " Retern val countinc=" << nmult << " for candidate, eigvl=" << eigvl[i] << std::endl;
for(j=0; j<nmult; j++){
nexist ++; eigint[nexist]=abtemp[j]; included[nexist]=1;
eigvr[nexist]=(abtemp[j].lower()+abtemp[j].upper())/rational::TWO;
}
countinc 関数はまず a と b での関数値 p(a) と p(b)，傾き p (a) と p (b)，曲率 p (a) と p (b) を Horner3 関
数で求める．図 3 の左に，区間 [a, b] に 2 根が存在する場合を示した55 ．プログラムでは mult=1 で把握されて
いる場合である．
p(a)p(b) > 0 で，p (a)p (b) < 0 の場合が 2 根ある場合で，このとき p (x) の零点を含む区間を bisect2 関数
で，反復回数を十分大きくとって探す．
55 左図のグラフは
PostScript のスタック型（後置記法）の言語で “x x mul x sub” で記述している．関数は x2 − x である．右図のグ
1
1
ラフは “x x x mul mul 3 div x x mul 2 div sub 0.05 add” で記述している．関数は x3 − x2 + 0.05 である．
3
2
50
countinc 関数の導関数の零点を含み両端で関数値が異符号の区間探索
if(fa*fb > rational::ZERO && fda*fdb <= rational::ZERO){
rtnval = bisect2(p,n,maxit+100,a,fa,fda,b,fb,fdb);
if(rtnval<=0) {
cout << "countinc return From bisect2 rtnval=" << rtnval << " return rtnval" << endl;
return rtnval;
}
}
bisect2 関数は，区間両端での導関数の値が異符号の状態を保ち，関数値が異符号になった時点で収束とした．図
示したように，傾き p (x) = 0 を含み，p (a1 )p (b1 ) < 0 でかつ p(a1 )p(b1 ) < 0 を満たす区間 [a1 , b1 ] が bisect2
によって得られる．
導関数の零点を含み両端で関数値が異符号の区間が見つかったら，区間端点の値と関数値を保存してから，軽
く（数回） bisect 関数で反復して包含を確認する（1 次式の因子探しに必要な精度での精度改良は後に行う）．
この区間の外側と，初めに与えられた区間端の間に，もう 1 つの根がある．この根は [a, a1 ] か [b1 , b] にあるの
で，それを探して，2 個めの包含を達成する（符号を調べて bisect 関数で反復する）．
4.4.3
近接する固有値の分離（3 根の場合）
図 3 の右に，区間 [a, b] に 3 根が存在する場合を示した．プログラムでは mult=2 で把握されている場合であ
る．p (a)p (b) > 0 で，p (a)p (b) < 0 の場合が 3 根ある場合で，このとき p (x) の零点を bisect3 関数で，反
復回数を十分大きくとって探す．つまり，両端で曲率が異符号の状態を保ち，関数の変曲点を含み，両端で傾き
が異符号になる区間を探す．bisect3 では，両端で傾きが異なる区間を見つけているので，続けて bisect2 を
呼び出すことができ，残る 2 つの根は，区間に 2 根がある場合と同じ方法で求めることができる．
4.4.4
零を含むグループに対する包含
浮動小数点演算で得られる零近傍の固有値は，精度が 1 桁もない場合がある．また，同じプログラムでありな
がら，システム（コンパイラやアーキテクチャ）の違いによって，得られる固有値がかなり差異がある．ここで
述べる方法は，これらの差異を吸収する一般的な方法ではない．CT2D（32 ビット環境での実行）から抽出した
データを，64 ビットの Linux の環境と，32 ビットの Cygwin 環境で同様に包含が成立させる方法を，応急処置
的に開発したプログラム例にすぎない．浮動小数点演算で，本来，零になるべき解が，丸め誤差のために小さな
非零の値をとる場合は，その扱い方は大変難しくなる．
interval 型の区間は両端点の符号が異なる場合は扱えない．そこで，零を挟むグループは，負の端点から零
までと，零から正の端点までの 2 つの区間に分解する56 ．
if(eigvl[i]*eigvh[i] < 0){
double dzero=0;
for(int ijk=1; ijk<=2; ijk++){
if(ijk==1){
idummy = chkinc(nsize,p,eigvl[i],dzero,mult[i],am,tempa,tempb);
}else{
idummy = chkinc(nsize,p,dzero,eigvh[i],mult[i],am,tempa,tempb);
}
if(idummy > 0){
nmult = countinc(nsize,p,tempa,tempb,mult[i],epsint,maxit,abtemp);
for(j=0; j<nmult; j++){
nexist ++; eigint[nexist]=abtemp[j]; included[nexist]=1;
eigvr[nexist]=(abtemp[j].lower()+abtemp[j].upper())/rational::TWO;
}
}
56 4.4.2 節に示した「sehfv3int.cpp の包含カウント」プログラムが，区間両端点の符号が同じ場合で，ここに示したプログラムは，符号
が異なる場合である．
51
}
}
chkinc 関数は，固有値の平均 am と相対誤差判定値 eps から，固有値が零近傍であると判断すると，零では
ない端点を 2 倍する．2 倍すれば包含が成立するという保証はないが，前述したように，2 通りの環境で，この
応急処置で包含が成立した．正しくは，多項式の次数と最高次の係数，および定数項の符号などを利用して，零
近傍の根を包含する頑強（robuat）なアルゴリズムを考案すべきである．
chkinc 関数
int chkinc(const int n, const rational_vector& p, const double e, const double h, const int mult,
const double am, rational& a, rational& b) {
double ta,tb;
int rtnval=0;
rational fa,fb;
if(fabs(e/am) > eps){
ta=e-fabs(e)*eps; tb=h+fabs(h)*eps;
if(mult==0) tb=e+fabs(e)*eps;
a = Rdset(&ta); fa=Horner(n,p,a);
b = Rdset(&tb); fb=Horner(n,p,b);
if(fa*fb<rational::ZERO) rtnval=1;
}else{
ta=e-eps; tb=h+eps;
// consider eigenvalues close to zero
if(mult==0) tb=e+eps;
a = Rdset(&ta); fa=Horner(n,p,a);
b = Rdset(&tb); fb=Horner(n,p,b);
if(fa*fb<rational::ZERO){
rtnval=1;
}else{
if(mult>0 && h==0){
// (a,0) 区間を
ta=2*ta;
// (2a,0) 区間にする
a = Rdset(&ta); fa=Horner(n,p,a);
std::cout << "chkink interval inflation ta=" << ta << " fa=" << (double)fa << std::endl;
if(fa*fb<rational::ZERO) rtnval=1;
}else if(e ==0){
// (0,b) 区間を
tb=2*tb;
// (0,2b) 区間にする
b = Rdset(&tb); fb=Horner(n,p,b);
std::cout << "chkink interval inflation tb=" << tb << " fb=" << (double)fb << std::endl;
if(fa*fb<rational::ZERO) rtnval=1;
}
}
}
return rtnval;
}
4.4.5
精度改良の sehfv2.cpp に対する変更点
sehfv2 では，判定式 (27) が使用する近似固有値は，そのグループに属する固有値のおよその値を使用した．
sehfv2 で用いた参照固有値
log2ev=log2(eigvl[i])+1;
// その近似固有値の小数点よりも上のビット数
int reqndeig=((int)log2ev+2)*d+accdouble; // 必要桁（ビット）数
int niteration=reqndeig-mxitr[i];
略
double dzero=0;
bisect(p,nsize,dzero,niteration,eiglow[i],faeig,eigup[i],fbeig,ceig);
n 桁の精度がある数同士の積は，約 n 桁の精度がある．図 4 の左図は，被乗数も乗数も真ん中に小数点があり，
積の真ん中に小数点がくる場合を示した．被乗数と乗数がともに同じ桁数の精度がある場合，積もその精度まで
はほぼ正しく，下半分は誤差が混入している．したがって，精度が n 桁あれば（小数点の上の n ÷ 2 桁と下に
52
.
.
×
.
×
?
.
?
.
図 4: 乗算の精度
n ÷ 2 桁が正しければ），積の小数点よりも上の n 桁はぎりぎりである．整数性判定は小数点以下の数桁を見る
ので，この場合は精度の改良が必要になる．
もし一方の数値の精度が悪いと，積の精度はどうなるだろうか．図 4 の右図は，乗数の真ん中の数字が怪しい
（精度が n ÷ 2 桁しかない）場合である．積の精度は短いほうの桁数しか保証されない．倍精度浮動小数点演算
で零固有値を求めても，行列そのものに丸め誤差が混入していて，零にならず，絶対値の小さな固有値が得られ
るが，この固有値は精度がない．零固有値の小数点よりも上の桁数を参照して必要精度を決め，さらに精度のな
い零固有値に対する近似解を，この必要精度から得られた反復回数だけ反復する sehfv2 の方法では，整数性判
定ができない．そこで sehfv3int では，要求精度の参照固有値は，包含を作成のループ反復で最大固有値を変数
evlarge に保存して，この桁数（ビット数）を判定式 (27) に代入し必要精度を決定した．
sehfv3int で用いた参照固有値と bisect 用区間幅
log2ev=log2(evlarge)+1;
// 最大近似固有値の小数点よりも上のビット数
reqndeig=((int)log2ev+2)*d+accdouble; // 必要桁（ビット）数
double stmp=pow(2,log2ev*d);
// 最大近似固有値の小数点よりも上のビット数を数値に変換
double epstmp=1/stmp;
// その逆数で区間幅に変換
略
int irtn=bisect(p,nsize,epstmp,reqndeig,eigint[i],faeig,fbeig,eigvr[i]);
さらに反復回数 k を，区間幅 2−k に変換することで，2 分法以外の区間縮小アルゴリズムでも使用できるよう
に改良した．
4.4.6
その他の sehfv2.cpp に対する変更点
固有値が零近傍の小さな値をとるため，上記の変更点のほかに，selectev 関数に，相対誤差を使用する目的
で，近似固有値の中間的な値を参照した．この値は変数 am に保存して chkinc 関数でも使用した．sehfv2 では
特性多項式が 1 次式の場合は goto Loopend によって因子探しを行わなかったが，sehfv3int では，1 次式でも
対応する近似固有値があるかどうかを見極める目的で，因子探しを行う．
因子多項式の符号は，sehfv2 ではこだわらなかったが，sehfv3int では問題の種類が増え，計算結果を diff コ
マンドでチェックするために，因子多項式の最高次の項は 1 とすることで符号を統一した．
4.5
計算例
はじめに解析内容を指定する入力データの与え方について述べる．プリプロセッサの変数 GENFTHERM を指定
した場合は，入力データの与え方は sehfv2 と同じである．指定しない場合は CT2D の行列を扱う．このとき “0
S13Free 0” のように，“数字，ファイル名，数字” を指定する．最初の数字はチェックポイント・ランかリスター
ト・ランかを指定する．0 のときは初期ラン，1 のときはリスタートで，ヘッセンベルグ変換とフロベニウス変換
をチェックポイント・ファイルから読み込む．ファイル名は CT2D の行列を格納したファイル名から拡張子 “.txt”
を除いた文字列で，3 つ目の数字が解析する行列のタイプである．0，1，2 で，S −1 KS −1 ，K ，M −1 K のいず
れかを指定する．初期ランでは常にチェックポイントをとる（ファイル名は 4.3 項）．
53
図 5: 13 節点トラス構造の振動モード
13 節点の剛体モードをもつ構造モデル S13Free.txt を例に説明する．図 5 に，最低次の非零モードを示した．
赤が振動モードで，この変形形状と直角方向に，同じ振動数のモードが存在し λ4 = λ5 = 1.5542 . . . である．上
辺と下辺が伸び縮みするモードで，これを 90 度回転したモードが同じ固有振動数で存在する．このような多重
度 2 の非零固有値が 6 ペア存在する．零固有値になる並進モードと剛体回転モードは，力学的には 3 重根であ
るが，誤差の影響で並進モードは重根，回転モードは異なる固有値をとる．表 6 にこれを示した．3 つの理論的
には零になってほしい固有値は，CT2D では零としては得られない．グラフ・ラプラシアン行列の零固有値は，
特性多項式の定数項が零になるので，浮動小数点演算で得た近似固有値が非零であっても解けた．しかし，浮動
小数点演算で作成された行列の特性多項式を正確に求めても，定数項は零にならない．これは行列そのものがす
でに丸め誤差で汚染されているからである．さらにこの絶対値の小さな近似固有値の有効桁数がない（正確に計
算してみると，符号すらもあっていない）うえに，システムによって異なる値を与える．これが問題を厄介にす
る．これが，精度改良の参照固有値を最大の値（evlarge）に変更しなくてはならなくなった理由である．本章
の冒頭に目標として掲げた改良項目の 1 つである，多倍精度演算のサポートについては，この零近傍の解に対す
る浮動小数点演算の弱点が理由で，4 倍精度，8 倍精度と多重精度算術演算を有理数計算プログラミング環境が
サポートしても意味がないと判断した．
4.5.1
最適化オプションなしで CT2D を作成した場合
熱伝導行列の場合のように，各小行列に属する固有値を表 7 にまとめる．LCM は 4,503,599,627,370,496 で，
これによって係数行列を整数行列とする．H1 から得られる固有値は有理数で λ15 = 22719829205299720 であっ
た．つまり H1 から得られる特性多項式は −λ + λ15 である．
”Sub Matrix 2” の表示後に包含のループが回り，H9 に対して “FOUND 9 intervals. nsize=9” が表示され
る．1 次の因子を見つけるために近似固有値の精度改良を行うが，零固有値に対しては，24 回の反復によって
(−0.72695829, −0.72695820) の区間に狭めるが，それ以外の近似固有値には，単根には反復 1 回，2 根のグルー
プには 10 回程度の反復を行う．1 次の因子はみつからない．
2 次の因子を見つけるための精度改良は，90 回ほどの反復によって区間幅を 4 から 7 × 10−35 に狭める．零
固有値に対しては，次のように表示される．
eig=[
-0.72695826351001392823598353775519177966215988907840:
-0.72695826351001392823598353775519170679139841389762]
|
なおこの表示は formatprn を呼ぶことで得られる．
std::cout << "eig="; formatprn(eigint[i]);
54
表 6: S13 モデルの固有値
group
i
λi
下位 3 桁の差
1
2
3
1,2,3
4,5
6
0.
1.554 . . .
1.776 . . .
4
5
7
8
2.175 . . .
2.955 . . .
6
7
9,10
11
3.380 . . .
3.456 . . .
3.38 . . . 410, 3.38 . . . 446
3.38026893907544201597
3.45573401676622660333
8
9
10
12,13
14
15
4.488 . . .
4.753 . . .
5.045 . . .
4.48 . . . 162, 4.48 . . . 189
4.48755519945811977008
4.75318573732433953243
5.04481549985496435795
11
12
16
17,18
6.139 . . .
8.161 . . .
13
14
19,20
21
10.71 . . .
10.87 . . .
15
16
22
23
11.25 . . .
11.68 . . .
17
18
24
25,26
12.26 . . .
12.93 . . .
−0.00 . . . 007, −0.00 . . . 006, 0.00 . . . 003
1.55 . . . 331, 1.55 . . . 335
−1.6150E − 16, −1.5117E − 16
1.55421534213563386560
1.77640027632889171915
2.17495713431908715053
2.95518450014503430978
6.13854338429956082551
8.16130634735060237338
8.16 . . . 960, 8.16 . . . 978
10.7 . . . 541, 10.7 . . . 648
10.70846185912325619640
10.86985836553587454288
11.24681426267565644372
11.67871248279866691460
12.9 . . . 590, 12.9 . . . 767
12.26431833879140533122
12.93226981213176417701
特性多項式は 2 次式 λ2 − 58748626224263676λ + 479506099890433391927153378920472 で割り切れて，7 次
式になる．この 2 次式は表 7 の λ7 と λ21 の単根から構成される．
3 次式のための精度改良は，区間幅を 10−52 程度にする．3 次の因子は見つからない（因子探しの時間は 0.1
秒以下である）．
この後，4 次，5 次，6 次の因子探しを行うが，この探索はスキップできる（チューニングのアイテムとして残
した）．7 次式のための精度改良は，区間幅を 10−120 程度にする．7 次の因子
-1.*x^7
+
185656740064426520.*x^6
-13355185 (20 digits) 9311196.*x^5
+47021772 (35 digits) 58533088.*x^4
-84245215 (52 digits) 1060672.*x^3
+71844297 (67 digits) 97565184.*x^2
-22232711 (84 digits) 2703744.*x^1
-16162253 (84 digits) 6976128.*x^0
は 9 つの包含された区間のうちの残りの近似固有値から構成されることを確かめる57 ．この 7 つの固有値は，7
つの重根の片方である．定数項と 1 次の項の桁数が 99 桁あるので，近似固有値の区間幅も 10−120 程度の精度
になっている．H9 の特性多項式は次のように分解された．
(λ2 + a1 λ + a0 )(λ7 + b6 λ6 + · · · + b0 )
57 x7
の係数が −1 で，定数項が負なので，負の根をもつことが分かる．
55
(31)
表 7: 小行列に属する固有値と因子多項式
λi
value
LCM 倍での区間（最適化なし）
H9
H16
λ 2 , λ3
λ1
±0.00 · · ·
(−1.513466502033021, 4.344940969805278)
7
7,3
λ4
λ6
λ5
1.554 · · ·
1.776 · · ·
(6999563635695548, 6999563635695552)
8000195622535650
7
7
3
λ7
9795136139666244
1.330896781366425e+16
(1.52233779344322e+16, 1.522337793443223e+16)
2
λ10
2.175 · · ·
2.955 · · ·
3.380 · · ·
λ13
3.456 · · ·
4.488 · · ·
1.556324243019992e+16
(2.02101519240841e+16, 2.021015192408413e+16)
4.753 · · ·
5.045 · · ·
2.140644551543668e+16
2.271982920529966e+16
3
6.139 · · ·
8.161 · · ·
10.71 · · ·
2.764554169812914e+16
(3.675525622478457e+16, 3.675525622478465e+16)
(4.822662483845869e+16, 4.822662483845872e+16)
3
7
7
λ22
10.87 · · ·
11.25 · · ·
4.895349008459746e+16
5.065114852249124e+16
λ23
λ24
11.68 · · ·
12.26 · · ·
5.25962451856991e+16
5.523357950053409e+16
12.93 · · ·
(5.824176550697125e+16, 5.824176550697141e+16)
λ8
λ9
λ11
λ12
λ14
λ15
λ16
λ17
λ19
λ18
λ20
λ21
λ25
λ26
H16
H9
H1
56
7
1
7
7
2
7
7
7
2
3
2
3
7
7
a1 と a0 は λ7 と λ21 の単根から構成され，b0 から b6 は λ1 ，λ5 ，λ10 ，λ13 ，λ18 ，λ20 ，λ26 の重根から構成さ
れる．
H16 から得られる特性多項式は，式 (31) の 7 次式 (λ7 + b6 λ6 + · · · + b0 ) で割り切れて 9 次式になる58 ．な
お，H9 から得られた 7 次式で割り切れる理由は，この 7 つの固有値が重根の残った片方なので，もとの行列の
特性多項式が平方の形 (−λ7 + · · · )2 になるからである．9 次式は 1 次式 = 13308967813664242 − λ で割り切れ
る（λ8 ）．つまり λ8 は有理数である．残りの 8 次式は 2 次式と 2 つの 3 次式に分解される．つまり，H16 の特
性多項式は次のように分解された．
(λ7 + b6 λ6 + · · · + b0 )(λ − λ8 )(λ2 + c1 λ + c0 )(λ3 + d2 λ2 + d1 λ + d0 )(λ3 + e2 λ2 + e1 λ + e0 )
(32)
c1 と c0 は λ11 と λ23 の単根から構成され，d0 から d2 は λ3 ，λ14 ，λ22 の単根から構成され，e0 から e2 は
λ6 ，λ16 ，λ24 の単根から構成される．
H16 の特性多項式に関係する 2 つの零固有値は，重根のペアで H9 と H16 の 7 次式に含まれ，負根である（2
つの並進モード）．残ったひとつは λ14 と λ22 と組み合わされた 3 次の因子になる．これは剛体の回転モードで，
固有値は正である．表 7 の最後の 2 列に，因子多項式の次数を示した．H16 は 2 つの 3 次の因子を持つが，3 と
イタリック体の 3 で識別した．
整数性判定のために精度改良される固有値の区間幅は 10−120 と狭く，小数点の上に 17 桁の数字があるので
140 桁近い精度が要求される．この精度は，多重精度算術演算では 8 倍精度でも不足する．
4.5.2
A = K を選択した場合
質量行列を考慮しないので，全節点に質量 1 を想定した場合の振動と考えられる．質量が変わるので固有値の
値に変化があり，重根値と単根の順序関係が変わる．各小行列に属する固有値を表 8 にまとめる．ヘッセンベルグ
行列は H16 と H10 に分かれる．零になる 3 つの固有値は，3 つとも小さな正の値をとった（λ1 = λ2 = 0.61 · · · ，
λ3 = 1.14 · · · ）．
H10 の特性多項式は次のように分解された．
(λ − λ16 )(λ2 + a1 λ + a0 )(λ7 + b6 λ6 + · · · + b0 )
(33)
a1 と a0 は λ8 と λ23 の単根から構成され，b0 から b6 は λ1 ，λ6 ，λ11 ，λ14 ，λ19 ，λ22 ，λ26 の重根から構成さ
れる．
H16 の特性多項式は次のように分解された．
(λ − λ7 )(λ7 + b6 λ6 + · · · + b0 )(λ2 + c1 λ + c0 )(λ3 + d2 λ2 + d1 λ + d0 )(λ3 + e2 λ2 + e1 λ + e0 )
(34)
c1 と c0 は λ9 と λ17 の単根から構成され，d0 から d2 は λ3 ，λ12 ，λ24 の単根から構成され，e0 から e2 は λ4 ，
λ15 ，λ20 の単根から構成される．
4.5.3
A = M −1 K を選択した場合
S −1 KS −1 と M −1 K は，演算の丸め誤差がなければ同じ物理モデルになるので，固有値の順序は S −1 KS −1
の場合と同じになる．ヘッセンベルグ行列は H16 と H10 に分かれる．各小行列に属する固有値を表 8 の形式で
まとめるならば，表 7 の第 1，2，3 列の λ の並びのうち λ15 を第 2 列に移動する（表は省略する）．
質量行列の対角項で正確に割り算を計算するので，行列の桁数は増えて，LCM は 6.76 × 1046 となり，もと
は零固有値になる 3 つの近似固有値も，倍精度浮動小数点演算では −4.69 × 1031 から 2.10 × 1031 に分布する．
最大の近似固有値は 8.74 × 1047 である．
58 この
9 次式の定数項は負で，最高次の係数は 1 なので，負根はもたない．
57
表 8: A = K の場合の小行列に属する固有値と因子多項式
λi
λ 2 , λ3
λ4
λ5
H10
H16
λ1
7
λ6
7
7,3
3
7
λ7
λ9
λ10
λ12
λ13
λ15
1
λ8
2
λ11
7
λ14
7
λ16
1
λ17
2
7
3
7
3
2
λ18
λ20
λ19
7
7
3
λ21
λ22
λ23
7
2
7
λ26
7
λ24
λ25
3
7
H10 の特性多項式は次のように分解された．
(λ − λ15 )(λ2 + a1 λ + a0 )(λ7 + b6 λ6 + · · · + b0 )
(35)
a1 と a0 は λ7 と λ21 の単根から構成され，b0 から b6 は λ1 ，λ5 ，λ10 ，λ13 ，λ18 ，λ20 ，λ26 の重根から構成さ
れる（S −1 KS −1 の場合と同じ）．
H16 の特性多項式は次のように分解された．
(λ7 + b6 λ6 + · · · + b0 )(λ − λ8 )(λ2 + c1 λ + c0 )(λ3 + d2 λ2 + d1 λ + d0 )(λ3 + e2 λ2 + e1 λ + e0 )
(36)
c1 と c0 は λ11 と λ23 の単根から構成され，d0 から d2 は λ3 ，λ14 ，λ22 の単根から構成され，e0 から e2 は
λ6 ，λ16 ，λ24 の単根から構成される（S −1 KS −1 の場合と同じ）．
最後の 7 次の因子を 7 つの近似固有値から構成するときの，近似固有値の精度は，500 桁を超える．
Refined accuracy of eigval eigvr=8.795420292271488e+30 Width of interval=0
eig=[
8795420292271488564653283220829.969762110309793004 (468 digits) 190808880683222488:
8795420292271488564653283220829.969762110309793004 (468 digits) 190808880683262954]
eiglow_n=61 eiglow_d=58 eigup_n=61 eigup_d=58
sehfv2.cpp ですでに，整数性判定に必要な精度まで近似固有値の精度を改良すれば，因子多項式を見つけられる
ことを確かめていたが，この程度のデータでも，必要な精度は，多倍精度算術演算では 32 倍精度を超える長さ
に達している．
58
a
b
図 6: [−4, 1] に 3 根が存在する場合（2 根は重根に近い状態）
4.5.4
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = S −1 KS −1 とした場合
ヘッセンベルグ行列の副対角項はすべて非零であり，26 次の特性多項式が得られ（定数項は 379 桁），多
重度の高い固有値はなくなる．この問題では，零固有値のグループに 3 根が存在するので，包含が困難にな
る．浮動小数点演算で近似固有値を求めるので，システムソフトや計算機の命令セットの違いで，得られる固
有値に違いが生ずる．ここでは 64 ビットアーキテクチャの Linux システムの結果で説明する．3 根を含む区
間 (a0 , b0 ) = (−5.7, −0.22) に対して bisect3 で 2 階導関数（曲率）がゼロを含み両端で傾きが異なる区間
(a1 , b1 ) = (−2.96, −0.22) が見つかる．この区間に対して bisect2 で導関数（傾き）がゼロを含み両端で関数
値が異なる区間 (a2 , b2 ) = (−2.340375387850572, −2.340375387850571) が見つかる．この区間を bisect で絞
り込み，区間 (a, b) が
a=
b=
-2.34037538785057119945927902383957475866260145381082
-2.34037538785057119945927902383943958258973208643083
|
に固有値が 1 つ存在することが分かる．図 3 の右で，区間に 3 根が存在するグラフを示したが，ここで解析して
いる問題は，本来は重根であったものが，倍精度浮動小数点演算の丸め誤差のために分離した 2 根なので，図示
すると図 6 のように描ける．
2 つめの包含をプログラムのコメントの図
//
//
------|----|----|---------|-----|------|-----------a0
a1
a2
b2
b1
b0
//
//
setinit
f’’(x)=0
//
setinit
1
2
second root
1
2
third
root
に従って説明する．(a1 , b1 ) 間には 2 根存在するはずなので，次に (a1 , a2 ) か (b2 , b1 ) を探す（変数 setinit
が 1 か 2）．ここでは setinit = 2 が選択されて，bisect で，区間 (a, b) が
a= -2.34037538785057115846081188823544830251686213817924:
b= -2.34037538785057068919606547918916441250341339773654]
|
に固有値が 1 つ存在することが分かる．
3 つめの包含は，(a0 , a1 ) か (b1 , b0 ) を探す．ここでは setinit = 1 が選択されて，bisect で，区間 (a, b) が
a= -3.65411544897164871498196347943121509160846471786499
b= -2.96869132933425497467005982343835057690739631652832
に固有値が 1 つ存在することが分かる．
他のグループについては，グループに 2 つ以下の近似固有値なので，bisect2 の使用で，FOUND 26 intervals.
となって包含が成立する．
59
有理数の固有値 λ8 は 1 次式として分解される．
2 次の因子探しでは近似固有値の精度は，最小固有値に対しては，
Refined accuracy of eigval eigvr=-2.340375387850571 Width of interval=6.600394183074579e-35
eig=[
-2.34037538785057119945927902383953132806887682308034:
-2.34037538785057119945927902383953126206493499233455]
eiglow_n=6 eiglow_d=6 eigup_n=6 eigup_d=6
の精度で多項式係数の整数性判定を行う．最後の行は，区間の下限の有理数の分子が 232 進数で 6 桁，分母が 6
桁，分母の有理数の分子が 6 桁，分母が 6 桁を表示している．
3 次から 11 次まで因子は見つからない．12 次の因子探しには 200 桁を超える精度が要求される．
Refined accuracy of eigval eigvr=-2.340375387850571 Width of interval=1.366351164922065e-205
eig=[
-2.340375387850571199 (173 digits) 626977603588361562:
-2.340375387850571199 (173 digits) 626977603588347899]
eiglow_n=24 eiglow_d=23 eigup_n=24 eigup_d=24
25 次の特性多項式は 12 次式で割り切れる．
13 次の因子探しには 230 桁程度の精度が要求される．12 次の因子と値が近いほうの零固有値は次の区間に絞
り込まれた．
Refined accuracy of eigval eigvr=-2.340375387850571 Width of interval=9.480965768874258e-223
eig=[
-2.340375387850571199 (191 digits) 060678401700349592:
-2.340375387850571199 (191 digits) 060678401700254782]
eiglow_n=25 eiglow_d=25 eigup_n=25 eigup_d=25
もう１つの理論的に零になる固有値は次の区間に絞り込まれた．
Refined accuracy of eigval eigvr=-3.457446882946081 Width of interval=9.480965768874258e-223
eig=[
-3.457446882946080503 (191 digits) 951572206194082491:
-3.457446882946080503 (191 digits) 951572206193987682]
eiglow_n=25 eiglow_d=25 eigup_n=25 eigup_d=25
13 次多項式が，2 つの零固有値を含む 13 個の近似固有値から構成された．
構造の形状対称性から重根になるペアは 12 次の因子多項式と 13 次の因子多項式に分かれた．また 3 つの理
論的に零固有値になる 3 根は，最適化なしの場合と同様に，2 つの並進モードが 12 次と 13 次に分かれ，剛体の
回転モードは 13 次の因子に含まれた．これは最適化なしの場合の H16 に剛体の回転モードが回ったことと関係
があるかもしれない（理由はわからない）．
H26 の特性多項式は次のように分解された．
(λ − λ8 )(λ12 + a11 λ11 + · · · + a0 )(λ13 + b12 λ12 + · · · + b0 )
(37)
λ8 は，最適化なしでは H1 で分離していた有理数の固有値である．最適化した場合の行列では，ヘッセンベル
グ変換で H1 として分離されなくなった．
a11 から a0 は，7 つの重根のペアと λ7 ，λ11 ，λ15 ，λ21 ，λ23 から構成される．表 7 では， λ7 と λ21 は H9 の
特性多項式の 2 次の因子，λ11 と λ23 は H16 の特性多項式の 2 次の因子，λ15 は有理数で H1 に分離していた．
b12 から b0 は，7 つの重根のペアと H16 で 2 つの 3 次多項式の因子を構成した固有値から構成される．この
ように，最適化なしでの因数分解の因子多項式を構成する固有値が，最適化のための誤差に影響されても因子多
項式のペアリングを崩さないことは，もとの行列の係数の作られ方と，丸め誤差の載り方に起因していると考え
られるが，細かいことは分からない59 ．
59 固有ベクトルと合わせて考えるとヒントが得られるかもしれない．
60
4.5.5
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = K とした場合
ヘッセンベルグ行列は H16 と H10 に分かれる．H10 と H16 の特性多項式は，最適化なしの場合の式 (33) お
よび式 (34) と同じであった．
4.5.6
Cygwin 環境で最適化オプション -O3 で CT2D を作成し A = M −1 K とした場合
ヘッセンベルグ行列は H16 と H10 に分かれる．H10 と H16 の特性多項式は，最適化なしの場合の式 (35) お
よび式 (36) と同じであった．
4.5.7
多重精度算術演算のサポートに関して
浮動小数点演算のヤコビ法などで求める固有値の精度は，通常は 10 数桁あるが，零固有値の近傍では，全く
精度がない場合がある．今回の計算では，符号すら違っていた．これを 4 倍精度にしても，改善の保証はない．
したがって，精度を出す必要がある場合は，多重精度算術演算よりも有理算術演算を用いるべきである．さらに
それを区間算術演算の形で使用するのがよい．
4.5.8
デバッグ例
有理算術演算は数値の 1 ビットの違いにも敏感に結果が反応する．浮動小数点演算では，有効桁の下位のビッ
トが違っても，結果に大差ないのが普通なので，著しい特長である．CT2D の行列で最適化ありとなしでは，行
列要素の違いは浮動小数点データの下位のわずかな（数値を目で眺めても分からない）ビットにしか存在しない
が，この差に敏感に反応して，ヘッセンベルグ変換で副対角項が零にならず，計算時間は数倍になる．プログラム
のバグで変換が正しく行われない場合も同様の結果になることが多い．このような問題を解決するために，ヘッ
センベルグ変換では，軸選択，副対角項の値（倍精度変換），その桁数の情報を表示している60 ．行列 K を選択
した場合次のように表示される．
剛性行列のヘッセンベルグ変換中の表示
elmhes start n=26
m=2 i=3 (double)x=-9.0072e+15 3
m=3 i=5 (double)x=-9.0072e+15 3
m=4 i=13 (double)x=4.72215e+16 10
m=5 i=7 (double)x=-2.11837e+16 13
m=6 i=12 (double)x=8.85753e+15 30
m=7 i=26 (double)x=-1.06297e+16 53
中略
m=15 i=19 (double)x=-2.89904e+15 340
m=16 i=20 (double)x=4216.64 364
m=17 i=17 (double)x=0 1
m=18 i=25 (double)x=-2.36378e+16 30
m=19 i=23 (double)x=-2.91652e+16 59
中略
m=24 i=26 (double)x=5.33917e+15 117
m=25 i=25 (double)x=1.46586e+16 137
入力行列を間違えたり，プログラムの下位の階層のバグに触れたりすると，副対角項 k17 が零にならず，計算時
間が極端に伸びる（1 分程度の分解の時間が 8 時間近くかかる）．
60 ここでのデバッグ例は，プログラムの更新のミスでバグがある longint.cpp を使用したことで発生したので，正しい計算（正解）と比
較することで問題解決を行うことができた．このような場合での問題追及の方法を紹介する．
61
剛性行列のヘッセンベルグ変換中の表示
* elmhes start n=26
m=2 i=3 x=-9.0072e+15 3
m=3 i=5 x=-9.0072e+15 3
m=4 i=13 x=4.72215e+16 10
m=5 i=7 x=-2.11837e+16 13
m=6 i=12 x=8.85753e+15 30
m=7 i=26 x=-1.06297e+16 70
<=== ここで違っている
m=8 i=8 x=2.02643e+16 120
中略
m=24 i=25 x=-4.50491e+15 55238
m=25 i=26 x=-2.51654e+15 59828
この例では正解と比較できるので，
（計算は 7 列目で停止するようにして）更新された行列要素をすべて出力する
プログラムを 2 通り実行して，正しい結果と比較する（diff をとる）ことで，問題が最初に現れた演算を捜した．
次に示す問題解決用のコードは，ヘッセンベルグ変換の関数 elmhes にはコメントで残した．
elmhes に加えたデバッグ用プリント
for(int j=m; j<=n; ++j) {
rational saved=a[i-1][j-1];
a[i-1][j-1] = a[i-1][j-1] - y * a[m-1][j-1];
if(saved != a[i-1][j-1]) std::cout << "change a[i-1][j-1]=" << a[i-1][j-1] << " i=" << i << " j=" << j
<< " a[m-1][j-1]=" << a[m-1][j-1] << " m-1=" << m-1 << " saved=" << saved << " y=" << y << std::endl;
}
乗加算 “saved - y * a[m-1][j-1]” に着目する．次に示すペアの 2 行は，上が正しくない計算の表示，下が正しい
計算の表示である．被乗数 y も a[m-1][j-1] も 被加数 seved も同じであるが結果 a[i-1][j-1] が異なるの
であるから，乗加算 “saved - y * a[m-1][j-1]” が正しく計算されていない．
elmhes に加えたデバッグ用プリント
change a[i-1][j-1]=766654595212970732358095103346411/47221509289895611 i=6 j=6
change a[i-1][j-1]=766654595212996819993745768917245/47221509289895611 i=6 j=6
a[m-1][j-1]=258359429617980736260547150317435671248383889011/40564819207303340847894502572032 m-1=3
a[m-1][j-1]=258359429617980736260547150317435671248383889011/40564819207303340847894502572032 m-1=3
saved=15376250927266764/1
saved=15376250927266764/1
y=-40564819207303340847894502572032/300756232722000874105869297291081
y=-40564819207303340847894502572032/300756232722000874105869297291081
この部分をテキストエディタのコピー&ペースト機能で，デバッグ用の mpyad.cpp プログラムに埋め込んで確
かめた．
デバッグ用 mpyad.cpp プログラム
int main(int argc, char** argv) {
longint an("15376250927266764");
rational a=RRset(an,longint::ONE);
longint bn("258359429617980736260547150317435671248383889011");
longint bd("40564819207303340847894502572032");
rational b=RRset(bn,bd);
longint cn("40564819207303340847894502572032");
longint cd("300756232722000874105869297291081");
rational c=RRset(cn,cd);
rational d=a+b*c;
std::cout << "d=" << d << std::endl;
}
これによって下位の階層（rational クラスか longint クラス）に問題があることがはっきりした．“g++ mpyad.cpp
62
rational.o longint.o mempool.o” でコンパイルでき，稼働するので，デバッグ作業は軽くなる．
63
5
高速化チューニング
本章で HPC プログラミングのテーマである sehfv3int.cpp のチューニングを行う．グラフ・ラプラシアン行
列や CT2D の最適化オプションを指定した場合の行列のように，特性多項式の係数の桁数が増え，次数も大き
くなり，それが低次の因子多項式で分解されない問題では，近似固有値の組合せの数が膨大になる（指数時間）．
また因子多項式の次数が大きくなると，多項式係数の作成にも多大な時間を要するようになる．組合せの数が指
数関数的に増大するので，計算量も指数時間になるが，それに加えてその計算のために有理数の分母と分子の桁
数も増加するので，計算時間は急激に増大する．そこでここでは数行から数 10 行で行える簡単なチューニング
を 4 つ試してみる．前章で触れた，特性多項式の次数は，nsize/2+1 から nsize-1 までは因子はないので，因
子探しを行わないでよい．はじめにこれを加える．次に整数性判定を高速化する．さらに，根と係数の関係を使
用する係数作成の乗算回数の最適化を行う．最後に，2 分法で近似固有値の精度を改良する際の分母・分子の桁
数の削減について考える．有理数を構成する分母と分子の桁数がチューニング対象になる点が，浮動小数点演算
とは異なる．チューニングを通して有理算術演算の特徴を学ぶ．
5.1
因子探しループでの反復のスキップ
次数 d を 2 から nsize まで行う因子探しのループでは，次数の半分まで探索して，最後に d = nsize で近
似固有値から構成される多項式で割り切れることを確かめる．因数分解だけが目的なら最後の nsize も不要であ
る61 ．この場合は while(d<=nsize/2) にする．
プログラムの変更は while ループ中に if ブロックを追加した．
while ループ内の反復スキップ
d=2;
#ifdef FACTORIZEONLY
while(d <= nsize/2){
// search for polynomial of degree 2 to nsize/2
#else
while(d <= nsize){
// search for polynomial of degree 2 to nsize
#endif
#ifndef NSIZEALL
if(d <= nsize/2 || d == nsize){ 次数 d は nsize の半分までと最後 d=nsize を計算
#endif
中 略
#ifndef NSIZEALL
}
#endif
d++;
}
//
while(d <= nsize){
効果を計測するために，NSIZEALL を #define し，これが指定されていない場合は d が nsize/2+1 から nsize−1
の反復をスキップする．この変更が処理時間の短縮につながる問題は，これまで扱った行列ではグラフ・ラプ
ラシアン行列だけである．“G 10” の場合，2 番目の小行列 H23 は分解されないので，NSIZEALL を指定しない
と nsize = 12 から 22 までは処理をスキップする．このスキップで近似固有値の精度改良はスキップされるが，
d = 23 で，その次数に合わせた精度要求に改善（bisect 反復）されるので，精度の改良に要する時間は変わら
ない．しかし計測結果は 2% 程度改善される．これは d = 12 や 13 で，先行判定をパスして vietatermschk で
d 次の係数を計算する時間である．プログラムには，先行判定をパスした回数を数えるカウンタ nok （Number
of OK）を入れたので，実行結果のファイルに対して “grep nok ファイル名” でこれを見ることができる．
61 プリプロセッサの変数
FACTORIZEONLY を #define したので，効果を測定可能とした．
64
.
.
+
.
.
図 7: 整数性判定で必要な桁
5.2
整数性判定の高速化
グラフ・ラプラシアン行列のように，次数の大きな特性多項式が低次の因子多項式で分解されない問題では，
近似固有値の組合せの数が膨大になる（指数時間）．また CT2D の最適化オプションを指定してコンパイルして
実行したときの行列では，近似固有値の精度改良によって，小数点以下の桁数が数 100 桁に及ぶ．このような問
題に対しては，整数性判定の高速化に効果がある．ここでは精度が十分に改良された固有値の小数部分だけを取
り出して倍精度浮動小数点数に戻し，整数性判定を短時間で処理する．図 7 に 3 つの近似固有値を加えて 1 次の
係数を求める場合に，整数性判定に必要となる桁の位置を示した．小数点よりも上の桁は外してよく，また小数
点以下数 10 桁の小さな桁も判定に関係しない．そこで setecanint 関数に倍精度型の dfcand 配列を引数に追加
して，小数点以下を 53 ビットの精度で取り出す．これには有理数型の固有値から，固有値を floor() で小数点
以下を切り捨てた有理数を引くことで得られ，これを倍精度浮動小数点数にキャストする．
有理数の小数点以下を倍精度で取り出す
#ifdef INTYDOUBLE
rational tmp = ercan[j]-ercan[j].floor(); dfcand[j]=(double)tmp;
#endif
効果を計測するために，INTYDOUBLE を #define した．
これを次のように 1 次の係数の総和計算で使用する．厳密な有理数の区間端点を計算するには，桁数の多い有
理算術演算を用いるが，ここでは 53 ビットの加算で間に合い，この判定を通過した場合だけ多項式係数を有理
算術演算で正確に求める．
65
整数性判定の倍精度演算による高速化
#ifdef INTYDOUBLE
for(j=1; j<=d; j++) figdaf[j]=dfcand[idx[j]];
double ud=figdaf[1];
for(j=2; j<=d; j++) ud=ud+figdaf[j];
#else
for(j=1; j<=d; j++){
figvr[j]=ercan[idx[j]]; figva[j]=ercana[idx[j]]; figvb[j]=ercanb[idx[j]];
figvaf[j]=efcana[idx[j]]; figvbf[j]=efcanb[idx[j]];
}
ura=figvaf[1]; urb=figvbf[1];
for(j=2; j<=d; j++){
ura=ura+figvaf[j]; urb=urb+figvbf[j];
}
#endif
#ifdef INTYDOUBLE
rational dummy;
if(askinty(ud,eps,dummy)){;
for(j=1; j<=d; j++){
figvr[j]=ercan[idx[j]]; figva[j]=ercana[idx[j]]; figvb[j]=ercanb[idx[j]];
}
#else
if(urb.floor() == ura.ceil()){
#endif
√
1 次の係数で行う判定は不完全である．10 ページの表 1 の例で説明すると，± 3 を 3 組の固有値が持つので，
和が整数になる組合せは 9 通りある．これは 2 次式の係数作成の例であるが，高次の係数でも同様の現象は現れ
る．したがってこの枝切りで落とされなかった組合せに対して，全係数を計算するときも，係数が確定した時点
で整数性判定を行い，不要な係数計算を避ける必要がある．CT2D の最適化オプションを指定してコンパイルし
て実行したときの行列では，判定に用いる閾値が大きい（たとえば 0.01 を使う）と，先行判定を通過するもの
が非常に多くなり，性能を低下させる（“grep nok ファイル名” で調べられる）．ここでは sehfv1 や sehfv2 で
使用した epsint よりも小さな，ファイル全体にスコープが及ぶ変数 eps = 10−9 を使用する．
5.3
根と係数の関係公式
先行判定を通過した組合せは，根と係数の関係を用いた多項式係数を求める．多項式が高次になると，この計
算量が爆発する．本節でこの計算量を減らす計算順序を，3 段階に分けて求め，プログラミングする．
5.3.1
インライン展開
17 ページに例示した次数 5 の場合の公式は，p2 から p5 の係数を作成するのに乗算を，p2 を求めるのに
= 10 回，p3 を求めるのに 5 C3 · 2 = 20 回，p4 を求めるのに 5 C4 · 3 = 15 回，p5 を求めるのに 5 C5 · 4 = 4
n
回の合計 49 回行っている62 ．最少回数は項数の
n Cr 回なので 26 回である．この数を求めるプログラムを
5 C2
r=2
十進 BASIC で示す．
62
n
X
n Cr
· (r − 1) 回．
r=2
66
表 9: vietaterm の乗算回数と VIETAY の一時配列長
i
k
m
k/m
p
2
1
1
1.000
1
3
4
5
4
11
26
5
17
49
.800
.647
.531
2
3
6
6
7
57
120
129
321
.442
.374
10
20
8
9
247
502
769
1793
.321
.280
35
70
10
11
12
1013
2036
4083
4097
9217
20481
.247
.221
.199
126
252
462
13
14
8178
16369
45057
98305
.182
.167
924
1716
15
16
32752
65519
212993
458753
.154
.143
3432
6435
17
18
131054
262125
983041
2097153
.133
.125
12870
24310
19
20
21
524268
1048555
2097130
4456449
9437185
19922945
.118
.111
.105
48620
92378
184756
22
23
4194281
8388584
41943041
88080385
.100
.095
352716
705432
24
25
16777191
33554406
184549377
385875969
.091
.087
1352078
2704156
complexityvieta.BAS
INPUT n
LET m=0
LET k=0
FOR r=2 TO n
LET ncr=comb(n,r)
LET m=m+ncr*(r-1)
LET k=k+ncr
NEXT r
LET d=INT((n-1)/2)
LET p=comb(n-1,d)
PRINT k;m;p
END
次数 6 でこの回数と（最少回数）は 129（57）回，次数 7 では 321（120） 回，次数 8 では 769（247） 回，次
数 9 では 1,793（502） 回，次数 10 では 4,097（1,013） 回，. . .，次数 15 では 212,9937（32,752） 回，. . .，次
数 20 では 9,437,185（1,048,555）回，. . .，次数 25 では 385,875,969（33,554,406） 回と n が 2 桁になると急
激に増加し，最少回数との差も開く．表 9 にこれを示した63 ．次数 10 で 4 倍，22 で 10 倍になる．
63 表は
complexityvietatab.BAS で LATEX の表の書式を指定して作成した．
67
vietaterm では nxtcmb で n 個のうちから r 個を選ぶ組合せを配列に作成すると，これを間接参照に使用し
て r 個の固有値の総積を作り，総和項に足し込む．この計算順序では，配列要素の最後だけが変わっても，毎回
1 つめから (r − 1) 番目までの積を計算する．
vietaterm 関数の総積ループ
for(int j=1; j<=r; j++) s=s*e[ith[j]];
この乗算の冗長さを p3 を例に示すと，オーバーブレースした項は，その後に出てくるアンダーブレースした項
と同一であることから分かる．
z }| {
z }| {
p3 = α1 α2 α3 + α1 α2 α4 + α1 α2 α5 + α1 α3 α4 + α1 α3 α5 + α1 α4 α5 + α2 α3 α4 + α2 α3 α5 + α2 α4 α5 + α3 α4 α5
| {z }
| {z }
| {z }
したがって，この計算順序でもアンダーブレースした項の乗算は省略できる．
これには，総積の最後の項だけが変わるところを別扱いすればよい．まず，上記のループ計算の r − 1 番目ま
での結果を残す．そして，nxtcmb 関数をインライン展開して，最後の項だけが置き換わるケースでは，残され
た中間結果に r 番目の項を掛ける．
vietatermx 関数
int vietatermx(const vector<interval>& e, const int n, const int r,
interval s,u,ss(rational::ONE,rational::ONE);
int ith[r+1];
for(int i=1; i<=r;i++){ith[i]=i;}
ith[r]=ith[r]-1;
interval t=interval(rational::ZERO,rational::ZERO);
long long nn=comb(n,r);
for(long long ii=1; ii<=nn; ii++){
//
nxtcmb(n,r,ith); // inline expansion
if(ith[r]==n){
for(int i=r-1; i>0; i--){ // search position to add
if(ith[i]<=n-(r-i)-1){
ith[i]++;
for(int j=i+1; j<=r; j++) ith[j]=ith[i]+j-i;
break;
}
}
s=interval(rational::ONE,rational::ONE);
for(int j=1; j<r; j++) s = s * e[ith[j]];
ss=s;
//
save prod excepting the last term
} else {
if(ii==1){
s=interval(rational::ONE,rational::ONE);
for(int j=1; j<r; j++) { s = s * e[ith[j]]; }
ss=s; //
save prod excepting the last term
}
ith[r]++; // differnce appears last term only
}
//
multiply the last term
s = ss * e[ith[r]];
t=t+s;
}
if(r%2==1) t=t.neg();
ta=t.lower(); tb=t.upper();
return 1;
}
rational& ta, rational& tb){
効果を計測するために，VIETAX を #define した．このような簡単な変更でも，最少回数にはまだ差があるが，
68
グラフ・ラプラシアン行列で “G 11” を指定すると，1 日かかった計算が 3 時間 40 分と約 7 倍の計算速度に向上
した．
5.3.2
一時配列の利用
17 ページに例示した次数 5 の場合の公式は，p2 から p5 の係数を作成するのに必要な乗算回数は，最低では
項数の 26 回である．これは p3 の計算では p2 の途中結果を，p4 の計算では p3 の途中結果を，p5 の計算では
p4 の途中結果を使用した場合である．ここでは全体のループ構成を変更することなく簡単にこのようなプログ
ラムへの変更を行う．
次数 5 の場合について，削減できる乗算を示す．
p2 = −(α1 α2 + α1 α3 + α1 α4 +α1 α5 + α2 α3 + α2 α4 +α2 α5 + α3 α4 +α3 α5 + α4 α5 )
p3 = α1 α2 α3 + α1 α2 α4 +α1 α2 α5 + α1 α3 α4 +α1 α3 α5 + α1 α4 α5 + α2 α3 α4 +α2 α3 α5 + α2 α4 α5 + α3 α4 α5
p4 = −(α1 α2 α3 α4 +α1 α2 α3 α5 + α1 α2 α4 α5 + α1 α3 α4 α5 + α2 α3 α4 α5 )
p5 = α1 α2 α3 α4 α5
p2 のアンダーブレースした 6 項を，p3 のオーバーブレースした項に渡せば，乗算が 6 回節約される（p3 のオー
バーブレースしていない項はのうちのはじめの 2 項の積は，その前のオーバーブレースした項を引き続き使用す
る）．p3 のアンダーブレースした 4 項を，p4 のオーバーブレースした項に渡せば，乗算が 4 × 2 回節約される．
p4 のアンダーブレースした項を，p5 のオーバーブレースした項に渡せば，乗算が 3 × 1 回節約される．
各次数において，添字列を数として並べた数は，昇順に計算される．次の次数の計算に引き渡す項は，添字に
“5” がない（一般的には d がない）項の積である．この条件に合致した中間項を，引数の配列 ss に格納する
（ss は変数 s を save するの意）．次に次数が高くなったときは，これを配列 sp から参照する（sp は「前の」
変数 s の値の意で previoous を付けた）．ループ反復中に，添字配列のはじめの (r − 1) 個が一致していれば，
同じ sp の要素を参照し，変化があれば次の要素を参照する．
69
重複計算を省く vietaterm 関数
int vietaterm(const vector<interval>& e, const int n, const int r, rational& ta, rational& tb,
const rational_vector& spa, const rational_vector& spb, rational_vector& ssa, rational_vector& ssb){
interval s;
int ith[r+1], kth[r+1];
for(int i=1; i<=r;i++) ith[i]=i;
ith[r]--;
for(int i=1; i<=r;i++) kth[i]=ith[i];
interval t=interval(rational::ZERO,rational::ZERO);
long long nn=comb(n,r);
int k=1, l=1;
for(long long i=1; i<=nn; i++){
nxtcmb(n,r,ith);
if(cmpary(ith,kth,r-1) != 1) l++;
if(r == 1){
s=e[ith[1]];
}else if( r == 2){
s = e[ith[1]]*e[ith[2]];
}else{
interval sp=interval(spa[l],spb[l]);
s = sp*e[ith[r]];
}
if(r > 1 && r < n){
if(ith[r] != n){
ssa[k]=s.lower(); ssb[k]=s.upper(); k++;
}
}
t=t+s;
for(int i=1; i<=r;i++) kth[i]=ith[i];
}
if(r%2==1) { t=t.neg(); }
ta=t.lower(); tb=t.upper();
return k-1;
}
2 つの配列 ith と kth の比較を cmpary 関数（rational クラス）で行う．
vietaterm 関数は引数（シグネチャ）が違うので，sehfv3int.cpp に含めた64 ．
一時配列のサイズは
d−1 C(d−1)/2
として，次のように配列を割当て，呼び出した．
64 関数の名称とそれらの型情報の組を合わせたものをシグネチャと呼ぶ．プログラム内でシグネチャが唯一に決まれば，関数名や演算子の
記号が重複していても呼び出すべき対象を唯一に決定することが出来るので，コンパイラは呼び分けることができ，多重定義が実現される．
70
vietaterm 関数の呼び出し
d=2;
while(d <= nsize){
中略
long long nnm=comb(d-1,(d-1)/2);
int sslen=nnm;
rational_vector spa(sslen+1), spb(sslen+1), ssa(sslen+1), ssb(sslen+1); // 配列を割当て
for(int m=1; m<=sslen; m++){ spa[m]=0; spb[m]=0;}
中略
int nterms;
if((r-1)%2 ==0){
nterms = vietaterm(figv,d,r,ca,cb,spa,spb,ssa,ssb);
}else{
nterms = vietaterm(figv,d,r,ca,cb,ssa,ssb,spa,spb);
}
}
}
ssa, ssb と spa, spb の割り当ては，
（各 d について）rational vector spa(sslen+1), spb(sslen+1), ssa(sslen+1),
ssb(sslen+1); で行われ，ブロックの最後（“}” の部分）でフリーされる（デストラクタがフリーする）65 ．配
列サイズは 次数 28 では
27 C13
=20,058,300 と大くなり，これが 4 本で，配列要素は有理数であるため，メモリ
領域を圧迫する．
効果を計測するために，VIETAY を #define した．
5.3.3
再帰呼出し
計算を次に示すように，アンダーブレースに付加した数字の順序にする．α1 α2 を計算したら，下（次の次数）
に行き，
「これ」に α3 を掛けると，α1 α2 α3 は乗算 1 回でできる（「これ」は α1 α2 を指す）．同様に「これ」に
α4 を掛けると，α1 α2 α3 α4 は乗算 1 回でできる（「これ」は α1 α2 α3 を指す）．同様に「これ」に α5 を掛ける
と，α1 α2 α3 α4 α5 は乗算 1 回でできる（「これ」は α1 α2 α3 α4 を指す）．プログラムは再帰呼び出しで，添字列
1 や 12 などは引数 num に，また「これ」も引数に渡される．また，次数ごとに「使用した添字列の最後」を記
憶しておく．
p2 = −(α1 α2 + α1 α3 + α1 α4 + α1 α5 + α2 α3 + α2 α4 + α2 α5 + α3 α4 + α3 α5 + α4 α5 )
1
9
13
15
16
20
22
23
25
26
p3 = α1 α2 α3 + α1 α2 α4 + α1 α2 α5 + α1 α3 α4 + α1 α3 α5 + α1 α4 α5 + α2 α3 α4 + α2 α3 α5 + α2 α4 α5 + α3 α4 α5
2
6
8
10
12
14
17
p4 = −(α1 α2 α3 α4 + α1 α2 α3 α5 + α1 α2 α4 α5 + α1 α3 α4 α5 + α2 α3 α4 α5 )
3
p5 = α1 α2 α3 α4 α5
5
7
11
19
21
24
18
4
ここまでで，次数の階層を 2 から 5 まで下り p5 が完成した．以上の操作は「下へ」の移動になる．
65 C では，配列を動的に獲得するには malloc 関数を用い，使用が終わった時点でこれを free 関数で解放した．C が流行った時代には，
ループ反復中に動的に獲得して解放し忘れるバグ（メモリリーク）が，システムに障害を与えることもあった．C++（や C++ のサブセッ
トとして定義された Java）では動的なメモリの獲得は new 「演算子」によって行う（ガーベッジ・コレクションつき）．管理はコンパイラ
に委ねられるので，解放は，制御が獲得したメモリのスコープの外に出た時点でコンパイラが用意した解放のコードによって行われる．メ
モリリークが起こらないという意味では安全になったが，オーバーヘッドは大きくなる．本「プログラミング環境」では longint 型の変数
に対して，下位の階層（longint クラス）でメモリの獲得を管理して，コンストラクション，デストラクションをプログラミングした．ガー
ベッジコレクションに対応する部分である．
71
添字列の最後が “5” になったので，上に戻る．
「リターン」すると，α1 α2 α3 α4 を作るときに渡された引数に，
添字列 123 と，呼出されたときの「これ」である α1 α2 α3 が入手できる．次数 4 の「使用した添字列の最後」で
ある 4 に 1 を加えた 5 を用いて α5 を掛けると α1 α2 α3 α5 は 1 回の乗算でできる．この操作は「横方向」の移
動になる．このように，次数の階層を下と上（リターン）と横に移動しながら計算すると，項数と乗算回数が一
致する最適な順序が得られる．
アルゴリズムは再帰呼び出し（recursive call）を使う．
• 下への呼出しは次数 r を 1 増やして，添字の最後が 5 が出てくるまで，呼ばれたときの添字列の最後の数
に 1 を加えた添字を 1 桁加えて呼び出す．
• 添字の最後が 5 だったらリターンする（上への移動）．
• リターンしてきて，その次数の最後に使用した添字に 1 を加えた添字を 1 桁加えて次数 r のまま呼び出す．
これによって前項の方法の一時配列を不要にできる．
十進 BASIC で再帰サブルーチン mulvieta を書いてみよう．下の（次数を増やす）呼出し（dir= 1）か横の
呼出し（dir= 0）かは引数 dir で制御する．num に添字列を 2 桁でコード化した数字が渡される（α1 α2 なら
“12”）．sp に α1 α2 などの値が渡され，全 αk , k = 1, . . . , n は配列 e() に渡される．n = 5 の場合を説明する．
CALL mulvieta(1,e(1),5,2,e,p,last,1) で呼ぶ．呼ばれた mulvieta は，num は 1，sp は α1 ，r= 2，
dir= 1 である．dgt= 1，dgt1= 2 より，引数の sp= α1 に α2 を掛けて p2 に足し込む．最後の桁が 5 より小
さいので，nump を 1 桁増やして，次数 3 の計算のために CALL mulvieta(12,e(1)*e(2),5,3,e,p,last,1) で
呼び出し，階層を下へ行く（last(r) に dgt1 を保存する）．このようにして α1 α2 α3 α4 を渡され，それに α5
を掛けて p5 に足しこんだら，最後の桁が 5 になっているので，サブルーチンはリターンする．
mulvietata 関数
EXTERNAL SUB mulvieta(num,sp,n,r,e(),p(),last(),dir)
OPTION ARITHMETIC RATIONAL
LET dgt=num-INT(num/10)*10
IF dir = 1 THEN
LET dgt1=dgt+1
ELSE
LET dgt1=last(r)+1
END IF
LET p(r)=p(r)+sp*e(dgt1)
LET last(r)=dgt1
! 添字列の最後を記憶
PRINT "mul r=";r;" num=";num;" dgt1=";dgt1;p(r)
IF dgt1 < n THEN
LET nump=num*10+dgt1
CALL mulvieta(nump,sp*e(dgt1),n,r+1,e,p,last,1)
IF dgt1 <n THEN
CALL mulvieta(num,sp,n,r,e,p,last,0)
END IF
END IF
END SUB
再帰呼び出しではプロシージャフレーム（アクティベーション・レコードとも呼ばれる）はスタック領域に積ま
れてゆくので，最後の呼びだしから 1 回戻り，mulvieta(123,e(1)*e(2)*e(3),5,4,e,p,last,1) の呼び出し
パラメータで，IF ステートメントに制御が来る66 ．
66 再帰呼び出しでプロシージャフレームに引数情報が置かれ，プロシージャフレームがスタックに積まれてゆく構造について，授業などで
習っていない人は教科書を読むことをお勧めする [5]．スタックにプロシージャフレームを置くことで，メモリ用が動的に管理しやすくなっ
「ステップ実行」で
た．また，再帰呼び出しとスタックに積まれた引数の状態に不慣れな人は，付録の HanoiStackMR.BAS プログラムを，
実行してみることをお勧めする．
72
そこで dgt1= 5 になっているので，CALL mulvieta(123,e(1)*e(2)*e(3),4,3,e,p,last,0) で呼び出し，
階層を横へ行く．最後の引数 dir= 0 で呼び出された場合は，dgt1=last(r)+1 で作成し， p4 に α1 α2 α3 に α5
を掛けて p4 に足し込む．
この再帰サブルーチンを CALL mulvieta(1,e(1),5,2,e,p,last,1) で呼ぶと，上記の α1 を含む 15 項の乗算
を行い，係数に加えることができる．したがって，CALL mulvieta(2,e(2),5,2,e,p,last,1) で呼ぶと，上記
の α2 を含む 7 項の乗算を行い，係数に加えることができる．同様に CALL mulvieta(3,e(3),5,2,e,p,last,1)
と CALL mulvieta(4,e(4),5,2,e,p,last,1) で上記の全 26 項の乗算を行い，係数に加えることができる．
multvieta.BAS
OPTION ARITHMETIC RATIONAL
LET n=5
DIM e(n),p(n),last(n)
LET e(1)=583767772594036.5
LET e(2)=4973981023983810
LET e(3)=1.236054583921164e+16
LET e(4)=2.039828008041568e+16
LET e(5)=2.653525991792998e+16
MAT p=ZER
FOR i=1 TO n-1
MAT last=ZER
CALL mulvieta(i,e(i),n,2,e,p,last,1)
NEXT i
FOR i=2 TO n
IF MOD(i,2)=1 THEN
PRINT -p(i)
ELSE
PRINT p(i)
END IF
NEXT i
END
正解はプログラムのコメントに記入してある．正解には一致しない．固有値の精度が倍精度計算のものをとって
きたので，精度が不足しているからである．
C++ でこの関数を書く．十進 BASIC で添字列を用いたが，C++では配列にする（第 1 引数）．また数値は
区間演算を使用するので，それぞれ a と b の文字を付加して上限と下限を抑える．
73
mulvieta
int mulvieta(vector<int>& num,const interval& sp,const int n,const int r,
const vector<interval>& e, vector<interval>& p,vector<int>& last,int dir) {
int dgt, dgt1;
interval t;
int rtnval=1;
dgt=num[r];
if(dir == 1){
dgt1=dgt+1;
}else{
dgt1=last[r]+1;
}
t=sp*e[dgt1];
p[r]=p[r]+t;
last[r]=dgt1;
if(r == n){
// check if the final term 123...n case
if(p[r].upper().floor() != p[r].lower().ceil()){
//interval
rtnval=0;
}
}
if(dgt1 < n){
num[r+1]=dgt1;
int inty=mulvieta(num,t,n,r+1,e,p,last,1);
if(inty == 1){
if(dgt1 <n){
mulvieta(num,sp,n,r,e,p,last,0);
}
}else{
rtnval=0;
}
}
return rtnval;
}
この関数は，全係数を作成するので，整数性判定も含める．整数性を pn が完成したときに判定して，その先の
計算を行わないようにした．
これを呼び出す関数 vietatermschk は，全係数を作成するので，インターフェースは係数ごとに計算する
vietaterm とは異なる．
74
vietatermschk
int vietatermschk(const vector<interval>& e, const int n, rational_vector& cor) {
int i,r,inty,ok=1;
vector<int> last(n+1);
vector<int> num(n+2);
vector<interval> p(n+1);
for(i=0; i<n+1; i++){ p[i]=interval(rational::ZERO,rational::ZERO);}
for(i=1; i<n; i++){
long long tstart0=gettime();
num[2]=i;
for(int j=1; j<n; j++) last[j]=0;
inty=mulvieta(num,e[i],n,2,e,p,last,1);
if(inty==0){ ok=0; break;}
std::cout << "loop in vietatermschk i=" << i << " time=" << (double)(gettime()-tstart0)/1000000.0
}
if(ok==1){
interval t=e[1];
for(r=2; r<=n; r++){
interval er=(e[r]);
t=t+er;
}
if((t.upper()).floor() == (t.lower()).ceil()){
cor[1]=-(t.upper().floor());
}else{
ok=0;
}
}
if(ok==1){
for(r=2; r<=n; r++){
if(p[r].upper().floor() == p[r].lower().ceil()){
if(r%2 == 1){
cor[r]=-p[r].upper().floor();
}else{
cor[r]=p[r].upper().floor();
}
}else{
ok=0; break;
}
}
}
return ok;
}
C や C++ではローカル変数や引数はスタックに積むので，再帰呼び出しが稼働する．有理数計算プログラミン
グ環境では，rational 型の変数は，桁数の増加とともに動的に配列を拡張するが，この配列そのものは 2 重にポ
イントされ，スタックに積まれるのは最初のポインタだけで，配列の本体はヒープに獲得されるので，再帰呼び
出し可能である．
「有理数計算プログラミング環境」を g++ 以外の C++ コンパイラで作成しても，再帰呼び
出しは稼働するはずである．また，多項式係数を格納する配列 for(i=0; i<n+1; i++){ pa[i]=0; pb[i]=0;}
が不要の理由は，ローカル変数として rational vector pa(n+1),pb(n+1); と定義したとき，コンストラクタ
が零で初期化しているからである67 ．
67 同様に
for(int j=1; j<n; j++) last[j]=0; も必要ない．
75
表 10: 性能の比較
s9
s 10
s 11
G9
G 10
INTY
INTY+X
6.6
6.5
6.5
33.3
28.6
28.8
7.7
7.5
7.6
144.9
20.2
14.4
INTY+Y
INTY+Z
6.5
6.6
28.7
28.8
7.5
7.6
Fonly
6.5
?8.7
?.5
なし
G 11
T13free
T13freeO3
1152
886
574
0.66
0.66
0.64
1257.7
51.0
50.6
7.74
8.25
210
227
0.??
0.64
49.5
50.0
3.06
8.6
0.??
49.5
15372
rational.cpp 内のコンストラクタ
rational::rational() {
this->s=0;
this->N= longint::ZERO;
this->D = longint::ONE;
}
ローカル変数は，プログラムの制御がローカル変数のスコープの外に出るとき，それらの変数のインスタンスを
取り除くデストラクタが起動して，取り除かれる．
効果を計測するために，VIETAZ を #define した．
5.3.4
計測結果
ここまでの改訂に対する性能測定の結果を表 10 に示す．計算時間が大きく変化する理由を，浮動小数点演算
との違いに着目して述べる．
• 熱伝導行列は次数の低い因子多項式の積に因数分解されるので，計算時間が短く，チューニングの効果は
現れない．
• グラフ・ラプラシアン行列
– m = 9, 10, 11 に対する計算時間が 8，200，15000 秒と急激に増加するのは，指数時間アルゴリズム
を，桁数増加が追い打ちするからである．
– “G 10” では H52 , H23 , H18 と 3 つの H1 に分解され，H18 は 2 つの 9 次式に分解される．H23 は分
解されない．したがって H23 の因子探しに大部分の時間を費やす．FACTORIZEONLY を指定すると 23
次の特性多項式を 11 次までしか因子探しをしないので，8 秒で計算が終了する．
• CT2D 行列は，小数点以下の桁数が多くなる．このため近似固有値の精度改良の 2 分法反復計算に多くの
時間を費やしている．多項式係数を求める vietatermschk にはあまり時間を要していない．このため，効
果の大きなチューニングアイテムは，整数性判定の INTYDOUBLE である．なお，先行判定の閾値を大きく
すると，これを通過する組合わせが多くなり，VIETAY が最速になる．これは vietatermschk 関数は d 次
の係数を計算するまで，整数性判定ができないからである．
“G 10” の問題では n = 96 の行列が H52 ，H23 ，H18 と 3 つの H1 に分割され，H18 の特性多項式が 2 つの 9
次の因子に分解される．H23 は分解されない．H52 の特性多項式は既知の因子のうち 1 次式，2 つの 9 次式，23
次式で割り切れて 10 次式になるが，1 次式（零固有値）と 4 次，5 次式に分解される．したがって計算時間は
23 次式まで次数を 1 次ずつ上げて行う整数性判定に大部分の時間を費やす．このため， INTYDOUBLE の効果
76
表 11: グラフ・ラプラシアン行列（11 分割）
INTYDOUBLE あり
VIETAX
VIETAY
VIETAZ
13,052(0.14)
?(?)
11,381(0.12)
は大きい．vietaterm の計算時間は 23 次式の係数を求めるところに約 200 秒（VIETAY の場合）時間がかかり，
大部分の計算時間をここで使用している（ヘッセンベルグ変換に約 4 秒，フロベニウス変換に約 3.7 秒）．この
データには，1 次の係数を通過して，他の係数ので整数性判定で除外されるケースが皆無に近く少ない特徴ある．
“G 11” を指定した場合の計算時間を測定した．チューニングなしの状態では 94,379 秒かかった計算が，チュー
ニングの変数の指定ごと（比率を括弧内）に表 11 に示した．
この問題では n = 117 の行列が H63 ，H42 ，H9 と 3 つの H1 に分割され，H9 の特性多項式は分解されない．
H42 の特性多項式は既知の因子で割れず，14 次，と 28 次式に分解される．したがって計算時間は 14 次式まで
次数を 1 次ずつ上げて行う整数性判定と 28 次式の係数を作成するのに多くの時間を費やす．H63 の特性多項式
は既知の 1 次，9 次，14 次，28 次の因子で割り切れて 11 次式になり，2 つの 1 次式（ひとつは零固有値）と 4
次と 5 次式に分解される．10 分割から 11 分割に上げると，因子多項式の最大の次数が 23 次から 28 次に上が
る．この計算は指数時間アルゴリズム （exponential time algorithm）であるため，因子多項式の次数がこれだ
け増えると爆発的に計算量が増える68 ．最速のオプション指定で 220 秒が 11381 秒（3 分 40 秒が 3 時間 9 分 41
秒）と 51 倍に増えている．また，一時配列を使用する VIETAY はジョブが終わらず kill したが，この理由は
メモリの不足にあると考えられる．使用する一時配列のサイズも O(n!) なので，再帰呼出しの順序で計算するこ
とは意味がある．
Cygwin 環境で最適化オプション “-O3” を指定して CT2D をコンパイルして実行した行列では，表 12 の計算
時間を要した（リスタート・ジョブで計測しており，ヘッセン変換に約 25 秒，フロベニウス変換に約 78 秒ほど
かかる）．この例では，チューニングなしの状態では 1,168 秒かかった．INTYDOUBLE の効果は大きい．12 次
と 13 次の vietatermschk の計算時間を表の下段に加えた．このデータには，1 次の係数を通過して，他の係数
ので整数性判定で除外されるケースが多い．この時間は全体の中では短く，大部分の計算時間が整数性判定と，
1 次項をパスして 2 次項（VIETAX と VIETAY）または d 次項（VIETAZ）までを計算する時間で占められて
いる．計算量にほぼ比例していると考えられる（桁数の違いがあるので乗算回数には比例しない）．VIETAZ で
最高次数の係数を求めてから除外するよりも，VIETAY で 2 次の係数を求めるほうが速いことが分る．
5.4
2 分法の桁数削減
a+b
an bd + ad bn
を用いるので，有理数として記述すると，c =
となっ
2
2ad bd
て桁数が増加していく．厳密な中点である必要はないので，ある範囲の中で，分母・分子を少し変更して桁数の
少ない有理数を選ぶ方法を考えてみる．
2 分法では反復ごとに区間中点 c =
表 12: S13freeO3 行列
VIETAX
VIETAY
VIETAZ
INTYDOUBLE あり
93(0.080)
69(0.059)
165(0.14)
vietatermschk
0.59 + 1.38
0.32 + 0.74
0.33 + 0.75
INTYDOUBLE なし
68 O(2n )
や O(n!) や O(nn ) は指数時間，O(n3 ) などを多項式時間（polynomial time algorithm）という．
77
5.4.1
最良近似分数
2 分法 bisect 関数の中で，中点を求めたら，収束判定値 tol を参照して，この範囲で中点の近傍で中点を最
良近似分数に置換えると，有理数の桁数を少なくすることができる．
#ifdef BISECTSHIFT
rational d = (a + b) / rational::TWO;
c = bestappfrac(d,tol);
#else
c = (a + b) / rational::TWO;
#endif
BISECTSHIFT を rational.cpp に #define したので，これによって中点をその近傍の最良近似分数に置換える69 ．
アルゴリズムは，いったん連分数展開を要求精度を満たすまで，配列 a に保存し，その連分数を（後半のループ
で）有理数に戻す．bestappfrac を示す．
rational.cpp の bestappfrac
rational bestappfrac(const rational& c, const double tol){
rational_vector a(1002);
rational p,q,x,y,z,t;
int i, n, m;
n=-(int)log10(tol)+10;
x=c;
p=x.ceil(); a[0]=p;
for(i=1;i<=n; i++){
y=x-p;
if(y == rational::ZERO) break;
x=rational::ONE/y; p=x.ceil(); a[i]=p;
}
m=i-1; p=rational::ONE; q=a[m];
for(i=m-1; i>=1; i--){
t=p+q*a[i]; p=q; q=t;
}
z=a[0]+p/q;
return z;
}
この置換えによって，ほぼ同じ区間幅の近似固有値が，桁数が半分程度で得られる．“G 10” の場合の 23 次式
の係数を作る vietatermschk の時間に着目する．e1 = 0.08 · · · から e51 = 7.45 · · · を精度改良して区間幅が
5 × 10−28 から 8 × 10−28 程度の区間を得る．e51 の区間両端点の有理数は，2 分法では分子 36 桁，分母 35 桁で
あるが，最良近似分数では分子 19 桁，分母 19 桁と半分近くになっている．しかし，計算時間は大幅に伸びる．
2 分法は（プロファイル取得の時間を含めて） 234 秒に対し，最良近似分数は 3038 秒と 13 倍にも増加する．
5.4.2
プロファイルの取得
プロファイルの取得用の関数 eprof は longint.h に用意されている．コンパイル時（Makefile）に EPROF2 を
指定することで longint クラスの基本的な関数で使用した時間が得られる．
測定の結果は次表のように得られた．
69 最良近似分数は，rational.cpp ではコメントアウトしてあり，BISECTSHIFT 変数は後述する「中点の有理数の分子の移動による桁数削
減」用にセットされている．
78
中点
最良近似分数
lgcd2
27
2792
ldivide
lmul
121
28
143
56
ladd
3
3
全体
244
3056
最良近似分数を用いると，GCD を求める時間が 100 倍にも増加している．
5.4.3
素因数分解
vietatermschk 関数に与える近似固有値の上限と下限の有理数の分母と分子を素因数分解してみる．UNIX
のコマンドには古くから factor があるが，近年このコマンドの実装が改良され，大きな自然数もかなり高速に素
an bn
,
を次の 6 行によっ
因数分解してくれるようになったのでこれを利用する．固有値を挟む区間 (a, b) =
ad bd
て標準出力にコマンドとして表示する．これを FACTORCOMP を#define された場合に入れた．
sehfv3int.cpp の素因数分解用 factor コマンドの表示出力
for(j=1; j<=d; j++){
figvr[j]=ercan[idx[j]]; figva[j]=ercana[idx[j]]; figvb[j]=ercanb[idx[j]];
longint ad=figva[j].denominator(); longint an=figva[j].numerator();
longint bd=figvb[j].denominator(); longint bn=figvb[j].numerator();
cout << "system’factor " << ad << "’;" << endl;
cout << "system’factor " << an << "’;" << endl;
cout << "system’factor " << bd << "’;" << endl;
cout << "system’factor " << bn << "’;" << endl;
}
“G 8” で実行し，Hb の 14 次多項式に対する因子探しの部分で，14 個の包含の最初の区間の下限の分母と分子，
および上限の分母と分子について次のような出力を得る（最良近似分数の場合）．
G 8 で最良近似分数の場合の出力
%comb(14,14)=1
Vieta’s formula check for d=14 ngrp1=33 ncandi=14
system’factor 20421520906721290’;
system’factor 2552853291422741’;
system’factor 286526009741918’;
system’factor 35818040702695’;
標準出力を hoge に tee してあれば，“grep system hoge > hoge.pl” としてから “ perl hoge.pl” とするか，system
コマンド部分をテキストエディタで取り出し，hoge.pl というファイル名で保存して “perl hoge.pl” をコマンド
行実行すると次の出力を得る．
G 8 で最良近似分数の場合の出力
20421520906721290: 2 5 41 349 743 192083467
2552853291422741: 293 7001 1244509337
286526009741918: 2 41 53 251 262664033
35818040702695: 5 2713 28493 92671
下限の有理数の 17 桁の分母は 9 桁の素数をもち，16 桁の分子は 10 桁の素数をもつ．上限は桁数は少ないが，や
はりかなり大きな素数をもつ．
一方 2 分法の場合の素因数分解は次のようになる．
79
G 8 で 2 分法の場合の出力
18889465931478580854784: 2 2 2 2 2
略
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2361334177638141317547: 3 13 21412799 2827609328227
略
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
9007199254740992: 2 2 2 2 2 2 2 2
1125971878832885: 5 7 7 160981 28548733
略
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
18889465931478580854784: 2 2 2 2 2 2
a+b
a + 3b
分母は 2 のべき乗である．これは区間 (a, b) の中点は
であり，中点と上限の中点は
であり，その
2
4
a + 7b
点と上限の中点は
であり，· · · ，と続く．中点は区間の 1 : 1 の内分点，次の点は 3 : 1 の内分点，次の
8
点は 7 : 1 の内分点である．結局，n 回の 2 分法反復で選択された点は
ma + (2n − m)b
,
2n
で表され，a =
1 ≤ m ≤ 2n − 1
(38)
an
bn
とb=
を代入すると
ad
bd
man bd + (2n − m)bn ad
2n ad bd
(39)
が得られる．近似固有値は倍精度浮動小数点数を有理数変換したので，ad も bd も 2 のべきであるから，2 分法
で精度改良を続けても，分母は 2 のべきのままである70 ．
2 進 gcd アルゴリズムは，分母が 2 のべき乗である有理数の四則演算の結果を既約にするために実行されると，
除算をシフト演算で行うので高速である．結局，倍精度浮動小数点数を有理数変換した数を，区間両端にもつ区
間演算で，2 分法は 2 進 gcd と相性がよい．
“G 10” を，vietatermschk 関数だけを gcd2 関数を gcd 関数に置き換えると，Hb に対する多項式係数を決
める時間が，212.9 秒から 33,531 秒と 157 倍にも長くなった．このように極端な計算時間の違いが現れるのは，
計算量のオーダーが変わるからである．
5.4.4
中点の有理数の分子の移動による桁数削減
分母の 2 べきを保ったまま，分子の下位のビットを零にすることで，桁数を削減してみる．中点に「分母 2 冪」
の中点を選択しさえすれば，区間両端は「分母 2 冪」を守れる．さらにこの有理数を，加減算と乗算だけによる
根と係数の関係公式で使用するかぎり，計算される有理数は「分母 2 冪」を保つ．
中点 c の分子は奇数である．この最後のビットを零にすると，中点が「分母分の 1」だけ移動する．厳密に中
点で 2 分法を実施する必要はないので，この移動で桁数が削減されて高速化されるなら，導入してみたい．中点
の座標値の分母，分子が 100 桁を超える大きな数のとき，1 ビット移動するよりも，数 10 ビット移動したい．そ
こで分子の桁数の 4 分の 1 を取り去るようにした．
70 倍精度浮動小数点数
a を，e を指数，di を 0 または 1 として 2 進数で次のように表す．
!
52
X
a= ±
di 2−i · 2e = A · 2e
(40)
i=0
B = A · 252 は 254 よりも小さい整数なので，10 進数では 17 桁以下で表すことができる．B を用いると式 (40) の数は，a = B · 2e−52
である．C = 52 − e とおくと a = B ÷ 2C である．B の最下位ビットが 1 であれば，分子は B ，分母は 2C と有理数表現される．B の
下位 k ビットが 0 であれば，分子は B · 2−k ，分母は 2C−k と有理数表現される．
80
rational.cpp の bisect に加えた分子のビットクリアー
#ifdef BISECTSHIFT
rational d = (a + b) / rational::TWO;
longint dd=d.denominator();
longint dn=d.numerator();
int nbts=numdgts(dd)*32;
int nzero=nbts*0.25;
int msk=1 << (nzero-1);
longint en=(dn >> nzero) << nzero; // nzero bits clear
luint32_t off = (luint32_t)(dn-en); // offed number
if(off >= msk){
luint32_t ofc = (1 << nzero) - off; // compliment of off
longint fn(ofc);
en=dn+fn;
}
rational e = RRset(en,dd);
if(e > a && e < b) d=e*rational::ONE; // prevend loop
c = d;
#else
c = (a + b) / rational::TWO;
#endif
CT2C を –O3 でコンパイルして実行したときの行列に対しては，リスタートで計測した場合， 173.7 秒が 163.4
秒と 10 秒程度短縮された．桁数の削減は，桁数の多い近似固有値で
<
>
eiglow_n=28 eiglow_d=27 eigup_n=28 eigup_d=27
eiglow_n=26 eiglow_d=24 eigup_n=26 eigup_d=24
であった．このデータでは，先行判定を通過する組合せが多く，vietatermschk で総積まで計算して，整数性の
判定を落とされるものが多い．25 次多項式の係数を
25 C11
= 4, 457, 400 通りの組合せに対して，先行判定を通
過するものが 89,637 通りあり，これに 49.7 秒かかっていたところが 47.2 秒に短縮されている．そのほか，精
度改良の反復計算もすこしずつ短縮され，全体で 10 秒の短縮になっている．
81
付
録
82
A
ヤコビ法
ヤコビ法（1846 年）は直交行列として平面回転を用いる．n 次元空間（座標軸を x1 , x2 , . . . , xn とする）内の
特定の平面（xi xj 平面）を選び，その面内で θ だけ回転を与える行列は c = cos θ および s = sin θ とすると次
のように表わすことができる．
⎛
H −1
⎜
⎜
⎜
⎜
= RT = ⎜
⎜
⎜
⎜
⎝
..
⎞
.
c
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
s
..
−s
.
c
..
(41)
.
この行列で c や s などを記入した以外の成分は，対角項は 1 で非対角項はゼロである．第 (i, i) 成分と (j, j) 成
分は c，(i, j) 成分と (j, i) 成分は s および −s になっている．
アルゴリズム
ヤコビ法は平面回転 RT AR を繰り返すことによって，A の非対角項をすこしずつ小さくしてゆ
き，最終的に対角行列に導く．そこで A(0) = A として，回転操作を A(k) = RkT A(k−1) Rk のように記述する．最
初の相似変換 A(1) = R1T AR1 を詳しく書き下してみよう．
⎛
⎜
⎜ ···
⎜
⎜
T
R1 AR1 = ⎜
⎜
⎜
⎜ ···
⎝
..
.
..
.
(1)
aii
..
.
(1)
aji
(1)
···
aij
..
.
···
(1)
ajj
..
.
..
.
⎞
⎟
··· ⎟
⎟
⎟
⎟
⎟
⎟
··· ⎟
⎠
(42)
更新を受ける行列要素は，第 i 行，第 j 行， 第 i 列，第 j 列に限られるので，更新される要素の位置が分かるよ
うに示した．なお A も A(k) も対称行列である．
(1)
aii = c2 aii + 2csaij + s2 ajj
(1)
(1)
aij = (c2 − s2 )aij − cs(aii − ajj ) = aji
(43)
(1)
ajj = c2 aii − 2csaij + s2 ajj
(1)
aim = caim + sajm ,
(1)
ajm = caim − sajm
m には i，j 以外の 1 から n が入る．元祖ヤコビ法では，Rk 行列として A(k−1) の非対角成分の絶対値が最大
(k−1)
(k−1)
を探し xi xj 平面を選ぶ．aij
をゼロにするためには，式 (43) より，
の aij
(1)
aij = (c2 − s2 )aij − cs(aii − ajj ) = 0
(44)
を満たせばよいことになる．回転角 θ は
φ = cot 2θ =
aii − ajj
c2 − s2
=
2cs
2aij
(45)
より得られる．t ≡ s/c と置くと，この φ の定義式は次の 2 次方程式の根として得られる．
t2 + 2φt − 1 = 0
83
(46)
この 2 次方程式の小さいほうの根は，絶対値が π/4 より小さい回転角に対応している．この根は次式で与えら
れる．
t=
sgnφ
|φ| + φ2 + 1
(47)
この t から次のように c と s が決まる．
c= √
1
,
+1
t2
s = tc
このようにして決めた xi xj 平面内での θ の回転を行列 A に繰り返し繰り返し掛ける A(k) = RkT A(k−1) Rk ので
ある．連立 1 次方程式のガウスの消去法とは異なって，一度ゼロになった要素が後続の変換で非ゼロに戻るので，
非対角項の項数 n(n − 1)/2 回の Rk の演算では対角化できない．したがってヤコビ法が収束することを確かめ
て，収束することを証明しておく必要がある．
Sk を A(k) の非対角項の平方和とする．最初の反復では
S0 − S1 =
a2rs −
r,s,r=s
2
2
(a(1)
rs ) = 2aij
(48)
r,s,r=s
なので，R1 による平面回転で行列 A(1) は A(0) よりも対角行列に近い．一方 RkT を左から掛ける操
作も Rk を右から掛ける操作も直交変換であるから，対角項の平方和は 2a2ij だけ増加する．このよ
うに，1 回の反復で非対角項の平方和 Sk は必ず小さくなるので，反復回数を大きくすると対角行列
に収束する（証明終）．
m 回の変換で非対角項が十分小さな値になり71 ，対角化が完了したとすると，アルゴリズムは次のように記述
される．
T T
Rm
Rm−1 · · · R3T R2T R1T A R1 R2 R3 · · · Rm−1 Rm = Λ
H −1
(49)
H
式 (49) の左から H を掛ければ AH = HΛ が得られるので H は固有ベクトルである．そこで相似変換と同時に
I を初期値として (· · · (((IR1 )R2 )R3 ) · · · Rm−1 )Rm によって固有ベクトルを求めることができる．
次数 5 のフランク行列の固有値と固有ベクトル
次数 5 のフランク行列の固有値と固有ベクトルをヤコビ法で求
めてみよう．
⎛
5
4
3
⎜
⎜
⎜
A=⎜
⎜
⎜
⎝ 2
1
4
4
3
3
3
3
2
2
2
2
1
2
1
2
1
1
1
1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
1 ⎠
1
(50)
元祖ヤコビ法は非対角要素の中から絶対値最大の要素を探したが，計算機で計算する場合は探す手間を省いて，
非対角項を 1 行 2 列目の要素から順番にゼロにしてゆく巡回ヤコビ法（cyclic Jacobi method）を用いる．最初
a11 − a22
1
= − である．この演算によって a12 = 0 となる．値が更新される
の平面回転は x1 x2 平面内で φ =
2a12
8
71 計算機では実数は浮動小数点演算という有限のビット数で近似される．したがって数学的な意味でのゼロとはすこし異なり，計算機が
与えられたビット数で表現可能な値よりもすべての非対角項の絶対値が小さくなったところで「ゼロとなった」と判定して反復を終える．
84
のは 1,2 行と 1,2 列である．
⎛
8.531
⎜
⎜ 0.000
⎜
⎜ 4.234
⎜
⎜
⎝ 2.823
1.411
0.000
0.469
4.234
0.264
2.823
0.176
0.264
0.176
0.088
3.000
2.000
1.000
2.000
2.000
1.000
しかし次の a13 をゼロにする演算で a12 は非ゼロ 0.125 になる．
⎛
10.823
0.125
0.000
3.435
⎜
⎜ 0.125
0.469
0.232
0.176
⎜
⎜ 0.000
0.232
0.708
0.415
⎜
⎜
0.176
0.415
2.000
⎝ 3.435
1.717
0.088
0.208
1.000
⎞
1.411
⎟
0.088 ⎟
⎟
1.000 ⎟
⎟
⎟
1.000 ⎠
1.000
(51)
⎞
1.717
⎟
0.088 ⎟
⎟
0.208 ⎟
⎟
⎟
1.000 ⎠
1.000
(52)
しかし大局的に見ると，a12 から a45 までをゼロにする 9 回の演算が完了した時点で行列 A(1) は次のようになる．
⎞
⎛
12.338
0.084 0.187 −0.033 −0.158
⎟
⎜
0.307 0.021 −0.037 −0.052 ⎟
⎜ 0.084
⎟
⎜
(53)
A=⎜
0.021 1.447 0.046
0.045 ⎟
⎟
⎜ 0.187
⎟
⎜
0.000 ⎠
⎝ −0.033 −0.037 0.046 0.330
−0.158
−0.052
0.045
0.000
0.578
非対角項の平方和は最初は S0 = 100 であったが，S1 = 0.153 になる．反対に対角項の平方和は 55 から 154.847
に増加する．
この操作を繰り返すと非対角項の平方和は S2 = 0.00226，S3 = 6.38 × 10−9 ，S4 = 1.15 × 10−24 と急速に小
さくなり，6 回の反復で計算機の 32 ビットの単精度浮動小数点数のゼロ認識の精度内で対角行列になる72 ．1 回
の走査（sweep）は 10 項の非対角項に対して演算するが，ゼロになった項はスキップするので，実際には 60 回
ではなく 52 回の平面回転を行った．
⎛
Λ ≈ A(6)
12.344
0.000
0.000
⎜
⎜
⎜
=⎜
⎜
⎜
⎝ 0.000
0.000
0.000
0.272
0.000
0.000
0.000
1.449
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.353
0.000
0.000
0.000
0.000
⎞
⎟
⎟
⎟
⎟
⎟
⎟
0.000 ⎠
0.583
(54)
固有値は大きい順には並ばないので，対角項を大きい順に並べて λ1 = 12.344，λ2 = 1.449，λ3 = 0.583，
λ4 = 0.353，λ5 = 0.272 を得る．
X T AX = Λ の両辺の逆（inverse）をとると (X T AX)−1 = Λ−1 であり，左辺は X T A−1 X だから，逆行列の
固有ベクトルは元の行列の固有ベクトルに一致する．実践的ではないが，固有値と固有ベクトルが分かっていれ
ば，A の逆行列を求めるのに式 (54) の逆行列数を求めて（対角項の逆数を求めて）A−1 = XΛ−1 X T と計算し
てもよい．
ヤコビ法によるフランク行列の固有値を求める symeig.BAS プログラムを BASIC プログラム例に含めた．
72 ヤコビ法は固有値に重なりがない場合には
2 次収束する．
85
B
シェルプログラム
sehfv3int.cpp では，プリプロセッサに与えるディレクティブ変数（GENFTHERM など）を指定して条件つきコ
ンパイルを行うことで，複数のチューニングアイテムから選択してプログラムを生成し，その効果を比較できる
ようにした．現状では INTYDOUBLE，VIETAX，VIETAY，VIETAZ，BISECTSHIFT，指定なし，の 6 通
りを GENFTHERM の有無で合計 12 通りの性能に関する計測が可能である．また，因子探しループでの反復を
FACTORIZEONLY，NSIZEALL なども計算時間に関係する．このような多くの選択肢や，さらに新しいチュー
ニングのアイデアをコードとして組み込んだとき，いくつかのデータに対してプログラムが正しく稼働している
かどうかを，表示出力を眼で追って確かめるのは困難である．ここでは簡単なシェルプログラムを介して実行す
ることで，この作業を補助する方法を紹介する．
B.1
検証用データの作成用シェルプログラム
熱伝導行列（T，t，s）とグラフ・ラプラシアン行列（G）を指定すると，“../veri/” ディレクトリに検証用の
ファイルを作成するシェルプログラムを示す．
genveriFTH.sh
#!/bin/bash
# generate run result of sehfv3int GENFTERM file for verification
for n in 2 3 4 5 6 7 8 9 10 11 12
do
./sehfv3int.exe $1 $n | grep ^% > ../veri/sehfv3intveri$1$n.txt
done
exit 0
sehfv3int.cpp プログラムは，ディレクティブ変数の指定にかかわらず最終的に見つかった因子多項式と，途中
のキーとなる経過に対しては同じ表示を出力する．これらの表示の 1 カラム目には % を付した．したがって
“./genveriFTH.sh T” のコマンド入力によって，sehfv3intveriT2.txt から sehfv3intveriT12.txt までの 11 個の
検証用ファイルが，カレントディレクトリと並列の位置に置かれた veri ディレクトリに作成される．
CT2D の行列データに対しては次のシェルプログラムを使用する．
genveriCT.sh
B.2
#!/bin/bash
# generate run result of sehfv3int CT2D file for error chk
#
# ./genveriCT.sh S13Free
#
for n in 0 1 2
do
# ./sehfv3int.exe 0 $1 $n| grep ^% > ../veri/sehfv3intveri0$1$n.txt
./sehfv3int.exe 1 $1 $n| grep ^% > ../veri/sehfv3intveri1$1$n.txt
done
exit 0
実行して検証するためのシェルプログラム
veri ディレクトリが存在すれば，次の実行シェルを通して実行すると，因子多項式が正しく生成されているか
どうかが，diff コマンドで作成されたファイル hoge のサイズが零で調べられる．
86
genveriFTH.sh
#!/bin/bash
# number of arguments must be two or three, otherwise error stop
if [ $# -eq 3 ]; then
./sehfv3int.exe $1 $2 $3 > sehfv3int$1$2$3.txt
grep ^% sehfv3int$1$2$3.txt > sehfv3int$1$2$3per.txt
diff sehfv3int$1$2$3per.txt ../veri/sehfv3intveri$1$2$3.txt
elif [ $# -eq 2 ]; then
./sehfv3int.exe $1 $2 > sehfv3int$1$2.txt
grep ^% sehfv3int$1$2.txt > sehfv3int$1$2per.txt
diff sehfv3int$1$2per.txt ../veri/sehfv3intveri$1$2.txt
else
echo "You must specify two or three arguments."
exit 1
fi
#cat <<__EOT__
exit 0
“./runsehfv3int.sh T 2” のコマンドで実行すると，実行結果は sehfv3intT2.txt に残るが，先頭が % の行だけを
抽出して，検証用ファイルと diff コマンドで比較する．FACTORIZEONLY を指定すると，因子多項式は同じ
になるが，小行列の最後に得られる因子多項式は，近似固有値から構成しないので，“nok=XX” のところに差異
が出る．
参考文献
[1] W. H. Press, Teukolsky S. A., Vetterling W. T., and B. P. Flannery. Numerical Recipes in Fortran, second edition.
Cambridge University Press, 1992.
[2] J. H. Wilkinson. Algebraic Eigenvalue Problem. Oxford University Press, 1965.
[3] 寒川 光, 藤野清次, 長嶋利夫, and 高橋大介. HPC プログラミング—ITText シリーズ. オーム社, 2009.
[4] D. Goldberg. What every computer scientist should know about floating–point arithmetic. ACM Computing Survey,
23(1), 1991.
[5] D. Patterson and J. Hennessy. Computer Organization & Design: The Hardware/Software Interface, fourth edition.
Morgan Kaufmann Publishers, Inc., 2011. 成田光彰訳: コンピュータの構成と設計–ハードウェアとソフトウェアのイ
ンターフェース–第 4 版,（上）（下），日経 BP 社 (2011).
87

「ラショナル・プランナーイベント」を6/18開催(PDF:515KB)

RATIONAL メーカー保証書 I. 購入契約と国内法の関係 1.Rational AG

アルプス電気株式会社 - IBM

ÿþ2 0 a RA

2015.2 春号

ソフトウェアのサポート終了: IBM Rational の一部製品

熱物理学演習問題(6)

宝印刷、シナルマス証券（インドネシア）と業務提携契約を締結