第 1 章 行列演算
ここでは, コンピュータを利用して行列の演算を行なう時の基本的なアルゴリズムとその実装方法について 解説する .
1.1 行列の基本演算
1.1.1 準備
以下では Z, Q, R, C などの環または体の元を係数とする , n × m 行列（すなわち , n 行 m 列の行列）全
体を M (K, n, m) と書く. 特に n = m の時, すなわち, n 次正方行列の時には M (K, n) と書く. ここで, K
は Z, Q, R, C のいずれかをあらわす. したがって, A ∈ M (K, n, m) は
A = (a
ij) =
a
11a
12· · · a
1ma
21a
22· · · a
2m... ... · · · ...
a
n1a
n2· · · a
nm
, a
ij∈ K と書くことができる. 以下では, 簡単のため,
A = (a
ij) = (a
1, · · · , a
m), a
i=
a
1i...
a
ni
または,
A = (a
ij) =
a
1...
a
n
, a
i= (a
i1, · · · , a
im)
と書くこともある. これら２つの書き方の間では a
iは異なるものを表していることに注意しよう. すなわ ち, 上の書き方では a
iは n 個の要素からなる列（縦）ベクトルをあらわし, 下の書き方では a
iは m 個の 要素からなる行（横）ベクトルをあらわしている.
また , A ∈ M (K, n
a, m
a), B ∈ M (K, n
b, m
b), C ∈ M (K, n
c, m
c), D ∈ M (K, n
d, m
d) に対して ,
• n
a= n
bならば, (A B) = (a
1, · · · , a
ma, b
1, · · · b
mb) ∈ M (K, n
a, m
a+ m
b),
• m
a= m
bならば A
B
=
a
1...
a
nab
1...
b
nb
∈ M (K, n
a+ n
b, m
a)
• n
a= n
b, n
c= n
d, m
a= m
c, m
b= m
dならば, A B
C D
=
(A B)
(C D)
= A
C B
D
∈ M (K, n
a+ n
b, m
a+ m
c)
と表しておくことにする. また, Z ⊂ Q ⊂ R ⊂ C に注意すると, 自然に M (Z, n, m) ⊂ M (Q, n, m) ⊂
M (R, n, m) ⊂ M (C, n, m) が成り立ち , 暗黙のうちにこの埋め込みを利用して計算を行うことがある .
1.1.2 基本演算
行列の基本演算には, スカラー倍・和（または差） ・積がある. ここでは, K, K
は Z, Q, R, C のいずれ かとし, K ⊂ K
と仮定する.
1.1.2.1 スカラー倍
行列のスカラー倍とは , 行列の各要素を定数倍することである . すなわち , A = (a
ij) ∈ M (K, n, m), λ ∈ K
に対して, λA ∈ M (K
, n, m) を λA = (λa
ij) によって定義する. 行列のスカラー倍を計算するた めには, nm 回の乗算が必要がある. 特に, n 次正方行列のスカラー倍には n
2回の乗算が必要となる.
1.1.2.2 和
行列の和とは, A = (a
ij) ∈ M (K, n, m), B = (b
ij) ∈ M (K
, n, m) に対して, A + B = B + A ∈ M (K
, n, m) を A + B = B + A = (a
ij+ b
ij) によって定義する. 行列の和を計算するためには, nm 回の 加算が必要である. 特に, n 次正方行列どうしの和には n
2回の加算が必要となる.
1.1.2.3 積
行列の積とは, A = (a
ij) ∈ M (K, n, k), B = (b
ij) ∈ M (K
, k, m) に対して, AB ∈ M (K
, n, m) を
AB = (c
ij), c
ij=
k
=1
a
ib
jによって定義する.
ここで , 行列の積を計算するための計算量を求めてみよう . AB の一つの要素 c
ijを求めるためには , k − 1 回の加算と k 回の乗算が必要である. これを nm 個の要素全てに対して計算するので, nm(k − 1) 回の加 算と nmk 回の乗算が必要となる. 特に, n 次正方行列どうしの積には n
2(n − 1) 回の加算と, n
3回の乗算 が必要となる.
1.1.3 行列を表現するデータ構造とプログラム開発
行列をプログラム内で表現するためのデータ構造として , 最も標準的な方法は , それを二重配列として表
すことである. ここでは, いくつかの言語に対して, 二重配列の定義を見直し, 行列をどのように二重配列
として格納するかを考えてみよう. さらに, Ｃ言語に関して, より便利な行列の格納方法について考察して
みる .
1.1.3.1 二重配列の復習
多くのプログラム言語において, 二重配列は「配列の配列」として定義されている. ここでは, Ｃ, Pascal について, 二重配列の定義, メモリアロケーションの方法を復習してみよう. なお, 簡単のため, 行列要素の 型は標準的な「整数型」としておく .
1.1.3.1.1 Ｃ Ｃ言語では , int 型の要素 N 個を持つ配列は int a[N] ;
と定義され, int 型の二重配列は 特定の要素数を持つ配列 の配列として定義される. すなわち, int a[N][M] ;
は , M 個の要素を持つ int 型の配列 の N 個の配列である .
この定義は二重配列のメモリ内での配置に直接影響を与え, 「配列はその要素の並び順にメモリ内で連続 した領域に割り当てられる」というＣの定義により, a[N] という M 個の要素をもつ配列がメモリ内に連 続して割り当てられる. すなわち, int 型 N 個の要素をもつ配列 a[0] が割り当てられ, その直後に, 同じ く int 型 N 個の要素をもつ配列 a[1] が割り当てられる. 最終的に, 二重配列 a[N][M] は
a[0][0] a[0][1] · · · a[0][M-1]
a[1][0] a[1][1] · · · a[1][M-1]
· · ·
a[N-1][0] a[N-1][1] · · · a[N-1][M-1]
とメモリ内に格納される. すなわち, Ｃでは二重配列は行優先 (row-major) で割り当てられる. また, int a[N][M], b[N][M] ;
と定義している場合には , b = a ;
b[0] = a[0] ;
といった一括代入も可能である. しかし, 特定の行の一括代入は許されていない.
1.1.3.1.2 Pascal Pascal では , integer 型の要素を N 個持つ配列は A : array [1..N] of integer ;
と定義され
1, integer 型の二重配列は 特定の要素数を持つ配列 の配列として定義される. すなわち, A : array [1..N] of array [1..M] of integer ;
と定義されるが, これを
A : array [1..N, 1..M] of integer ;
と書くこともできる. この場合においても, 配列の意味は上の定義と同じである. すなわち, ここで定義し た A は M 個の要素を持つ integer 型の配列 の N 個の配列である. Pascal ではこのような配列をメモリ 内にどのように配置するかは未定義（処理系依存）とされているが , この定義から行優先であることがわか り, 文法的には,
1Pascalでは配列の添字範囲は型定義において自由に設定することが可能である. ここでは簡単のため, 1からNとおいた.
A, B : array [1..N, 1..M] of integer ; と定義したとき,
B = A ; B[1] = A[1] ;
といった一括代入が可能である. しかし, 特定の行の一括代入は許されていない.
1.1.3.1.3 その他の言語 Ｃ, Pascal では, 固定サイズの二重配列は, 行優先となっていることがわかる.
同様に Fortran でも二重配列は行優先となっているが, PL/I では列優先 (column-major), すなわち,
A[1,1] A[1,2] A[1,3]
A[2,1] A[2,2] A[2,3]
を
A[1,1] A[2,1] A[1,2] A[2,2] A[1,3] A[2,3]
と割り当てる言語も存在する.
1.1.3.2 行列の表現方法
ここでは, Ｃ言語を用いて, 行列をどのように表現するかを考えてみよう. 目標は, 行列演算（スカラー 倍・和・差・積）を行列サイズや係数環を意識せずに計算する関数を作ることであり, 実行速度の速い関数 を作ることではない
2.
1.1.3.2.1 最も簡単な実現 簡単のため, A ∈ M (Z, 2, 3) を格納することを考えてみる. Ｃでは二重配列 は行優先であることを考慮すると,
int a[2][3] ;
という変数を用意し, 要素 a[i][j] に a
i+1j+1を代入すれば良いことがわかる
3. 実際の初期化代入方法は, int a[2][3] = {{1,2,3},{2,3,4}} ;
とすれば,
A =
1 2 3
2 3 4
(1.1) と代入したことに等しくなる.
1.1.3.2.1.1 【この方法の問題点】 このような定義をしたときに, 最も単純なスカラー倍を例にして,
行列の演算の実現がどうなるかを考えてみよう. int 型の要素をもつ行列のスカラー倍を実現する関数 int_matrix_scalar_multiple は以下のように定義する必要がある.
void int_matrix_scalar_multiple(int b[][3], int a[][3], int k) なぜなら,
2実行速度の速い関数を作りたいときには,係数環は何か一つに固定した方が良い.
3Ｃの配列添字は0から始まることに注意.
• Ｃでは配列を戻り値とする関数は定義できない
4.
• 関数仮引数に a[][] という記述は出来ない.
このような定義をしてしまうと, 二重配列の第一添字に関するポインタのインクリメント量が確定し ない.
という２つの理由による.
したがって , 異なったサイズの行列に関して , スカラー倍の関数さえも異なる関数を用意しなければなら ない. 例えば,
int a[2][3], b[3][3], c[2][4] ;
と定義された３つの行列（を表していると考える二重配列）を考察してみよう. 上の int_matrix_scalar_multiple 関数でスカラー倍を計算できるものは, a, b の２つに限られる. a, b は第二添字のサイズが 3 であるた め , int_matrix_scalar_multiple 関数の定義に合致する . しかし , c は第二添字のサイズが 4 であるた め, int_matrix_scalar_multiple 関数の定義に合致しない. さらに, a, b のスカラー倍を計算するには, int_matrix_scalar_multiple 関数内部において, 第一添字の動く範囲を知る必要がある. そのため, 関 数定義のレベルでは a, b ともに int_matrix_scalar_multiple 関数を利用することが可能なように思え るが, 関数の実装のために必要な第一添字の動く範囲を指定する必要がある.
この問題は , 関数に対して第一添字のサイズを実引数として渡すことで回避することが可能である . すな わち,
void int_matrix_scalar_multiple(int b[][3], int a[][3], int k, int size)
と定義して, 仮引数 size は a の第一添字のサイズを表すことにすればよい. だが, これでも常に各 a, b の 第一添字サイズをプログラム内部で管理する必要があることに注意しよう.
1.1.3.2.2 行列の表現方法（改良その１） 上の表現方法の基本的な問題点は , たとえスカラー倍とい
う単純な関数（手続き）であっても, 行列のサイズに対して, それぞれ別の関数を用意しなければならな いことにあった. この問題点を回避する単純な方法は, プログラム内部で表現する行列のサイズに最大値
(MAX_ROW_SIZE, MAX_COL_SIZE または, 行・列ともに同じ値 MAX_SIZE) を設定し,
#define MAX_SIZE 16 int a[MAX_SIZE][MAX_SIZE]
としてしまう方法である. この場合, (1.1) の 2 × 3 行列 A は
int a[MAX_SIZE][MAX_SIZE] = {{1,2,3},{2,3,4}} ;
と定義及び初期化すれば, a の 意味のある部分 に関しては, 正しい値で初期化が行われる. この時, 上述の スカラー倍の関数 int_matrix_scalar_multiple は,
void int_matrix_scalar_multiple(int b[][MAX_SIZE], int a[][MAX_SIZE], int k, int row, int col)
と定義すればよい. ここで, 仮引数 row は意味のある行のサイズ, col は意味のある列のサイズを表す. さ らに, 実際の関数の実装としては,
4もちろん,配列へのポインタを利用して,むりやり返すことも不可能ではないが...
void int_matrix_scalar_multiple(int b[][MAX_SIZE], int a[][MAX_SIZE], int k, int row, int col)
{
int i,j ;
for(i=0;i<row;i++) for(j=0;j<col;j++)
b[i][j] = k*a[i][j] ; return ;
}
とすればよい.
1.1.3.2.2.1 【この方法の問題点】 この方法での問題点の第一はすぐに見つかる.
• 行列サイズの最大値は可変には出来ない .
もし, より大きなサイズの行列を扱いたくなった場合には, MAX_SIZE を変更して, プログラムを再コ ンパイルする必要がある.
Ｃ , Pascal, Fortran をはじめとする多くの言語では , 配列サイズを可変にとることが出来ないため
5, 二重
配列を使う限りは, 扱える行列のサイズに上限が出てくることは避けることが出来ない
6. 問題はこれだけでは終わらない.
• 配列サイズが大きいため, 実引数の値を渡すためのオーバヘッドが大きい.
関数を呼び出す際には, 二重配列の第２添字の数だけのポインタが値として渡される.
• プログラム中で各行列を表す二重配列の 有効な 範囲を管理しなければならない .
すなわち, 各行列（を表す二重配列）の行・列それぞれのサイズは, プログラム中で管理しなければ ならない.
本質的にここに関わる問題は, 次のような場面にもあらわれる.
• 和・積を計算する際に, 演算を定義できない状況を関数内でチェックするのが難しい.
二重配列そのものに行列サイズが入っていないことが原因である .
• 内部自動変数として二重配列を定義したとき, 有効でない部分の要素の初期化が出来ていない.
演算関数を正しくサイズチェックをするように定義すれば問題は回避できるが, 無効な要素の判断が 難しい.
1.1.3.2.3 行列の表現方法（改良その２） 「改良その１」での問題点である, 「行列のサイズを保持し
ていない」という問題を解消しよう. そのためには, 行列を保持するデータ構造として, 行列のサイズを付 加すればよい. すなわち, M (Z, n, m) の要素を表すデータ構造として,
struct int_matrix {
int coeff[MAX_SIZE][MAX_SIZE] ;
5PL/Iでは実行時に配列サイズを指定することが可能である.
6本来は,メモリが許す限り大きなサイズの行列を扱えることが望ましい.この方法については後述する.
unsigned int row ; unsigned int col ; } ;
という構造体を採用すればよい. ここで, 構造体メンバー row は, その行列の「行サイズ」をあらわし, col は「列サイズ」を表すことする. つまり, (1.1) の行列を定義及び初期化するには,
struct int_matrix a = {{{1,2,3},{2,3,4}}, 2, 3} ;
と定義すれば良い . この定義であれば , 構造体内部に行列のサイズを持っているので , 関数に対して行列の サイズを明示的に渡す必要が無い.
この場合, スカラー倍を求める関数は以下のようにすれば良い.
struct int_matrix int_matrix_scalar_multiple(struct int_matrix a, int k) {
int i,j ;
struct int_matrix b ;
b.row = a.row ; b.col = a.col ; for(i=0;i<a.row;i++)
for(j=0;j<a.col;j++)
b.coeff[i][j] = k*a.coeff[i][j] ; return b ;
}
この関数の呼出しは ,
struct int_matrix a = {{{1,2,3},{2,3,4}}, 2, 3} ; struct int_matrix b ;
b = int_matrix_scalar_multiple(a,4) ; とすれば良いことがわかる .
1.1.3.2.3.1 【この方法の問題点】 しかし依然として,
• 行列サイズの最大値は可変には出来ない.
• 配列サイズが大きいため , 実引数の値を渡すためのオーバヘッドが大きい .
という２つの問題はクリア出来ていない. この問題はひとまず棚上げにして, 次の問題を考えてみよう.
• 要素の型が異なる行列の演算が行えない.
例えば,
– Z 係数の行列に Q の元をスカラー倍して, Q 係数の行列を得る.
– Q 係数の行列と R 係数の行列の和を計算して, R 係数の行列を得る.
– Z 係数と Q 係数の行列の積をとり, 全ての係数が Z に入るならば Z 係数の行列を得る.
という関数を作成する必要があるかもしれない. もちろん, （たとえば「和」の関数を, 係数の組合 わせごとに書いていたら大変なことになってしまうので）これらの関数をできる限り要素の環を意識 すること無く使えるようにしたい.
次の改良では, この問題をクリアすることにしよう.
1.1.3.2.4 行列の表現方法（改良その３） 上の問題をクリアするために , 次のような構造体を定義しよ
う
7.
struct matrix { union {
int int_c[MAX_SIZE][MAX_SIZE] ; /* 整数 */
struct fractional frac_c[MAX_SIZE][MAX_SIZE] ; /* 有理数 */
double double_c[MAX_SIZE][MAX_SIZE] ; /* 実数 */
} coeff ;
int type ; /* 0: 整数, 1: 有理数, 2: 実数 */
unsigned int row ; unsigned int col ; } ;
ここで , 構造体メンバー type は , どの係数を持っているかを表すデータとする . また , 係数は共用体で保 持することとして, （わずかばかりだが）メモリを節約することとする. ここで, struct fractional は
struct fractional {
int n ; /* 分子 */
unsigned int d ; /* 分母 */
} ;
によって定義される, それぞれ, 有理数, 複素数を表す構造体であり, 有理数を表す構造体では, 常に既約分 数を表し, 符号は分子に与えることとする. この構造体に (1.1) の行列を定義及び初期化をしたいのだが, 直接は初期化は出来ない. この構造体に要素を代入するために,
struct int_matrix {
int coeff[MAX_SIZE][MAX_SIZE] ; unsigned int row ;
unsigned int col ; } ;
struct frac_matrix {
struct fractional coeff[MAX_SIZE][MAX_SIZE] ; unsigned int row ;
unsigned int col ; } ;
struct real_matrix {
double coeff[MAX_SIZE][MAX_SIZE] ; unsigned int row ;
7本来はさらにC を考慮すべきだが,余りに組合わせが多くなることや,そこまで考えるのならば,C だけではなく,Z[√
−1]や
Q(√
−1)も考慮したくなる.そうなると,もはや手に負えない...
unsigned int col ; } ;
という構造体を補助的に定義しておき,
struct int_matrix ai = {{{1,2,3},{2,3,4}}, 2, 3} ;
struct frac_matrix fb = {{{{1,2},{2,3},{3,4}},{{2,3},{3,4},{4,5}}}, 2, 3} ; struct real_matrix rc = {{{1.2,2.3,3.4},{2.3,3.4,4.5}}, 2, 3} ;
として補助的な構造体に初期化をしておく. この後,
a.row = 2 ; a.col = 3 ; a.type = 0 ; for(i=0;i<a.row;i++)
for(j=0;j<a.col;j++) a.coeff.int_c[i][j] = ai.coeff[i][j] ; と matrix 構造体に代入する.
Remark 1.1.1 ここで,
a.coeff.int_c = ai.coeff ;
この左辺は配列であり, 右辺は配列へのポインタとなっているため, このような代入は出来ない.
このように定義した matrix 構造体が表す行列の演算を行う関数を作ってみよう. 試しにスカラー倍を 計算する関数を作りたいのだが, 「スカラー」も Z, Q などという区別をしなければならない. そのために,
「スカラー」を表す構造体を定義しておこう . struct scalar {
union {
int int_s ;
struct fractional frac_s ;
double real_s ;
} s ; int type ; } ;
この時,
struct scalar lambda ; lambda.s.int_s = 2 ; lambda.type = 0 ;
と定義すれば, λ = 2 ∈ Z を代入したと考えることが出来る.
この時, 行列のスカラー倍を実現する関数は, 以下のような形式にすればよい.
struct matrix matrix_scalar_multiple(struct matrix a, struct scalar k) そして , この関数内部のアルゴリズムは ,
struct matrix matrix_scalar_multiple(struct matrix a, struct scalar k) {
int i,j ;
struct matrix b ;
b.row = a.row ; b.col = a.col ; for(i=0;i<b.row;i++)
for(j=0;j<b.col;j++)
b.coeff.XXX_c[i][j] = k.s.XXX_s * a.coeff.XXX_c[i][j] ; return b ;
}
しかし, ここでの積 * はそれぞれのオブジェクトの型が異なるため, 単純に演算子 * で処理することが出 来ないため, 実際には以下のように書く必要がある
8.
struct matrix matrix_scalar_multiple(struct matrix a, struct scalar k) {
int i,j ;
struct matrix b ;
b.row = a.row ; b.col = a.col ;
if (a.type == k.type) { /* 両方のタイプが一致している場合 */
b.type = a.type ;
if (a.type == 0) { /* 行列, スカラー：整数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.int_c[i][j] = k.s.int_s * a.coeff.int_c[i][j] ; }
else if (a.type == 1) { /* 行列, スカラー：有理数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.frac_c[i][j] = fractional_mul(k.s.frac_s,a.coeff.frac_c[i][j]) ; }
else { /* 行列, スカラー：実数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.double_c[i][j] = k.s.real_s * a.coeff.double_c[i][j] ; }
}
else if (a.type < k.type) { /* タイプが一致していない場合 */
if (k.type == 2) { /* 結果が実数になる場合 */
b.type = 2 ;
if (a.type == 0) { /* 行列：整数, スカラー：実数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.double_c[i][j] = k.s.real_s * a.coeff.int_c[i][j] ; }
else if (a.type == 1) { /* 行列：有理数, スカラー：実数 */
8もちろん,器用に書けばもう少し簡潔に書ける可能性はある.
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.double_c[i][j]
= k.s.real_s * a.coeff.frac_c[i][j].n / a.coeff.frac_c[i][j].d ; }
}
else if (k.type == 1) { /* 結果が有理数になる場合 */
b.type = 1 ;
for(i=0;i<b.row;i++) /* 行列：整数, スカラー：有理数 */
for(j=0;j<b.col;j++) {
b.coeff.frac_c[i][j].n = k.s.frac_s.n * a.coeff.int_c[i][j] ; b.coeff.frac_c[i][j].d = k.s.frac_s.d ;
} } }
else { /* (a.type > k.type) */
if (a.type == 2) { /* 結果が実数になる場合 */
b.type = 2 ;
if (k.type == 0) { /* 行列：実数, スカラー：整数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.double_c[i][j] = k.s.int_s * a.coeff.double_c[i][j] ; }
else if (k.type == 1) { /* 行列：実数, スカラー：有理数 */
for(i=0;i<b.row;i++) for(j=0;j<b.col;j++)
b.coeff.double_c[i][j]
= k.s.frac_s.n * a.coeff.double_c[i][j] / k.s.frac_s.d ; }
}
else if (k.type == 1) { /* 結果が有理数の場合 */
b.type = 1 ;
for(i=0;i<b.row;i++) /* 行列：有理数, スカラー：整数 */
for(j=0;j<b.col;j++) {
b.coeff.frac_c[i][j].n = k.s.int_s * a.coeff.frac_c[i][j].n ; b.coeff.frac_c[i][j].d = a.coeff.frac_c[i][j].d ;
} } }
/* 有理数の場合, 既約化して, 整数かどうかをチェック */
if ((b.type == 1)&&(check_fractional(&b))) b.type = 1 ; return b ;
}
ここで, check_fractional, fractional_gcd 関数は以下のように記述される.
/* 成分の既約化と整数かどうかのチェック.
* もし, 全ての成分が整数なら 1 を返し, そうでなければ 0 を返す */
int check_fractional(struct matrix *b) {
int i,j,flag=1 ;
if (b->type != 1) return 0 ; for(i=0;i<b->row;i++)
for(j=0;j<b->col;j++) {
fractional_gcd(&(b->coeff.frac_c[i][j])) ; if (b->coeff.frac_c[i][j].d != 1) flag = 0 ; }
return flag ; }
void fractional_gcd(struct fractional *c) {
int d, n, r ; n = c->n ; d = c->d ; if (c->n < 0) n = -(c->n) ; while(r = n%d) {
n = d ; d = r ; }
c->d /= d ; c->n /= d ; return ;
}
1.1.3.2.4.1 【この方法の問題点】 この方法を用いることにより, 要素の型が異なる行列でさえも演算
が行えるようになった. その代償として,
• 関数内部での演算が複雑
という欠点を持つことになる. 当然, 関数の演算速度はかなり低下することは避けられないが, 一旦関数を 作ってしまえば, 要素の型を意識することなく関数を利用できる利点を得ることになる.
1.1.3.2.5 行列の表現方法（改良その４） 次に, 要素の型を int に固定して,
• 行列サイズの最大値は可変には出来ない.
• 配列サイズが大きいため , 実引数の値を渡すためのオーバヘッドが大きい . という２つの問題をクリアしてみよう. この場合, 行列を表す構造体として struct int_matrix {
int **coeff ;
unsigned int row ;
unsigned int col ;
} ;
と定義する. この構造体のメンバー coeff は, 行列の要素を表す配列へのポインタであるので, メンバー coeff 自体のサイズはポインタのサイズ
9に過ぎない . この構造体を用いると , (1.1) の行列は次のように初 期化すればよい.
struct int_matrix a ;
int a1[2][3] = {{1,2,3},{2,3,4}} ; int *a2[2] = {a1[0], a1[1]} ; a.col = 2 ; a.row = 3 ;
a.coeff = a2 ;
この場合には, a.coeff は 「int 型のポインタの配列」a2 の先頭のアドレスを表し, a2 には, int 型の配 列 a1[0] と a1[1] の先頭のアドレスが格納される . この時 , ２つのオブジェクト a1, a2 のメモリ領域は , それぞれを定義及び初期化した時点で確保されていることに注意しよう.
一方 , A = (a
ij) ∈ M (Z, 2, 3) を a
ij= i+j −1 と
10プログラム中で値を代入したいときには , 値を格納する メモリ領域を malloc 関数を用いて確保する必要がある. メモリ領域確保のための関数を alloc_int_matrix とすると, 値の代入は以下のように行えばよい.
struct int_matrix a ; int i,j ;
a.col = 2 ; a.row = 3 ;
a.coeff = alloc_int_matrix(a.col, a.row) ; for(i=0;i<a.col;i++)
for(j=0;j<a.row;j++) a.coeff[i][j] = i+j+1 ; メモリ確保のための関数は
/* row * col の int のメモリ領域を確保する
* 失敗したら NULL を返し, 成功したらその領域へのポインタを返す */
int **alloc_int_matrix(unsigned int col, unsigned int row) {
int **ppc ; int *pc, i ;
if ((ppc = (int **)malloc(sizeof(int *)*col)) == NULL) return NULL ; for(i=0;i<col;i++) {
if ((pc = (int *)malloc(sizeof(int)*row)) == NULL) return NULL ;
*(ppc+i) = pc ; }
return ppc ; }
とすればよい. この関数の最初の malloc では, col 個の要素からなる int へのポインタの配列のメモリ 領域を確保し, 次の malloc で int 型 row 個の要素からなる配列のメモリ領域を確保しながら, その先頭 アドレスを, 最初に確保した col 個のポインタに代入している. この操作により, int 型のオブジェクトを row * col 個格納するメモリ領域を確保したことになる.
この方法を用いて, 整数係数の行列のスカラー倍を計算する関数を作ってみよう. ここでは,
9３２ビットのシステムであれば,３２ビット,６４ビットのシステムであれば,６４ビットで記述される.
10この行列Aは(1.1)そのものである.
struct int_matrix *int_matrix_scalar_multiple(struct int_matrix a, int k)
という関数を作ることにしよう . すなわち , 戻り値としては struct int_matrix へのポインタを返す
11. この関数の中身は以下のようにすれば良い.
/* 行列のスカラー倍 k a を返す.
* メモリ確保に失敗したら NULL を返す */
struct int_matrix *int_matrix_scalar_multiple(struct int_matrix a, int k) {
static struct int_matrix c ; int i,j ;
c.row = a.row ; c.col = a.col ;
if ((c.coeff = alloc_int_matrix(c.col,c.row)) == NULL) return NULL ; for(i=0;i<c.col;i++)
for(j=0;j<c.row;j++)
*(*(c.coeff+i)+j) = *(*(a.coeff+i)+j) * k ; return &c ;
}
この関数の注意事項は以下の通りである.
• 正常終了した場合の戻り値は c へのポインタを返している . 関数内部で , 計算結果を格納する変数 c
は, static に定義してあることに注意しよう. これにより, c のメモリ領域は関数外部からでも参照
可能になっている. もし, static に定義しなかったら, わざわざポインタを返しても, その領域への アクセスが出来ないことに注意.
• 実際にスカラー倍した要素を格納するメモリ領域は, alloc_int_matrix 関数を用いて確保している ことに注意 . メモリ確保に失敗した場合には NULL ポインタを返すことで , 関数を呼出した部分で関 数が正常に実行できたかどうかを判断することが可能になる.
この関数を呼出す手順は以下のようにすれば良い . struct int_matrix a, *pc ;
pc = int_matrix_scalar_multiple(a,2) ;
ただし, pc が指し示すメモリ領域は, int_matrix_scalar_multiple 関数内部で static に定義された領 域であるので, もしも,
struct int_matrix a, *pc, *pd ;
pc = int_matrix_scalar_multiple(a,2) ; pd = int_matrix_scalar_multiple(a,3) ;
のように, pc が指し示す内容を保存せずに再び関数を呼出すと, pc の内容は２回目の呼出しで上書きされ てしまう. これを防ぐためには,
struct int_matrix a, c, *pc, *pd ; pc = int_matrix_scalar_multiple(a,2) ; c = *pc ;
pd = int_matrix_scalar_multiple(a,3) ;
11戻り値は構造体へのポインタではなく,構造体を返すことも出来るが,わざわざポインタを返す理由はあとでわかる.
と, １回目の呼出しの後に結果を保存しておく必要がある.
Remark 1.1.2 このように値を保存するのが面倒だからといって ,
struct int_matrix a, c ;
&c = int_matrix_scalar_multiple(a,2) ; という方法は通用しない.
Remark 1.1.3 このように malloc 関数でメモリ領域を確保しながら計算を実行する場合には, もう一つ
面倒なことがある. 次のような２種類の呼出しを考えてみよう.
struct int_matrix a, *pc, *pd ;
pc = int_matrix_scalar_multiple(a,2) ; pd = int_matrix_scalar_multiple(a,3) ;
struct int_matrix a, c, *pc, *pd ; pc = int_matrix_scalar_multiple(a,2) ; c = *pc ;
pd = int_matrix_scalar_multiple(a,3) ; c = *pd ;
という呼出しを考えてみよう. このいずれの場合であっても, このコードの実行終了後には, １回目の関数 呼出しの内部で確保したメモリ領域を指し示すポインタは存在しない. すなわち, メモリリークが発生して いることになる. このようなメモリリークが大量に発生すると, プログラムのメモリ使用量が極めて大きく なり, 実行速度の低下や実行時エラーの発生といった問題が生じる.
このようなメモリリークの発生を防ぐことを考えよう. 一般に, malloc などの「メモリ取得関数」を用い て確保したメモリ領域は, その利用終了時に free 関数によりメモリ領域の開放を行わなければならない
12.
今回の場合は , ２重ポインタにメモリを割り当てているので , 単純に free(pc->coeff) ;
としただけでは完全にメモリを開放したことにはならない . alloc_int_matrix 関数で取得したメモリを 開放するには,
void free_int_matrix(struct int_matrix *a) {
int i ;
for(i=0;i<a->row;i++) {
if (a->type == 0) free((int *)*(a->coeff+i)) ;
else if (a->type == 1) free((struct fractional *)*(a->coeff+i)) ; else if (a->type == 2) free((double *)*(a->coeff+i)) ;
}
free(a->coeff) ; return ;
}
という関数を作り, free_int_matrix(pc) ; と呼出すのが簡単である
1314.
12つまり,freeにより,そのメモリ領域を再び他の用途に利用できるようにすることである. malloc関数はシステムコールによ りＯＳからメモリを取得するが,free関数はＯＳにメモリを返却するわけではなく,プログラム中で再利用可能にするだけである.
13一般に,オブジェクトのメモリ領域をmallocで確保した後にfreeで開放するタイミングがわからなくなることがあれば,何か マズイ書き方をしていると思えば良い.
14この関数内で,ポインタをキャストしているのは,境界整合をとるためである.
1.1.3.2.5.1 【この方法の問題点】 このように malloc 関数でメモリ操作をする方法は, 各種のプログ ラムに応用可能で , 最も標準的な方法と考えて良い . ただし , 以下のような問題点は存在する .
• 計算速度は必ずしも速くなるとは限らない. 計算速度を求めるのであれば, 簡単な２重配列の方が良 いだろう.
• メモリ割り当てと開放がちょっと面倒で, 注意深く書かないとメモリリークが発生しやすい.
1.1.3.2.6 行列の表現方法（改良その５） ここでは , 「改良その３」 (Section 1.1.3.2.4) と「改良その
４」 (Section 1.1.3.2.5) を合せてみよう. すなわち, 「その４」で使った２重ポインタを用いてメモリ確保
をする方法で, 「その３」のように Z, Q, R のいずれを係数にすることが出来るようなデータ構造を定義 し, さらに, 基本演算関数を構成してみよう.
このような演算を実現させるデータ構造としては, 「その４」で用いたものを改良して, struct int_matrix {
void **coeff ;
unsigned char type ; unsigned int row ; unsigned int col ; } ;
と定義する. ここで, メンバー type の意味は「その３」と同じとする. 「その３」では共用体を用いて行 列の成分を格納する２重配列を構成したが, ここでは「汎用データポインタ」を用いることが出来る
15. 以 下では , このデータ構造を用いて ,
• メモリ領域の確保,
• データの代入,
• スカラー倍の関数の作成 ,
• 結果の表示,
• メモリ領域の開放
という基本的な５つのパートを順に構成していこう .
1.1.3.2.6.1 メモリ領域の確保 メモリ領域の確保の段階であっても, 行列成分の違いにより確保するメ
モリ量が異なる . したがって , メモリ確保関数の引数として , 行列サイズと成分の型を渡さなければならな い. ここでは, 以下のような関数を構成することにしよう.
void **alloc_matrix(unsigned int col, unsigned int row, unsigned char type)
「その４」と同様に , メモリ確保の失敗を検出するために , 確保したメモリ領域を指し示すポインタを返す ことにする.
この関数の中身は, 基本的には「その４」で構成した alloc_int_matrix 関数と同じであるが, 確保す るポインタの指し示すオブジェクトの型の違いによって分岐をつくる必要がある. 具体的な関数の構成は 以下の通りである.
15汎用データポインタとは,どのような型のオブジェクトをも指し示すことが出来るポインタである.オブジェクトの型に関わら ず,そのアドレスを表す「数値」の型は同じであるため,このような「便利」な（でも,厄介な）概念を持ち出すことが出来る.
汎用データポインタが指し示すオブジェクトの参照のためには,「キャスト」を用いて参照の解読(dereference)を行わなければ ならない.
/* row * col のメモリ領域を確保する
* 失敗したら NULL を返し, 成功したらその領域へのポインタを返す */
void **alloc_matrix(unsigned int col, unsigned int row, unsigned char type) {
void **ppc ; void *pc ; int i ;
if (type == 0) {
if ((ppc = (void **)malloc(sizeof(int *)*col)) == NULL) return NULL ; for(i=0;i<col;i++) {
if ((pc = (void *)malloc(sizeof(int)*row)) == NULL) return NULL ;
*(ppc+i) = pc ; }
}
else if (type == 1) {
if ((ppc = (void **)malloc(sizeof(struct fractional *)*col)) == NULL) return NULL ; for(i=0;i<col;i++) {
if ((pc = (void *)malloc(sizeof(struct fractional)*row)) == NULL) return NULL ;
*(ppc+i) = pc ; }
}
else if (type == 2) {
if ((ppc = (void **)malloc(sizeof(double *)*col)) == NULL) return NULL ; for(i=0;i<col;i++) {
if ((pc = (void *)malloc(sizeof(double)*row)) == NULL) return NULL ;
*(ppc+i) = pc ; }
}
return ppc ; }
1.1.3.2.6.2 データの代入方法 データを格納するメモリ領域を確保したら, その領域に実際に値を代入
することにしよう. この部分は関数として実現するのではなく, 例えば main 関数の一部として実現すれば よい.
最も簡単なのは, 「単純型」である int または double の値を代入する方法である.
struct matrix a ; int i,j ;
a.col = 2 ; a.row = 3 ; a.type = 0 ;
a.coeff = alloc_matrix(a.col, a.row, a.type) ; for(i=0;i<a.col;i++)
for(j=0;j<a.row;j++)
(int)*((int *)*(a.coeff+i)+j) = i+j+1 ;
とすれば, A = (a
ij), a
i+1j+1= i + j − 1 を代入したことになる. ここで, (int *) を使って, int へのポ インタであることを明示的に示していることに注意しよう . これを行わないと , 元々 a.coeff は void へ のポインタ（汎用データポインタ）であるため, 参照を解決できなくなる. また, 代入式において, 汎用デー タポインタとして参照された ((int *)*(a.coeff+i)+j) に単純代入を行う際に,
*((int *)*(a.coeff+i)+j) = i+j+1 ;
としても, 右辺値が int 型の値を持つので, 正しく代入できるように思えるのだが, 代入式においては, 「式 の値の型は左オペランドの型」であり, 単純代入では, 「右オペランドの値を代入式の型に型変換し値を置 き換える」と定義されている (cf. [1, X3010, 6.3.16.1]) ため, 汎用データポインタとして参照された ((int
*)*(a.coeff+i)+j) を int へのポインタとして明示的にキャストをした方が安全である.
したがって, double の値の代入も全く同様であり, a.col = 2 ; a.row = 3 ; a.type = 2 ;
a.coeff = alloc_matrix(a.col, a.row, a.type) ; for(i=0;i<a.col;i++)
for(j=0;j<a.row;j++)
(double)*((double *)*(a.coeff+i)+j) = (double)(i+j+1) ; とすればよい.
次に , 有理数を成分に持つ行列の代入を考えよう . 有理数を表すデータ構造は , もはや単純型ではなく
「構造体」 struct fractional で表されていたことを思い出そう. したがって, 上のように単純に代入が 出来るわけではない. その代入方法には,
a.col = 2 ; a.row = 3 ; a.type = 1 ;
a.coeff = alloc_matrix(a.col, a.row, a.type) ; for(i=0;i<a.col;i++)
for(j=0;j<a.row;j++) {
((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).n = 1 ; (struct fractional *)((struct fractional *)*(a.coeff+i)+j)->d = i+j+1 ; }
とすれば, a
i+1j+1= 1/(i + j − 1) を代入したことになる. ここでは, 分子・分母で異なる代入方法を用い ている
16.
1.1.3.2.6.3 行列の表示 次に, 成分に値を代入した行列を表示する関数を作成しよう. この関数は本質
的には行列を表す構造体のみを渡せば良いが, ここでは多少汎用性を持たせるために, 出力する FILE 構造 体も渡すことにする. つまり, 以下のような関数を構成しよう.
int fprint_matrix(FILE *file, struct matrix a)
基本的には, a の中の各成分を順に fprintf 関数で出力していけば良いが, この場合も, 成分の型により出 力を切り分けなければならない.
int fprint_matrix(FILE *file, struct matrix a) {
int i,j ;
16異なる代入方法を用いる必要はなく,単に,どちらでも正しい代入式となっていることを意味する.
for(i=0;i<a.col;i++) { fprintf(file,"[") ; for(j=0;j<a.row;j++) {
if (a.type == 0) fprintf(file," %d", (int)*((int *)*(a.coeff+i)+j)) ; else if (a.type == 1) fprintf(file," %d/%u",
((struct fractional *)((struct fractional *)*(a.coeff+i)+j))->n, ((struct fractional *)((struct fractional *)*(a.coeff+i)+j))->d
) ;
else if (a.type == 2) fprintf(file," %f", (double)*((double *)*(a.coeff+i)+j)) ; }
fprintf(file,"]\n") ; }
fprintf(file,"\n") ; return 0 ;
}
1.1.3.2.6.4 メモリ領域の開放 free 関数の引数は void へのポインタ（汎用データポインタ）となっ
ているので, 「その４」で用いたものをほぼそのまま利用することが出来る.
void free_matrix(struct matrix *a) {
int i ;
for(i=0;i<a->row;i++) {
if (a->type == 0) free((int *)*(a->coeff+i)) ;
else if (a->type == 1) free((struct fractional *)*(a->coeff+i)) ; else if (a->type == 2) free((double *)*(a->coeff+i)) ;
}
free(a->coeff) ; return ;
}
1.1.3.2.6.5 スカラー倍の関数の作成 最も面倒なのが実際の演算関数の作成である . ここでは , 次のよ
うなスカラー倍の関数を作成してみよう.
struct matrix *matrix_scalar_multiple(struct matrix a, void *k, unsigned char s_type)
「その３」では「スカラー」を表す構造体を定義したが , ここでは , （手抜きをして）スカラー値を格納し たオブジェクトへのポインタと, その意味を表す type を渡すことにしよう.
基本的な関数のアイデアは「その３」で作ったように, 行列成分とスカラーの型ごとに「要素の掛け算」
を書いていけば良い.
struct matrix *matrix_scalar_multiple(struct matrix a, void *k, unsigned char s_type) {
static struct matrix c ;
int i,j ;
if (a.type == s_type) c.type = a.type ; else if (a.type > s_type) c.type = a.type ; else c.type = s_type ;
c.row = a.row ; c.col = a.col ;
if ((c.coeff = alloc_matrix(c.col,c.row,c.type)) == NULL) return NULL ; for(i=0;i<c.col;i++)
for(j=0;j<c.row;j++) {
/* 行列：整数, スカラー：整数 */
if (c.type == 0)
(int)*((int *)*(c.coeff+i)+j)
= (int)*((int *)*(a.coeff+i)+j) * (*(int *)k) ; else if (c.type == 1) {
/* 行列：整数, スカラー：有理数 */
if (a.type == 0) {
((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).n
= (int)*((int *)*(a.coeff+i)+j) * ((struct fractional *)k)->n ; ((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).d
= ((struct fractional *)k)->d ; }
/* 行列：有理数, スカラー：整数 */
else if (s_type == 0) {
((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).n
= ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).n
* (*(int *)k) ;
((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).d
= ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).d ; }
/* 行列：有理数, スカラー：有理数 */
else {
((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).n
= ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).n
* ((struct fractional *)k)->n ;
((struct fractional)*((struct fractional *)*(c.coeff+i)+j)).d
= ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).d
* ((struct fractional *)k)->d ; }
} else {
/* 行列：実数, スカラー：実数 */
if (a.type == s_type)
(double)*((double *)*(c.coeff+i)+j)
= (double)*((double *)*(a.coeff+i)+j) * (*(double *)k) ;
/* 行列：整数, スカラー：実数 */
else if (a.type == 0)
(double)*((double *)*(c.coeff+i)+j)
= (int)*((int *)*(a.coeff+i)+j) * (*(double *)k) ; /* 行列：有理数, スカラー：実数 */
else if (a.type == 1)
(double)*((double *)*(c.coeff+i)+j)
= ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).n
* (*(double *)k)
/ ((struct fractional)*((struct fractional *)*(a.coeff+i)+j)).d ; /* 行列：実数, スカラー：整数 */
else if (s_type == 0)
(double)*((double *)*(c.coeff+i)+j)
= (double)*((double *)*(a.coeff+i)+j) * (*(int *)k) ; /* 行列：実数, スカラー：有理数 */
else if (s_type == 1)
(double)*((double *)*(c.coeff+i)+j)
= (double)*((double *)*(a.coeff+i)+j)
* ((struct fractional *)k)->n / ((struct fractional *)k)->d ; }
}
return &c ; }
この関数では, ２重の for ループの中で条件分岐を行っているため, 実際の実行速度はかなり遅いと想像で きる . 条件分岐を for ループの外に出せば（すなわち , 各条件ごとに for ループを書けば）条件判断の回 数が減るため, 実行速度をあげることが出来る.
1.1.3.2.6.6 実際のプログラム ここまでで, 行列のスカラー倍を計算するプログラムを書く準備は出来
た. これを用いて, 実際にスカラー倍を計算するには,
struct matrix a, *pc ; int i,j ;
struct fractional f ;
/* 行列：実数, スカラー：有理数 */
a.col = 2 ; a.row = 3 ; a.type = 2 ; f.n = 1 ; f.d = 2 ;
a.coeff = alloc_matrix(a.col, a.row, a.type) ; for(i=0;i<a.col;i++)
for(j=0;j<a.row;j++)
(double)*((double *)*(a.coeff+i)+j) = (double)(i+j+1) ; fprint_matrix(stdout,a) ;
pc = matrix_scalar_multiple(a,&k,1) ; fprint_matrix(stdout,*pc) ;
free_matrix(&a) ; free_matrix(pc) ;
とすれば, 実数を成分に持つ行列の有理数倍を, 実数を成分に持つ行列として計算し, その結果を表示させ ることが出来る .
より安全なプログラムにしたいときには,
a.coeff = alloc_matrix(a.col, a.row, a.type) ; pc = matrix_scalar_multiple(a,&k,1) ;
のかわりに
if ((a.coeff = alloc_matrix(a.col, a.row, a.type)) == NULL) ....
if ((pc = matrix_scalar_multiple(a,&k,1)) == NULL) ....
というスタイルで, メモリ確保に失敗したときの処理を書く
17のが望ましい.
1.1.3.2.6.7 【この方法の問題点】 この方法はかなり汎用性が高い方法であるが ,
• 実行速度がかなり遅いと思われる.
高々 2 × 3 程度の行列であれば問題にはならないが, 非常に大きなサイズの行列になると, 単純なデー タ構造の方が遥かに実行速度が速いだろう.
• 演算の関数が極めて複雑になる.
という欠点を持つ . ただし , プログラムの「美しさ」の観点からいえば , 具体的な演算（成分の積や和）に 関しては, かなり工夫の余地が残されている.
1.1.3.2.6.8 演習問題
Exercise 1.1.1 この構造体を用いて , 行列の和（差）を計算する関数を書け . ２つの行列の間に和（差）
が定義できないときには, NULL を返すことにせよ.
Exercise 1.1.2 この構造体を用いて, 行列の積を計算する関数を書け. ２つの行列の間に積が定義できな
いときには, NULL を返すことにせよ.
1.1.3.2.7 最後の改良 「改良その５」では, 演算関数の中身が極めて複雑になるという欠点を持ってい
た . これを改良することで , 行列演算ライブラリの最終版を作ってみよう . 演算が複雑になっている原因は , b = k * a という演算において, k, a の型が異なり, 構造体の場合も存在することと, その結果の型が多様 であることに起因する
18. C ではこの問題を解決するために, 汎用データポインタを用いることができる が, 演算結果をポインタとして返してしまうと, 演算の記述が困難になる.
Example 1.1.1 a = b*c を計算するために, 以下のような関数を考えてみよう.
void _mul(void *a, unsigned char *a_type, void *b, unsigned char b_type, void *c, unsigned char c_type) ;
17具体的に失敗した後の処理をどのように書くかは,前後の文脈に強く依存するので,ここでは明示的には述べないこととする.
例えば,プログラムの終了の手続きに入っても良いし,シグナルハンドラを定義しておき,自分自身のプロセスに対して適当なシグ ナルを送って,ハンドラに処理をまかせても良いだろう.
18C++, Javaなどのオブジェクト指向言語の中には,「派生型」であっても単純演算子*で演算を定義できるものがある. この
場合,各クラスの間の自動型変換まで定義できるのであれば,このような複雑な問題は生じない.
