計算機基礎B - 04. 繰り返し

目次

リンク

一定回数の繰り返し

次のプログラムは 1 から 10 までの和を求めて結果を表示します。

program work0401
  implicit none
  integer :: i, tmp = 0

  do i = 1, 10
     tmp = tmp + i
  end do

  write(*,*) tmp

  stop
end program work0401

変数の初期化

Fortran 90 以降では型宣言と同時に変数に初期値を与えることができます。

integer :: i, tmp = 0

では整数型の変数 i と tmp が宣言されていますが、 tmp だけに初期値 0 が与えられています。 変数に初期値を与えることを「初期化」といいます。 i は初期化されず、値は不定（どのような値になるかわからない）となります。

FORTRAN77 では変数の初期化は data 文で行い、型宣言文の中で初期値を与えることはできません。 data 文についてはあとで説明します。

do 構文の基本形

特定の処理をある回数だけ繰り返す場合には do 構文（do ループ）が用いられます。 do 文の一般的な書式は次のようになります。

do var = e1, e2 [, e3]
   block
end do

var は整数型の変数で「do 変数」と呼ばれます。 e1 は初期値パラメータ、 e2 は終値パラメータ、 e3 は増分パラメータと呼ばれ、いずれも整数型です。 増分パラメータ e3 は省略可能で、省略時は 1 が指定されたものとみなされます。 do 文が実行されると、変数 var は e1 にセットされ、 do と end do で囲まれたブロックが実行されます。 次に var の値が e3 だけ増やされ、再度ブロックが実行されます。 var の値が e2 より大きくなるまでこれが繰り返されます。

プログラムの流れ

プログラム中の do ループに注目してみましょう。

do i = 1, 10
   tmp = tmp + i
end do

初期値パラメータは 1、 終値パラメータは 10 です。 増分パラメータは省略されているので 1 と見なされます。 したがって、整数型の変数 i の値は 1, 2, 3, …, 10 と変化することがわかります。

do ループの処理が開始されると、変数 i は 1 にセットされます。 tmp は 0 に初期化されていることに注意しましょう。

   tmp = tmp + i

の右辺は 0 + 1 となるので、 tmp の値は 1 に更新されます。 これで1回目の処理が完了します。

次に i は 2 にセットされます。 tmp の値は 1 になっているので

   tmp = tmp + i

の右辺は 1 + 2 となり、 tmp の値は 3 に更新されます。

この処理は i が 10 になるまで繰り返されます。 結果として tmp には 1 から 10 の和が入ることになります。

do ループの実行回数

do 変数の初期値が終値を超えている場合、do ループの繰り返し回数は 0 となり、ループの内側のブロックは1回も実行されません。

do i = 1, 0
   :               (1回も実行されない)
end do

初期値と終値が等しい場合は1回だけ実行されます。

do i = 1, 1
   :               (i = 1 として1回だけ実行される)
end do

ループ終了後の do 変数の値

do ループの実行が終わったとき、do 変数の値はどうなるでしょうか？ 実はループの実行が終了したときの値が保持されます。 例えば

do i = 1, 10
   :
end do

では i は 11 となります（10 ではないことに注意しましょう）。 exit 文などで do ループの処理が中断された場合は、その時点での i の値となります。

練習問題

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

program work0402
  implicit none
  integer :: i, n, tmp = 0

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

  do i = 1, n
     tmp = tmp + i
  end do

  write(*,*) tmp

  stop
end program work0402

減少ループ

キーボードから整数 n を入力し、n の2重階乗 n!! を求めるプログラムを考えてみましょう。 2重階乗 n!! は次のように定義されます。

 n!! = 1 * 3 * 5 * … * n    (n が奇数の場合)
       2 * 4 * 6 * … * n    (n が偶数の場合)
 0!! = 1
-1!! = 1

定義をみると n の偶奇で場合分けする必要があるように思えますが、do ループをうまく使うと不要です。

program work0403
  implicit none
  integer :: i, n
  real(8) :: tmp = 1

  write(*,*) 'Input an integer (>= -1):'
  read(*,*) n

  do i = n, 1, -2
     tmp = tmp * i
  end do

  write(*,*) tmp

  stop
