BOJ 31939 멀티버스를 여행하는 성재를 위한 안내서
sorohue가 PS하는 블로그
문제 링크입니다.
문제 요약
\[ {1 \over \pi R^2} \int_{0}^{R} \int_{0}^{2\pi} r \prod\limits_{k=1}^{N}{(r \cos \theta -x_k)^2 + (r\sin \theta-y_k)^2} \, d \theta \, dr \] 을 구하세요.
풀이
복소평면을 들고 옵시다. $(x_k, y_k) = x_k + y_ki = p_k$로, $(r\cos\theta,r\sin\theta) = re^{i\theta}$로 표현할 수 있습니다. 따라서 주어진 식을 다음과 같이 변형할 수 있습니다.
\[ {1 \over \pi R^2} \int_{0}^{R} \int_{0}^{2\pi} r\prod\limits_{k=1}^{N}{\vert re^{i\theta}-p_k\vert ^2} \, d\theta \, dr \] \[ ={1 \over \pi R^2} \int_{0}^{R} \int_{0}^{2\pi} r\prod\limits_{k=1}^{N}{re^{i\theta}-p_k}\prod\limits_{k=1}^{N}\overline{re^{i\theta}-p_k} \, d \theta \, dr \]
이때 각 $\prod$들을 $re^{i\theta}$에 대한 $N$차 복소다항식으로 여길 수 있습니다. 따라서, 주어진 식은
\[ {1 \over \pi R^2} \int _0 ^ R \int_0 ^ {2\pi} r\sum a_k r^k e^{ik\theta} \sum \overline{a_k} r^k e^{-ik\theta}\ d\theta dr \] 로 나타낼 수 있습니다.
$e^{i\theta}$는 $0 \le \theta \le 2\pi$에서 단위원을 그립니다. 이를 $\theta$에 대해 적분한 결과가 0이 됨은 쉽게 알 수 있습니다.
즉, $\sum a_k r^k e^{ik\theta} \sum a_k r^k e^{-ik\theta}$에서 양 $\sum$의 $k$ 값이 서로 다른 항끼리 곱해지는 경우는 적분값이 모두 0이 됩니다. 이로부터 우리가 주목해야 하는 항들만 모아 계산하면, 주어진 식은
\[ {1 \over \pi R^2} \int_0 ^ R \int_0 ^ {2\pi} \sum \vert a_k\vert ^2 r^{2k+1}\ d\theta dr \] \[ = {2 \over R^2} \int_0 ^ R \sum \vert a_k\vert ^2 r^{2k+1}\ dr = \sum { \vert a_k\vert ^2 R^{2k} \over k+1} \] 로 정리됩니다.
이제 $\prod (x-p_k)$의 계수를 구하면 됩니다. 복소수를 실수부와 허수부로 나누고, 실수부를 담을 공간을 한 칸 더 마련해 $3N$차 다항식처럼 여기고 FFT를 진행합니다. 분할 정복을 통해 다항식을 합쳐주면 시간 복잡도 ${\cal O} (N \lg^2 N)$에 모든 $a_k$를 구할 수 있습니다.
코드
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include<bits/stdc++.h>
#define mid (l+r>>1)
using namespace std;
using ll = long long;
using pll = pair<ll, ll>;
const ll mod = 998244353; // =(119<<23)+1
const ll w = 3;
ll x[101010], y[101010];
vector<ll> p[101010];
ll pw(ll n, ll r){
ll ret = 1;
while(r){
if(r&1) ret = ret*n%mod;
n = n*n%mod;
r >>= 1;
}
return ret;
}
inline ll inv(ll n){return pw(n, mod-2);}
void ntt(vector<ll>& f, bool flag = 0){
int n = f.size(), j = 0;
vector<ll> root(n>>1);
for(int i = 1; i < n; i++){
int bit = n>>1;
while(j >= bit){
j -= bit; bit >>= 1;
}
j += bit;
if(i < j) swap(f[i], f[j]);
}
ll ang = pw(w, (mod - 1) / n); if(flag) ang = inv(ang);
root[0] = 1; for(int i=1; i<(n >> 1); i++) root[i] = root[i-1] * ang % mod;
for(int i=2; i<=n; i<<=1){
int step = n / i;
for(int j=0; j<n; j+=i){
for(int k=0; k<(i >> 1); k++){
ll u = f[j | k], v = f[j | k | i >> 1] * root[step * k] % mod;
f[j | k] = (u + v) % mod;
f[j | k | i >> 1] = (u - v) % mod;
if(f[j | k | i >> 1] < 0) f[j | k | i >> 1] += mod;
}
}
}
ll t = inv(n);
if(flag) for(int i=0; i<n; i++) f[i] = f[i] * t % mod;
}
void mult(int l, int r){
vector<ll>& a = p[l]; vector<ll>& b = p[r];
int sz = a.size()+b.size()-3;
int n = 2; while(n < a.size()+b.size()) n <<= 1;
a.resize(n); b.resize(n); ntt(a); ntt(b);
for(int i = 0; i < n; i++) a[i] = a[i]*b[i]%mod;
ntt(a, 1);
a.resize(sz);
for(int i = 2; i < n; i+=3){
a[i-2] += (mod-a[i])%mod;
a[i-2] %= mod;
a[i] = 0;
}
return;
}
void solve(int l, int r){
if(l == r){
p[l] = {(mod-x[l])%mod, (mod-y[l])%mod, 0, 1, 0, 0};
return;
}
solve(l, mid);
solve(mid+1, r);
return mult(l, mid+1);
}
int main(){
cin.tie(0); cout.tie(0); ios::sync_with_stdio(false);
ll n, R; cin >> n >> R; R = R*R%mod;
for(int i = 1; i <= n; i++) cin >> x[i] >> y[i];
solve(1, n); vector<ll>& ans = p[1];
ll ret = 0, r = 1;
for(int i = 0; i <= n; i++){
ret += inv(i+1)*r%mod*(ans[3*i]*ans[3*i]%mod+ans[3*i+1]*ans[3*i+1]%mod)%mod;
ret %= mod;
r = r*R%mod;
}
cout << ret;
}