計算機基礎B - 06. 配列 (2)

この文書の URL は http://www.cc.kyoto-su.ac.jp/~mtkg/lecture/comp_B/2012/06.html です。

内積の計算

キーボードから2つのベクトルを入力し、その内積を計算するプログラムを考えてみましょう。 あらかじめ大きい配列（要素数 100）を用意しておき、その一部を利用することによって、任意の要素数のベクトルの内積が計算できるようにします。

!
! 内積の計算
!
program work0601
  implicit none
  integer :: dim
  integer, parameter :: dimmax = 100
  real(8) :: v1(dimmax), v2(dimmax)

  write(*,*) 'Input dimension:'
  read(*,*) dim                        ! 次元数を入力

  write(*,*) 'Input vector 1:'
  read(*,*) v1(1:dim)

  write(*,*) 'Input vector 2:'
  read(*,*) v2(1:dim)

  write(*,*) 'v1 . v2 = ', dot_product(v1(1:dim), v2(1:dim))

  stop
end program work0601

キーボードからベクトルを入力するには各要素を空白で区切って与えます。

% ./a.out
 Input dimension:
3
 Input vector 1:
1 2 3
 Input vector 2:
3 4 5
 v1 . v2 =    26.000000000000000

(1, 2, 3) . (3, 4, 5) = 1*3 + 2*4 + 3*5 = 26 ということで、正しく計算できているようです。

組込み関数 dot_product

組込み関数 dot_product は2つの1次元配列 vector_a, vector_b を引数にとり、それらの内積を計算します。

dot_product(vector_a, vector_b)

ただし、引数の型によって次のような違いがあるので注意しましょう。

vector_a が整数型または実数型の配列の場合

sum(vector_a * vector_b) を返します。 この場合、 dot_product(v1(1:dim), v2(1:dim)) は次の do ループと等価です。

real(8) :: tmp
  :
tmp = 0
do i = 1, dim
   tmp = tmp + v1(i) * v2(i)
end do

vector_a が複素数型の配列の場合

sum(conjg(vector_a) * vector_b) を返します。 組込み関数 conjg は複素共役を返す関数です。

行列の和

次のような行列 A, B の和 C = A + B を考えてみましょう。

A = ( 1 2 3 ),   B = ( 4 6 2 )
    ( 4 5 6 )        ( 3 1 5 )

2行3列の行列は (2, 3) 配列によって表すことができます。 2次元配列を初期化するには data 文を使うのが簡単です。 要素を並べる順番に注意して下さい（次の「配列要素順序」を参照）。

program work0602
  implicit none
  integer :: ma(2, 3), mb(2, 3), mc(2, 3), i
  data ma / 1, 4, 2, 5, 3, 6 /
  data mb / 4, 3, 6, 1, 2, 5 /

  mc = ma + mb

  do i = 1, 2
     write(*,*) mc(i, 1:3)
  end do

  stop
end program work0602

練習問題

1. mc = ma + 1 の結果を確認せよ。
2. mc = ma - mb の結果を確認せよ。
3. mc = ma * mb の結果を確認せよ。

!-----------------------------------------------------------------
! 行列の計算
! A = ( 1 2 3 ),   B = ( 4 6 2 )
!     ( 4 5 6 )        ( 3 1 5 )
!-----------------------------------------------------------------
program work0603
  implicit none
  integer :: ma(2, 3), mb(2, 3), mc(2, 3), i
  data ma / 1, 4, 2, 5, 3, 6 /
  data mb / 4, 3, 6, 1, 2, 5 /

  ! mc = ma + 1
  mc = ma + 1
  write(*,*) 'mc = ma + 1'
  do i = 1, 2
     write(*,*) mc(i, 1:3)
  end do

  ! mc = ma - mb
  mc = ma - mb
  write(*,*) 'mc = ma - mb'
  do i = 1, 2
     write(*,*) mc(i, 1:3)
  end do

  ! mc = ma * mb
  mc = ma * mb
  write(*,*) 'mc = ma * mb'
  do i = 1, 2
     write(*,*) mc(i, 1:3)
  end do

  stop
end program work0603

配列要素順序