end program work0403

プログラムの流れ

変数 i の値がどのように変化するか考えてみましょう。

n = 5 （奇数）の場合

増分パラメータが -2 なので i の値は 5 から 2 ずつ減少していきます。 終値パラメータが 1 なので i の値は 5, 3, 1 と変化することがわかります。

n = 6 （偶数）の場合

i の値は 6, 4, 2 と変化します。 終値パラメータが 1 なので do ループの中で i の値は 0 にはなりません。

tmp の値

実数型の変数 tmp には初期値 1.0 が与えられています。 n = 5 の場合、do ループの処理が開始されると変数 i は 5 にセットされます。

tmp = tmp * i

の右辺は 1.0d+0 * 5 となるので、 tmp は 5.0d+0 に更新されます。

次に i には 3 がセットされます。

tmp = tmp * i

の右辺は 5.0d+0 * 3 となるので、 tmp は 15.0d+0 に更新されます。

この処理は i が 1 になるまで繰り返されます。 結果として tmp には 5 * 3 * 1 の値が入ります。

整数の範囲

このプログラムでは途中結果を記憶する作業変数として real(8) 型の変数 tmp を利用しています。 これは整数型では計算できる値の範囲が狭いためです。

整数型 integer が表すことのできる値の範囲は次のようになっています。

最小値: -2147483648
最大値:  2147483647

もっと大きな整数が必要な場合は integer(8) 型を使います。 integer(8) 型の範囲は

最小値: -9223372036854775808
最大値:  9223372036854775807

となります。

多重ループ

do ループを入れ子にすることで複数の do 変数を独立に変化させることができます。 例として、かけ算九九の表を作成するプログラムを示します。 1 から 9 までのすべての i に対して j の do ループが実行されます。

program work0404
  implicit none
  integer :: i, j

  do i = 1, 9
     do j = 1, 9
        write(*,*) i, ' * ', j, ' = ', i * j
     end do
  end do

  stop
end program work0404

字下げ

変数 j のループは行頭に空白を入れて右にずらして書いてあります。 このように書くことで、変数 j の do ループが変数 i の do ループの中に包含されていることが視覚的にはっきりします。 このような記法を「字下げ」（インデント）といいます。

Emacs では TAB キーを押すことによって適切な字下げを行うことができます。

漸化式の計算

do ループを利用して、フィボナッチ数列を第 20 項まで計算するプログラムを考えてみましょう。 フィボナッチ数列 Fn は次のように定義されます。

Fn = Fn-1 + Fn-2    (n >= 3)

初項 F1, F2 はキーボードから与えます。 プログラムは次のように書けます。

program work0405
  implicit none
  integer :: i, a, b, c
  integer, parameter :: imax = 20

  write(*,*) 'Input F1 and F2:'
  read(*,*) a, b

  do i = 3, imax
     c = a + b
     write(*,*) i, c
     ! 繰り返しの準備
     a = b
     b = c
  end do

  stop
end program work0405

プログラムの流れ

c = a + b で漸化式を計算します。 続く a = b, b = c は次の繰り返しのための準備です。 結果として c には F3, F4, F5, …, F20 の値が入ります。

漸化式を 20 項計算するからといって 20 個の変数が必要になるわけではないことに注意しましょう。 一般に、n 項間漸化式の計算には n 個の変数が必要です。

名前付き定数

プログラム中に定数を直接埋め込んでしまうと、定数の値を変更するときに変更漏れなどの間違いが生じやすくなります（プログラム中に埋め込まれた定数のことを「マジックナンバー」といいます）。 変数を名前付き定数として宣言しておけば変更は宣言部分だけで済みます。 名前付き定数を定義するには、型宣言に parameter 属性を付加します。

integer, parameter :: imax = 100
real(8), parameter :: gravity = 9.8d+0

もちろん定数の値をプログラム中で変更することはできません。

練習問題

次の漸化式で定義される数列 An を 10 項まで計算せよ。

An = 2 * An-1 - 1      (n >= 2)
A1 = 2

2項間漸化式なので数列を表現する変数が2つ必要である。

