Fortranのanyとallで条件判定を簡略化する

# 概要
`all`と`any`関数を用いることで，複数の条件の真偽を簡潔に考慮できます．例えば，配列の形状や値の等価性を簡単に判別できるようになります．

# 環境
- Windows 10 64 bit
- gfortran 10.3
- Intel Fortran classic 2021.5.0
- NAG Fortran 7.1.0
- fypp 3.1

# Fortranの関係演算子
`all`関数，`any`関数の使い方を説明する前に，Fortranにおける関係演算子（`==`, `/=`, `<`, `<=`, `>`, `>=`）について言及しておきます．

関係演算子は，演算子右と左の値を比較して，演算子が示す関係を満たしていれば`.true.`，そうでなければ`.false.`を返します．

```Fortran
print *, (1 == 2), (1 /= 2), (1 < 2), (1 > 2) ! F T T F
```

言語によっては，関係演算子の左右にはスカラしか置けない場合もありますが，Fortranは配列を置くことができます．その場合，関係演算子の戻値も配列になります．

```Fortran
integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))     ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b)) ! [[6, 5], [4, 3], [2, 1]]
print *, a < b ! [[T, T], [T, F], [F, F]]
```

配列要素の関係を，関係演算子を用いてまとめて評価できるのは非常に便利なのですが，その結果に基づいて`if`分岐をしようとすると，問題が生じます．`if`文の条件式は，スカラの論理値にしか対応していないからです．

```Fortran
program main
    use, intrinsic :: iso_fortran_env
    implicit none

    integer(int32) :: a(2, 3), b(2, 3), i
    a = reshape([(i, i=1, product(shape(a)))], shape(a))
    b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b))
    if (a < b) print *, "a is smaller than b"
end program main
```

このプログラムをコンパイルしようとすると，エラーが発生します．
```console
    8 |     if (a < b) print *, "a is smaller than b"
      |        1
Error: IF clause at (1) requires a scalar LOGICAL expression
```

こういう場合に，`any`関数や`all`関数が役に立ちます．

# `all`関数, `any`関数
Fortranにおける`all`関数，`any`関数は，配列のreductionを行う関数の一つです．

|関数|戻値|機能|
|:--|:--|:--|
|`all(mask [, dim])`|`logiacl`|論理型配列`mask`が全要素が`.true.`の場合，真を返す．`dim`で検査する配列の次元を指定できる．|
|`any(mask [, dim])`|`logiacl`|論理型配列`mask`のいずれかの要素が`.true.`の場合，真を返す．`dim`で検査する配列の次元を指定できる．|

インタフェースとしては非常に簡単で，論理型配列を渡せば，その内容に応じて`.true.`または`.false.`が返ってきます．

前節の例のように，関係演算子による比較結果に基づいて`if`文で判別を行う場合に，比較結果の論理型配列からさらに一つの結果（全要素が小さい，全要素が等しい等）を得る場合に利用します．

```Fortran
integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b))
if (all(a < b)) print *, "a is smaller than b" ! aはbより大きい要素も小さい要素も持っているのでこれは表示されない
```

配列の全要素の値を比較する場合は，下記のように書けるでしょう．

```Fortran
integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))     ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b)) ! [[6, 5], [4, 3], [2, 1]]

if (all(a == b)) then
    print *, "a is equal to b"
else
    print *, "a is not equal to b."
    print *, "maximum absolute difference = ", maxval(abs(a - b)), "at ", maxloc(abs(a - b))
    print *, "minimum absolute difference = ", minval(abs(a - b)), "at ", minloc(abs(a - b))
    ! abs(a - b) == [[5, 3], [1, 1], [3, 5]]
end if
```
実行結果は下記のようになります．
```console
 a is not equal to b.
 maximum absolute difference =            5 at            1           1
 minimum absolute difference =            1 at            1           2
```

配列要素のどれか一つでも一致していればよい場合については，下記のような例を挙げることができます．

```Fortran
integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))      ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=1, product(shape(b))*2, 2)], shape(b)) ! [[1, 3], [5, 7], [9, 11]]

if (any(a == b)) then
    print *, "a has", count((a - b) == 0), " value(s) equal to b"
else
    print *, "a does not have the same value as b."
end if
```

実行結果は下記のようになります．
```console
 a has           1  value(s) equal to b
```

## その他のreductionを行う関数
既にいくつか現れていますが，`all`, `any`以外に，配列のreductionを行う関数を列挙しておきます．

