FortranのPureなFunction

# はじめに

## pure

　モダンなFortran90になって、言語仕様として`pure`属性が定義されました。

`pure function`は日本語で「純粋関数」です。

```math
   f(x) = y
```

という計算をするときに、純粋関数$f$は、

* 副作用がない：関係ない変数 $z$ に影響されないし、影響しない
* 参照透過性：$x$ の値が同じなら $y$ の値は必ず同じ

といった性質があります。ちなみに、プログラミング言語での「関数」はサブルーチン的な役割があるので純粋でないことも多いですが、数学で出てくる「関数」は普通は純粋関数じゃないかなと思います。

　副作用がない関数というのはプログラミングにおいて利点がいろいろあります。これについては関数型言語（LISPやHaskellなど）での議論が詳しく、あまり深入りするとFortranから離れてしまいますので、ここではFortranにおいての利点に注目します。
　Fortrannerにとっての純粋関数の利点は次の２つかと思います。

　* 関数の外の変数を意識しなくて良いのでデバッグがやりやすい
　* 関数の実行が独立なので並列化しやすい

　２つ目の利点は「`pure`は並列計算の時にさらに活きる属性」ということです。どういうことかという説明の前に、`elemental`について触れておきます。

## elemental

Fortranは配列を扱うのが得意な言語です。
特化しすぎてリストとか辞書型とか木構造を使うと辛い

このとき、配列の各成分に関数を適用して別の配列を作るというのはよくあると思いますが、そのたびに`do`ループを回すのは冗長です。そこでfortranでは`elemental`という属性が定義されています。
例えば三角関数の`sin`がありますが、これに配列を渡すことを考えます。

```fortran
double precision, allocatable :: v(:), w(:)
double precision, allocatable :: A(:, :), B(:, :)
! 中略
w = sin(v)
B = sin(A)
```

結果の`w`と`B`は、`v`や`A`と同じ形状の配列で、各成分ごとに`sin`を適用したものだと想像できます。Fortranの`elemental`関数は、上記の例のようにその直感をそのまま実現しています。
　なんか「`W = v.map(sin)`としたほうが直感的だ」という声が聞こえた気がしましたが、Fortranタグの前には無力です。[^1]

ここで注目したいこととして、各成分に`sin`関数を適用する際、`sin`は成分を超えて情報をやりとりしないということがあります。`sin(v(1))`の計算時に`sin(v(2))`を考えなくても大丈夫です。このため、`sin(v(1))`、`sin(v(2))`、・・・、`sin(v(i))`を独立に計算できます。これは並列計算にもってこいな性質です。

　逆に言うとある成分での計算が他の成分に依存する場合、色々な問題が発生します。なので、この`sin`のような`elemental`関数は引数だけに依存するように、つまり`pure`関数としなければなりません。(追記：最近のFortran2008ではimpureなelementalが可能になりました。記事の最後にまとめておきます。)
 例えば、配列の要素をどのような順番で適用するかとかが問題になりますが、順番というのは関数の外にあって、配列ループを回すときのインデックスと同値です。これでは「引数が決まれば結果が決まる」という純粋関数の性質を満たしません。また、コンパイラにどうやって順番を教えるかという問題もあります。

　別の利点として、上記の例のように引数の配列の形状を限定しません。もし`elemental`でない関数で上記のように複数の形状の配列に対応しようとすると、スカラー値用の`sin`、1次元配列用の`sin1d`、２次元配列用の`sin2d`、・・・と関数の定義がとんでもなく面倒になります。

　まとめますと、`pure`や`elemental`というのは配列の取り扱いや並列計算において重要なキーワードです。

この記事では主に`pure`な関数の振る舞いを紹介したいと思います。

# 環境

以下の環境で確認しました。

```
$ gfortran --version
GNU Fortran (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609
Copyright (C) 2015 Free Software Foundation, Inc.

$ uname -rosi
Linux 4.4.0-139-generic x86_64 GNU/Linux
```

申し訳ありませんがintel fortranは現在所持していないので、確認できません。

# 基本

## PUREの条件

ある関数が`pure`であるための条件は、導入の話を反映させると以下のようになります。

