8. 副プログラム 目次へ
8.1. サブルーチン
8.1.1. 実引数について
8.1.2. 変数の有効範囲
8.1.3. 引数の入出力特性 [90]
8.2. 外部関数
8.2.1. 関数副プログラム
8.2.2. EXTERNAL宣言
8.2.3. 文関数
8.3. 内部手続き [90]
8.3.1. 内部手続きに特有の約束
8.4. 配列の扱い
8.4.1. 整合配列（形状明示仮配列）［（付）乱数の発生］
8.4.2. 形状引継ぎ配列(INTERFACE) [90]
8.4.3. 内部手続き中の配列 [90]
8.5. キーワード引数とOPTIONAL属性 [90]
8.6. 再帰的呼び出し [90]
まとまったルーチンワークを何度も使ったり，多数の作業の組合わせから構成される
複雑なプログラムを書くときには，サブルーチンや関数を副プログラムとして定義して
呼び出す形に書く方がプログラムの流れが理解しやすく，また，副プログラムを分業し
て作成したり，切り取って他のプログラムで利用したりすることもできる。
8.1. サブルーチン
例題 8_1 「３つの整数を読み込み，大きい順に並べ替えて出力する。」
整数をi,j,kとして
(1)i,jを比べて，jがiより大きければi,jの中身を入れ替える
(2)j,kを比べて，kがjより大きければj,kの中身を入れ替える
(3)もう一度i,jを比べる
とすればよい。３回繰返される「中身を入れ替える単純作業」をサブルーチンにする。
(プログラム例） 正規のプログラムは 8.1.2. へ
! --- とりあえずのプログラム ---
INTEGER :: i,j,k
READ*, i,j,k
IF(i < j) CALL swap(i,j)
IF(j < k) CALL swap(j,k)
IF(i < j) CALL swap(i,j)
PRINT *, i,j,k
END
!
SUBROUTINE swap(p,q)
INTEGER :: p,q,r
r=p; p=q; q=r
END
サブルーチンは，SUBROURTINE文で始まり，END文で終わる副プログラムで定義する。
SUBROUTINE サブルーチン名 [(仮引数のリスト)]
..........
END [ SUBROUTINE [サブルーチン名 ]]
主プログラム，あるいは別のプログラム単位から引用するときには
CALL サブルーチン名 [(実引数のリスト)]
とする。（注：引数のないサブルーチンも可能である。）副プログラムは，8.6.の
「再帰的呼び出し」を指定しない限り，自分自身を引用することはできない。
副プログラムのEND文の実行により，CALL文の次の文へ実行がもどる。何らかの
条件の成立により副プログラムの途中で元のプログラムへ復帰したい場合には，
RETURN文
RETURN
を用いる。
8.1.1. 実引数について
実引数には，仮引数と同じ型をもつ値の定まった変数や配列要素，配列（→8.4.）
および定数（またはPARAMETER属性を指定された定数名）が許される。ただし，定数の
場合には次の例のように「定数がすり替えられる」という予期せぬ結果を来すことも
あるので，できるだけ変数に代入してからサブルーチンにわたすようにした方がよい：
例
SUBROUTINE bai(x)
REAL :: x
x = 2.0 * x
END
!
CALL bai(1.0)
s = 1.0
PRINT*, s ! → なんと，2.000000 が出力されることがある。
END
このようなミステリアスな混乱を防止するために，[90]では引数の入出力特性指定
(8.1.3.)が新設された。この例では，少なくとも
s=1.0
CALL bai(s)
としておけば，定数1.0が２倍されて返されることはない。
8.1.2. 変数の有効範囲
プログラムで使われる変数名は，引数として使われている場合に引用によって結合
する以外には，他のすべての変数はプログラム単位内でのみ有効である。したがって
主プログラムと副プログラムで同じ名前が使われていても，（あとで述べる内部手続き
の場合を除いて）値は全く共有されない。プログラム単位間で共有したい場合には，
モジュール，または COMMON文 (説明省略）を用いる。
このことは，同一のプログラムが繰返し呼ばれる場合も同様である。副プログラムで
用いられているある変数に対して代入された値は，同じ副プログラムが次に呼び出され
たときに，前回の値が保存されているとは限らない。これを保証したい場合は，変数の
宣言文で「SAVE属性」を指定するか，SAVE宣言をする：
型, SAVE :: 変数名のリスト [90]
SAVE 型宣言ずみの変数名のリスト
文のラベルや文番号の有効範囲についても同様であるが，変数のように共有を指定する
方法はない。
8.1.3. 引数の入出力特性 [90]
引数による結合：変数は，CALL文で呼び出されたときにはじめて，仮引数と実引数が
リストに並べられた順に対応づけられ，値を共有する。上の例の CALL swap(i,j) では，
i,jの値がそれぞれp,qに送られ，値を入れ替えた後，pがiに，qがjに返される。
仮引数の入出力特性：仮引数には，副プログラムに対する入力に用いられるもの，
出力に用いられるもの，入出力両方に用いられるもの（上の例のp,q）の３種類がある
ため，副プログラムの中で変えられた値が元のプログラムに戻されていることに気が
つかないことが起こりうる（例 8.1.1.のミステリ）。この危険を避けるために変数の
型宣言文で，仮引数になっている変数の「INTENT属性」を指定する：
型，INTENT(IN, OUT, INOUT のいずれか) :: 変数のリスト [90]
この指定を行った場合，たとえば IN属性 が指定された変数の値を副プログラム中で
変えるような文があると，コンパイルエラーとなる。
[ex8_1.f90]
! --- 簡単なサブルーチンの例 ---
PROGRAM main
INTEGER :: i,j,k
READ*, i,j,k
IF(i < j) CALL swap(i,j)
IF(j < k) CALL swap(j,k)
IF(i < j) CALL swap(i,j)
PRINT *, i,j,k
END
! --- Swapping ---
SUBROUTINE swap(p, q)
INTEGER, INTENT(INOUT) :: p,q
INTEGER :: r
r=p; p=q; q=r
END SUBROUTINE
例題 8_2 「与えられた10点幅の人数分布を，ヒストグラムに描く。」
[ex8_2.f90]
! --- Drawing a Histogram ---
INTEGER :: nn(0:9), j
CHARACTER(LEN=80):: cn
CHARACTER(LEN=30):: fmt
nn = (/ 1, 5, 12, 22, 28, 35, 48, 28, 10, 3 /)
fmt='(1X, I2, "-", I2, 1X, A50)'
DO j = 0, 9
CALL str( nn(j), cn )
PRINT fmt, 10*j, 10*j+9, cn
END DO
END
! --- 文字列だけ作成する下請け工場 ---
SUBROUTINE str( n, c )
INTEGER, INTENT(IN) :: n
CHARACTER(LEN=80), INTENT(OUT) :: c
INTEGER :: k
c = ' '; c(1:1)='I'
DO k = 2, n+1
c(k:k) = 'o'
END DO
END
（実行）
0- 9 Io
10-19 Iooooo
20-29 Ioooooooooooo
30-39 Ioooooooooooooooooooooo
40-49 Ioooooooooooooooooooooooooooo
50-59 Iooooooooooooooooooooooooooooooooooo
60-69 Ioooooooooooooooooooooooooooooooooooooooooooooooo
70-79 Ioooooooooooooooooooooooooooo
80-89 Ioooooooooo
90-99 Iooo
注）サブルーチンの中の c=' ' を削除して実行してみよ。次に宣言文の中に c=' '
と書いて実行してみよ。SAVE指定がなくても，普通は前回に呼ばれたときの値が
保持される。また，宣言文やDATA文による初期値設定は最初に呼ばれたときにし
か有効ではなく，実行後は値がSAVEされてしまうので，これも注意がいる。
例：
! --- 小心！ 初期値設定 ---
SUBROUTINE count(n)
INTEGER :: sum = 0
DO i= 1, 10
sum = sum + i
END DO
n = sum
END
!
PROGRAM main
INTEGER :: n !（実行結果）
CALL count(n); PRINT*, n ! → 55
CALL count(n); PRINT*, n ! → 110：2回目のCALLではsum=0に設定されない。
END
この場合，宣言文で初期値設定ではなく，DO文の前で実行文として sum=0 を書いて
おけば，サブルーチンがCALLされるたびに実行されて0になるので問題は起こらない。
8.2. 外部関数
8.2.1. 関数副プログラム
自分で定義する関数は FUCTION文で始まり，END文で終わる副プログラムで定義する：
[型] FUNCTION 関数名(仮引数のリスト) [ RESULT(変数名)]
...................
END [FUNCTION [ 関数名]]
注）関数の型は，副プログラムの中で宣言してもよい：
型 :: 関数名
引数の有効範囲などの約束はサブルーチンの場合と同様である。また，関数の引用は
組込み関数の場合と同じである。
例題 8_3 「４バイトの正の整数の16進数表現を返す関数」
[ex8_3.f90]
! --- Generating a Hexadecimal Expression ---
FUNCTION hex(n)
CHARACTER(LEN=8) :: hex
CHARACTER(LEN=1) :: h(0:15) = (/ '0', '1', '2', '3', '4', '5', '6', &
'7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F' /)
INTEGER :: n, j, nn
hex = ' '
DO j = 8, 1, -1
nn = n/16
hex(j:j) = h(n-16*nn)
IF(nn == 0) EXIT
n = nn ! n に INTENT(IN) を指定してあるとエラー
END DO
END
! --- Main ---
CHARACTER(LEN=8) :: hex ! 注
INTEGER :: i
PRINT *, 'Input a positive Integer or negative one to stop:'
DO
READ *, i; IF( i < 0 ) EXIT
PRINT *, hex(i)
END DO
END
注）関数名については，これを引用する側でも型宣言が必要であることに注意。
この例のように，関数プログラム中で関数の値は，( )のない「関数名＝ 」という
形の代入文で与えられる。これが気持ち悪いときには，RESULTオプションで関数値を
与える変数を指定することができる。上の例では
FUNCTION hex(n) RESULT(hx) [90]
CHARACTER(LEN=8) :: hx
...............
hx(j:j) =
...............
となる。しかし，こうなるとわざわざ関数にすることもなく，サブルーチンでも大差は
ない。
SUBROUTINE(n, hx)
CHARACTER(LEN=8), INTENT(OUT) :: hx
...............
8.2.2. EXTERNAL宣言
副プログラムとして作った関数名が，組込み関数と同じであった場合には，副プログ
ラムとして定義された方が優先されるが，これを明示するには，呼び出す方のプログラ
ム中で EXTERNAL属性 または EXTERNAL文 を用いる：
型，EXTERNAL :: 関数名 [90]
EXTERNAL 関数名
これに対して，組込関数であることを明示するには，INTRINSIC属性 または INTRINSIC文
を用いる：
例 INTRINSIC SIN, COS
8.2.3. 文関数
一つの文で書けるような簡単な関数は「文関数」で定義できる。文関数は実行文では
ないので，プログラム単位中では全ての宣言文のあと，全ての実行文の前に書く。
例１ REAL :: f, x, s
REAL, PARAMETER :: a0=1.0, a1=2.0, a2=2.5, a3=1.5, a4=3.0, a5=1.0
f(x) = ((((a5*x + a4)*x + a3)*x + a2)*x + a1)*x + a0
READ*, s; PRINT*, f(s)
END
例２ REAL:: f, x, y
INTRINSIC SIN, COS
f(x, y) = 2.0*SIN(x)+ COS(y)
文関数定義文の右辺には，左辺の（ ）に書かれた仮引数以外の変数名（PARAMETER
属性を指定された定数名を除く）は書くことができない。仮引数は複数あってもよい。
また右辺には算術式のほか，組込み関数も許される。
注）例１の５次式の書き方に注意。すなおに
f(x) = a5*x**5 + a4*x**4 + a3*x**3 + a2*x**2 + a1*x + a0
と書く方がわかりやすいかもしれないが，これが24回の演算を要するのに対し，上の
書き方は10回の演算ですむ。また，高次の多項式になると上の方が誤差が少なくなる。
8.3. 内部手続き [90]
サブルーチンや関数は副プログラムとしてではなく，CONTAINS文 [90] で区切ること
により主プログラム（または別の副プログラム）中に含めてしまうことができる。
例題 8_4 「年月日を入れて曜日を知る。」
[ex8_4.f90]
! --- All Round Calendar ---
PROGRAM main
CHARACTER(LEN=2) :: name(0:6) = (/'日','月','火','水','木','金','土'/)
INTEGER :: md(1:12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
INTEGER :: y, m, d, ld
DO
PRINT *,'何年? (0 で終了)'
READ *, y; IF( y <= 0 ) EXIT
md(2) = 28 + leap(y) ! Add 1 for Leap.
PRINT *,'何月?'
DO
READ *, m ; IF( 0 < m .AND. m <= 12 ) EXIT
PRINT *, '12以下の正の数を入れ直して下さい。'
END DO
PRINT *,'Day?'
READ *, d
IF( d <= 0 .OR. d > md(m) ) THEN
PRINT*, 'そんな日は存在しません。'
ELSE
CALL count(ld)
PRINT "( 1X, I4, '年', I2, '月', I2, '日は ', A2, A/ )", &
y, m, d, name(ld), '曜日 です。'
END IF
END DO
!
CONTAINS ! ここからが内部手続き
! --- その日が西暦１年１月１日から何日目かを数え，曜日を算出する
SUBROUTINE count(n)
INTEGER, INTENT(OUT) :: n
n = yday(y) + SUM( md(1:m-1) ) + d
n = MOD(n, 7)
END SUBROUTINE
! --- 前年の12月31日が，西暦１年１月１日から何日目かを数える関数
FUNCTION yday(x) RESULT(days)
INTEGER :: x, days
days = 365*(x-1) + (x-1)/4 - (x-1)/100 + (x-1)/400
END FUNCTION yday
! --- うるう年なら１，平年なら０を返す関数
INTEGER FUNCTION leap(x)
INTEGER :: x
leap = yday(x+1) - yday(x) - 365
END FUNCTION leap
!
END PROGRAM
8.3.1. 内部手続きに特有な約束
1)内部サブルーチンや内部関数を呼び出すプログラム（親手続き）中で，サブルー
チン名や関数名の型宣言，およびEXTERNAL指定は不要（してはならない）。
2)親手続き中で使われている変数名，配列名などは内部手続き中でも有効で，共通
の値を持つ。
3)もし，同じ名前を内部手続き中で型宣言すると，その手続き中でのみ有効な独立
した変数とみなされ，共通の値は持たない。
4)内部手続きの中にCONTAINS文（すなわち内部手続き）を書くことはできない。
5)内部手続きから別の内部手続きを引用することは可能である。
8.4. 配列の扱い
自動割付配列 副プログラム中で使われている配列で，仮引数になっていない局所
的なものは，その寸法を仮引数になっている整変数（内部手続きの場合には親プログラ
ム中で使われている整変数でもよい）を用いて指定することができる。また，もちろん
ALLOCATABLE属性 を指定しておき， ALLOCATE文 で仮引数を含めて一般の整変数を用い
て寸法を指定することができる。
例 SUBROUTINE sub(x, n)
REAL :: x
INTEGER :: n, a(1:n), j
REAL, ALLOCATABLE :: temp(:)
.....
ALLOCATE( temp(0:a(j)) )
.....
配列引数 配列要素ではなく，配列そのものを引数にすることができる。ただし，
仮引数になっている配列を，ALLOCATABLE属性を指定して動的に大きさを割り付けるこ
とはできない。宣言文で普通どおりに定数の上下限を指定するか，あるいは場合に応じ
て大きさを変えたいときには，次の２通りのいずれかの方法による。
8.4.1. 整合配列（形状明示仮配列）
仮引数になっている配列は，やはり仮引数になっている整変数（内部手続きの場合に
は親プログラム中で使われている整変数でもよい）を用いて大きさ（寸法）の宣言をす
ることができる。配列も整変数も両方とも仮引数になっている場合に限る。 （文字型
仮引数の長さ指定にも類似の方法がある。→例題 9_6の注）
例題 8_5 「応募者数nと，当選者数mを与えて抽選をするサブルーチン」
現在(i-1)番目の当選者まで決まっているとして，残りの(n-i+1)枚のカードから
任意抽出し，これをi番目の当選者とする。
[ex8_5.f90]
! --- 抽選を行うサブルーチン ---
! --- n1 = 応募者数， n2= 当選者数， kk : 当選者番号 ---
SUBROUTINE chusen(kk, n1, n2)
INTEGER, INTENT(IN) :: n1, n2
INTEGER, INTENT(OUT) :: kk(n1)
INTEGER :: i, ir, j
REAL :: x
kk = (/ ( i, i = 1, n1) /)
PRINT*,'乱数発生のシード（なるべく大きい整数）を入れてください：'
READ*, ir
DO i = 1, MIN(n1 ,n2)
CALL ran(ir, x) ! 0 < x < 1 の一様乱数 x が返ってくる
j = i + INT( x*(n1 - i + 1) ) ! 残りのカードから任意抽出
CALL swap( kk(j), kk(i) ) ! これを i 番目としておく
END DO
END
!
! --- 入れ替え ---
SUBROUTINE swap(i, j)
INTEGER, INTENT(INOUT) :: i, j
INTEGER :: k
k = i; i = j; j = k
END
!
! --- 乱数発生 ---
SUBROUTINE ran(i, r)
INTEGER, INTENT(INOUT):: i
REAL, INTENT(OUT) :: r
INTEGER, PARAMETER :: mask = 2147483647, a = 48828125
i = IAND( a*i, mask )
r = REAL(i)/REAL(mask)
END
!
! --- 入出力 ---
PROGRAM main
INTEGER, ALLOCATABLE :: number(:)
INTEGER :: n, m, i
PRINT*, '応募者数は？'; READ*, n
PRINT*, '当選者数は？'; READ*, m
ALLOCATE( number(n) )
CALL chusen( number, n, m )
! --- 出力 ---
PRINT *,' 順位 番号'
DO i = 1, MIN(n, m)
PRINT '(2X, I4, I8 )', i, number(i)
END DO
END
==============================================================================
乱数の発生−−−乗算合同法
漸化公式
n(i) = A×n(i-1) mod M （A倍し，Mで割った剰り）
は，整数AとMおよび初期値(シード）n(0)を上手に選ぶことにより，実用に耐えうる
疑似整数乱数列を発生するのに用いることができる。この方法では，数列の周期は長
くてもM以下であるので，上の例ではMとして４バイト整数の最大値 231 を採用して
ある。（mask = 2147483647 は，２進数では 01111111 11111111 11111111 11111111
であり，ビット毎のAND演算により下31桁を取り出している。）係数Aとしては，最大
周期が実現可能な条件をみたす整数のうち，乱数列の性質が比較的よく調べられてい
る例であるA = 511 を採用してある。この整数乱数をMで割って規格化することにより，
(0,1) 間の実数乱数Rとする。これを (1,2,...K) の整数乱数Jにするには，(0,1) 間
を K等分 した区間に分け，J=1+INT(R*K) とすればよい。 シードは自分で入力する
か，システムクロックなどを利用すればよい。
なお，(0,1)間の一様乱数の発生には，組込みサブルーチンも用意されている：
CALL RANDOM_NUMBER(x)
================================================================================
8.4.2. 形状引継ぎ配列 [90]
外部手続きの仮引数になっている配列の大きさを具体的に指定しないで，引用による
結合が発生した段階で自動的に割り付ける方法である（形状引継ぎ配列）。この場合，
手続きを引用する側が仮引数の属性を知っておかなければならない。これをINTERFACE文
で指定する。内部手続きに対しては，この指定を行う必要はない。
例題 8_6「正方形行列のトレース（対角成分の和）を返す関数」
[ex8_6.f90]
! --- INTERFACEの例
PROGRAM main
INTERFACE
REAL FUNCTION trace(sqr)
REAL, INTENT(IN) :: sqr(:,:) ! 形状引継ぎ配列ですよ。
END FUNCTION
END INTERFACE
!
REAL, ALLOCATABLE :: a(:,:)
INTEGER :: n, i
PRINT*, '正方行列の行(または列）の数は？'
READ*, n
ALLOCATE( a(n,n) )
DO i = 1, n
PRINT '(A, I2, A)', '第', i, '行の要素？'
READ *, a(i, 1:n)
END DO
PRINT '(A, F7.4)', 'Trace = ', trace(a)
END
!
! --- Trace of a Square Matrix ---
FUNCTION trace(x) RESULT(tr)
REAL, INTENT(IN) :: x(:,:) ! ２次元配列であることだけ宣言する
REAL :: tr
INTEGER :: n1, n2, i
n1 = SIZE(x, 1); n2 = SIZE(x, 2)
IF( n1 /= n2 ) THEN
PRINT*,'正方行列ではありません。' ; RETURN
ELSE
tr = 0.0
DO i = 1, n1
tr = tr + x(i, i)
END DO
END IF
END FUNCTION
INTERFACE文は，このほかに「ユーザ定義演算」や「ユーザ定義代入文」を定義する
のに用いられる。（略）
8.4.3. 内部手続き中の配列 [90]
すでに説明したように，内部手続き中の配列は，局所的なものであれ仮引数になって
いるものであれ，親プログラム中で定義されている整数変数を用いて寸法を指定するこ
とができる。内部手続きを利用すると，以下のように配列を返す関数なども定義するこ
とができる。これは９章のモジュール副プログラム中で，内部手続きとしてモジュール
関数を作成するときに有用となる。
例題 8_7「２つの配列の対応する要素の大きい方を要素とする配列を返す
関数。」
[ex8_7.f90]
!----- Maximum Array
PROGRAM main
INTEGER :: n = 3
INTEGER, ALLOCATABLE :: a(:), b(:)
ALLOCATE( a(1:n), b(1:n) )
a = (/ 1, 2, 3 /); b = (/ 2, 1, 4 /)
PRINT *, maxa(a, b)
CONTAINS
FUNCTION maxa(x, y)
INTEGER :: maxa(1:n), x(1:n), y(1:n), i
DO i = 1, n
maxa(i) = MAX(x(i), y(i))
END DO
END FUNCTION maxa
END
この例は，もちろん親プログラムに配列演算として直接 MAX(a, b) と書いてもよい
のであるが，内部関数としてではなく，独立したプログラム単位の外部関数として定義
しようとすると困ったことが起きる例である。
8.5. キーワード引数とOPTIONAL指定 [90]
実引数は普通，仮引数に書かれた順に，書かれた数だけ書かねばならないが，引数の
名前を指定することで書く順序を変えたり，OPTIONAL属性を指定しておくことで場合に
応じて省略したりすることも可能である。ただし，外部手続きの引数にこの方法を適用
するためには，引用する側で仮引数名とその属性がわかっていなければならないから，
INTERFACE文を用いて教えておく必要がある。
例題 8_8 「点(x,y,z)の原点からの距離（の２乗）を返す。もし時間tも与えられ
た場合は，４次元時空距離 x**2 + y**2 + z**2 - t**2 を返す。」
[ex8_8.f90]
! --- 4D Space-Time ---
FUNCTION space_time(x, y, z, time) RESULT(ss)
REAL, INTENT(IN) :: x, y, z
REAL, INTENT(IN), OPTIONAL :: time
REAL :: ss
ss = x*x + y*y + z*z
IF( PRESENT(time) ) ss = ss - time*time
END FUNCTION
!
PROGRAM main
INTERFACE
REAL FUNCTION space_time(x, y, z, t)
REAL, INTENT(IN) :: x, y, z
REAL, INTENT(IN), OPTIONAL :: t
END FUNCTION
END INTERFACE
!
REAL :: a, b, c, t
a = 1.0; b= 2.0; c = 1.0
PRINT*, space_time(x=a, y=b, z=c) ! time は省略
t = 1.5
PRINT*, space_time(a, b, c, t)
END
実引数を全て仮引数どおりに並べるときは，通常の書き方でよい。順序を変えたり，
OPTIONAL属性をもつ引数を省略したりするときには，実引数を
（ INTERFACEで指定された仮引数名＝実引数，．．．）
の形で書く。なお，副プログラム中でOPTIONAL属性指定された引数が，実際にわたされたか
どうかを確認するには，引数問い合わせ関数 PRESENT を用いる：
PRESENT( OPTIONAL属性を持つ引数名 ）
該当する引数が実際に送られて来た場合には .TRUE.，そうでなければ .FALSE. となる。
8.6. 再帰的呼び出し [90]
副プログラムは，その内部で自分自身を引用することを許す場合には，RECURSIVE指定
をする：
RECURSIVE SUBROUTINE サブルーチン名[(...)] [90]
RECURSIVE [型] FUNCTION 関数名(...) RESULT(変数名) [90]
注）再帰関数の場合，必ず RESULT 形式にしなければならない。
例題 8_9 「次の漸化式で与えられるルジャンドル多項式のグラフを描く。
n P(n,x) = (2n-1) x P(n-1,x) - (n-1) P(n-2,x)
P(0,x)= 1, P(1,x) = x ( -1 <= x <= 1 ) 」
[ex8_9.f90]
! --- ルジャンドル多項式 ---
RECURSIVE FUNCTION P(n,x) RESULT(pnx) ! 必ず RESULT を使うこと
REAL :: pnx
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: x
IF(n == 0) THEN
pnx = 1.0
ELSE IF(n == 1) THEN
pnx = x
ELSE
pnx = (2.0-1.0/n) * x * P(n-1,x) - (1.0-1.0/n) * P(n-2,x)
END IF
END FUNCTION
!
PROGRAM main
INTEGER :: n, i
CHARACTER(LEN=1) :: c(-25:25)
REAL :: x
PRINT*,'何次の多項式にしますか？'
READ*,n
PRINT '( "Legendre Polynomial P(", I2, ", x)" )' , n
DO i = 0, 20
x = 0.05*i
y = p(n, x)
CALL presentation
PRINT '(F5.2, 2X, F10.6, 1X, 51A1)', x, y, c
END DO
CONTAINS
SUBROUTINE presentation ! 内部手続き故，変数を運ぶ必要なし
INTEGER :: j
IF(i /= 0) THEN
c = ' '
ELSE
c = (/ ('-', j = -25, 25) /)
END IF
c(-25) = 'I'; c(0) = 'I'; c(25) = 'I'
c(NINT(25*y)) = '*'
END SUBROUTINE
END
（出力）
何次の多項式にしますか？
6
Legendre Polynomial P( 6, x)
0.00 -0.312500 I----------------*-------I------------------------I
0.05 -0.296217 I * I I
0.10 -0.248829 I * I I
0.15 -0.174646 I * I I
0.20 -0.080576 I * I I
0.25 0.024277 I I* I
0.30 0.129181 I I * I
0.35 0.222511 I I * I
0.40 0.292636 I I * I
0.45 0.328981 I I * I
0.50 0.323242 I I * I
0.55 0.270766 I I * I
0.60 0.172096 I I * I
0.65 0.034675 I I* I
0.70 -0.125286 I * I I
0.75 -0.280777 I * I I
0.80 -0.391796 I * I I
0.85 -0.402996 I * I I
0.90 -0.241164 I * I I
0.95 0.187454 I I * I
1.00 1.000000 I I *
再帰手続きを使用する場合，必ず有限回の再帰で終了することに確信がないと使わ
ない方が無難である。そうでないと，あっという間に「スタックがあふれました。」
と警告される。上の例なら，繰り返し回数が予め与えられるから，DOループで書いて
も大差はないであろう。
［DOループで書いた例］
[ex8_9_1.f90]
REAL FUNCTION P(n,x)
REAL :: pnx(0:n)
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: x
INTEGER :: i
pnx(0) = 1.0 ; pnx(1) = x
DO i = 2, n
pnx(i) = (2.0-1.0/i)*x*pnx(i-1) - (1.0-1.0/i)*pnx(i-2)
END DO
P = pnx(n) ! 左辺には，引数のない関数名だけを書く
END FUNCTION
例題 8_10 「整数nを正整数の和に分ける組み合わせを全て求めよ。 6
あわせて，異なる正整数に分ける場合も求められるよう 5+1
にせよ。」 4+2
4+1+1
簡単そうに見えるが，母関数の方法以外には組み合わせ数に対する 3+3
一般公式が与えられておらず，数え上げるしかない問題である。右図 3+2+1
はn=6の例であるが，こうやって手で数えあげてみると，再帰的な手順 3+1+1+1
が含まれていることがすぐにわかる。これを多重のDOループで書くの 2+2+2
は不可能であろう。 2+2+1+1
2+1+1+1+1
1+1+1+1+1+1
[ex8_a.f90]
! --- Partition of Integer ---
INTEGER :: n, i, kind, id
INTEGER, ALLOCATABLE :: parts(:)
DO
PRINT*, " Select Case(0): N into arbitrary positive integers"
PRINT*, " Case(1): N into distinct positive integers"
PRINT*, " Default: Exit"
WRITE(*,'(A)', ADVANCE='no') "Case? "
READ *, id; IF(id < 0 .OR. id > 1) EXIT
WRITE(*,'(A)', ADVANCE='no') " n ? "; READ*, n
ALLOCATE( parts(n) ); i = 0; kind = 0; parts = 0; PRINT*
CALL Partition(n, n)
DEALLOCATE(parts)
PRINT "(' Problem (', I1,') for n =', I3 /)", id, n
END DO
!
CONTAINS
RECURSIVE SUBROUTINE Partition(p, px)
INTEGER, INTENT(IN) :: p, px
INTEGER :: k, ist
IF(p == 0) THEN
kind = kind + 1;
PRINT "( I6,') ', 80I3)", kind, parts(1:i)
ELSE
DO k = px, 1, -1
ist = i; i = i + 1; parts(i) = k
CALL Partition( p-k, min(p-k, k-id) )
i = ist
END DO
END IF
END SUBROUTINE Partition
END
整数pを分けるとしよう。まず先頭の数kを決め，残りのp-kを分ける分け方を調べれ
ばよい。ただし，構成整数は大きさの順に並べるとしていいから，残りの数の分け方の
うち先頭の数が「残りの数p-k以下」かつ「直前の数k（問題後半の場合はk-1）」以下
のものだけ考えればよい。これを残りが0となるまで繰り返せばよい。iは構成整数に
つける番号であるが，再帰手続きの最後（残りが0）まで行きついて途中まで戻った
ときに，どこから始めるかを覚えておかなければならないため，再帰手続きに行く前に
iの値をistという変数に一時保管してある。
例題 8_11 「ハノイの塔の手順を出力せよ。」
これは再帰手続きの例題の定番で「３本のポールＡ，Ｂ，Ｃがあって，最初はＡに
ｎ個の輪が，大きいものから順に積み上げられている。これを１個ずつ移動して山を
そっくりＢに移しなさい。ただし，どの時点でも小さい輪が大きい輪の下にあっては
いけません。」．．．プログラムを解読するより，クイズを解く方がやさしいかもし
れない。
| | |
(1) | |
(_2_) | |
(__3__) | |
(___4___) | |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Ａ Ｂ Ｃ
４つをＢに移すためには，まず３つをＣに移しておいてから，４をＢに移し，Ｃの
３つをＢに移せばよい。このためには３つを移す方法がわかればよい。３つを移す方法
を考えるためには，２つを移す方法がわかれば．．．。
[ex8_b.f90]
! --- ハノイの塔の解答 ---
INTEGER :: n = 0
CHARACTER(LEN=1) :: a, b, c
a = 'A'; b = 'B'; c = 'C'
DO WHILE( n <= 0 )
PRINT*, "最初の輪の個数は？"
READ*, n
END DO
CALL hanoi(n, a, b, c)
END
!
RECURSIVE SUBROUTINE hanoi(i, x, y, z)
INTEGER, INTENT(IN) :: i
CHARACTER(LEN=1), INTENT(IN) :: x, y, z
IF( i > 0 ) THEN
CALL hanoi(i-1, x, z, y)
PRINT '("第", I2, "番目の輪を ", A1," から ", A1, &
" へ移動する")', i, x, y
CALL hanoi(i-1, z, y, x)
END IF
END
⇒８章先頭へ ⇒９章へ