!
! 数列 An を 20 項まで求める
! An = 2 * An-1 - 1      (n >= 2)
! A1 = 2
!
program work0406
  implicit none
  integer :: i, a, b

  a = 2                         ! A1 をセット
  do i = 2, 10
     b = 2 * a - 1              ! 漸化式 Ai = 2 * Ai-1 - 1 の計算
     write(*,*) i, b
     a = b                      ! 次の繰り返しのための準備
  end do

  stop
end program work0406

反復回数を指定しない do 文

あらかじめ繰り返し回数がわからない場合にも do 文が使えます。 その場合、do 変数や初期値パラメータ、終値パラメータはすべて省略されます。

do
   :
   if (ループからの脱出条件) exit
   :
end do

exit 文

明示的に繰り返しを中断しないとループ内の処理がいつまでも繰り返されます（無限ループ）。 そのため、通常は do ループの中に脱出のための条件判定が必要になります。 do ループの中で exit 文が実行されると処理が end do 以降に移され do ループから脱出することができます。

例として、Newton 法で正数 a の平方根を求めるプログラムを考えましょう。

Newton 法

方程式 f(x) = 0 の根は次のように求めることができます。 x0 を根の第1近似とします。 曲線 y = f(x) 上の点 (x0, f(x0)) における接線が x 軸と交わる点を x1 とすると、 x1 は x0 よりも f(x) = 0 の根に近い値となります。 第 n 近似 xn が求まったとすると次の近似 xn+1 は次の式で求めることができます。

xn+1 = xn - f(xn) / f'(xn)

これを繰り返すことによって方程式の根を得ることができます。 これを Newton 法といいます。 【Wikipedia の説明をみる】

解の収束は xn と xn+1 の差が十分小さくなったかどうかで判定します。

平方根の計算

f(x) = x*x - a とすると f(x) = 0 の根が a の平方根を与えます。 Newton 法の一般論より

xn+1 = (xn + a / xn) / 2

で作られる数列が √a に収束します。

!
! Newton 法で a の平方根を求める
!
program work0407
  implicit none
  real(8) :: a, xold, xnew
  real(8), parameter :: eps = 1.0d-5

  write(*,*) 'Input a positive real:'
  read(*,*) a

  ! 第 0 近似を a とする
  xold = a

  ! 漸化式の計算
  do
     xnew = (xold + a / xold) / 2
     if (abs(xnew - xold) < eps) exit
     xold = xnew
  end do

  ! 結果表示
  write(*,*) xnew

  stop
end program work0407

2.0 の平方根を求めてみましょう。

% ./a.out
 Input a positive real:
2.0
   1.4142135623746899

正しく計算できているようです。

収束判定

if (abs(xnew - xold) < eps) で解が収束したか判定しています。 組込み関数 abs(x) は実数型変数 x の絶対値を返す関数です。 xnew と xold の差の絶対値がしきい値 eps より小さくなると exit が実行され do ループから脱出します。

本日の課題

2乗和

キーボードから自然数 n, m を入力し、1 から n までの m 乗和

1**m + 2**m + ... + n**m

を求めるプログラムを作成せよ。

!
! 1**m + 2**m + ... + n**m を求める
!
program work0408
  implicit none
  integer :: i, n, m, tmp = 0

  write(*,*) 'Input N and M:'
  read(*,*) n, m

  do i = 1, n
     tmp = tmp + i**m
  end do

  write(*,*) tmp

  stop
end program work0408

階乗の計算

キーボードから整数 n を入力し、n の階乗 n! を求めるプログラムを作成せよ。 0! = 1 も正しく計算できるようにすること。

!
! n の階乗 n! を計算する
!
program work0409
  implicit none
  integer :: n, i
  real(8) :: tmp

  write(*,*) 'Input an integer (N >= 0):'
  read(*,*) n

  ! n が負の場合は stop 文でプログラムを途中終了
  if (n < 0) then
     write(*,*) 'N must be >= 0'
     stop
  end if

  tmp = 1.0d+0
  do i = 1, n
     tmp = tmp * i
  end do

  write(*,*) 'N! = ', tmp

  stop
end program work0409