* 副作用がない
  - 引数が`intent(in)`である[^2]
  - 外部の変数への代入がない
  - ファイル・端末への出力がない
* 同じ引数なら必ず同じ結果を返す
  - 中で`random_number`のような呼ぶたびに結果が違う関数は呼べない
  - ファイル・端末からの`read`はできない

端末への出力が副作用だというのは、副作用についての議論に馴染みがないとわかりにくいかもしれませんが、ここでは詳しくは割愛します。

## 具体例

早速`pure`や`elemental`な関数の定義例を見てみます。

```fortran
module func
    implicit none
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        r = n*n
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        r = n*10 + m
    end function addmul
end module func
```

`square(n)`、`addmul(n, m)`は引数の`n`や`m`が決まれば結果が一意に決まり、また関数外の変数にも手を出しておらず、副作用もありません。また、`n`や`m`も`intent(in)`の効果で変更されません。`pure`としての資格は十分です。

これを以下のプログラムのように呼び出します。

```fortran
program test
    use func
    implicit none
    integer :: n, a(10), b(10), i

    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)
end program test
```

これはコンパイルして出力すると以下のようになります。

```plaintext:出力
         100
           9          18          27          36          45          54         63          72          81          90
```

　ひとつ目の`print`では`square(10)`の結果、つまり100が出力されます。挙動自体は普通の関数ですね。
　ふたつ目の`print`では、09, 18, 27, ..., 90 が出力されます。というのは`addmul`はひとつ目の引数を10倍して、ふたつ目の引数を足した整数を出力しますが、`a`は`0, 1, 2, ..., 9`で、`b`は`9, 8, 7, ..., 0`が代入されていますので、`addmul`の結果は`9, 18, 27, ..., 90`と９の倍数が並んだものになります。`elemental`なので配列を直接引数においてやると、同じ形の配列が返って来ました。

## バリエーション例

### elementalではpure宣言は無くてもよい

```fortran
    elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        r = n*10 + m
    end function addmul
```

`pure`がなくても`pure`としての性質は満たさないといけません。

### 文字列への変換のためのwrite

`addmul`について、結果を整数型ではなく文字列で返してみます。

```fortran
    elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        character(len=4) :: r
        write(r, '(i4)') n*10 + m
    end function addmul
```

この`write`文は副作用のない使用法なので、`pure`関数内で使えます。ただし、文字列長は固定です。

### 「値を変えない」 共有変数

たとえばモジュールの変数やCOMMONブロックでも、左辺に現れなければ使用可能です。

```fortran
module func
    implicit none
    integer :: i_local
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        r = i_local*n*n
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        integer :: c1, c2
        common /hoge/ c1, c2
        r = c1*n*10 + c2*m
    end function addmul
end module func
 
program test
 ! 中略
    common /hoge/ c1, c2
    i_local = 10
    c1 = 2
    c2 = 4
    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)
end program test
```

```:出力
        1000
          36          52          68          84         100         116         132         148         164         180
```

これは注意してやらないと参照透過性が保てないかもしれません。

## NG例

### 左辺に現れた共有変数

`pure`関数の中で外部の変数を変更することはできません。

```fortran
module func
    implicit none
    integer :: i_local
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        i_local = i_local + 1
        r = n*n
    end function square
```

```:コンパイルエラー
         i_local = i_local + 1
        1
Error: Variable ‘i_local’ can not appear in a variable definition context (assignment) at (1) in PURE procedure
```

`COMMON`ブロックの変数についても、同様のエラーが出ます。

### SAVE

```fortran
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        integer, save :: j = 0
        r = n*n
    end function square
```

```:コンパイルエラー
 integer, save :: j = 0
                     1
Error: SAVE attribute at (1) cannot be specified in a PURE procedure
```
使っていようといまいと、現れた時点でOUTです。

### Write文

IOは副作用の代表格です。

```fortran
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        write(*, *) n
        r = n*n
    end function square
```

```:コンパイルエラー
         write(*, *) n
                     1
Error: IO UNIT in WRITE statement at (1) must be an internal file in a PURE procedure
```

### INTENT忘れ

`INTENT(IN)`もしくは`VALUE`を厳密に要求します。