配列の要素は「配列要素順序」 (array element order) という特別な順序で扱われることがあります。 ma(1:2, 1:3) という2次元配列では

ma(1, 1)
ma(2, 1)
ma(1, 2)
ma(2, 2)
ma(1, 3)
ma(2, 3)

が配列要素順序となります。 これは早い次元（左の添字）ほど速く変化するような順番 (column-major order) になっています。

Fortran の規格として定められているわけではありませんが、 多くの処理系ではこの順序で連続したメモリが配列に割り当てられます。 このため、左の添字に関するループを内側にすると、メモリアクセスが連続になり、高速な動作が期待できるといわれています。

data 文で配列を初期化する場合、要素の並べ方は配列要素順序となります。 例えば、配列 ma(1:2, 1:3) を

A = ( 1 2 3 )
    ( 4 5 6 )

という行列として初期化するには、次のように要素を並べます。

integer :: ma(2, 3)
data ma / 1, 4, 2, 5, 3, 6 /

行列の積

次の例はキーボードまたはファイルから n 次元の行列 A とベクトル x を読み込み、 両者の積 y = A x を計算するプログラムです。 内積の例と同様、あらかじめ大きい配列（要素数 nmax=100）を用意しておき、その一部を利用することによって、任意の次元の計算ができるようになっています。

program work0604
  implicit none
  integer :: i, j, n
  integer, parameter :: nmax = 100
  real(8) :: a(nmax, nmax), x(nmax), y(nmax)

  write(*,*) 'Input dimension:'
  read(*,*) n

  write(*,*) 'Input a vector:'
  read(*,*) x(1:n)

  write(*,*) 'Input a matrix:'
  do i = 1, n
     read(*,*) a(i, 1:n)
  end do

  y(1:n) = matmul(a(1:n, 1:n), x(1:n))

  write(*,*) y(1:n)

  stop
end program work0604

ファイルからの入力

ベクトルや行列の要素ををキーボードから毎回入力するのは面倒です。 そこで、シェルのリダイレクト機能を利用して、ファイルから入力してみましょう。 以下のような内容のファイルを作成しておきます。 ファイル名は input.dat としましょう。

4
 1.0  -2.0   3.0   4.0
 2.0   1.0  -3.0   2.0
 3.0   4.0   5.0   6.0
-1.0   2.0   0.0   3.0
 2.0  -6.0  -7.0   5.0

各行の意味は次のようになります。

• 1行目は次元 n の値 (read(*,*) n で読み込まれる)
• 2行目はベクトル x(1:n) の内容 (read(*,*) x(1:n) で読み込まれる)
• 3-6行目は行列 a(1:n, 1:n) の内容
• 3行目が行列の第1行、4行目が第2行、…、6行目が第4行に対応

キーボードの代わりにファイルから入力するにはリダイレクト記号 < を使います。

% ./a.out < input.dat
 Input dimension:
 Input a vector:
 Input a matrix:
  -1.0000000000000000  34.000000000000000   7.0000000000000000   13.000000000000000

cat というプログラムを利用して、次のようにすることもできます。

% cat input.dat | ./a.out

組込み関数 mutmul(matrix_a, matrix_b)

組込み関数 mutmul を使うと行列と行列、行列と列ベクトル、行ベクトルと行列の積が計算できます。 matrix_a, matrix_b は2次元または1次元の配列で、2次元配列の場合は行列、1次元配列の場合はベクトルと解釈されます。 列ベクトルと行ベクトルはどちらも同じ1次元配列で表されるので、文脈に応じて適切に解釈しましょう。

• matrix_a の形状が (n, m), matrix_b の形状が (m, k) の場合、 行列と行列の積が計算され、結果の形状は (n, k)、つまり (n, k) 行列となる。
• matrix_a の形状が (m), matrix_b の形状が (m, k) の場合、 m 次元行ベクトルと (m, k) 行列の積が計算され、結果の形状は (k)、つまり k 次元行ベクトルとなる。
• matrix_a の形状が (n, m), matrix_b の形状が (m) の場合、 (n, m) 行列と m 次元列ベクトルの積が計算され、結果の形状は (n)、つまり n 次元列ベクトルとなる。