|関数|戻値|機能|
|:--|:--|:--|
|`count(mask [, dim])`|`integer`|論理型配列`mask`における`.true.`の要素数を返す．`dim`で検査する配列の次元を指定できる．|
|`sum(array, [, dim] [, mask])`|`array`と同じ型|`array`の各要素の和を返す．`dim`で検査する配列の次元を指定できる．`mask`で計算対象となる要素を指定できる．|
|`product(array, [, dim] [, mask])`|`array`と同じ型|`array`の各要素の積を返す．`dim`で検査する配列の次元を指定できる．`mask`で計算対象となる要素を指定できる．|
|`maxval(array [, dim] [, mask])`|`array`と同じ型（整数型もしくは実数型）|配列`array`内の最大値を返す．`dim`で検査する配列の次元を指定できる．`mask`で検査する要素を指定できる．|
|`minval(array [, dim] [, mask])`|`array`と同じ型（整数型もしくは実数型）|配列`array`内の最小値を返す．`dim`で検査する配列の次元を指定できる．`mask`で検査する要素を指定できる．|
|`maxloc(array [, dim] [, mask] [, kind])`|整数型配列（次元は`array`と同じ．型種別は`kind`で指定）|配列`array`内で最初に最大値が現れる配列添字を返す．`dim`で検査する配列の次元を指定できる．`mask`で検査する要素を指定できる．|
|`minloc(array [, dim] [, mask] [, kind])`|整数型配列（次元は`array`と同じ．型種別は`kind`で指定）|配列`array`内で最初に最小値が現れる配列添字を返す．`dim`で検査する配列の次元を指定できる．`mask`で検査する要素を指定できる．|
|`findloc(array, value [, dim] [, mask] [, kind] [, back])`|整数型配列（次元は`array`と同じ．型種別は`kind`で指定）|配列`array`内で最初に`value`が現れる配列添字を返す．`dim`で検査する配列の次元を指定できる．`mask`で検査する要素を指定できる．`back=.true.`にすると，配列を後ろから検査する．|

# `all`と`any`の使いどころ
## 配列の値の等価性比較
既に例示していますが，二つの配列の等価性を比較する際に，`all`や`any`が活用できます．

二つの配列`a`と`b`の全要素が同じであることを判定したいのであれば，

```Fortran
all(a==b)
```
と書けばよく，どれか一つの要素でも値が等しくないことを判定したいのであれば，

```Fortran
any(a/=b)
```
と書けるでしょう．

実数の場合は，`==`で等価性を比較してはいけないと教えられたと思います．（私も教えられました）
実数型配列の等価性の比較は，二つの配列の差が，全要素においてあるしきい値（一般的にはマシンイプシロン）以下であれば等しいと判別することが多いと思います．
その場合でも，`all`を使えば

```Fortran
all(abs(a - b) < epsilon(a))
```

といった書き方ができます．

ユーザ定義演算子を用いれば，実数に対する $\approx$ 演算子が定義できます．

```Fortran
    interface operator(.approx.)
        procedure :: approximately_equal_real32
        procedure :: approximately_equal_real32_rank1
    end interface
contains
    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a
        real(real32), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real32

    pure function approximately_equal_real32_rank1(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: approx_eq

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank1
```

## 配列の形状チェック
二つの配列`a`と`b`を比較する場合，前提として配列の形状が同じでなければなりません．そのため，処理を安全側に振るのであれば，配列の値を比較する前に，配列の形状が一致することを確認しておく方がよいでしょう．

その場合も，`all`を使えば簡単です．

```Fortran
all(shape(a)==shape(b))
```

`shape`手続は，引数の配列形状を配列で返すので，`all`を使えば配列の各次元の要素数を個別に比較できます．

```Fortran
integer(int32) :: a(2, 3), b(3, 2)
print *, shape(a) ! [2, 3]
print *, shape(b) ! [3, 2]
print *, all(shape(a) == shape(b))
```
```console
           2           3
           3           2
 F
```

配列の次元が異なる場合はコンパイルエラーが出て比較自体ができないので，注意してください．

## 複数の条件が`.and.`あるいは`.or.`で評価されている場合の論理演算の簡略化
色々な条件が関係する場合，`if`文で判定しようとすると，複数の条件を`.and.`または`.or.`で評価する必要がでてきて，条件式が非常に長くなり，見通しも悪くなります．

全ての条件が`.and.`で評価されている場合，あるいは`.or.`で評価されている場合は，`all`あるいは`any`を用いて簡略化できます．

