Fortranにおける余り/剰余の計算

数学的に厳密な定義を知らないので，情報をご存じの方はご教示ください．

# 余り
整数同士の除算において，被除数が除数で割り切れない場合，除算で割り切れなかった量を余り，または剰余と呼びます．つまり

```math
被除数=除数 \times 商 + 余り
```

と表せます．余りは，除数の絶対値より小さくなければなりません．

被除数と除数の符号が同じ時は何も考えなくてよいのですが，被除数と除数の符号が異なるとき，余りをどちらの符号に合わせるかで，2通りの計算が考えられます．

例えば，被除数を27，除数を4とした場合，余りは3です．同様に，被除数を-27，除数を-4とした場合，余りは-3です．どちらも，上の定義式に当てはめてみれば簡単に確認できます．

```math
27=4 \times 6 + 3
```

```math
-27=-4 \times 6 -3
```

次に，被除数と除数の符号が異なる場合を考えてみましょう．被除数を-27，除数を4とした場合，余りは除数の絶対値より小さければよいので，-3を選ぶことも，1を選ぶこともできます．このとき，商も異なります．余りが-3のとき商は-6であり，余りが1のとき，商は-7です．

```math
-27=4 \times -6 -3
```

```math
-27=4 \times -7 + 1
```

前者は，余りの符号が被除数と同じで，後者は除数と同じです．厳密な定義と一致しているか分かりませんが，前者をremainder，後者をmodulo，両者を総称して余りと呼ぶことにします．

# Fortranにおける余りの演算

プログラミング言語のいくつかは，remainderもしくはmoduloを計算する演算子が定義されており，そうでない言語では，remainderもしくはmoduloを計算する関数が定義されています．

Fortranには余りを計算する演算子は定義されておらず，関数が定義されています．プログラミング言語によっては，remainderもしくはmoduloのどちらかしか定義されていない事もありますが，Fortranでは両方定義されています．

|関数|処理内容|
|:--:|:--|
|mod|remainderを計算する（余りの符号が被除数と同じ）|
|modulo|moduloを計算する(余りの符号が除数と同じ)|

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

    integer(int32) :: a,b
    a = 27
    b = -4
    print '(4I4)',a,b,mod(a,b),modulo(a,b)
    !  27   4   3   3
    ! -27  -4  -3  -3
    !  27  -4   3  -1
    ! -27   4  -3   1
end program main
```

「おいおい，remainderを計算するのにmodかよwさすが絶滅危惧種だなw」という声が聞こえてきそうですが，Fortranはまだましな方です．Wikipediaの剰余演算[^1]の項目には多くのプログラミング言語における余りの演算子あるいは関数が載っていますが，この統一性の無さは壮観です．

[^1]: https://ja.wikipedia.org/wiki/剰余演算

# moduloの使い道
## 差分法
差分法は，微分を数値的に計算する手法の一つです．

$x$の関数$f(x)$がなめらかであれば，$f$の微分は次式で定義されます．

```math
\frac{df}{dx} = \lim_{\Delta x \rightarrow0}\frac{f(x+\Delta x)-f(x)}{\Delta x}
```

関数$f$に対して，$x$方向に$\Delta x$だけ離れた2点を考えて，その2点の間隔を無限小に近づけた時の極限が$f$の傾きです．しかし，この無限小というのが曲者で，計算機では表現できません．このとき，$\Delta x$に対する制約を緩め，「無限小」ではなく「関数$f$の変化に対して十分小さい」と考えると，微分の定義式は減算と除算を組み合わせた差分で表現できるようになります．

つまり，微分を

```math
\frac{df}{dx} \approx \frac{f(x+\Delta x)-f(x)}{\Delta x}
```

として計算します．他にも，

```math
\frac{df}{dx} \approx \frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}
```

とも計算できます．これは，中心差分などと呼ばれます．

微分は局所的な計算なので，微分値は位置$x$の1点のみで求まります．1次元空間全体の微分を計算するには，下記表のように，計算したい関数の値を$\Delta x$間隔で記録しておき，繰り返し処理を用いて異なる配列の要素を参照しながら，各点における微分値を計算することになります．

|i|1|2|3|4|5|6|7|8|9|
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
|x|0.|0.785|1.57|2.36|3.14|3.93|4.71|5.50|6.28|
|f|0.|0.707|1.0|0.707|0.|-0.707|-1.0|-0.707|0.|

プログラムで書くと，次のようになります．

```Fortran
integer(int32) :: i
do i = 1, 9
    d_f_dx(i) = (f(i+1)-f(i-1))/(2d0*dx)