立方根の計算

Newton 法で正数 a の立方根を求めるプログラムを作成せよ。 f(x) = x^3 - a とすると f(x) = 0 の解が立方根を与える。 Newton 法の一般論から、漸化式は

xn+1 = xn - (xn**3 - a) / (3 * xn**2)
     = (2 * xn + a / xn**2) / 3

で与えられることを確認せよ。

!
! a の立方根を計算する
!
program work0410
  implicit none
  略
  stop
end program work0410

チャレンジ問題

コラッツの問題 (Collatz problem)

任意の正整数 n に対して次の操作を繰り返すと有限回で 1 に到達すると予想されている。

• n が偶数なら 2 で割る
• n が奇数なら n に 3 をかけて 1 を足す

これを「コラッツの問題」という。 例えば、n が 3 のときは7回の操作で 1 となる。

3 → 10 → 5 → 16 → 8 → 4 → 2 → 1

キーボードから n を入力し 1 になるまでこの操作を行うプログラムを作成せよ。

!
! コラッツの問題
!
program work0411
  implicit none
  integer :: n, nmax            ! nmax は n の最大値を記憶
  integer :: c = 0              ! n が 1 になるまでの操作の回数

  write(*,*) 'Input an integer (>= 2):'
  read(*,*) n

  nmax = n                      ! nmax を初期化

  do
     write(*,*) n
     if (n == 1) exit           ! 1 になったら do ループから脱出
     if (mod(n, 2) == 0) then
        n = n / 2
     else
        n = 3 * n + 1
     end if
     if (n > nmax) nmax = n     ! n > nmax なら nmax を更新
     c = c + 1                  ! 操作の回数を数える
  end do

  write(*,*) 'Nmax:  ', nmax
  write(*,*) 'Total: ', c

  stop
end program work0411

素数の判定（単純版）

奇数 n が与えられたとき、それが素数かどうか、次のように判定するプログラムを作成せよ。

• k = 3, 5, 7, ..., n-2 に対して n が k で割り切れるか調べる。
• 途中で割り切れれば素数でないと判定する。
• すべての k に対して割り切れなければ素数と判定する。

!
! 素数の判定（単純版）
!
program work0412
  implicit none
  integer :: n, k

  ! n を入力する
  write(*,*) 'Input an odd integer:'
  read(*,*) n

  ! k = 3, 5, 7, ..., n-2 で n が割りきれるかチェック
  do k = 3, n-2, 2
     ! n が k で割りきれたら素数ではない。
     ! その場合は stop 文でプログラムを終了する。
     if (mod(n, k) == 0) then
        write(*,*) n, ' is not prime'
        stop
     end if
  end do

  ! 割りきれなかったら素数
  write(*,*) n, ' is prime'

  stop
end program work0412

約数を求める

キーボードから整数 n を読み込んで、n の約数をすべて求めるプログラムを作成せよ。

!
! 約数を求める
!
program work0413
  implicit none
  integer :: n, i

  write(*,*) 'Input an integer:'
  read(*,*) n

  ! n を i で割って割りきれたら約数
  do i = 2, n
     if (mod(n, i) == 0) write(*,*) i
  end do

  stop
end program work0413

数表の作成

0.1 から 0.2 きざみで 3.5 までの x に対する平方根の表を作成せよ。 平方根の計算には組み込み関数 sqrt(x) を用いてよい。

注意: do 変数は整数型、 x は real(8) 型にすること。 整数 i から 0.1, 0.3, 0.5, ..., 3.5 を作るには次のようにすればよい。

do i = 0, 17
   x = 0.1d+0 + 0.2d+0 * i
     :
end do

!
! 数表の作成
!
program work0414
  implicit none
  integer :: i
  real(8) :: x

  do i = 0, 17
     x = 0.1d+0 + 0.2d+0 * i
     write(*,*) x, sqrt(x)
  end do

  stop
end program work0414

等比級数

どんな高さから落としても、もとの高さの 0.65 倍の高さまで跳ね返るボールがある。 このボールを h メートルの高さから落としたら、何回目でボールが 1 メートル以下になるか求めよ。 （はじめの高さ h はキーボードから入力する。）

