AtCoder に登録したら解くべき精選過去問 10 問を Fortran90以降で解いてみた

[AtCoder に登録したら次にやること ～ これだけ解けば十分闘える！過去問精選 10 問 ～](https://qiita.com/drken/items/fd4e5e3630d0f5859067)を見かけました。

AtCoderは昔1度参加したのですが、昔過ぎてコンテスト履歴が見つかりませんでした。
それはともかく、アカウント情報を引っ張り出してきてやってみることにしました。

[先駆者の方](https://yukiutaai.wordpress.com/2018/03/21/atcoder%E3%81%AB%E7%99%BB%E9%8C%B2%E3%81%97%E3%81%9F%E3%82%89%E8%A7%A3%E3%81%8F%E3%81%B9%E3%81%8D%E7%B2%BE%E9%81%B8%E9%81%8E%E5%8E%BB%E5%95%8F10%E3%82%92fortran%E3%81%A7%E8%A7%A3%E3%81%84%E3%81%A6/)とは異なる解法もしくはこの方が使っていない応用的な関数を使って解いていこうと思います。

その過程で、基本的にはFortran90/95に基づいていますが、場合によってはFortran2003、2008の機能も使っています。


# 環境

手元のPCは以下の環境で、これで動作チェックしています。

* gcc 5.4.0
* Ubuntu 16.04
* コンパイルオプション：　`gfortran hoge.f90`

また、AtCoderでも提出して動くことを確認しています。

* gcc 4.8.4
* コンパイルオプション：　`gfortran -O2 -o a.out Main.f08`


# 個人的なFortran利点

* 頑張ればPython並＋おまじない程度のコード量
* 参照透過性を意識できる`pure`、`elemental`キーワード
  - なおやり過ぎてprintデバッグができなくなってすごい苦労した記憶も
* 数学関数に限れば`map`、`reduce`の発想に馴染みやすい

# 問題と回答

## 001: [Product](https://beta.atcoder.jp/contests/abc086/tasks/abc086_a)

偶奇判定問題です。

```fortran
program main
    implicit double precision (a-h,o-z)
    read(*, *)i, j
    k = ibits(i, 0, 1) * ibits(j, 0, 1)
    if(k == 0)then
        write(*, *) "Even"
    else
        write(*, *) "Odd"
    endif
    stop
end program main
```

`ibits(x, i, len)`は`x`の2bit表記時の`i`桁目から`len`桁取り出した値を取り出し、`x`と同じ変数型で返します。
偶数かどうかを見るには最後のビットだけ見れば十分です。
掛けてから取り出すか、取り出してから掛けるかですが、後者の方が関数呼び出しが１回減りますね。

なお、以降は`implicit double precision (a-h,o-z)`をサボります。これは`implicit none`を書かなければデフォルトで有効になります。
また、最後の`stop`も要りません。

[IBITS — Bit extraction](https://gcc.gnu.org/onlinedocs/gfortran/IBITS.html#IBITS)

## 002: [Placing Marbles](https://beta.atcoder.jp/contests/abc081/tasks/abc081_a)

ビット数え上げ問題です。

```fortran
program main
   read(*,'(b3)') i
   write(*,*) sum([(ibits(i, j, 1), j=0,2)])
end
```
`read`のフォーマット指定子で３桁の２進数だとして読み込ませ、そのまま`ibits`で数えてしまいます。なんかビット版`count_nonzero`みたいな関数がなくて苦労した。

なお`i`が普通の`integer`なので、この方法だと拡張できるのは32桁までとなります。
もし問題が64桁までの入力となった場合、

```
integer(8) :: i
```
の行を最初に足すか、もしくはコンパイルオプションに`-fdefault-integer-8`をつけたら良いかと思います。
６５桁以上となる場合、根本的にアルゴリズムを修正する必要があります。

(追記)

コメントで指摘いただきました。`popcnt`関数で数え上げができます。Fortran2008〜です。

```fortran
read(*,'(b3)')i
print *,popcnt(i)
end
```

AtCoder提出で39byte。

[POPCNT — Number of bits set](https://gcc.gnu.org/onlinedocs/gfortran/POPCNT.html#POPCNT)

## 003: [Shift only](https://beta.atcoder.jp/contests/abc081/tasks/abc081_b)

ビット数え上げ問題。

```fortran
program main
    integer, allocatable :: a(:)
    read(*,*) N
    allocate(a(N))
    read(*, *)(a(i), i = 1,N)
    write(*, *) minval(trailz(a))
    deallocate(a)
end
```
２で何回わり切れるかは、２進数表記だと０が末尾にいくつ付いているかということになります。
Fortran2008から`trailz`というそのものの関数が追加されました。
誰が使うんだろう。


[TRAILZ — Number of trailing zero bits of an integer](https://gcc.gnu.org/onlinedocs/gfortran/TRAILZ.html#TRAILZ)

そしてこの関数は`elemental`だったので、配列を渡せば要素ごとに適用して同shapeの配列を返します(つまりmap)。なので答えに必要な数字を`minval(trailz(a))`だけで得ることができます。

## 004: [Coins](https://beta.atcoder.jp/contests/abc087/tasks/abc087_b)

場合の数問題です。

```fortran
program main
   integer, allocatable :: s(:, :, :)
   read(*, *)ia, jb, kc, ix
   allocate(s(0:ia, 0:jb, 0:kc))
   forall(i=0:ia, j=0:jb, k=0:kc) s(i, j, k) = i*500 + j*100 + k*50
   write(*, *) count(s == ix)
   deallocate(s)
end
```

`ix`がターゲットの金額です。与えられた枚数から可能な支払額をすべて`s`に入れて`count(s == ix)`とまとめて比較、数え上げることにしました(`==`をmapして`count`でreduce)。
`s`の構成に`forall`を使い、実質３重ループを１行に収めました。

## 005: [Some Sums](https://beta.atcoder.jp/contests/abc083/tasks/abc083_b)

１０進数表記の任意の桁を抽出できるかどうかです。

```fortran
program main
   read(*, *) N, ia, ib
   j = 0
   do i = 1, N
      k = sum((/( mod(i/(10**m), 10), m=0,4 )/))
      if (k < ia .or. k > ib) cycle
      j = j + i
   enddo
   write(*, *)j
end
```

`mod(i/(10**m), 10)`は`m`桁目の数をとりだします。これをそのままArray Constructorをつかって `(/(m桁目の数, m=0, 4)/)`のように配列をその場で作り、そのまま`sum()`に渡します。
あとは定義通り。

ちなみにArray Constructor `(//)`ですが、DOループを使う場合には`(/(f(i), i=1,N)/)`と内側にもうひとつ`()`が必要なようです。

## 006: [Card Game for Two](https://beta.atcoder.jp/contests/abc088/tasks/abc088_b)

ソート問題です。適当にクイックソートを実装します。

```fortran
program main
   integer, allocatable :: j(:), k(:)
   read(*, *) N
   allocate(j(N))
   read(*, *)(j(i), i=1,N)
   k = iqsort(j)
   write(*, *)sum(k(1:N:2)) - sum(k(2:N:2))
   deallocate(j)
   deallocate(k)
   stop
contains
recursive function iqsort(a) result(r)
  integer, intent(in) :: a(:)
  integer, allocatable :: r(:)
  logical, allocatable :: q(:)
  integer:: p
  if(size(a) < 2)then
    r = a
  else
    p = a(1)
    q = a(2:) < p
    r = [iqsort(pack(a(2:), .not. q)), p, iqsort(pack(a(2:), q))]
  endif
end function iqsort
end
```

ロジックの行数とそれ以外の言語のお約束行数が同じくらいになる、この感じ。
`iqsort`では`deallocate`がありませんが、これはFortran95からの仕様で、「スコープを抜けた`allocatable`アレイは自動的な`deallocate`が保証される」というのがあります。
一方で`program`本体の終了時には`deallocate`されないようです。

ついでにもうひとつ、ここで初出の`allocate`周りの挙動。
`iqsort`内の`r`と`q`が`allocate`されていませんが、問題ありません。
`r = a`などの代入時に左辺のALLOCATABLEアレイが未割当だった場合、あるいは`shape`が一致せず代入できない場合、自動的に代入される行列の`shape`に合わせて`allocate`しなおし、代入までやってくれます。
一方で`read`は代入ではないので、前もって`allocate`を明示的にやっておく必要があります。

あと、FortranのArray Constructorの`[]`ですが、`[A, B, C, ...]`とすると中に並べた変数の型さえ一致しておけば、ベクトルかスカラーかの区別なくつなげたベクトルを返してくれます。

```fortran
    r = [iqsort(pack(a(2:), .not. q)), p, iqsort(pack(a(2:), q))]
```
を見てみると、`iqsort()`は長さ０の可能性もあるアレイを返す一方、`p`はスカラーです。
何気に便利です。今回のようなクイックソートではすごいコード量が減りました。
クイックソート自体の工夫はしておりません。ピボットも先頭の要素を使います。

ソートし終わったあとは、普通に１つ飛びに要素を取り出して合計し、差を取るだけです。

```fortran
write(*, *)sum(k(1:N:2)) - sum(k(2:N:2))
```

## 007: [Kagami Mochi](https://beta.atcoder.jp/contests/abc085/tasks/abc085_b)

重複検出、もしくはバケツ法。

```fortran
program main
   integer, allocatable :: d(:), c(:)
   read(*, *)N
   allocate(d(N))
   read(*, *) d
   c = (/d(1)/)
   do i = 2, N
      if(any(c == d(i))) cycle
      c = [c, d(i)]
   enddo
   write(*, *) size(c)
   deallocate(d)
   deallocate(c)
end
```
頭から順番に見ていって、見つかった数字を格納する`c`アレイを用いることにしました。
` if(any(c == d(i))) `で既に出現したかどうかの確認を行います。
`      c = [c, d(i)]`で新しく出現した数字を追加します。
最後は`c`の配列の長さを出して終わり。

(追記)

$N$も$d_i$もどうせ１００以下だったので、もう必要な配列はあらかじめ宣言してしまってバケツ法を使うことにしました。

```fortran
integer::d(100),a(100)=0
read(*,*)N
read(*,*)d(:N)
forall(i=1:N)a(d(i))=1
print *,count(a>0)
end
```
`a(d)`が1であれば、半径`d`が少なくとも一回は出現したことになります。

## 008: [Otoshidama](https://beta.atcoder.jp/contests/abc085/tasks/abc085_c)

素直に考えて$O(N^2)$としてもいいのですが、みんなやっていると思うので、$O(1)$のアルゴリズムで解きます。

まず問題の状況を式で書くと、

```math
x + y + z - N = 0\\\
10000x + 5000y + 1000z - Y = 0
```
の2つの式を満たす整数の探索となります。ここで$Y = 1000Y'$として、2番めの式を

```math
10x + 5y + z - Y' = 0
```
としておきます。$Y'$はすべて1000円札としたときの1000円札の枚数です。
この連立方程式は三次元空間中の２つの平面の式となります。そしてこれらの法線ベクトルは異なるので、交線があるはずです。求める答えは、この交線上ですべての成分が整数となる点です。この交線上で答えを探索すれば$O(N)$アルゴリズムになるかと思います。

ふたつの式の差を取ります。

```math
9x + 4y - (Y' - N) = 0
```
交線上の2点$r_1$、$r_2$を適当に取ります。計算しやすいように$x$, $y$をどちらか０にしてやります。

```math
\vec{r_1} = \left(0, \frac{Y' - N}{4}, N - y \right) = \left(0, \frac{Y' - N}{4}, \frac{5N - Y'}{5}\right)\\ 
\vec{r_2} = \left(\frac{Y' - N}{9}, 0, N - x \right) = \left(\frac{Y' - N}{9}, 0, \frac{10N-Y'}{9}\right)
```
この2点間の方向ベクトルを得ます。

```math
\vec{t_N} = \vec{r_1} - \vec{r_2} = \left(\frac{N - Y'}{9},\frac{Y' - N}{4},\frac{5N - 5Y'}{36} \right)
```
共通項を消して整数化します。

```math
\vec{t} = \frac{36}{N-Y'}\vec{t_N} = (4, -9, 5)
```

（あるいは2つの平面の法線ベクトルのクロス積が $\vec{t}$ となります。この場合 $\vec{r_2}$ を計算する必要はありません。）

交線上の点は適当なパラメータ$p$を用いて$\vec{r_1} + p\vec{t}$ と書けます。
ここで$p = (N - Y')\vec{t}/4$のケースを考えると

```math
\vec{r} = \vec{r_1} - \frac{N - Y'}{4}\vec{t} = (Y' - N, 2N - 2Y', Y')
```
となり、$\vec{r}$の要素はすべて整数となります。ここで$k \in \mathbb{Z}$として

```math
\vec{r} + k\vec{t}
```
の成分すべてが０以上$N$以下となる$k$が存在すればそれが答えとなります。

あとは$k$の条件を絞ります。

* $Y' < N$ の場合、条件を満たすものは枚数が多すぎて存在しない：$N - Y' \le 0$
* 1000円札枚数 $Y' + 5k$ は $Y'$ 以下のはず：$k \le 0$
* 5000円札枚数 $2N - 2Y' - 9k$ は0以上のはず ：$k \le \frac{2(N - Y')}{9}$
* 万札枚数 $Y' - N + 4k$ も０以上のはず：$k \ge \frac{N - Y'}{4}$

よって、

```math
\frac{N - Y'}{4} \le k \le \frac{2(N - Y')}{9}
```

の範囲で整数 $k$ を調べれば良いことがわかります。そして問題は存在すれば１つだけ挙げたらよいので、最大値直下の整数

```math
k_0 = \lfloor \frac{2(N - Y')}{9} \rfloor
```

として、$\vec{r} + k_0\vec{t}$ が条件を満たすかどうかをチェックしたら終了です。（これが条件を満たさないと、解となる整数は存在しないとして良いかなと思います。）

以上の話をFortranコードにします。

```fortran
program main
   integer :: x(3),t(3) = (/ 4, -9, 5 /)
   read(*, *)n,m
   i = n - m/1000
   dmax = i*2.d0/9.d0
   k = floor(dmax)
   x = (/-i, 2*i, m/1000/) + k*t
   if(any(x < 0)) x = -1
   write(*, *) x
end
```

`   k = floor(dmax)`ですが、`int(dmax)`とすると、`dmax`が必ず負であり、また負の場合は`int(-15.2) = -15`と切り上げとして働いてしまうため、`floor`を使いました。

答えをえる処理を行う `x = (/-i, 2*i, m/1000/) + k*t` ですが、$\vec{r}$ に当たる変数を宣言するのが面倒だったので、Array Constructor `(/ /)`を使って式の中で$\vec{r}$を定義して使いました。

とりあえず得た答えが条件（0以上$N$以下）を満たすかどうかは、全ての成分が０以上かどうかだけをチェックしています。

正直アルゴリズム自体は自信ないですが、AtCoderのテストを通ったからいいかなと。

## 009: [Daydream](https://beta.atcoder.jp/contests/abc049/tasks/arc065_a)

文字列の操作とGreedy法。

```fortran
program main
   character(len=7) :: s(4) = (/"dream  ", "dreamer", "erase  ", "eraser "/)
   character(len=100000) :: b
   integer :: L(4) = (/5, 7, 5, 6/)
   read(*, *) b
   N = len_trim(b)
   a: do while(N > 0)
      do i = 1, 4
         if(N < L(i)) cycle
         if(trim(s(i)) == b(N-L(i)+1:N))then
            N = N - L(i)
            cycle a
         endif
      enddo
      write(*, *)"NO"
      stop
   enddo a
   write(*, *) "YES"
end
```
ここで、`read`が自動で`allocate`してくれない問題が出てきてしまいました。
Fortranにはdeferred length characterという文字列の長さを可変にできる機能があります。宣言や使い方は以下のようにします。

```fortran
   character(len=:), allocatable :: b
   b = "hogehoge"
   b = b // "fuga"
```

しかし以下はエラーです。

```fortran
   character(len=:), allocatable :: b
   read(*, *) b  ! error : bが未割当的なメッセージ 
```

https://stackoverflow.com/questions/31084087/using-a-deferred-length-character-string-to-read-user-input

このため、`b`を問題の入力の最大長の固定長文字列にする必要があります。

`len_trim`で実際の文字列長を取得。
文字列の比較自体は`==`で可能です。スペースがあると良くないので、`s`の方は`trim`してから比較します。もし一致していたら文字列を取り除くのと同等の操作、文字列長`N`を減らしていきます。
このとき内側のループから外側のループの先頭に戻る必要があるので、外側のループに名前`a`をつけています。
`      write(*, *)"NO"`の行ですが、ここに来たということは４つの文字列すべてと異なる、つまり与えられた文字列を作ることができないということになります。なので内側のループが終わってから"NO"を出力してプログラムを停止させます。

## 010: [Traveling](https://beta.atcoder.jp/contests/abc086/tasks/arc089_a)

結局偶奇判定。

```fortran
program main
   integer, allocatable :: i(:, :)
   read(*, *)n
   allocate(i(3, 0:n))
   i(:, 0) = (/0, 0, 0/)
   read(*, *)i(:, 1:n)
   do k = 1, n
      jx = abs(i(2, k) - i(2, k - 1)) + abs(i(3, k) - i(3, k - 1)) 
      jt = i(1, k) - i(1, k - 1)
      if(jx > jt .or. ibits(jx, 0, 1) /= ibits(jt, 0, 1))then
         write(*, *) "No"
         stop
      endif
   enddo
   write(*, *)"Yes"
end
```

入力の$t$と$x, y$をまとめて`i`に入れてしまいます。また、`i`の最初は$0, 0, 0$ですが、これは入力に与えてくれませんので、自分で前もってつけておく必要があります。
あとは入力を読み込み、順番に行程をチェックしていきます。所要時間を`jt`、移動距離を`jx`として、これを元に判断します。１つでも不可能なコースがあれば`"No"`を出力して停止、問題なくDOループが終わった時に限り`"yes"`を出力させます。
パリティのチェックは偶奇の一致、つまり最終bitの一致を見ます。
またパリティは成分ごとに`abs`を適用しても変わらないので、`jx`を距離判定にもパリティ判定にも使いまわします。

あまり変なことはしていないと思います。

あ、`deallocate(i)`を忘れてた。

# その他

全体的に`deallocate`をサボる、というか問題の最大長配列をあらかじめ確保して宣言する、`enddo` `endif`を`end`とする、`program main`をサボる、スペースを削るなどすれば文字数はもっと減るはずですが、ゴルフ方向のチャレンジはやめておきました。

他の人の提出履歴見たら、Fortranでやってる人は僕の他に２，３人しか見かけませんね・・・こんなに有名な言語なのに・・・。

# 参考

* [AtCoder に登録したら次にやること ～ これだけ解けば十分闘える！過去問精選 10 問 ～](https://qiita.com/drken/items/fd4e5e3630d0f5859067)
* 先駆者様：[AtCoder に登録したら解くべき精選過去問 10 問を FORTRAN で解いてみた](https://yukiutaai.wordpress.com/2018/03/21/atcoder%E3%81%AB%E7%99%BB%E9%8C%B2%E3%81%97%E3%81%9F%E3%82%89%E8%A7%A3%E3%81%8F%E3%81%B9%E3%81%8D%E7%B2%BE%E9%81%B8%E9%81%8E%E5%8E%BB%E5%95%8F10%E3%82%92fortran%E3%81%A7%E8%A7%A3%E3%81%84%E3%81%A6/)
* [GNU Fortranの組み込み関数一覧](https://gcc.gnu.org/onlinedocs/gfortran/Intrinsic-Procedures.html)