`all`や`any`は，引数に論理型配列を渡せばよいので，配列構成子`[]`と組み合わせれば簡単に複数の条件を考慮できます．

例えば，二つの配列の添字上下限を比較する場合，`all`を使わないと，条件がかなり長くなります．

```Fortran
if (lbound(a, dim=1) == lbound(b, dim=1) .and. &
    lbound(a, dim=2) == lbound(b, dim=2) .and. &
    ubound(a, dim=1) == ubound(b, dim=1) .and. &
    ubound(a, dim=2) == ubound(b, dim=2)) then
    ! 添字の上下限が全て一致した場合の処理
end if
```

これは2次元配列の例ですが，Fortran規格では15次元配列（Coarrayを使う場合は15-corank）まで認められているので，最大で30個の条件を`.and.`で評価しなければなりません．

`all`を用いれば，非常に簡潔になります．

```Fortran
if (all([lbound(a) == lbound(b), &
         ubound(a) == ubound(b)])) then
    ! 添字の上下限が全て一致した場合の処理
end if
```

複数の条件が`.or.`で評価されている場合は，`any`を用います．

ただし，簡単になるからと全ての条件を`all`や`any`でまとめればよいかというと，そうは考えていません．条件の意味のまとまりごとに考えて，`all`でまとめる場合と`.and.`を用いる場合を使い分けた方がよいでしょう．

## 配列内に`NaN`が存在しているかの確認
Fortranの面白い機能に`elemental`属性があります．`elemental`属性を付与された手続は，定義上はスカラ値のみを引数に取りますが，配列を受け取れるようになります．配列の全要素に対して，手続内で定義された演算を実行し，配列を返します．

例えば，組込手続の`sin`はelemental functionに分類されており，リファレンス上は実数（あるいは複素数）のスカラ値を引数に取るようになっていますが，配列を渡すことも可能です．

```Fortran
real(real32), parameter :: PI = acos(-1.)
print *, sin(0.), sin(PI)               ! 0.0000000  -8.7422777E-08
print *, sin([0., 0.5, 1., 1.5, 2.]*PI) ! 0.0000000   1.0000000  -8.7422777E-08  -1.0000000   1.7484555E-07
```

実数が`NaN`かどうかをFortran標準に従って判別するには，`ieee_arithmetic`モジュールで定義されている`ieee_is_nan`関数を用います．`ieee_is_nan`もelemental functionなので，配列を渡すと全要素について`NaN`かどうかを調べ，結果を配列で返します．

返ってきた配列のなかに一つでも`.true.`があれば`NaN`が生じているということですがら，`any`関数を用いてやれば簡単に検証できます．

```Fortran
use, intrinsic :: ieee_arithmetic
real(real32) :: a(10, 10)
a(:, :) = 0d0

a(2, 3) = ieee_value(0._real32, ieee_quiet_nan)
if (any(ieee_is_nan(a))) then
    print *, "a has ", count(ieee_is_nan(a)), " NaN at ", findloc(ieee_is_nan(a), .true.)
end if
```

```console
 a has            1  NaN at            2           3
```

下記のような，配列を引数にとり，その中に一つでも`NaN`があれば`.true.`を返すような手続があると，数値計算では役に立つかも知れません．

```Fortran
    interface has_nan
        procedure :: has_nan_real32_rank1
        procedure :: has_nan_real32_rank2
    end interface
contains
    function has_nan_real32_rank1(array) result(has_nan)
        implicit none
        real(real32), intent(in) :: array(:)
        logical :: has_nan

        has_nan = any(ieee_is_nan(array))
    end function has_nan_real32_rank1

    function has_nan_real32_rank2(array) result(has_nan)
        implicit none
        real(real32), intent(in) :: array(:, :)
        logical :: has_nan

        has_nan = any(ieee_is_nan(array))
    end function has_nan_real32_rank2
```

# `.approx.`演算子および`has_nan`手続のためのメタプログラミング
`all`と`any`の使い方をいくつか紹介してきました．その中で，演算子や関数をいくつか示しました．引数は，実数のスカラ値か配列を取ります．Fortranでは，異なる引数を取る同じような手続を，総称名を利用してまとめることができます．しかし，個別の引数の型に対する手続は書く必要があります．

実行する処理自体はFortranの機能によって，配列次元を問わず全く同一に書くことができます．一方で，Fortranにはメタプログラミングの機能が無いため，同じ処理で引数の異なる手続を何個も書く必要があります．