end do
```

計算には配列に記録された値を用いるので，`i=1`, `i=9`では配列外参照が起こります．関数$f$に周期性があるとき，周期関数を利用することで，この問題を回避できます．例えば，$f(x)=\sin ⁡x$は周期$2\pi$の周期関数なので，$f(2\pi)=f(0)$, $f(2\pi+i\Delta x)=f(i\Delta x)$,$f(0-i\Delta x)=f(2\pi-i\Delta x)$という関係が得られます（$i$は任意の自然数）．つまり，周期関数の1周期だけを取り出して計算します．$x=2\pi$の関数値は$x=0$と等しいので，上の表のうち，`i=9`の値は配列に保持する必要はなくなります．そして，`i=1`のときは，$f(0-\Delta x)=f(2\pi-\Delta x)$なので，`i-1(=0)`のときは`i=8`の値を参照すればよいことがわかります．同様に，`i=8`のときは，$f([2\pi-\Delta x]+\Delta x)=f(0)$なので，`i+1(=9)`のときは`i=1`の値を参照すればよいわけです．

### moduloの利用
`i-1`や`i+1`を`modulo`関数によってたたみ込んでしまえば，配列外参照を回避しつつ，周期境界条件の反映を1行で行えます．

ところで，Fortranは配列の添字を任意の整数の範囲で定められますが，特に指定しない場合には1から始まります．ところが，`modulo`は結果を0から除数-1の範囲で返すので，結局配列外参照が起こります．そのため，被除数から1を引いてからmoduloを計算し，結果に1を足す必要があります．

|i|1|2|3|4|5|6|7|8|
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
|`modulo(i,8)`|1|2|3|4|5|6|7|0|
|`modulo(i-1,8)`|0|1|2|3|4|5|6|7|
|`modulo(i-1,8)+1`|1|2|3|4|5|6|7|8|
|`modulo((i-1)-1,8)+1`|8|1|2|3|4|5|6|7|
|`modulo((i+1)-1,8)+1`|2|3|4|5|6|7|8|1|

プログラムでは，`i+1`, `i-1`をそれぞれ変数`ip1`, `im1`に代入して利用します．

```Fortran
integer(int32) :: i, ip1, pm1
do i = 1, 8
    ip1 = modulo((i+1)-1,8) + 1
    im1 = modulo((i-1)-1,8) + 1
    d_f_dx(i) = (f(ip1)-f(im1))/(2d0*dx)
end do
```

性能的には不利になるでしょうが，かなり簡潔に書けるようになります．この書き方は，差分の計算に使う点が増えても（ステンシルが広がっても）同じ書き方でいけるのが利点です．

## 符号なし整数
Fortranは符号なし整数を持っていませんが，moduloを使うことで擬似的に再現できます．
例えば，`integer(int8)`は-128~127の値を取ります．除数を256としてmoduloを計算すると，0から255までの値が得られます．

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

    integer(int8)  :: c
    integer(int16) :: i
    do i = -128_int16, 127_int16
        c = int(i,int8) ! int8の整数
        print *,i,modulo(int(c,int16),256_int16) ! 256はint8で表現できないので，int16へ変換
    end do
end program main
```

