ABC275をPythonとC++で

# 0 はじめに
## 0-1 記事について
AtCoder Beginner Contest 275の解説です。
実装はPythonとC++で書きます。
公式解説と違いがあるかも知れませんがご了承。
問題の本質が簡素過ぎる場合は、問題文の欄には問題の要約を載せています。(面倒なだけ)
ミス等発見されましたらコメント欄にお願いします。

# 1 ABC275 解説
## 1-1 個人的な感想
たのしい回でした。
Diffは、Aが灰前半、Bが灰後半、Cが緑前半、Dが茶後半、Eが水前半、Fが水後半といった感じです。
Cで沼りましたがEをササっと通したので割とよいPerfが取れ、14増えました。
Perf:1361 Rating:1234→1248(+14)

## 1-2  A問題 Find Takahashi

### 問題
長さが$N$の整数列$H$が与えられるので、$H$の最大の要素が$H$の何項目なのか求めて下さい。
尚、$H$のうちどの異なる$2$要素も相異なります。

#### 制約
・$1 \leq N \leq 100$
・$1 \leq H_i \leq 10^9$

### 解説

$H$の最大値を求め、`for`文を回しそれが何項目にあるかを探索するだけでよいです。
(もっと短い時間で求める方法もありますがこれが簡単ではと。)
出力すべき値は`1-indexed`にする必要があることに注意しましょう。

言語によっては、最大値や、要素が配列の何項目にあるかのライブラリが用意されている場合があります。


<details>
<summary> Pythonでの実装例 </summary>

```py
N=int(input())
H=list(map(int,input().split()))

print(H.index(max(H))+1)
```
`<配列>.index`や`max(<配列>)`は$O(N)$です。
容易に`for`文の中で使わないようにしましょう。
</details>
<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
#define rep(i,N,M) for(int i=N; i<M; i++)

int main(){
  int N; cin>>N;
  vector<int> H(N);
  
  int m=0;
  rep(i,0,N){
    cin>>H[i];
    m=max(m,H[i]);
  }
  
  rep(i,0,H){
    if(H[i]==m){
      cout<<i+1<<endl;
    }
  }
}
```
`std::max_element`なるものも存在するらしいですがこれも$O(N)$です。
`for`文の中で容易に使わないようにしましょう。
</details>

はい。

## 1-3  B問題 ABC-DEF

### 問題
$6$整数$A,B,C,D,E,F$が与えられるので、$A×B×C-D×E×F$を$998244353$で割った余りを求めて下さい。

#### 制約
・$0 \leq A,B,C,D,E,F \leq 10^{18}$
・$A×B×C \geq D×E×F$


### 解説

($mod(n,m)$を、$n$を$m$で割った余りと定義します。)

PythonやRubyのように、多倍長整数がある言語であれば工夫無しで求められますが、C++など他の言語ではオーバーフローが起こってしまいます。
$A×B×C$は最大で$10^{54}$にもなり得るからです。

ここで、剰余演算の特徴を使いましょう。
$mod((S×T),M)=mod(mod(S,M)×mod(T,M),M)$です。

これにより、掛け算を$998244353$未満の数同士まで小さくできました。
$998244353 \leq 10^9$より、扱う数は$10^{18}$未満になりました。
C++では`long long`型を使えば$10^{18}$程度まで扱えるのでオーバーフローせずに済みます。

<details>
<summary> Pythonでの実装例 </summary>

```py
A,B,C,D,E,F=map(int,input().split())
print((A*B*C-D*E*F)%998244353)
```
単純ですね。
</details>

<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
using ll=long long;

int main(){
  ll A,B,C,D,E,F; cin>>A>>B>>C>>D>>E>>F;
  
  const ll mod=998244353;
  ll x=((A%mod)*(B%mod)%mod);
  x*=(C%mod); x%=mod;
  
  ll y=((D%mod)*(E%mod)%mod);
  y*=(F%mod); y%=mod;
  
  ll ans=(x+mod-y)%mod;
  
  cout<<ans<<endl;
}
```
C++の場合、処理環境によっては、剰余の値が負になることもあります。
その場合の処理に気を付けましょう。
因みに、ACLのmint型というものを使えば自動で998244353で割った余りを出してくれるらしいです。
</details>

