ABC191をHaskellで

# A - Vanishing Pitch

[問題 ABC191A](https://atcoder.jp/contests/abc191/tasks/abc191_a)

シグネチャを決める。

```haskell
abc191a :: Int  -- V
        -> Int  -- T
        -> Int  -- S
        -> Int  -- D
        -> Bool -- 答え
abc191a v t s d = d < v * t || v * s < d
```

# B - Remove It

[問題 ABC191B](https://atcoder.jp/contests/abc191/tasks/abc191_b)

シグネチャを決める。

```haskell
abc191b :: Int    -- N
        -> Int    -- X
        -> [Int]  -- Ai
        -> [Int]  -- 答え
abc191b _n x as = filter (x /=) as
```

# C - Digital Graffiti

[問題 ABC191C](https://atcoder.jp/contests/abc191/tasks/abc191_c)

シグネチャを決める。

```haskell
abc191c :: Int       -- H
        -> Int       -- W
        -> [String]  -- Sij
        -> Int       -- 答え
```

[フレンズさんのヒント](https://twitter.com/kyopro_friends/status/1358048374079909890)にあるとおり、角を数えればいい。
つまり、2x2 の領域だけを見て、

こういうところには角はなく
```
.. | ## | ## | .. | #. | .#
.. | ## | .. | ## | #. | .#
```

これとこれの回転3とおりの場所には角があり、
```
.. | .#
.# | ##
```

自己交差がないという保証があるのでこういう箇所はない。
```
#. | .#
.# | #.
```

つまり、2x2で見て `#` が奇数個な箇所を数えればよい。

## 結果

```haskell
import Data.Array

abc191c :: Int -> Int -> [String] -> Int
abc191c h w ss = length
  [ ()
  | i <- [1 .. pred h], j <- [1 .. pred w]
  , odd $ length [() | p <- [i, succ i], q <- [j, succ j], cs ! (p,q)]
  ]
  where
    cs = listArray ((1,1),(h,w)) $ map ('#' ==) $ concat ss
```

# D - Circle Lattice Points

[問題 ABC191D](https://atcoder.jp/contests/abc191/tasks/abc191_d)

AtCoder Problemsで見ると、妙に難易度が高く評価されている。
問題文から、浮動小数点数の誤差でギリギリのところを突かれそうな匂いがプンプンする。

## こまけぇこたぁいいんだよ！

その気配を無視して `Double` でやってみる。
シグネチャを決める。

```haskell
abc191d :: Double  -- X
        -> Double  -- Y
        -> Double  -- R
        -> Int     -- 答え
```

$\lceil X - R \rceil \leq x \leq \lfloor X + R \rfloor$ の範囲で整数座標 $x$ を考える。この垂直な線の上で、円の中にある格子点の一番下の $y$ 座標は、$e = \sqrt{R^2 - (X - x)^2}$ として $\lceil Y - e \rceil$ 一番上は $\lfloor Y + e \rfloor$ である。

```haskell
abc191d x y r = sum
    [ floor (y + e) - ceiling (y - e) + 1
    | x1 <- [ceiling $ x - r .. floor $ x + r]
    , let e = sqrt $ r * r - (x - fromIntegral x1) ^ 2
    ]
```

ほとんと正解だが3つほどWAになる。後はその計算誤差を克服する方法だけだが…

## Rationalで無限精度

有理数型を使えば精度の心配はしなくてよくなるが、平方根は無理数になるので上の手順そのままでは　できない。

まず、引数は文字列のまま受け取り、`Double`にすることなく読み込む。
小数点がないときは、整数として読む。
小数点があるときは、小数点を除いた列を読んだ整数を分子、小数点以降の文字列を、`.`を`1`、数字を`0`に書き換えた列を読んだ整数を分母とすればよい。

```haskell
readDec :: String -> Rational
readDec xs
  | null bs   = read as % 1
  | otherwise = read (as ++ bs1) % read ('1' : map (const '0') bs1)
  where
    (as,bs) = span ('.' /=) xs
    bs1 = tail bs
```

`sqrt`が使えないので、$(X - x)^2 + (Y - y)^2 \leq R^2$ を満たす整数座標 $y$ の最大値と最小値を整数の二分探索で求めることにする。

```haskell
abc191d :: String -> Int
abc191d args = sum
    [ y2 - y1 + 1
    | x1 <- map fromIntegral [ceiling $ x - r .. floor $ x + r]
    , let (_,y1) = binarySearch (prop x1) (pred $ floor   $ y - r) (succ $ floor   y)
    , let (_,y2) = binarySearch (prop x1) (succ $ ceiling $ y + r) (pred $ ceiling y)
    , y1 <= y2
    ]
  where
    [x,y,r] = map readDec $ words args
    prop x1 y1 = sq (y - fromIntegral y1) <= sq r - sq (x - x1)
    sq x = x * x

binarySearch :: (Int -> Bool) -> Int -> Int -> (Int, Int)
binarySearch prop unsat sat = loop unsat sat
  where
    loop a b
      | ende   = (a, b)
      | prop m = loop a m
      | True   = loop m b
      where
        ende = a == m || b == m
        m = div (a + b) 2
```

これだとTLEする。
`prop`の計算が何度も行われる中で、`prop x1` という部分適用の結果、`y1`以外の変数は確定するが、しかし飽和するまで計算はされないので、毎回全ての計算が無駄に行われている。これを

```haskell
    prop x1 = let d = sq r - sq (x - x1) in \y1 -> sq (y - fromIntegral y1) <= d
```

と変更すると、`x1`が与えられたとき（正確にはさらに最初に`y1`が与えられて実行されたとき）に、内部変数`d`の値が確定し、以後の呼び出しで共有される。
結果：[AC, 1071ms](https://atcoder.jp/contests/abc191/submissions/51954147) かなりかかる。

## 公式解説のやり方

小数点以下4桁までと決まっているので、$10^4$のゲタを履かせた固定小数点で考えよ、と。
なるほどそれならRationalのようなデータ型を持たない言語でも普通に書ける。

まず、小数点以下4桁を補って整数で読み込む。
小数点がある場合は、それより後ろ4文字を取り出す、足りない分は`0`を追加する。
小数点がない場合は、`0`を4つ追加する。
という計算を少々トリッキーに書いてみた。

```haskell
readDec :: String -> Int
readDec xs = read $ (as ++) $ take 4 $ tail $ bs ++ "00000"
  where
    (as,bs) = span ('.' /=) xs
```

本体では、`floor`を使う代わりに$10^4$で切り捨て除算、`ceiling`を使う代わりに切り上げ除算をすればよい。

```haskell
abc191d args = sum
    [ div y2 base - divrup y1 base + 1  -- ⌊y2/base⌋ - ⌈y1/base⌉ + 1
    | x1 <- [x10, x10 + base .. x + r]
    , let dy2 = sq r - sq (x - x1)
    , let y1 = snd $ binarySearch (prop dy2) (pred y - r) y
    , let y2 = snd $ binarySearch (prop dy2) (succ y + r) y
--    y1 <= y2
    ]
  where
    [x,y,r] = map readDec $ words args
    x10 = base * quot (x - r) base
    sq x = x * x
    prop dy2 y1 = sq (y - y1) <= dy2

base :: Int
base = 10^4

divrup :: Integral a => a -> a -> a
divrup x y = negate $ div (negate x) y
```

結果 [AC, 51ms](https://atcoder.jp/contests/abc191/submissions/51967170)

コンテスト時間中にここにたどり着ける気がしない。

# E - Come Back Quickly

[問題 ABC191E](https://atcoder.jp/contests/abc191/tasks/abc191_e)

シグネチャを決める。

```haskell
abc191e :: Int      -- N
        -> Int      -- M
        -> [[Int]]  -- Ai, Bi, Ci
        -> [Int]    -- 答え
```

$N \leq 2000$ は、ワーシャルフロイド法で距離を数えるには大きすぎる。
ダイクストラ法でそれぞれ数える方が間に合う、というのも変な感じだが、そういうものだと納得するしかない。

それぞれの町 $i$ について、$i$ の複製を $N + 1$ として作り、$N+1$から$i$への距離を数えることを行えばよい。

## 結果

```haskell
import Data.Array
import qualified Data.Heap as H
import Control.Monad.ST
import Data.Array.ST

abc191e :: Int -> Int -> [[Int]] -> [Int]
abc191e n m abcs = map f [1 .. n]
  where
    g = accumArray (flip (:)) [] (1,n) [(a, (b,c)) | a:b:c:_ <- abcs]
    n1 = succ n
    f s = runST $ do
      d <- dijkstra n1 (e s) n1                 -- ダイクストラ法を実行して
      v <- readArray d s                        -- N+1からsまでの距離を取り出す
      return $ if v == maxBound then -1 else v
    e s i                                       -- スタートsを複製したグラフでのiからの出辺を渡す関数
      | i == n1   = g ! s
      | otherwise = g ! i

-- @gotoki_no_joe
dijkstra :: Int                  -- 頂点数N (1～N)
         -> (Int -> [(Int,Int)]) -- 隣接頂点とその辺の重み、グラフの情報
         -> Int                  -- 開始点
         -> ST s (STArray s Int Int)
dijkstra n edges start =
  do
    dist <- newArray (1,n) maxBound
    writeArray dist start 0
    loop dist $ H.singleton (H.Entry 0 start)
  where
    loop dist queue | H.null queue = return dist
    loop dist queue = do
      let Just (H.Entry cost u, queue1) = H.uncons queue
      du <- readArray dist u
      if du < cost then loop dist queue1 else do
        queue2 <- foldM (\q (v, we) -> do
            let d1 = du + we
            dv <- readArray dist v
            if d1 >= dv then return q else do
              writeArray dist v d1
              return $ H.insert (H.Entry d1 v) q
            ) queue1 (edges u)
        loop dist queue2
```

特定のノードの距離が確定した時点で打ち切るような変種を作れば、実行時間をより節約することもできるだろう。

## 解説のやり方

公式解説、ユーザ解説ともに、グラフはオリジナルのままでダイクストラ法を実行することが前提で、町$a$から$b$への距離を$d(a,b)$として、

各町$i$について

- $i$ から $i$ への自己辺があるならその長さは候補の一つ
- ダイクストラ法の結果を使って、全ての $j \neq i$ について、$d(i,j) + d(j,i)$ は最大 $N-1$ 個の候補

から最小値を選ぶという方針を説明している。

どうせ各$i$について別々にダイクストラ法を実行するのだから、上のアプローチの方が簡潔なのでは…

[フレンズさんによると](https://twitter.com/kyopro_friends/status/1358050084236058627)、ダイクストラ法で間に合うのは、辺の数が少なめだからという要素があるらしい。

# F - GCD or MIN

[問題 ABC191F](https://atcoder.jp/contests/abc191/tasks/abc191_f)

シグネチャを決める。

```haskell
abc191f :: Int    -- N
        -> [Int]  -- Ai
        -> Int    -- 答え
```

わからなくて[フレンズさんのヒント](https://twitter.com/kyopro_friends/status/1358050768570306561)を見た。けど、「残せるか」がちょっと理解できなかった。その先はわかった。

- $\min(x,y)$を選ぶとは、つまり、大きい数を任意に消せるということ。こればかり繰り返すと、書かれている最小値が最終的に残る。
- $\gcd(x,y)$を選ぶことで、既存の数よりさらに小さい数が作り出せる。

つまり、

- $A_i$ の最小値は、常に$\min$を使うことで残せる。この値を$Amin$と呼ぶことにする。
- $Amin$ より大きな数は決して残せない。
- $A_i$ より小さい数は、$\gcd$の組みあわせで作ることができれば、残せる。

ということになる。そして「$\gcd$の組みあわせで作ることができる$Amin$以下の数」を全て見つけるには、
まず全ての $A_i$ について、その（$Amin$以下の）全ての約数 $f$ に対して、$A_i$ は $f$ の倍数であると登録した辞書を作る。
辞書に登録されたそれぞれの約数 $f$ について、関連付けられた倍数 $A_i$ は、それらの任意の組みあわせについて $\gcd$ をとると $f$ の倍数だけが得られ、可能な最小値は $f$ である。それより小さい値は作れない。よって、全体の $\gcd$ が $f$ になるなら $f$ は作ることができるとわかる。

実際には、辞書に登録するとき、これまでに登録された $A_i$ との $\gcd$ をとるようにすればよい。

## 結果

```haskell
abc191f :: Int -> [Int] -> Int
abc191f n as = length $ filter (uncurry (==)) $ IM.assocs im
  where
    amin = minimum as
    im = IM.fromListWith gcd [(f, a) | a <- as, f <- takeWhile (amin >=) $ factors a]

-- @gotoki_no_joe
factors :: Int -> [Int]
factors 1 = [1]
factors n = 1 : loop 2 [n]
  where
    loop k us
      | k2 >  n =     us
      | k2 == n = k : us
      | r  == 0 = k : next (q:us)
      | True    =     next    us
      where
        (q,r) = divMod n k
        next = loop (succ k)
        k2 = k * k
```

約数を求める `factors` を除いた本体 `abc191f` はたった 3 行。