実行結果を見ると，確かに負の数から正の数（256を法として合同な正の整数と言えばよいのか？）が得られています．

これを`integer(int8)`型の変数に保存すると，結局-128～127にたたまれてしまうので，演算や出力の時にしか使うことができません．

```console
   -128    128
   -127    129
   -126    130
   -125    131
   -124    132
   -123    133
   -122    134
   -121    135
   -120    136
   -119    137
   -118    138
   -117    139
   -116    140
   -115    141
   -114    142
   -113    143
   -112    144
   -111    145
   -110    146
   -109    147
   -108    148
   -107    149
   -106    150
   -105    151
   -104    152
   -103    153
   -102    154
   -101    155
   -100    156
    -99    157
    -98    158
    -97    159
    -96    160
    -95    161
    -94    162
    -93    163
    -92    164
    -91    165
    -90    166
    -89    167
    -88    168
    -87    169
    -86    170
    -85    171
    -84    172
    -83    173
    -82    174
    -81    175
    -80    176
    -79    177
    -78    178
    -77    179
    -76    180
    -75    181
    -74    182
    -73    183
    -72    184
    -71    185
    -70    186
    -69    187
    -68    188
    -67    189
    -66    190
    -65    191
    -64    192
    -63    193
    -62    194
    -61    195
    -60    196
    -59    197
    -58    198
    -57    199
    -56    200
    -55    201
    -54    202
    -53    203
    -52    204
    -51    205
    -50    206
    -49    207
    -48    208
    -47    209
    -46    210
    -45    211
    -44    212
    -43    213
    -42    214
    -41    215
    -40    216
    -39    217
    -38    218
    -37    219
    -36    220
    -35    221
    -34    222
    -33    223
    -32    224
    -31    225
    -30    226
    -29    227
    -28    228
    -27    229
    -26    230
    -25    231
    -24    232
    -23    233
    -22    234
    -21    235
    -20    236
    -19    237
    -18    238
    -17    239
    -16    240
    -15    241
    -14    242
    -13    243
    -12    244
    -11    245
    -10    246
     -9    247
     -8    248
     -7    249
     -6    250
     -5    251
     -4    252
     -3    253
     -2    254
     -1    255
      0      0
      1      1
      2      2
      3      3
      4      4
      5      5
      6      6
      7      7
      8      8
      9      9
     10     10
     11     11
     12     12
     13     13
     14     14
     15     15
     16     16
     17     17
     18     18
     19     19
     20     20
     21     21
     22     22
     23     23
     24     24
     25     25
     26     26
     27     27
     28     28
     29     29
     30     30
     31     31
     32     32
     33     33
     34     34
     35     35
     36     36
     37     37
     38     38
     39     39
     40     40
     41     41
     42     42
     43     43
     44     44
     45     45
     46     46
     47     47
     48     48
     49     49
     50     50
     51     51
     52     52
     53     53
     54     54
     55     55
     56     56
     57     57
     58     58
     59     59
     60     60
     61     61
     62     62
     63     63
     64     64
     65     65
     66     66
     67     67
     68     68
     69     69
     70     70
     71     71
     72     72
     73     73
     74     74
     75     75
     76     76
     77     77
     78     78
     79     79
     80     80
     81     81
     82     82
     83     83
     84     84
     85     85
     86     86
     87     87
     88     88
     89     89
     90     90
     91     91
     92     92
     93     93
     94     94
     95     95
     96     96
     97     97
     98     98
     99     99
    100    100
    101    101
    102    102
    103    103
    104    104
    105    105
    106    106
    107    107
    108    108
    109    109
    110    110
    111    111
    112    112
    113    113
    114    114
    115    115
    116    116
    117    117
    118    118
    119    119
    120    120
    121    121
    122    122
    123    123
    124    124
    125    125
    126    126
    127    127
```

# まとめ
Fortranの`mod`はremainder．