この問題であれば、多倍長整数を堂々と使っても`TLE`は出ませんが、勿論その分データ容量も消費するので、値が$10^{100000}$くらいになると余裕で`TLE`が出てしまうので無暗に扱うのは危険です。
[ABC233-E](https://atcoder.jp/contests/abc233/tasks/abc233_e)は、そんな多倍長整数の濫用の危険さを教えてくれる問題です。

## 1-4 C問題 Counting Squares
### 問題
$9×9$の二次元グリッド上に、ポーンが置いてあります。
$S$はグリッドの状態を表し、$S_{r,c}$が`#`であれば$(r,c)$にはポーンが置いてあり、`.`であれば何も置いていません。

この平面上の正方形のうち、四頂点全てにポーンが置いてあるものの個数を求めて下さい。

#### 制約
特になし

### 解説

#### 速めの言語を使用する場合
速い言語であれば、ポーンの置いてある四点の組合せを全探索し、それらが正方形を為し得るか判定していっても間に合います。

判定は、各点間の距離を全て探索し、小さいものに並べた時比が$1:1:1:1:\sqrt{2}:\sqrt{2}$
となっていれば真、そうでなければ偽とすればよいです。
計算する際は誤差による影響を鑑み、距離の二乗の比が$1:1:1:1:2:2$となっているかを判定すればよいです。

全探索は、四頂点の組合せを全探索するので、同じ組を二回以上見ないように`for`文に工夫を少し入れましょう。(実装を見ればわかります)

<details>
<summary> Pythonでの実装例 </summary>

```py
from itertools import combinations

S=[list(input()) for _ in range(9)]

N=9
pair=[]
for i in range(9):
  for j in range(9):
    if S[i][j]=="#":
      pair.append((i,j))

calc=lambda a,b:(a[0]-b[0])**2+(a[1]-b[1])**2

def hantei(L):
  le=[]
  for n in combinations(range(4),2):
    a,b=n
    le.append(calc(L[a],L[b]))
  le.sort()
  for i,n in enumerate([1,1,1,1,2,2]):
    if le[0]*n!=le[i]:
      return 0
  return 1

N=len(pair)
ans=0
for i in range(N):
  for j in range(i+1,N):
    for k in range(j+1,N):
      for l in range(k+1,N):
        a,b,c,d=pair[i],pair[j],pair[k],pair[l]
        if hantei([a,b,c,d]):
          ans+=1

print(ans)
```
PyPyで提出したら990msで通りました。
PythonだとTLEが出ました。
</details>

<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
#define rep(i,N,M) for(int i=N; i<M; i++)

int calc(pair<int,int> a,pair<int,int> b){
  int ax,ay,bx,by;
  tie(ax,ay)=a; tie(bx,by)=b;
  return (ax-bx)*(ax-bx)+(ay-by)*(ay-by);
}

bool judge(vector<pair<int,int>> L){
  vector<int> le;
  rep(i,0,4){
    rep(j,i+1,4){
      pair<int,int> a=L[i],b=L[j];
      le.push_back(calc(a,b));
    }
  }

  sort(le.begin(),le.end());
  vector<int> times={1,1,1,1,2,2};
  rep(i,0,6){
    if(le[0]*times[i]!=le[i]){
      return false;
    };
  }
  return true;
}
  
int main(){
  vector<vector<char>> S(9,vector<char>(9));
  
  vector<pair<int,int>> p;
  rep(i,0,9){
    rep(j,0,9){
      cin>>S[i][j];
      if(S[i][j]=='#'){
        p.push_back(make_pair(i,j));
      }
    }
  }
  
  int N=p.size();
  int ans=0;
  rep(i,0,N){
    rep(j,i+1,N){
      rep(k,j+1,N){
        rep(l,k+1,N){
          pair<int,int> a=p[i],b=p[j],c=p[k],d=p[l];
          if(judge({a,b,c,d})){
            ans++;
          }
        }
      }
    }
  }
  
  cout<<ans<<endl;
}
```
C++速いですね...
</details>

#### 遅い言語の場合 or ポーンが1000個まで増えた場合

遅い言語の場合、上の解法ではACできない場合もあり得ます。
又、この問題ではポーンは最大でも$81$個しかないので、$O(N^4)$($N$はポーンの数)の解法が通りますが、$N=1000$とかになるとC++でももう無理です。
どうすればよいでしょうか。

正方形は、向かい合う$2$頂点$A,B$が決まれば残りの$2$頂点$C,D$も自然に決まります。

<details>
<summary> どのように決まる？ </summary>
$A,B$の中点を$M$とすると、$AM=CM,BM=DM$であり$AM \perp CM,BM \perp DM$です。
$C,D$は$A,B$の座標から簡単な計算で求められます。
詳しくは実装をご覧ください。
</details>

なので、向かい合う$2$頂点$A,B$の組合せを全探索し、$C,D$が格子点であり$C,D$上にポーンが存在するか判定すればよいです。$A,B$の組合せを全探索するので重複がないようにしましょう。

又、同じ正方形について、$A,B$でも$C,D$でも数えるので、答えを$\tfrac{1}{2}$するのを忘れないようにしましょう。


$A,B$の組の全探索に$O(N^2)$掛かり、判定は$O(1)$でできるのでこの時の計算量は$O(N^2)$です。

<details>
<summary> Pythonでの実装例 </summary>

```py
S=[list(input()) for _ in range(9)]

pair=[]
for i in range(9):
  for j in range(9):
    if S[i][j]=="#":
      pair.append((i,j))

N=len(pair)

ans=0
judge=lambda n:n==int(n) and 0<=n<9 #値は整数でなくグリッド内にあるか
for i in range(N):
  for j in range(i+1,N):
    ax,ay=pair[i]; bx,by=pair[j]
    mx,my=(ax+bx)/2,(ay+by)/2 #Mの座標
    s,t=ax-mx,ay-my #MからAへ
    cx,cy=mx-t,my+s #C
    dx,dy=mx+t,my-s #D
    if judge(cx) and judge(cy) and judge(dx) and judge(dy):
      cx,cy,dx,dy=int(cx),int(cy),int(dx),int(dy)
      if S[cx][cy]=="#" and S[dx][dy]=="#":
        ans+=1

print(ans//2)
```
PyPyだと67msでした。はやい！
</details>
<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
#define rep(i,N,M) for(int i=N; i<M; i++)

bool judge(double n){
  return (n==(int)n) and 0<=n and n<9;
}

int main(){
  vector<vector<char>> S(9,vector<char>(9));
  
  vector<pair<int,int>> p;
  rep(i,0,9){
    rep(j,0,9){
      cin>>S[i][j];
      if(S[i][j]=='#'){
        p.push_back(make_pair(i,j));
      }
    }
  }
  
  int N=p.size();
  int ans=0;
  rep(i,0,N){
    rep(j,i+1,N){
      int ax,ay,bx,by;
      tie(ax,ay)=p[i]; tie(bx,by)=p[j];
      double mx=(ax+bx)/(double)2,my=(ay+by)/(double)2;
      double s=ax-mx,t=ay-my;
      double cx=mx-t,cy=my+s;
      double dx=mx+t,dy=my-s;
      if(judge(cx) and judge(cy) and judge(dx) and judge(dy)){
        cx=(int)cx; cy=(int)cy; dx=(int)dx; dy=(int)dy;
        if(S[cx][cy]=='#' and S[dx][dy]=='#'){
          ans++;
        }
      }
    }
  }
  
  cout<<ans/2<<endl;
}
```
7msでした。はやい。
</details>

## 1-5 D問題 Yet Another Recursive Function

### 問題

非負整数$x$に対し定義される関数$f(x)$を以下に定義します。

> $f(0)=1$
> $f(n)=f(\lfloor \tfrac{n}{2} \rfloor)+f(\lfloor \tfrac{n}{3} \rfloor ) \\　$

整数$N$が与えられるので、$f(N)$の値を求めて下さい。

#### 制約
・$1 \leq N \leq 10^{18}$

### 解説

愚直に、問題文にある式 $f(n)=f(\lfloor \tfrac{n}{2} \rfloor )+f(\lfloor \tfrac{n}{3} \rfloor ) \\　$を使って求めることを考えます。

これにはメモ化再帰が有効です。

#### メモ化再帰

メモ化再帰とは、既に呼び出した関数を再び呼び出されないようにし、指数オーダーの計算量を$O(N)$に圧縮するテクニックです。
DPみたいなものです。

#### メモ化再帰で間に合うか

メモ化再帰を使っても、関数の呼び出される回数が多ければ意味がありません。

$1$つの関数につき最大$2$種類の関数に分岐すると言っても、$f(x)$の引数としてあり得るものはそう多くない筈です。

$f(x)$の引数としてあり得るのは$ \lfloor \tfrac{N}{2^A3^B} \rfloor$($A,B$は非負整数)と表せる数です。
$A \leq \log_2 N$、$B \leq \log_3 N$なので、$f(x)$の引数としてあり得るのは$O(\log^2 N)$個程しかないことがわかります。

$N$は最大でも$10^{18}$ですので余裕で間に合います。

#### 最後に

メモ化再帰の実装についてですが、既に呼び出した関数であればその値を返し、そうでなければ値を求めメモするという風にすればよいです。

呼び出したかの判定には集合を扱う型を用いるとよいです。

$x=0$の場合を忘れないようにしましょう。

<details>
<summary> Pythonでの実装例 </summary>

```py
N=int(input())

memo={}
seen=set()
def calc(now):
  if now in seen:
    return memo[now]
  if now==0:
    memo[now]=1
    return 1
  seen.add(now)
  memo[now]=calc(now//2)+calc(now//3)
  return memo[now]

print(calc(N))
```
</details>
<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
using ll=long long;

map<ll,ll> memo;
set<ll> seen;

ll calc(ll now){
  if(seen.count(now)){
    return memo[now];
  }
  if(now==0LL){
    return memo[now]=1LL;
  }
  seen.insert(now);
  return memo[now]=calc(now/2)+calc(now/3);
}

int main(){
  ll N; cin>>N;
  cout<<calc(N)<<endl;
}
```
`long long`型を使わないと無理です。
</details>

このように、一見計算量が$O(N)$のように見えるものを$O(\log N)$だと気付くのは重要です。
頻出はしませんが覚えておきましょう。

## 1-6 E問題 Sugoroku 4

### 問題
$0$から$N$まで番号のついた$N+1$個のマスのある双六で遊んでいます。
マス$0$からスタートし、$N$を目指します。

この双六では$1$から$M$までの数が等確率で出るルーレットを使用します。
ルーレットを回して出た数だけ進みますが、マス$N$を越えて進むことになる場合は、その分マス$N$から引き返します。
マス$N$に到達したらゲームを終了します。

ルーレットを$K$回まで回せる時、ゴールできる確率を求めて下さい。

尚、確率は `mod 998244353`で出力して下さい。

#### 制約
・$M \leq N \leq 1000$
・$1 \leq M \leq 10$
・$1 \leq K \leq 1000$

### 解説

確率のDPの問題です。苦手な方も多いのではないでしょうか。私は期待値DPの方が苦手です。
(これがDPで解けるというのは説明していませんが大丈夫ですよね...極力詳しい説明をするよう心がけていますが限界はあるものなのです)

#### DPの定義

どのようにDPを定義すればいいか考えてみます。

今回のゲームに於いて、ルーレットを回す回数が制限となっています。
$i(1 \leq i \leq K)$回でゴールできる確率を求めたいです。
それだけでは無理なので、$i(1 \leq i \leq K)$回目にどこにいるかの情報も載せたいです。

ですのでDPはこのように定義すればよいです。

> $DP[i][j]=i回ルーレットを回した後にマスjにいる確率$

遷移は、$k$を$i$回目のルーレットで出た数として、$DP[i+1][j+k]$に$DP[i][j] ×\tfrac{1}{M}$を加算すればよいです。
マス$N$にゴールした以降はもうルーレットを回さない為、マス$N$にいる状態からの遷移は行わないことに注意して下さい。

求めるべき値は、$DP[i][N](1 \leq i \leq K)$の総和です。

#### mod 998244353

(は？小数の`mod`ってなんだよって方は[こちらの記事](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a)をお読みください。)

単に求めて出力するだけなら上を行うだけでよいのですが、答えを`mod 998244353`で出力せよとあります。

最後には`mod 998244353`の処理はできないので、一回毎の遷移で`mod`の処理をしましょう。
一回毎の遷移で$×\tfrac{1}{M}$をしますが、これは$M$の逆元を掛けるのと同じことです。$M$の逆元を既に求めておくと軽くなります。


<details>
<summary> mod 998244353の逆元の求め方</summary>
ここでは逆元の求め方を説明します。

$\bmod m$に於ける自然数$n$の逆元とは、$\bmod m$に於ける$1 ÷ n=n^{-1}$の値を指します。

以下に求め方を示します。
$998244353$は素数なので、フェルマーの小定理 $a^{p-1} \equiv 1 \bmod p$が成立ち、$a^{-1} \equiv a^{p-2} \bmod p$が成立ちます。

よって、$a^{p-2} \bmod p$を求めればよいです。

これは繰り返し二乗法を用いればよく、$O(\log mod)$で求められます。

以下に繰り返し二乗法のコードを貼っておきます。

<details>
<summary> Pythonでの実装例 </summary>

```py
def pow(x,n,m): #x^nをmで割った余り
  ret=1
  n%=m-1
  while n:
    if n%2:
      ret*=x
    x*=x; x%=m
    ret%=m
    n>>=1
  return ret
```

Pythonの場合標準ライブラリにありますが、PyPy3の場合AtCoderでのバージョンでは対応していないので、ライブラリとして持っていても良いのではないでしょうか。
</details>

(この逆元の求め方は`mod`が素数でないと動かないのでご注意。)

<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
using ll=long long;

ll pow(ll x,ll n,ll m){
  ll ret=1LL;
  n=(n+m-1)%m;
  while(n>0){
    if(n%2){
      ret*=x; ret%=m;
    }
    x*=x; x%=m;
    n>>=1;
  }
  return ret;
}
```

`long long`型で扱いましょう。
</details>
</details>

#### 最後に

`mod 998244353`の処理をしながら前述のDPをすればよいです。
DPは、$DP[i][j]$($1 \leq i \leq K,1 \leq j \leq N+1$)について$M$のループを回すので計算量は$O(NMK)$となります。
前述の通り、$M$の逆元の計算量については、$O(\log {mod})$となります。(今回の場合$mod=998244353$)
全体の計算量は$O(NMK+\log mod)$となり間に合います。

<details>
<summary> Pythonでの実装例 </summary>

```py
mo=998244353
def pow(x,n,m):
  ret=1
  n%=(m-1)
  while n:
    if n%2:
      ret*=x
    x*=x; x%=m
    ret%=m
    n>>=1
  return ret

N,M,K=map(int,input().split())

DP=[[0]*(N+1) for _ in range(K+1)]
DP[0][0]=1
po=pow(M,-1,mo)

for i in range(K):
  for j in range(N):
    for k in range(1,M+1):
      n=j+k
      if n>N:
        n=2*N-n
      DP[i+1][n]+=(DP[i][j]*po)%mo

ans=0
for i in range(K+1):
  ans+=DP[i][-1]
print(ans%mo)
```
</details>
<details>
<summary> C++での実装例 </summary>

```cpp
#include <bits/stdc++.h>
using namespace std;
using ll=long long;
#define rep(i,N,M) for(int i=N; i<M; i++)

const ll mo=998244353LL;

ll pow(ll x,ll n,ll m){
  ll ret=1LL;
  n=(n+m-1)%m;
  while(n>0){
    if(n%2){
      ret*=x; ret%=m;
    }
    x*=x; x%=m;
    n>>=1;
  }
  return ret;
}

int main(){
  int N,M,K; cin>>N>>M>>K;
  
  vector<vector<ll>> DP(K+1,vector<ll>(N+1,0LL));
  DP[0][0]=1LL;
  ll po=pow(M,-1,mo);
  
  rep(i,0,K){
    rep(j,0,N){
      rep(k,1,M+1){
        int n=j+k;
        if(n>N){
          n=2*N-n;
        }
        DP[i+1][n]+=DP[i][j]*po;
        DP[i+1][n]%=mo;
      }
    }
  }
  
  ll ans=0;
  rep(i,0,K+1){
    ans+=DP[i][N];
  }
  cout<<ans%mo<<endl;
}
```
</details>

ここまで書いて気付きましたが、確率を求めて一回毎にmodの処理をしなくても、全体の通り数に対するゴールできる通り数を求めればよいですね。

繰り返し二乗法のライブラリはよく使うのでライブラリに埋め込んでおいても良いでしょう。

# 2 さいごに
最後までお読み頂きありがとうございます。
今回は早解きをすれば良いPerfが取れた回だったように思います。
F問題の解説は時間があれば書きます。ご了承。
これで以上です。
ありがとうございました。
よい競プロライフを!!
~~いいねして頂ければモチベになります...~~