Relaxed Convolution (2)
개요
이전에 작성한 글에서는 Relaxed Convolution의 개념과 성질, 구현 방법에 대해서 다루었습니다. 이 글에서는 Relaxed Convolution을 Problem Solving에 활용하는 방법을 다룹니다.
연습 문제
AtCoder Beginner Contest 213 H. Stroll
문제
정점이 $N$개이고 간선이 매우 많은 무방향 그래프가 주어집니다.
간선 집합은 다음과 같은 방식으로 정의됩니다.
- $0 \leq i < M, 1 \leq d \leq T$를 만족하는 모든 $(i,d)$ 쌍에 대해, 두 정점 쌍 $(a_i, b_i)$를 잇는 길이가 $d$인 간선이 $p_{i,d}$개 존재합니다.
이 그래프의 $1$번 정점에서 출발해서 $1$번 정점으로 되돌아오는, 길이가 정확히 $T$인 walk의 개수를 $998244353$으로 나눈 나머지를 출력해야 합니다.
$(2 \leq N \leq 10; 1 \leq M \leq \min \left(10, \frac{N(N-1)}{2} \right); 1 \leq T \leq 40000; 0 \leq p_{i,j} < 998244353)$
풀이
$dp[u][k]$를, $1$번 정점에서 출발해서 $u$번 정점에 도착하는 길이가 정확히 $k$인 walk의 개수라고 정의합시다.
모든 간선의 길이가 1 이상이므로 상태 전이에는 사이클이 없고, 점화식을 아래와 같이 작성할 수 있습니다.
$dp[u][k] = \sum_{0 \leq i < M, {a_i, b_i} = {u, v}} \sum_{1 \leq d \leq T} dp[v][k-d] \cdot p_{i,d}$
위 점화식을 단순히 계산하면 $O(MT^2)$로 시간 초과가 발생하므로, $dp[u]$의 상태 전이가 $dp[v]$와 $p_{i}$의 Convolution임을 이용해야 합니다.
$dp[u][k]$를 구할 때 $dp[v][0 \dots k-1]$이 필요하고, $dp[v][k]$를 구할 때에도 $dp[u][0 \dots k-1]$이 필요하므로 Convolution을 바로 적용할 수는 없지만, $k=1,2,\cdots,T$ 순서대로 $dp[1 \dots N][k]$을 계산한다면 Semi-Relaxed Convolution을 적용할 수 있습니다.
최종 시간복잡도는 $O(MT \log^2 T)$가 됩니다.
코드
코드는 아래와 같습니다.
#include<bits/stdc++.h>
#include<atcoder/convolution>
#include<atcoder/modint>
using namespace std;
using namespace atcoder;
using modular = modint998244353;
#define INF 1234567890
#define ll long long
int N, M, T;
int A[10], B[10];
vector<modular> p[10];
vector<modular> dp[11];
void dnc(int l, int r)
{
if (l == r) return;
int mid = (l+r)/2;
dnc(l, mid);
for(int i=0; i<M; i++) // A[i] -> B[i]로 이동하는 경우
{
auto AA = vector<modular>(p[i].begin() + 0, p[i].begin() + r-l + 1);
auto BB = vector<modular>(dp[A[i]].begin() + l, dp[A[i]].begin() + mid + 1);
auto CC = convolution(AA, BB);
for(int j=mid+1; j<=r; j++)
dp[B[i]][j] += CC[j-l];
}
for(int i=0; i<M; i++) // B[i] -> A[i]로 이동하는 경우
{
auto AA = vector<modular>(p[i].begin() + 0, p[i].begin() + r-l + 1);
auto BB = vector<modular>(dp[B[i]].begin() + l, dp[B[i]].begin() + mid + 1);
auto CC = convolution(AA, BB);
for(int j=mid+1; j<=r; j++)
dp[A[i]][j] += CC[j-l];
}
dnc(mid+1, r);
}
int main()
{
ios::sync_with_stdio(0); cin.tie(0);
cin.exceptions(ios::badbit | ios::failbit);
cin >> N >> M >> T;
for(int i=0; i<M; i++)
{
cin >> A[i] >> B[i];
p[i].resize(T+1);
for(int d=1,x; d<=T; d++)
cin >> x, p[i][d] = x;
}
for(int i=1; i<=N; i++)
dp[i].resize(T+1);
dp[1][0] = 1; // base case
dnc(0, T);
cout << dp[1][T].val() << "\n";
return 0;
}
AtCoder Beginner Contest 315 Ex. Typical Convolution Problem
문제
수열 $(A_1, A_2, \cdots, A_N)$이 주어지면, 다음과 같이 정의되는 수열 $(F_0, F_1, \cdots, F_N)$의 각 원소를 $998244353$으로 나눈 나머지를 구해야 합니다.
-
$F_0 = 1$
-
$1 \leq n$에서 $F_n = A_n \sum_{i+j < n} F_i F_j$
$(1 \leq N \leq 200,000; 0 \leq A_i < 998244353)$
풀이
$F$의 누적합 수열인 $S_n = \sum_{0 \leq i \leq n} F_i$를 정의합시다.
수식을 다시 풀어 적으면, $F_n = A_n \sum_{0 \leq i \leq n-1} F_i S_{n-1-i}$이 됩니다.
$F$와 $S$의 Relaxed Convolution을 구하면 $O(N \log^2 N)$에 해결할 수 있습니다.
코드
$F_i$와 $S_j$를 곱한 결과물을 $F_{i+j}$에 더하는 것이 아니라 $F_{i+j+1}$에 더해야 한다는 디테일에 주의해서 구현해야 합니다.
코드는 아래와 같습니다.
#include<bits/stdc++.h>
#include<atcoder/convolution>
#include<atcoder/modint>
using namespace std;
using namespace atcoder;
using modular = modint998244353;
#define INF 1234567890
#define ll long long
int N;
int A[200001];
vector<modular> F, S;
void dnc(int l, int r)
{
if (l == r) // F[r+1], S[r+1]을 구한다.
{
F[r+1] += F[r] * S[0];
if (r) F[r+1] += F[0] * S[r];
F[r+1] *= A[r+1];
S[r+1] = S[r] + F[r+1];
return;
}
int mid = (l+r)/2;
dnc(l, mid);
if (l == 0)
{
auto CC = convolution(vector<modular>(F.begin() + l, F.begin() + mid + 1), vector<modular>(S.begin() + l, S.begin() + mid + 1));
for(int i=mid+1; i<=r; i++)
F[i+1] += CC[i-l];
}
else
{
auto AB = convolution(vector<modular>(F.begin() + 0, F.begin() + r-l + 1), vector<modular>(S.begin() + l, S.begin() + mid + 1));
auto BA = convolution(vector<modular>(S.begin() + 0, S.begin() + r-l + 1), vector<modular>(F.begin() + l, F.begin() + mid + 1));
for(int i=mid+1; i<=r; i++)
F[i+1] += AB[i-l] + BA[i-l];
}
dnc(mid+1, r);
}
int main()
{
ios::sync_with_stdio(0); cin.tie(0);
cin.exceptions(ios::badbit | ios::failbit);
cin >> N;
for(int i=1; i<=N; i++)
cin >> A[i];
F.resize(N+1), S.resize(N+1);
F[0] = 1, S[0] = 1; // base case
dnc(0, N-1);
for(int i=1; i<=N; i++)
cout << F[i].val() << " ";
cout << "\n";
return 0;
}
삼각분할
일반적으로 알려진 풀이와 다르게 Relaxed Convolution으로도 해결할 수 있습니다.
문제
$N$과 $M$이 주어집니다.
$i$각형을 삼각분할하는 경우의 수를 $T_i$라 했을 때, $T_3 + T_4 + \cdots + T_N$을 $M$으로 나눈 나머지를 구해야 합니다.
$(3 \leq N \leq 100,000; 2 \leq M \leq 10^9)$
풀이
$C_i$를 $i$번째 카탈란 수라고 하면 $T_i = C_{i-2}$임이 잘 알려져 있습니다. 따라서 $C_1 + C_2 + \cdots + C_{N-2} \pmod M$을 구하면 됩니다.
카탈란 수의 점화식은 아래와 같습니다.
-
$C_0 = 1$
-
$C_i = \sum_{0 \leq j \leq i-1} C_i C_{i-1-j}$
Relaxed Convolution을 적용할 수 있는 형태이므로 $O(N \log^2 N)$에 해결할 수 있습니다.
코드
$M$이 임의로 주어지기 때문에 Convolution 구현체를 잘 선택해야 하는데, 저는 Aeren님의 코드와 justiceHui님의 코드를 적당히 섞어서 사용했습니다.
코드는 아래와 같습니다.
#include<bits/stdc++.h>
using namespace std;
#define INF 1234567890
#define ll long long
using real_t = double; using cpx = complex<real_t>;
void FFT(vector<complex<double>> &a, bool invert = false){
int n = (int)a.size();
static vector<complex<double>> root(2, 1);
static vector<complex<long double>> root_ld(2, 1);
for(static int k = 2; k < n; k <<= 1){
root.resize(n), root_ld.resize(n);
auto theta = polar(1.0L, acos(-1.0L) / k);
for(auto i = k; i < k << 1; ++ i) root[i] = root_ld[i] = i & 1 ? root_ld[i >> 1] * theta : root_ld[i >> 1];
}
for(auto i = 1, j = 0; i < n; ++ i){
int bit = n >> 1;
for(; j & bit; bit >>= 1) j ^= bit;
j ^= bit;
if(i < j) swap(a[i], a[j]);
}
for(auto k = 1; k < n; k <<= 1) for(auto i = 0; i < n; i += k << 1) for(auto j = 0; j < k; ++ j){
auto x = (double *)&root[j + k], y = (double *)&a[i + j + k];
complex<double> z(x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]);
a[i + j + k] = a[i + j] - z, a[i + j] += z;
}
if(invert){
reverse(a.begin() + 1, a.end());
auto inv_n = 1.0 / n;
for(auto &x: a) x *= inv_n;
}
}
vector<ll> multiply_mod(const vector<ll> &a, const vector<ll> &b, const ll mod){
int N = 2; while(N < a.size() + b.size()) N <<= 1;
vector<cpx> v1(N), v2(N), r1(N), r2(N);
for(int i=0; i<a.size(); i++) v1[i] = cpx(a[i] >> 15, a[i] & 32767);
for(int i=0; i<b.size(); i++) v2[i] = cpx(b[i] >> 15, b[i] & 32767);
FFT(v1); FFT(v2);
for(int i=0; i<N; i++){
int j = i ? N-i : i;
cpx ans1 = (v1[i] + conj(v1[j])) * cpx(0.5, 0);
cpx ans2 = (v1[i] - conj(v1[j])) * cpx(0, -0.5);
cpx ans3 = (v2[i] + conj(v2[j])) * cpx(0.5, 0);
cpx ans4 = (v2[i] - conj(v2[j])) * cpx(0, -0.5);
r1[i] = (ans1 * ans3) + (ans1 * ans4) * cpx(0, 1);
r2[i] = (ans2 * ans3) + (ans2 * ans4) * cpx(0, 1);
}
vector<ll> ret(N); FFT(r1, true); FFT(r2, true);
for(int i=0; i<N; i++){
ll av = llround(r1[i].real()) % mod;
ll bv = ( llround(r1[i].imag()) + llround(r2[i].real()) ) % mod;
ll cv = llround(r2[i].imag()) % mod;
ret[i] = (av << 30) + (bv << 15) + cv;
ret[i] %= mod; ret[i] += mod; ret[i] %= mod;
}
while(ret.size() > 1 && ret.back() == 0) ret.pop_back();
return ret;
}
int N, M;
vector<ll> C;
void dnc(int l, int r)
{
if (l == r) // C[r+1]을 구한다.
{
(C[r+1] += C[r] * C[0]) %= M;
if (r) (C[r+1] += C[0] * C[r]) %= M;
return;
}
int mid = (l+r)/2;
dnc(l, mid);
if (l == 0)
{
auto CC = multiply_mod(vector<ll>(C.begin() + l, C.begin() + mid + 1), vector<ll>(C.begin() + l, C.begin() + mid + 1), M);
CC.resize(r-l+1);
for(int i=mid+1; i<=r; i++)
(C[i+1] += CC[i-l]) %= M;
}
else
{
auto CC = multiply_mod(vector<ll>(C.begin() + 0, C.begin() + r-l + 1), vector<ll>(C.begin() + l, C.begin() + mid + 1), M);
CC.resize(r-l+1);
for(int i=mid+1; i<=r; i++)
(C[i+1] += CC[i-l] * 2) %= M;
}
dnc(mid+1, r);
}
int main()
{
ios::sync_with_stdio(0); cin.tie(0);
cin.exceptions(ios::badbit | ios::failbit);
cin >> N >> M;
N -= 2;
C.resize(N+1);
C[0] = 1; // base case
dnc(0, N-1);
cout << accumulate(C.begin() + 1, C.end(), 0LL) % M << "\n";
return 0;
}
Codeforces Round 309 (Div. 1) E. Kyoya and Train
문제
정점이 $N$개, 간선이 $M$개인 단순 방향 그래프가 주어지고 두 정수 $T,X$가 주어집니다.
$0 \leq i < M$에 대해 간선 $i$는 정점 $a_i$에서 정점 $b_i$로 비용 $c_i$를 지불하고 이동할 수 있습니다. 이 때 이동에 걸리는 시간은 확률 분포로 주어지고 매 이동마다 랜덤하게 결정되는데, 정수 $1 \leq d \leq T$에 대해 $d$의 시간이 걸릴 확률은 $p_{i,d}$입니다. $(0 \leq p_{i,d} \leq 1; \sum_{1 \leq d \leq T} p_{i,d} = 1)$
이 그래프에서 매번 어떤 간선을 따라 이동할지를 최적으로 선택해서, $1$번 정점에서 출발해서 $N$번 정점에 도착하는 비용의 기댓값을 최소화해야 합니다. 단, 걸린 시간이 $T$를 초과하면 지각이므로 비용을 $X$만큼 추가로 지불해야 합니다.
$(2 \leq N \leq 50; 1 \leq M \leq 100; 1 \leq T \leq 20000; 0 \leq X \leq 10^6)$
풀이
위 연습 문제인 Stroll과 전반적으로 비슷한 문제입니다.
$dp[u][k]$를, 앞으로 $k$ 이상의 시간을 소모하면 지각일 때 $u$번 정점에서 출발해서 $N$번 정점에 도착하는 비용의 기댓값의 최솟값이라고 정의합시다.
$dp[u][k] = \min_{0 \leq i < M, (a_i,b_i) = (u,v)} \left( \left(\sum_{1 \leq d \leq T} dp[v][\max(0, k-d)] \cdot p_{i,d} \right) + c_i \right)$
위 점화식을 단순히 계산하면 $O(MT^2)$로 시간 초과가 발생하므로, $dp[u]$의 상태 전이가 $dp[v]$와 $p_i$의 Convolution임을 이용해야 합니다.
$ep[i][k] = \sum_{1 \leq d \leq k-1} dp[b_i][k-d] \cdot p_{i,d}$를 새로 정의하고 수식 정리를 합시다.
$dp[u][k] = \min_{0 \leq i < M, (a_i,b_i) = (u,v)} \left(ep[i][k] + dp[v][0] \cdot \left( \sum_{k \leq d \leq T} p_{i,d} \right) + c_i \right)$
이제 $ep$가 Semi-Relaxed Convolution이 가능하므로 $O(MT \log^2 T)$에 해결할 수 있습니다.
$k = 0$일 때에만 상태 전이에 사이클이 존재하므로, 이 경우는 다익스트라 등의 알고리즘으로 따로 전처리해야 합니다.
코드
$N$ 제한이 매우 작아서 구현의 편의를 위해 플로이드-워셜 알고리즘을 사용했습니다.
이 문제에서는 기댓값을 double로 구해야 하기 때문에, Convolution 구현체로 AtCoder Library 대신에 Aeren님의 코드와 justiceHui님의 코드를 적당히 섞어서 사용했습니다.
코드는 아래와 같습니다.
#include<bits/stdc++.h>
using namespace std;
#define INF 1e18
#define ll long long
using real_t = double; using cpx = complex<real_t>;
void FFT(vector<complex<double>> &a, bool invert = false){
int n = (int)a.size();
static vector<complex<double>> root(2, 1);
static vector<complex<long double>> root_ld(2, 1);
for(static int k = 2; k < n; k <<= 1){
root.resize(n), root_ld.resize(n);
auto theta = polar(1.0L, acos(-1.0L) / k);
for(auto i = k; i < k << 1; ++ i) root[i] = root_ld[i] = i & 1 ? root_ld[i >> 1] * theta : root_ld[i >> 1];
}
for(auto i = 1, j = 0; i < n; ++ i){
int bit = n >> 1;
for(; j & bit; bit >>= 1) j ^= bit;
j ^= bit;
if(i < j) swap(a[i], a[j]);
}
for(auto k = 1; k < n; k <<= 1) for(auto i = 0; i < n; i += k << 1) for(auto j = 0; j < k; ++ j){
auto x = (double *)&root[j + k], y = (double *)&a[i + j + k];
complex<double> z(x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]);
a[i + j + k] = a[i + j] - z, a[i + j] += z;
}
if(invert){
reverse(a.begin() + 1, a.end());
auto inv_n = 1.0 / n;
for(auto &x: a) x *= inv_n;
}
}
vector<double> multiply(const vector<double> &_a, const vector<double> &_b){
vector<cpx> a(_a.begin(), _a.end()), b(_b.begin(), _b.end());
int N = 2; while(N < a.size() + b.size()) N <<= 1;
a.resize(N); b.resize(N); FFT(a); FFT(b);
for(int i=0; i<N; i++) a[i] *= b[i];
vector<double> ret(N); FFT(a, 1);
for(int i=0; i<N; i++) ret[i] = a[i].real(); //
while(ret.size() > 1 && ret.back() == 0) ret.pop_back();
return ret;
}
int N, M, T, X;
int A[100], B[100], C[100];
vector<double> p[100], P[100], dp[51], ep[100];
void dnc(int l, int r)
{
if (l == r) // dp[u][r] 구하기
{
for(int i=0; i<M; i++)
dp[A[i]][r] = min(dp[A[i]][r], ep[i][r] + dp[B[i]][0] * (P[i][T] - P[i][r-1]) + C[i]);
return;
}
int mid = (l+r)/2;
dnc(l, mid);
for(int i=0; i<M; i++)
{
auto AA = vector<double>(p[i].begin() + 0, p[i].begin() + r-l + 1);
auto BB = vector<double>(dp[B[i]].begin() + l, dp[B[i]].begin() + mid + 1);
auto CC = multiply(AA, BB);
CC.resize(r-l+1);
for(int k=mid+1; k<=r; k++)
ep[i][k] += CC[k-l];
}
dnc(mid+1, r);
}
int main()
{
ios::sync_with_stdio(0); cin.tie(0);
cin.exceptions(ios::badbit | ios::failbit);
cout << setprecision(12) << fixed;
cin >> N >> M >> T >> X;
for(int i=0; i<M; i++)
{
cin >> A[i] >> B[i] >> C[i];
p[i].resize(T+2);
P[i].resize(T+1);
for(int d=1,x; d<=T; d++)
{
cin >> x, p[i][d] = x / 100000.0;
P[i][d] = P[i][d-1] + p[i][d]; // 누적 합
}
ep[i].resize(T+2);
}
for(int i=1; i<=N; i++)
dp[i].resize(T+2, INF);
for(int i=1; i<=T+1; i++)
dp[N][i] = 0;
// dp[u][0] 전처리하기 (플로이드)
vector<vector<double> > dist(N+1, vector<double>(N+1, INF));
for(int i=1; i<=N; i++)
dist[i][i] = 0;
for(int i=0; i<M; i++)
dist[A[i]][B[i]] = C[i];
for(int k=1; k<=N; k++)
for(int i=1; i<=N; i++)
for(int j=1; j<=N; j++)
dist[i][j] = min(dist[i][j], dist[i][k] + dist[k][j]);
for(int i=1; i<=N; i++)
dp[i][0] = dist[i][N] + X;
dnc(1, T+1);
cout << dp[1][T+1] << "\n";
return 0;
}
Gennady Korotkevich Contest 5. Bin
문제
$N,K$가 주어지면, 다음의 두 조건을 만족하는 리프가 $N$개인 이진 트리의 개수를 $998244353$으로 나눈 나머지를 구해야 합니다.
-
모든 정점은 0개 또는 2개의 자식을 갖는다.
-
2개의 자식을 갖는 모든 정점에 대해, (왼쪽 서브트리의 리프 개수) $\leq$ (오른쪽 서브트리의 리프 개수) + $K$를 만족한다.
$(2 \leq N \leq 10^6, 0 \leq K \leq 100)$
풀이
$D_n$을 위의 두 조건을 만족하는 리프가 $n$개인 이진 트리의 개수라고 정의하면, 간단하게 아래 점화식을 찾을 수 있습니다.
-
$D_0 = 0, D_1 = 1$
-
$i \geq 2$에서 $D_i = \sum_{1 \leq l+r \leq i-1; \ l \leq r+K} D_l D_r$
수식 정리를 해서 문제를 풀기 좋게 변형해 봅시다.
점화식을 다시 작성하면 $D_i = \sum_{1 \leq j \leq \min(i-1, \frac{i+K}{2})} D_j D_{i-j}$이 됩니다.
$i \leq K+1$에 대해 나이브하게 계산하면 $i > K+1$부터는 $D_i = \sum_{1 \leq j \leq \frac{i+K}{2} } D_j D_{i-j}$이 됩니다.
$j > \frac{i}{2}$를 나이브하게 계산한다고 치면, 결국 $D_i = \sum_{1 \leq j \leq \frac{i}{2} } D_j D_{i-j}$를 구하기만 해도 됩니다.
마지막으로 parity에 따라 케이스를 나누면,
-
$i$가 홀수일 때 $D_i = \left( \sum_{1 \leq j \leq i-1} D_j D_{i-j} \right) / 2$
-
$i$가 짝수일 때 $D_i = \left( (D_{i/2})^2 + \sum_{1 \leq j \leq i-1} D_j D_{i-j} \right) / 2$
이제 Relaxed Convolution으로 문제를 해결할 수 있습니다.
시간복잡도는 $O(K^2 + NK + N \log^2 N)$이 됩니다.
코드
코드는 아래와 같습니다.
#include<bits/stdc++.h>
#include<atcoder/convolution>
#include<atcoder/modint>
using namespace std;
using namespace atcoder;
using modular = modint998244353;
#define INF 1234567890
#define ll long long
int N, K;
vector<modular> D;
void dnc(int l, int r)
{
if (l == r)
{
int i = r;
if (i > K+1)
{
// parity에 따라 케이스 나누기
if (i % 2 == 0) D[i] += D[i/2] * D[i/2];
D[i] /= 2;
// j > i/2 나이브
for(int j=i/2+1; j<=(i+K)/2; j++)
D[i] += D[j] * D[i-j];
}
return;
}
int mid = (l+r)/2;
dnc(l, mid);
if (l == 0)
{
auto DD = convolution(vector<modular>(D.begin() + l, D.begin() + mid + 1), vector<modular>(D.begin() + l, D.begin() + mid + 1));
DD.resize(DD.size() + 1);
for(int i=mid+1; i<=r; i++)
if (i > K+1)
D[i] += DD[i-l];
}
else
{
auto DD = convolution(vector<modular>(D.begin() + 0, D.begin() + r-l + 1), vector<modular>(D.begin() + l, D.begin() + mid + 1));
for(int i=mid+1; i<=r; i++)
if (i > K+1)
D[i] += DD[i-l] * 2;
}
dnc(mid+1, r);
}
int main()
{
ios::sync_with_stdio(0); cin.tie(0);
cin.exceptions(ios::badbit | ios::failbit);
cin >> N >> K;
D.resize(max(N+1, K+2));
D[0] = 0;
D[1] = 1;
// 작은 케이스 전처리
for(int i=2; i<=K+1; i++)
for(int j=1; j<=min(i-1, (i+K)/2); j++)
D[i] += D[j] * D[i-j];
dnc(0, N);
cout << D[N].val() << "\n";
return 0;
}