書式

数値の全体の桁数や小数点以下の桁数といった出力形式のことを「書式」といいます。 Fortran では書式を詳しく指定することができます。

program work0605
  implicit none
  integer :: n
  real(8) :: x, y

  do n = 1, 5
     x = dble(n)
     y = sqrt(x)
     ! 書式を指定した write 文の例
     write(*,'(I5,5X,F15.5,5X,E15.7)') n, x, y
  end do

  stop
end program work0605

このプログラムを実行すると次のような出力が得られます。 （注意：1行目は桁数をみやすくするための目盛りで、プログラムの出力ではありません。）

12345678901234567890123456789012345678901234567890
% ./a.out
    1             1.00000       0.1000000E+01
    2             2.00000       0.1414214E+01
    3             3.00000       0.1732051E+01
    4             4.00000       0.2000000E+01
    5             5.00000       0.2236068E+01

書式は write 文の2つ目の要素に指定します。 上の例では次のような書式が使われています。

• I5 : 整数型変数 n を5桁
• 5X : 空白を5桁
• F15.5 : real(8) 型の変数 x を全体を 15 桁、小数点以下を 5 桁
• 5X : 空白を5桁
• E15.7 : real(8) 型の変数 y を全体を 15 桁、小数点以下を 7 桁

個々の数字の桁数などを指定する I5 や F15.5 などを「編集記述子」と呼びます。 書式は一般に

'(s1,s2,...,sn)'

と書きます。 ここで s1, s2, ..., sn は個々の編集記述子です。

編集記述子

編集記述子には非常に多くの種類があります。 ここでは代表的なものだけを示します。

整数型の編集記述子

整数型データの書式指定には Iw または Iw.m 形の編集記述子を用います。 Iw は整数を w 桁に出力することを表します。 Iw.m も同様ですが、符号なしの整数が少なくとも m 桁になるように 0 が付加されます。 したがって w >= m でなければなりません。

整数値     編集記述子          出力欄
-----------------------------------------------------------------
 12345     I5                  12345
-12345     I5                  *****
 12345     I10                      12345
 12345     I10.6                   012345
-12345     I10.7                 -0012345

指定した書式で出力できない場合は * が出力されます。

実数型の編集記述子

実数型データの書式指定には Fw.d または Ew.d 形の編集記述子を用います。 w は欄の幅、 d は小数部の桁数を表します。 E 形の編集記述子は次のように実数を仮数部と指数部に分けて表示します（基数は 10 です）。

[±][0].x1x2...xdE±z1z2

例えば 0.12E+03 は 0.12×10^3 を表します。

実数値     編集記述子          出力欄
-----------------------------------------------------------------
 123.456   F7.2                 123.456
-123.456   F10.5               -123.456
 12.3456   F5.3                *****
0.012345   E12.5                0.12345e-01
-123.456   E12.5               -0.12346e+03
 1.23456   E5.2                *****

以下の点に注意しましょう。

• Fw.d 形編集記述子では符号と小数点のための欄を考慮して w > d+2 でなければならない。
• Ew.d 形編集記述子では符号、整数部の 0、仮数部、指数部のための欄を考慮して w > d+7 でなければならない。

A 形編集記述子

文字型データの書式指定には Aw 形編集記述子を用います。 w は欄の幅です。

文字列     編集記述子          出力欄
-----------------------------------------------------------------
'ABCDE'    A                   ABCDE
'ABCDE'    A7                    ABCDE
'ABCDE'    A3                  ABC

X 形編集記述子

X 形編集記述子は空白を意味します。 一般に nX （n は整数）という形式で用い、 出力においては n 個の空白の出力、 入力においては n 桁だけ読み飛ばすことを意味します。

反復可能編集記述子

I 形、F 形、E 形、X 形などの編集記述子は 10I5 のように左側に反復回数を指定することができます。

練習問題

次のプログラムの出力を確認してみましょう。

program work0606
  implicit none
  integer :: i

  write(*,'(4I10)') (i, i = 1, 10)

  stop
end program work0606

素数を求める