```fortran
    pure function square(n) result(r)
        integer :: n
        integer :: r
        r = n*n
    end function square
```

```:コンパイルエラー
     pure function square(n) result(r)
                          1
Error: Argument ‘n’ of pure function ‘square’ at (1) must be INTENT(IN) or VALUE
```

### 引数をターゲットとしたポインタの割当

基本的に`pure`関数の中でも`allocatable`な配列を使ったり`pointer`を割り当てたりが可能ですが、使えないケースというのもあり、例えば引数が関係するようなポインタ割当はエラーとなります。

```fortran
    pure function square(n) result(r)
        integer, intent(in), target :: n
        integer :: r
        integer,pointer :: pn
        pn => n
        r = pn*pn
    end function square
```

```:コンパイルエラー
         pn => n
              1
Error: Bad target in pointer assignment in PURE procedure at (1)
```

### 乱数

`pure`の中で乱数は使えません。

```fortran
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        double precision :: x
        call random_number(x)
        r = int(x*n*n)
    end function square
```

```:コンパイルエラー
call random_number(x)
                    1
Error: Subroutine call to intrinsic ‘random_number’ at (1) is not PURE
```

### elementalの形状間違い

任意の配列形状を受け取れるelementalですが、複数の引数をとるときにその次元が一致しない場合はどう計算してよいかわかりませんので、エラーとなります。

```fortran
    a = (/(i, i=0, 4)/)
    b = (/(i, i=9, 0, -1)/)
    print *, addmul(a, b)
```

```:コンパイルエラー
print *, addmul(a, b)
                       1
Error: Different shape for elemental procedure at (1) on dimension 1 (10 and 5)
```

### 他にも

* 中で呼び出される関数・サブルーチンはすべて`pure`でなければならない。
* プログラムのSTOP文を入れられない。
* ファイルの読み書きなどができない。
* `external`による関数渡しはできない。

などといった決まりがあります。まあ、純粋関数の性質を考えたら妥当な仕様かなと思います。

もし何らかのコンパイラで、こういったことができるものがあったとしても、`pure`の目的を考えると使うべきではないでしょう。

# デバッグ

`pure`はメリットが多い属性ですが、IO文が使えないため、printデバッグができなくなってしまうデメリットがあります。
すべての関数にPUREをつけるという遊びをしている時にいかにデバッグするかを考えてみます。

## プリプロセッサ

```fortran:main.f90
module func
    implicit none
contains 
#ifdef DEBUG
#define PUREF function
#else
#define PUREF pure function
#endif

    PUREF square(n) result(r)
        integer, intent(in) :: n
        integer :: r
#ifdef DEBUG
        write(*, *) "DEBUG"
        write(*, *) "input: ", n
#endif
        r = n*n
    end function square
end module func

program test
    use func
    implicit none
    integer :: n

    n = 10
    print *, square(n)
end program test
```

要はプリプロセッサのマクロで`DEBUG`が定義されているとデバッグモードとして適切にソースコードが切り替わるようにしています。

* `DEBUG`あり：`PUREF`は`function`と定義、関数内の`write`文を有効に
* `DEBUG`なし：`PUREF`は`pure function`と定義、関数内の`write`文を取り除く

```:デバッグありコンパイル
$ gfortran main.f90 -cpp -DDEBUG
$ ./a.out
 DEBUG
 input:           10
         100
```

```:デバッグなしコンパイル
$ gfortran main.f90 -cpp
$ ./a.out
         100
```

なお、この方法では`elemental`を外すとインターフェースもしくは呼び出し側も大幅に変えなければならない可能性があり、工夫が必要です。


## デバッガ

正直慣れるとこれが一番かと思います。
基本のプログラムについて、デバッグのデモ目的で中間変数を追加します。

```fortran:main.f90
module func
    implicit none
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        integer :: s
        s = n*n
        r = s
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        integer :: t10
        integer :: t1
        t10 = n*10
        t1 = m
        r = t10 + t1
    end function addmul
end module func

program test
    use func
    implicit none
    integer :: n, a(10), b(10), i

    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)

end program test
```

デバッガ情報を持たせてコンパイルします。

```
gfortran -g main.f90
```

デバッガから起動します。