!
! 等比級数
!
program work0415
  implicit none
  real(8) :: h
  integer :: n                  ! 跳ねた回数

  write(*,*) 'Input initial height:'
  read(*,*) h

  n = 0                         ! 跳ねた回数を初期化
  do
     h = h * 0.65d+0
     n = n + 1
     write(*,*) n, h
     if (h <= 1.0d+0) exit      ! 高さが 1 m 以下になったら do ループを抜ける
  end do

  stop
end program work0415

整数の問題

正の整数 n をキーボードから読み込んで m**2 <= n < (m+1)**2 となるような整数 m を求めよ。

!
! 整数の問題
! m**2 <= n < (m+1)**2 となるような整数 m を求める。
!
program work0416
  implicit none
  integer :: n, m

  write(*,*) 'Input N (>= 1):'
  read(*,*) n

  do m = 1, n
     if (m**2 <= n .AND. n < (m+1)**2) exit
  end do

  write(*,*) 'M = ', m

  stop
end program work0416

整数の問題 (2)

3桁の整数で、各桁の数字の3乗和がもとの数と等しくなるようなものをすべて求めよ。

!
! 整数の問題 (2)
! 3桁の整数で、各桁の数字の3乗和がもとの数と等しくなるようなもの
!
program work0417
  implicit none
  integer :: m, n               ! 3桁の整数 (234 など)
  integer :: i, j, k            ! 各桁の数字 (i=2, j=3, k=4 など)

  do m = 100, 999
     i = m / 100                ! 百の位
     j = (m - i * 100) / 10     ! 十の位
     k = m - (i * 100 + j * 10) ! 一の位
     n = i**3 + j**3 + k**3     ! 3乗和
     if (m == n) then
        write(*,*) n
     end if
  end do

  stop
end program work0417

円の面積

半径 1 の円に内接する正 n 角形の面積 Sn, (n = 4, 8, 16, ..., 2**k, ...) を計算せよ。 |Sn+1 - Sn| の値が十分小さくなるまで（例えば 1.0d-5 など）計算すること。

!
! Sn : 半径 1 の円に内接する正 n 角形の面積 Sn
!      (n = 4, 8, 16, ..., 2**k, ...)
!
program work0417
  implicit none
  integer :: k, n
  real(8) :: sn, sold, pi, th
  real(8), parameter :: eps = 1.0d-5

  pi = acos(-1.0d+0)            ! 円周率をセット
  sold = 0                      ! ひとつ前の Sn
  k = 2                         ! k = 2, 3, 4, ...

  do
     n = 2**k
     th = 2.0d+0 * pi / n       ! (360度 / n) をラジアンで表す
     sn = sin(th) * n / 2       ! Sn を求める (半径は 1)
     write(*,*) k, n, sn
     if (abs(sold - sn) < eps) exit
     k = k + 1
     sold = sn
  end do

  stop
end program work0417

双子素数

1 から n までの間の双子素数（差が 2 であるような素数の組）を求めよ。 n はキーボードから入力できるようにすること。

!
! 双子素数
!
program work0419
  implicit none
  integer :: n, p, q, i
  logical :: fp, fq             ! p, q が素数かどうかを表す論理型変数（フラグ）
                                ! p が素数なら fp = .true., 素数でないなら fp = .false.

  ! n は 5 以上でないとダメ
  write(*,*) 'Input N (>= 5)'
  read(*,*) n

  ! p を小さい方、q を大きい方の素数候補とする
  do p = 3, n-2, 2
     q = p + 2

     ! p が素数であるかチェック
     fp = .true.                ! p は素数である、とする。
     do i = 3, p-2, 2
        ! p が i で割り切れたら p は素数ではない => fp = .false.
        if (mod(p, i) == 0) fp = .false.
     end do

     ! q が素数であるかチェック
     fq = .true.
     do i = 3, q-2, 2
        if (mod(q, i) == 0) fq = .false.
     end do

     ! p, q がともに素数（fp, fq がともに .true.）なら表示
     if (fp .AND. fq) then
        write(*,*) p, q
     end if
  end do

  stop
end program work0419