配列の形状に関しては，Fortran 2018から加わったassumed-rank（次元引継？）仮引数の機能である程度簡略化は可能です．assumed-rank仮引数を使うには仮引数の配列形状に`..`を用います．

```Fortran
real(real32), intent(in) :: a(..)
```

実引数の形状に応じた処理を記述するには，`select rank`文を用います．
```Fortran
select rank (a)
rank (0) ! 実引数がスカラの場合
    has_nan = ieee_is_nan(a)
end select
```

これは`select type`と類似で，なかなか融通が利きません．これまで上で見てきたように，`ieee_is_nan`を用いて配列内の`NaN`を検出するには，配列の次元がいくつであったとしても，`any(ieee_is_nan(a))`と書けました．assumed-rankの機能を利用する場合，`rank`を一つずつ書いていく必要があります．

```Fortran
select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1)
    has_nan = any(ieee_is_nan(a))
rank (2)
    has_nan = any(ieee_is_nan(a))
rank (3)
    has_nan = any(ieee_is_nan(a))
rank (4)
    has_nan = any(ieee_is_nan(a))
rank (5)
    has_nan = any(ieee_is_nan(a))
rank (6)
    has_nan = any(ieee_is_nan(a))
rank (7)
    has_nan = any(ieee_is_nan(a))
rank (8)
    has_nan = any(ieee_is_nan(a))
rank (9)
    has_nan = any(ieee_is_nan(a))
rank (10)
    has_nan = any(ieee_is_nan(a))
rank (11)
    has_nan = any(ieee_is_nan(a))
rank (12)
    has_nan = any(ieee_is_nan(a))
rank (13)
    has_nan = any(ieee_is_nan(a))
rank (14)
    has_nan = any(ieee_is_nan(a))
rank (15)
    has_nan = any(ieee_is_nan(a))
end select
```

おそらくコンパイル時に配列次元が静的に決まっていないといけないことが理由だと思いますが，もう少し簡潔に

```Fortran
select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1,2,3,4,...)
    has_nan = any(ieee_is_nan(a))
end select
```
とか
```Fortran
select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1:15)
    has_nan = any(ieee_is_nan(a))
end select
```
とか
```Fortran
select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank default
    has_nan = any(ieee_is_nan(a))
end select
```
とかが通ってくれればよいのですが，前二つは規格上対応しておらず，後者はコンパイラによっては`ieee_is_nan`がassumed-rank仮引数に対応していないようです．Intel Fortranではビルドができましたが，gfortranやNAG Fortranでは無理でした．

今までは，関数インタフェースも含めて同じ処理を何個も書かなければいけなかったことを思うと，大分便利になったと思います．

```Fortran
    pure function has_nan_real64_rank1(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(:)
        logical :: has_nan
        has_nan = any(ieee_is_nan(a))
    end function has_nan_real64_rank1
    pure function has_nan_real64_rank2(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(:, :)
        logical :: has_nan
        has_nan = any(ieee_is_nan(a))
    end function has_nan_real64_rank2
```