```
% gdb ./a.out 
GNU gdb (Ubuntu 7.11.1-0ubuntu1~16.5) 7.11.1

(中略)

For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./a.out...done.
(gdb) run
Starting program: /path/to/a.out 
         100
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 29438) exited normally]
(gdb) q
```

`q`でデバッガが終了します。

さて、中間の変数`s`にちゃんと平方数が入っているか確認したいとなったとき。

```fortran
        s = n*n
```

これは8行目にありました。

```
(gdb) break 8
Breakpoint 1 at 0x4007c5: file main.f90, line 8.
```

これで8行目にブレークポイントを設定できました。実行してみます。

```
(gdb) run
Starting program: /path/to/a.out 

Breakpoint 1, func::square (n=10) at main.f90:8
8               s = n*n
```

8行目のソースを表示しつつ止まりました。ここで変数の内容を出力します。

```
(gdb) print n
$1 = 10
(gdb) print s
$2 = 0
```
引数`n`には正しく10が入っていますね。ですが、まだ`s`には代入されていないようです。1行進めてみます。それには`step`を用います。

```
(gdb) step
9               r = s
(gdb) print s
$3 = 100
```
8行目をすぎて9行目に到達し、`s`にちゃんと代入できていることが確認できました。見たいものが終わったのでプログラムを次のブレークポイントもしくはプログラム最後まで実行します。

```
(gdb) continue
Continuing.
         100
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 30346) exited normally]
```

次は`t10`にちゃんと値が入っているか確認しようと思います。

```fortran:17行目
        t10 = n*10
```

　ここでもうひとつ、`gdb`のコマンドは`gnuplot`のように１文字まで省略できます。それを示しながらやってみます。

　１７行目のブレークポイントの設定と実行。ブレークポイントより前にあるので、`square(10)`の出力が出ていますね。

```
(gdb) b 17
Breakpoint 1 at 0x400792: file main.f90, line 17.
(gdb) r
Starting program: /path/to/a.out 
         100
```

まずは一回目のブレークポイント到着。`addmul`の関数内で、引数もちゃんと表示されています。

```
Breakpoint 1, func::addmul (n=0, m=9) at main.f90:17
17              t10 = n*10
(gdb) p n
$1 = 0
(gdb) s
18              t1 = m
(gdb) s
19              r = t10 + t1
(gdb) p t10
$2 = 0
(gdb) c
Continuing.
```

`p`は`print`、`s`は`step`、`c`は`continue`と、最初の1字で大体OKです。
`t10`には`0`が入っていますが、これは初期化されたままなのか、正しく`0`が代入されたのかわかりません。なので別の引数で確認したいところ。

今回`addmul(n, m)`は`elemental`で、引数に長さ10の配列`a`、`b`を渡しており、内部的にはループとなりますので、17行目のブレークポイントは何度も通ります。ということで2回めのブレークポイントです。

```
Breakpoint 1, func::addmul (n=1, m=8) at main.f90:17
17              t10 = n*10
(gdb) p n
$5 = 1
(gdb) p t10
$6 = 0
(gdb) s
18              t1 = m
(gdb) p t10
$7 = 10
(gdb) p m
$8 = 8
(gdb) s
19              r = t10 + t1
(gdb) c
Continuing.
```

引数も配列の次の成分のものとなっていること、また`t10`にも正しく`10`が入っていることが確認できました。

3回目のブレークポイント。

```
Breakpoint 1, func::addmul (n=2, m=7) at main.f90:17
17              t10 = n*10
(gdb) p n
$9 = 2
(gdb) s
18              t1 = m
(gdb) s
19              r = t10 + t1
(gdb) s
20          end function addmul
(gdb) p r
$10 = 27
```
`end function`の行ではまだ関数の中のようですので、`r`を出力することでこの関数がどんな値を返しているかも見られます。

`continue`コマンドに整数を渡してやるとその数だけブレークポイントをスキップします。

```
(gdb) c 10
Will ignore next 9 crossings of breakpoint 1.  Continuing.
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 30419) exited normally]
```

これで最低限の変数の確認ができるかと思います。

ほかにも関数名で指定してステップ実行など、便利な機能がいろいろあります。たとえば今いる行の前後のプログラムをその場で確認できるコマンド`list`があります。