今後は素数を求めるプログラムを考えてみましょう。 素数は1と自分自身以外の約数を持たない正整数のことです。 整数 n が素数か判定するには、n 未満のすべての素数で n が割り切れるかチェックします。 割り切れなければ n は素数、割り切れれば合成数となります。 以下のプログラムではみつかった素数を prime(:) という配列に記憶し、チェックのために利用しています。 何をやっているか、がんばって解読してみましょう。

!-----------------------------------------------------------------
! 1000 以下の素数を求める
!-----------------------------------------------------------------
program work0607
  implicit none
  integer, parameter :: nmax = 1000
  integer :: np = 0             ! みつかった素数の個数
  integer :: prime(nmax)        ! 素数を記憶する配列（最大 nmax 個）
  integer :: n, i

  ! 2 は素数
  np = 1
  prime(1) = 2

  ! 3 以上 nmax 以下の奇数 n だけをチェックする。
  do n = 3, nmax, 2
     ! n がこれまでにみつかった素数で割り切れるか調べる。
     do i = 1, np
        if (mod(n, prime(i)) == 0) exit
     end do
     ! n がそれまでに見つかった素数で一度も割り切れなければ素数である。
     ! その場合は exit が実行されないので、do ループが終了した時点で i
     ! の値は np + 1 になる。したがって、次の if 文の中身は n が素数の
     ! ときだけ実行される。
     if (i > np) then
        np = np + 1             ! みつかった素数の個数を 1 増やす
        prime(np) = n           ! n を素数表に登録する
     end if
  end do

  ! 結果表示
  write(*,'(I4,X,A)') np, 'primes found.'
  write(*,'(10I5)') prime(1:np)

  stop
end program work0607

exit 文

prime(i) にはそれまでにみつかった素数が記憶されています。

• n が prime(i) で割り切れると exit 文が実行され i に関する do ループから脱出します。 この場合、 i の値は np 以下なので i > np は偽となり if 文の中のブロックは実行されません。
• n が prime(i) で一度も割り切れない場合は do ループの中の exit 文は実行されません。 この場合、do ループの終了時に i の値は np + 1 となることに注意しましょう。 i > np が真となり if 文の中のブロックが実行されます。

書式の指定

write(*,'(10I5)') prime(1:np) では「各行に5桁の整数を10個ずつ出力する」という書式が指定されています。

本日の課題

1. 平均値

キーボードから5つの整数を読み込み、平均値を求めるプログラムを作成せよ。

!------------------------------------------------------------
! 平均値
!------------------------------------------------------------
program work0609
  implicit none
  integer :: x(5)

  write(*,*) 'Input 5 integers:'
  read(*,*) x(1:5)

  write(*,*) 'mean = ', sum(x) / 5.0d+0

  stop
end program work0609

2. 平均値 (2)

キーボードから n 個の整数を読み込み、平均値を求めるプログラムを作成せよ。

!------------------------------------------------------------
! 平均値 (2)
!------------------------------------------------------------
program work0610
  implicit none
  integer :: x(100), n

  write(*,*) 'Input N (<= 100):'
  read(*,*) n

  write(*,*) 'Input N integers:'
  read(*,*) x(1:n)

  ! 部分配列 x(1:n) の和を取り、要素数 n で割る。
  ! dble(n) は整数を real(8) 型の実数に変換する関数。
  write(*,*) 'mean = ', sum(x(1:n)) / dble(n)

  stop
end program work0610

3. 平均値 (3)

ファイルから次のような形式のデータを読み込み、身長と体重の平均値を求めるプログラムを作成せよ。 1行目は人数、2行目以降が身長と体重の数値である。

5
155.8   40.8
178.3   67.1
167.2   50.6
162.8   45.7
170.9   63.9

!------------------------------------------------------------
! 平均値 (3)
!------------------------------------------------------------
program work0611
  implicit none
  integer :: n, i, w(1:100), h(1:100)

  read(*,*) n                   ! 人数を n にセット
  do i = 1, n
     read(*,*) w(i), h(i)       ! 体重と身長を w(i), h(i) に読み込む
  end do

  write(*,*) 'W mean: ', sum(w(1:n)) / dble(n)
  write(*,*) 'H mean: ', sum(h(1:n)) / dble(n)

  stop