とはいえ，いちいち繰り返すのも大変なので，以前紹介した[fypp](https://github.com/aradi/fypp)を用います．

fyppは，Pythonによるプリプロセッサであり，（他の言語でも使えるものの，主に）Fortranのプリプロセッサとして開発されました．Pythonを使ったことがある人には，非常に馴染みやすいと思います．

変数の定義や処理については[以前紹介した記事](https://qiita.com/implicit_none/items/47680973fb98d85d6cbb)を参考にしてもらうとして，fyppを利用すると，assumed-rankを用いた`NaN`の検出手続が下記のように書けます．

```Fortran
    public :: has_nan

    interface has_nan
    #:for real_kind in REAL_KINDS
        procedure :: has_nan_${real_kind}$
    #:endfor
    end interface
contains
    #:for real_kind in REAL_KINDS
    pure function has_nan_${real_kind}$(a) result(has_nan)
        implicit none
        real(${real_kind}$), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        #:for rank in RANKS
        rank (${rank}$)
            has_nan = any(ieee_is_nan(a))
        #:endfor
        end select
    end function has_nan_${real_kind}$

    #:endfor
```

この中で出ている`REAL_KINDS`や`RANKS`はfyppの機能で定義できる変数であり，ソースの先頭で下記のように定義しています．

```python
#! nubmer of bits of real type
#:set REAL_BITS = [32, 64]

#:if defined('real128')
#:set REAL_BITS = REAL_BITS + [128]
#:endif

#! real kinds
#:set REAL_KINDS = ["real{}".format(bits) for bits in REAL_BITS]

#! target array rank
#:set MAX_RANK = 15

#! rank
#:set RANKS = [rank for rank in range(1, MAX_RANK+1)]
```

これをfyppで処理すると，下記のようなFortranのコードが生成されます．生成されたコードに適切なモジュール名をつけてやれば完成です．

<details>
<summary>生成されたコード</summary>

```Fortran
    public :: has_nan

    interface has_nan
        procedure :: has_nan_real32
        procedure :: has_nan_real64
    end interface

contains
    pure function has_nan_real32(a) result(has_nan)
        implicit none
        real(real32), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        rank (1)
            has_nan = any(ieee_is_nan(a))
        rank (2)
            has_nan = any(ieee_is_nan(a))
        rank (3)
            has_nan = any(ieee_is_nan(a))
        rank (4)
            has_nan = any(ieee_is_nan(a))
        rank (5)
            has_nan = any(ieee_is_nan(a))
        rank (6)
            has_nan = any(ieee_is_nan(a))
        rank (7)
            has_nan = any(ieee_is_nan(a))
        rank (8)
            has_nan = any(ieee_is_nan(a))
        rank (9)
            has_nan = any(ieee_is_nan(a))
        rank (10)
            has_nan = any(ieee_is_nan(a))
        rank (11)
            has_nan = any(ieee_is_nan(a))
        rank (12)
            has_nan = any(ieee_is_nan(a))
        rank (13)
            has_nan = any(ieee_is_nan(a))
        rank (14)
            has_nan = any(ieee_is_nan(a))
        rank (15)
            has_nan = any(ieee_is_nan(a))
        end select
    end function has_nan_real32

    pure function has_nan_real64(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        rank (1)
            has_nan = any(ieee_is_nan(a))
        rank (2)
            has_nan = any(ieee_is_nan(a))
        rank (3)
            has_nan = any(ieee_is_nan(a))
        rank (4)
            has_nan = any(ieee_is_nan(a))
        rank (5)
            has_nan = any(ieee_is_nan(a))
        rank (6)
            has_nan = any(ieee_is_nan(a))
        rank (7)
            has_nan = any(ieee_is_nan(a))
        rank (8)
            has_nan = any(ieee_is_nan(a))
        rank (9)
            has_nan = any(ieee_is_nan(a))
        rank (10)
            has_nan = any(ieee_is_nan(a))
        rank (11)
            has_nan = any(ieee_is_nan(a))
        rank (12)
            has_nan = any(ieee_is_nan(a))
        rank (13)
            has_nan = any(ieee_is_nan(a))
        rank (14)
            has_nan = any(ieee_is_nan(a))
        rank (15)
            has_nan = any(ieee_is_nan(a))
        end select
    end function has_nan_real64
```
</details>

`.approx.`演算子の実装としていた`approximately_equal_real32_rank1`等の手続については，assume-rank仮引数を用いるのは現状では止めた方がよいでしょう．
二つの仮引数のrankを同時に調べられないので，`a`のrankを調べた後`b`のrankを調べる必要があり，あまりスマートではありません．

```Fortran
    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(..)
        real(real32), intent(in) :: b(..)
        logical :: approx_eq

        select rank (a)
        rank (0); select rank (b); rank (0)
                approx_eq = abs(a - b) <= epsilon(0.); end select
        rank (1); select rank (b); rank (1)
                approx_eq = all(abs(a - b) <= epsilon(a)); end select
...
```

また，`select rank`内において，総称名でまとめられた手続を呼ぶと，引数の配列次元を正しく認識できず，コンパイルエラーになる事例も確認できました．そのため，assumed-rankの使用について，下記のように自身の基準を決めました．

- assumed-rank仮引数を用いてもよい
    - 引数が一つの手続
- assumed-rank仮引数は用いない方がよい
    - 引数が二つ以上
    - 引数が一つ以上で，`select rank`内で総称名でまとめられた手続を呼び出している

そのため，メタプログラミングでコードを生成する場合は，仮引数の配列次元が異なる手続を生成するようにします．

```Fortran
    public :: operator(.approx.)
    public :: are_same_shape

    interface operator(.approx.)
    #:for real_kind in REAL_KINDS
        procedure :: approximately_equal_${real_kind}$
        #:for rank in RANKS
        procedure :: approximately_equal_${real_kind}$_rank${rank}$
        #:endfor
    #:endfor
    end interface

    interface are_same_shape
    #:for real_kind in REAL_KINDS
    #:for rank in RANKS
        procedure :: are_same_shape_${real_kind}$_rank${rank}$
    #:endfor
    #:endfor
    end interface

contains
    #:for real_kind in REAL_KINDS
    pure function approximately_equal_${real_kind}$(a, b) result(approx_eq)
        implicit none
        real(${real_kind}$), intent(in) :: a
        real(${real_kind}$), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a-b)<=epsilon(a)
    end function approximately_equal_${real_kind}$

        #:for rank in RANKS
        pure function approximately_equal_${real_kind}$_rank${rank}$(a, b) result(approx_eq)
            implicit none
            real(${real_kind}$), intent(in) :: a${ranksuffix(rank)}$
            real(${real_kind}$), intent(in) :: b${ranksuffix(rank)}$
            logical :: approx_eq

            if(.not.are_same_shape(a, b))then
                approx_eq = .false.
                return
            end if

            approx_eq = all(abs(a-b)<=epsilon(a))
        end function approximately_equal_${real_kind}$_rank${rank}$

        #:endfor

        #:for rank in RANKS
        pure function are_same_shape_${real_kind}$_rank${rank}$(a, b) result(are_same_shape)
            implicit none
            real(${real_kind}$), intent(in) :: a${ranksuffix(rank)}$
            real(${real_kind}$), intent(in) :: b${ranksuffix(rank)}$
            logical :: are_same_shape

            are_same_shape = all(shape(a) == shape(b))
        end function are_same_shape_${real_kind}$_rank${rank}$

        #:endfor

    #:endfor
```

同様にfyppで処理すると，Fortranのコードが生成されます．異様に長くなるので，`MAX_RANK`を3として生成したコードを示します．

<details>
<summary>生成されたコード</summary>

```Fortran
    public :: are_same_shape
    public :: has_nan

    interface operator(.approx.)
        procedure :: approximately_equal_real32
        procedure :: approximately_equal_real32_rank1
        procedure :: approximately_equal_real32_rank2
        procedure :: approximately_equal_real32_rank3
        procedure :: approximately_equal_real64
        procedure :: approximately_equal_real64_rank1
        procedure :: approximately_equal_real64_rank2
        procedure :: approximately_equal_real64_rank3
    end interface

    interface are_same_shape
        procedure :: are_same_shape_real32_rank1
        procedure :: are_same_shape_real32_rank2
        procedure :: are_same_shape_real32_rank3
        procedure :: are_same_shape_real64_rank1
        procedure :: are_same_shape_real64_rank2
        procedure :: are_same_shape_real64_rank3
    end interface

contains
    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a
        real(real32), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real32

    pure function approximately_equal_real32_rank1(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank1

    pure function approximately_equal_real32_rank2(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:, :)
        real(real32), intent(in) :: b(:, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank2

    pure function approximately_equal_real32_rank3(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:, :, :)
        real(real32), intent(in) :: b(:, :, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank3

    pure function are_same_shape_real32_rank1(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank1

    pure function are_same_shape_real32_rank2(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:, :)
        real(real32), intent(in) :: b(:, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank2

    pure function are_same_shape_real32_rank3(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:, :, :)
        real(real32), intent(in) :: b(:, :, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank3

    pure function approximately_equal_real64(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a
        real(real64), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real64

    pure function approximately_equal_real64_rank1(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:)
        real(real64), intent(in) :: b(:)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank1

    pure function approximately_equal_real64_rank2(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:, :)
        real(real64), intent(in) :: b(:, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank2

    pure function approximately_equal_real64_rank3(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:, :, :)
        real(real64), intent(in) :: b(:, :, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank3

    pure function are_same_shape_real64_rank1(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:)
        real(real64), intent(in) :: b(:)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank1

    pure function are_same_shape_real64_rank2(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:, :)
        real(real64), intent(in) :: b(:, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank2

    pure function are_same_shape_real64_rank3(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:, :, :)
        real(real64), intent(in) :: b(:, :, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank3
```
</details>

# まとめ
非常に便利なのにあまり使われていない`all`と`any`の使い方について説明しました．

特に配列を多用して数値計算を行うFortranでは，配列の等価性の比較や配列形状の比較，配列内の`NaN`の検出といったところでの利用価値が高いと考え，それらの例も示しました．

Fortran 2018の機能であるassumed-rankについても言及しました．Fortranの強みである配列操作とは非常に相性が良さそうですが，コンパイラの対応状況なのか，実装の都合なのか，もう一息頑張って欲しいと思います．