```
(gdb) list
12          pure elemental function addmul(n, m) result(r)
13              integer, intent(in) :: n, m
14              integer :: r
15              integer :: t10
16              integer :: t1
17              t10 = n*10
18              t1 = m
19              r = t10 + t1
20          end function addmul
21      end module func
```

ブレークポイントの設定についても、現在位置からの相対行数、関数名などでの設定などもあります。
http://flex.phys.tohoku.ac.jp/texi/gdb-j/gdb-j_18.html

FortranのみならずCやC++でも使えるコマンドですので、Linuxで生きている人なら覚えて損はないと思います。

計算化学の人みてるとprintデバッグばっかでデバッガなんて使っちゃいない


# Cとの連携

Fortran2003から`ISO_C_BINDING`が導入され、C言語との連携が強化されました。C言語で定義された関数に対しPUREが使えるか試してみます。

なお、先に言っておきますと`elemental`属性は`BIND(C)`と衝突しますので、使えません。

```
     pure elemental function addmul(n, m) result(r) bind(C, name="addmul")
                                  1
Error: BIND(C) attribute conflicts with ELEMENTAL attribute at (1)
```

## 普通に

```C:func.c
int square(int n){
    return n*n;
}
```

```fortran:main.f90
module func
    use iso_c_binding
    implicit none
    interface
        pure function square(n) result(r) bind(C, name="square")
            import :: c_int
            integer(c_int), value :: n
            integer(c_int) :: r
        end function square
    end interface
end module func

program test
    use func
    implicit none
    integer :: n

    n = 10
    print *, square(n)
end program test
```

```:コンパイルと実行
gcc -c func.c
gfortran main.f90 func.o
./a.out
```

```:出力
         100
```

普通のC関数の呼び出しと同じように使えますね。

## ルール違反

さて、ここでFortranの目が届かないところで良くないことをしてみます。例えばPURE関数内で端末出力。

```C:func.c
#include <stdio.h>

int square(int n){
    printf("inside C: %d\n", n);
    return n*n;
}
```

先と同様にコンパイルをしてみますと、何の問題もなくリンクできます。
実行してみます。

```:出力
inside C: 10
         100
```

あれ、pure関数が画面出力してる（棒）。

うまく使えば`pure`関数用の画面出力関数を作ってprintデバッグできるかもしれません。

# まとめ

PURE関数をうまく使うといろいろ良いことがあります。

# その他1

組み込みの数学関数は大体 `pure` です。

# その他2

Visual StudioのFortran対応きぼん。ソリューション管理とProfilerマジ便利。

# 追記：impureなelemental

Fortran2008になって、`impure`キーワードが導入されました。

https://software.intel.com/en-us/articles/fortran-2008-impure-support-in-intel-fortran-compiler

これを用いると、一部のimpureな操作がelementalでできるようになります。

```fortran
    impure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        write(*, *) n  ! コンパイルエラーが出ない！実行できる！
        r = n*10 + m
    end function addmul
```

考えてみれば、`elemental`と`pure`は関連が深いとはいえ同値ではないですね。

|pure |elemental  | 例  |
|---|---|---|
|○  |○  |各種数学関数  |
|○  |×  |`sum`などのreduce操作関数  |
|×  |○  |（理想的な）乱数 |
|×  |×  | 状態を保持する擬似乱数、I/O |

なおこの表では数学的な定義に基づいて考えていますので、Fortranのimpure elementalではまた別の判断基準（要は仕様書）に基づいて判断したらいいかなと思います。

@cure_honey さん、情報ありがとうございました。

[^1]: `v.map(sin)`において、`map`は高階関数で引数に関数をとっていて、`sin`関数を引数に渡していますが、Fortranは関数が第一級オブジェクトではないので、`map`のような高階関数は使えません。まあ関数ポインタとか使えばなんとかできるんですが、正直めんどい。しかも実数型配列`v`にメンバ関数を追加することもできないので、まあFortranではできない書きかたかなと。
[^2]: subroutineでは`intent(inout)`とか`intent(out)`が可能なようです。確かに、これができないと`pure subroutine`の存在意義が・・・