end program work0611

4. 平均と分散

キーボードから n 個の整数 x1, x2, ..., xn を読み込み、平均 xmean と分散 sigma を求めよ。 分散の2乗は次のように定義される。

            Σ (xi - xmean)**2
sigma**2 = --------------------     （総和 Σ は i = 1, 2, .., n について取る）
                    n

!------------------------------------------------------------
! 分散 (variance)
!------------------------------------------------------------
program work0612
  implicit none
  integer :: x(100), n, i
  real(8) :: xmean, sigma

  write(*,*) 'Input N (<= 100):'
  read(*,*) n

  write(*,*) 'Input N integers:'
  read(*,*) x(1:n)

  ! 平均
  xmean = sum(x(1:n)) / dble(n)

  ! 分散
  sigma = 0
  do i = 1, n
     sigma = sigma + (x(i) - xmean)**2
  end do
  sigma = sqrt(sigma)

  ! 上の5行の代わりに次のようにしてもよい。
  ! 配列演算を使えば1行で済む。
  !!! sigma = sqrt( sum( (x(1:n) - xmean)**2 ) )

  write(*,*) 'mean:     ', xmean
  write(*,*) 'variance: ', sigma

  stop
end program work0612

チャレンジ問題

1. 覆面算

次の覆面算を解くプログラムを作成せよ。 5つの未知数バ、ナ、シ、モ、ンを do ループで回して総当たりでチェックすればよい。 覆面算を知らない人はこちらをチェック。 ちなみに、正解は 877+877=1754 である。

バナナ＋バナナ＝シナモン

!------------------------------------------------------------
! 覆面算：バナナ＋バナナ＝シナモン
! 本当は次のようなチェックも必要だが、やっていない。
!  - ba, si は 0 ではない。
!  - ba, na, si, mo, nn はすべて異なる。
!------------------------------------------------------------
program work0613
  implicit none
  integer :: ba, na, si, mo, nn
  integer :: banana, sinamonn

  do ba = 0, 9
     do na = 0, 9
        do si = 0, 9
           do mo = 0, 9
              do nn = 0, 9
                 banana   =             ba * 100 + na * 10 + na
                 sinamonn = si * 1000 + na * 100 + mo * 10 + nn
                 if (banana + banana == sinamonn) then
                    write(*,*) banana, sinamonn
                 end if
              end do
           end do
        end do
     end do
  end do

  stop
end program work0613

2. 公約数（ユークリッドの互除法）

ユークリッドの互除法を用いてキーボードから入力した2整数 a, b の最大公約数を求めるプログラムを作成せよ。 ユークリッドの互除法のアルゴリズムは次のようなものである。

• step 1: a, b をそれぞれ変数 m, n に代入する。
• step 2: m を n で割ったあまりを k とする（組み込み関数 mod を利用）。
• step 3: k が 0 なら n が求める最大公約数である。k が 0 でなければ step 4 へ。
• step 4: n の値を m に代入し、k の値を n に代入して step 2 へ戻る。

step 4 → step 2 のループには繰り返し回数を指定しない do 文を用いる。

!------------------------------------------------------------
! 最大公約数
!------------------------------------------------------------
program work0614
  implicit none
  integer :: a, b
  integer :: m, n, k

  write(*,*) 'Input two positive integer:'
  read(*,*) a, b

  ! step 1
  m = a
  n = b
  do
     ! step 2
     k = mod(m, n)
     ! step 3
     if (k == 0) then
        write(*,*) 'GCM: ', n
        exit
     end if
     ! step 4
     m = n
     n = k
  end do

  stop
end program work0614

3. 組み合わせ

n 個の整数 (1, 2, .., n) から r 個の整数を取り出す方法は nCr 通りある。 例えば、n = 4, r = 2 の場合は次の 6 通りである。

(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)

キーボードから n と r を入力し、すべての組み合わせを出力するプログラムを作成せよ。

!------------------------------------------------------------
! n 個の数字 (1, 2, .., n) から r 個の数字を取り出す
! --------------------------------------------------
! 1 から n までの数 i を、i 桁目が 1 で他の桁はすべて 0 であるような
! n 桁の2進数と対応させる。例えば n = 4 のときは
!     1 -> 0001   (10進数で表すと 2**0)
!     2 -> 0010   (10進数で表すと 2**1)
!     3 -> 0100   (10進数で表すと 2**2)
!     4 -> 1000   (10進数で表すと 2**3)
! と対応させる。i に対応する2進数の10進数表記は 2**(i-1) であること
! に注意。
!
! n 個の数字 (1, 2, .., n) から r 個の数字を取り出した組み合わせは、
! r 個の桁が 1 で他の桁が 0 であるような2進数として表される。例えば、
! n = 4, r = 2 の場合、2 個の桁だけが 1 であるような2進数と、数字の
! 組み合わせの対応は
!     0011 -> (1, 2)     (0011 の1桁目が 1, 2桁目が 2 に対応)
!     0101 -> (1, 3)     (0101 の1桁目が 1, 3桁目が 3 に対応)
!     1001 -> (1, 4)     (以下同様)
!     0110 -> (2, 3)
!     1010 -> (2, 4)
!     1100 -> (3, 4)
! となる。したがって、(1, 2, 3, 4) から 2 個の数字を取り出した組み
! 合わせを得るには、0000 から 1111 までのすべての2進数から、各桁の
! 和が 2 になるようなものを探し出せばよい。ここで 1111 を10進数で
! 表すと (2**4 - 1) であることに注意する。
!
! 以上の考察から、n 個の数字 (1, 2, .., n) から r 個の数字の組み合
! わせを取り出すには、0 から (2**n - 1) までの整数を2進数と対応させ、
! 各桁の和が r になるもの（つまり、r 個の桁が 1 で他の桁は 0 であ
! るもの）を探せばよいことが分かる。
!------------------------------------------------------------
program work0615
  implicit none
  integer :: n, r
  integer :: i, j, k, tmp
  integer :: b(100), d(100)     ! 100 は適当。大きめに取っておけばよい。

  write(*,*) 'Input N and R:'
  read(*,*) n, r

  do i = 0, 2**n - 1
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     ! i を n 桁の2進数で表現し、配列 b(1:n) に入れる。
     ! n = 4 の場合:
     !   i   ->  b(1:n)
     !   ------------------------
     !   0   ->  (/ 0, 0, 0, 0 /)
     !   1   ->  (/ 1, 0, 0, 0 /)
     !   2   ->  (/ 0, 1, 0, 0 /)
     !   3   ->  (/ 1, 1, 0, 0 /)
     !       …
     !   15  ->  (/ 1, 1, 1, 1 /)
     ! b(1:n) は通常の2進数表記と左右が逆になっていることに注意。
     tmp = i
     do j = 1, n
        b(j) = mod(tmp, 2)
        tmp = tmp / 2
     end do

     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     ! チェック用に i の値と b(1:n) の内容を出力
     ! 見やすいように b(n:1:-1) として通常の2進数表記と同じ順序に直す。
     ! チェックが終わったらコメントアウト。
     !!! write(*,'(I5,3X,30I1)') i, b(n:1:-1)

     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     ! b(1:n) の各桁を足して r にならなければ次の i へ移る。
     ! cycle 文は do ループ内の残りの部分を省略し、次の繰り返しに移る命令。
     if (sum(b(1:n)) /= r) cycle

     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     ! b(1:n) から数字の組み合わせに戻し、配列 d(1:r) に格納する。
     ! n = 4, r = 2 の場合:
     !   b(1:n)            ->  d(1:r)
     !   --------------------------------
     !   (/ 1, 1, 0, 0 /)  ->  (/ 1, 2 /)
     !   (/ 1, 0, 1, 0 /)  ->  (/ 1, 3 /)
     !      …
     !   (/ 0, 0, 1, 1 /)  ->  (/ 3, 4 /)
     k = 0
     do j = 1, n
        ! b(1:n) の j 桁目が 1 なら j は組み合わせの要素。
        if (b(j) == 1) then
           k = k + 1
           d(k) = j
        end if
     end do

     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     ! 結果表示
     write(*,'(10I5)') d(1:r)
  end do

  stop
end program work0615