快速数论变换(NTT)原理及实现
概述 上次写了一篇狗屎文章快速傅里叶变换(FFT)原理与实现 | KSkun’ …
May all the beauty be blessed.
题目地址:Codeforces:Problem – 364D – Codeforces、洛谷:【CF364D】Ghd – 洛谷
John Doe offered his sister Jane Doe find the gcd of some set of numbers a.
Gcd is a positive integer g, such that all number from the set are evenly divisible by g and there isn’t such g’ (g’ > g), that all numbers of the set are evenly divisible by g’.
Unfortunately Jane couldn’t cope with the task and John offered her to find the ghd of the same subset of numbers.
Ghd is a positive integer g, such that at least half of numbers from the set are evenly divisible by g and there isn’t such g’ (g’ > g) that at least half of the numbers from the set are evenly divisible by g’.
Jane coped with the task for two hours. Please try it, too.
有一个数列,求最大的数,使得该数是数列中一半以上的数的因数。
输入格式:
The first line contains an integer n (1 ≤ n ≤ 10^6) showing how many numbers are in set a. The second line contains space-separated integers a1, a2, …, an (1 ≤ ai ≤ 10^12). Please note, that given set can contain equal numbers.
Please, do not write the %lld specifier to read or write 64-bit integers in С++. It is preferred to use the %I64d specifier.
输出格式:
Print a single integer g — the Ghd of set a.
输入样例#1:
6 6 2 3 4 5 6
输出样例#1:
3
输入样例#2:
5 5 5 6 10 15
输出样例#2:
5
参考资料:CFR 364 D. Ghd ( Random, Math ) – 0w1
随机化的题目,直接求显然不好办,n的数据范围看着就吓人,直观感觉像O(n \log n)。
我们知道这个数字一定是这个数列中某个数的因数,因此我们需要选择一个数求它与其他数的最大公因数。这些公因数中的一个就是答案。求一遍公因数是O(n \log n)的。
对于求出来的公因数,我们去从大到小找一个会成为超过一半数的因数的数字。具体做法是,选择一个因数,去找比它大的因数,如果它能整除大因数,说明大因数对应的数字也可以被这个数整除,应当把加到这个数的计数上。这个过程直观看上去是O(n^2)的,但是实际上我们不会对比当前最优解还小的因数计算,并且当计数超过n的一半的时候就可以停止计数并计入答案,在这样的处理下我们可以让这个计数的时间变得可接受。
答案肯定是这些数字中超过一半数的因数,因此每一次随机找到正确解的概率超过0.5。我们考虑重复找解的过程多次,这里我使用10次(超过10次似乎会TLE 18)。
// Code by KSkun, 2018/3
#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <map>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF
: *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
const int MAXN = 1000005;
inline LL gcd(LL a, LL b) {
LL t;
while(b) {
t = a % b;
a = b;
b = t;
}
return a;
}
int n;
LL a[MAXN];
int main() {
srand(time(NULL));
n = readint();
for(int i = 1; i <= n; i++) {
a[i] = readint();
}
LL ans = 1;
for(int rep = 1; rep <= 10; rep++) { // bad rand
int rnd = (rand() * RAND_MAX + rand()) % n + 1;
std::map<LL, int> fact;
for(int i = 1; i <= n; i++) {
LL t = gcd(a[i], a[rnd]);
if(!fact.count(t)) fact[t] = 1;
else fact[t]++;
}
std::map<LL, int>::iterator it = fact.end();
do {
it--;
if((*it).first <= ans) continue;
int cnt = 0;
for(std::map<LL, int>::iterator it1 = it;
it1 != fact.end() && cnt << 1 < n; it1++) {
if(!((*it1).first % (*it).first)) {
cnt += (*it1).second;
}
}
if(cnt << 1 >= n) ans = (*it).first;
} while(it != fact.begin());
}
printf("%I64d", ans);
return 0;
}
题目地址:洛谷:【P3343】[ZJOI2015]地震后的幻想乡 – 洛谷、BZOJ:Problem 3925. — [Zjoi2015]地震后的幻想乡
傲娇少女幽香是一个很萌很萌的妹子,而且她非常非常地有爱心,很喜欢为幻想乡的人们做一些自己力所能及的事情来帮助他们。 这不,幻想乡突然发生了地震,所有的道路都崩塌了。现在的首要任务是尽快让幻想乡的交通体系重新建立起来。
幻想乡一共有n个地方,那么最快的方法当然是修复n-1条道路将这n个地方都连接起来。 幻想乡这n个地方本来是连通的,一共有m条边。现在这m条边由于地震的关系,全部都毁坏掉了。每条边都有一个修复它需要花费的时间,第i条边所需要的时间为ei。地震发生以后,由于幽香是一位人生经验丰富,见得多了的长者,她根据以前的经验,知道每次地震以后,每个ei会是一个0到1之间均匀分布的随机实数。并且所有ei都是完全独立的。
现在幽香要出发去帮忙修复道路了,她可以使用一个神奇的大魔法,能够选择需要的那n-1条边,同时开始修复,那么修复完成的时间就是这n-1条边的ei的最大值。当然幽香会先使用一个更加神奇的大魔法来观察出每条边ei的值,然后再选择完成时间最小的方案。 幽香在走之前,她想知道修复完成的时间的期望是多少呢?
输入格式:
第一行两个数n,m,表示地方的数量和边的数量。其中点从1到n标号。 接下来m行,每行两个数a,b,表示点a和点b之间原来有一条边。 这个图不会有重边和自环。
输出格式:
一行输出答案,四舍五入保留6位小数。
输入样例#1:
5 4 1 2 1 5 4 3 5 3
输出样例#1:
0.800000
提示:
(以下内容与题意无关,对于解题也不是必要的。)
对于n个[0,1]之间的随机变量x1,x2,…,xn,第k小的那个的期望值是k/(n+1)。
样例解释:
对于第一个样例,由于只有4条边,幽香显然只能选择这4条,那么答案就是4条边的ei中最大的数的期望,由提示中的内容,可知答案为0.8。
数据范围:
对于所有数据:n<=10, m<=n(n-1)/2, n,m>=1。
对于15%的数据:n<=3。
另有15%的数据:n<=10, m=n。
另有10%的数据:n<=10, m=n(n-1)/2。
另有20%的数据:n<=5。
另有20%的数据:n<=8。
题目要求最小生成树上的最大边权Y的期望,因此我们只需要关心这个最大边权在所有边权中的排名,我们发现有下面的式子
\mathrm{E}(Y) = \sum_{i=1}^n \mathrm{E}(i)\mathrm{P}(i) = \sum_{i=1}^n \frac{i}{m+1} \cdot \mathrm{P}(i) = \frac{\sum_{i=1}^n i\mathrm{P}(i)}{m+1} = \frac{\mathrm{E}(Z)}{m+1}
最后我们发现要求的就是最小生成树上最大边权排名的期望值Z。考虑怎么求这个Z,设L为同定义的随机变量,有
\mathrm{E}(Z) = \sum_{i=1}^m i\mathrm{P}(L=i) = \sum_{i=1}^m \mathrm{P}(L \geq i)
如果说最大边排名大于等于i不好求,我们就正难则反,求排名小于i的边集无法构成生成树的概率,由于所有的边权都是随机分布的变量,这个概率可以转化为从边集中随机选择i条边无法构成生成树的概率,我们考虑计算出选择i条边无法构成生成树的方案数,进而计算概率。
通过观察题目数据范围,我们发现n的范围很适合状压。考虑状压DP,用f[S][i]表示在S点集构成的子图中,选i条边使得该点集不连通的方案数,这个状态似乎没法转移,但是我们考虑与其意义互补的量:g[S][i]表示在S点集构成的子图中,选i条边使得该点集连通的方案数,显然有
f[S][i] + g[S][i] = \mathrm{C}_{ecnt[S]}^i
ecnt表示该点集构成子图内的边数。如果知道了其中一个,我们就可以知道另外一个。那么怎么转移成了问题,我们考虑固定S中的某个点,对于S这个点集的一个包含定点的真子集T,只要使T连通但\complement_S T与T之间没有连边就能保证S不连通了,而枚举每一个T可以实现对S不连通情况的遍历。最后我们的转移式就是
f[S][i] = \sum_{T \subsetneqq S, U \in T} \sum_{j=0}^{ecnt[T]} g[T][j] \cdot \mathrm{C}_{ecnt[\complement_S T]}^{i-j}
最后,我们获得了 \mathrm{E}(Z) = \sum_{i=0}^{n-1} f[\mathbb{U}][i] ,直接计算答案即可。
注意爆int。
需要证明的内容:对于n个[0,1]之间的随机变量x1,x2,…,xn,第k小的那个的期望值是k/(n+1)。
我们首先令x_1为第k小值,我们现在来求它的期望。
对于这种情况,我们显然有其中k-1个值比它小,n-k个值比它大,我们分别计算概率,乘起来,就得到了概率密度,对其积分就是x_1为第k小值时x_1的期望值,即\mathrm{E}(x_1, x_1是第k小数)
\begin{aligned} & \int_0^1 x_1 \cdot x_1^{k-1} (1-x_1)^{n-k} \mathrm{d}x_1 \\ = & \int_0^1 x_1^k (1-x_1)^{n-k} \mathrm{d}x_1 \\ = & \mathrm{B}(k+1, n-k+1) \\ = & \frac{k!(n-k)!}{(n+1)!} \end{aligned}
我们知道了这个,但是要求的是E(第k小数),即\mathrm{E}(x_1|x_1是第k小数)。我们记x_1是第k小数这一事件为A,根据下面的关系
\mathrm{E}(x_1|A) = \frac{\mathrm{E}(x_1, A)}{\mathrm{P}(A)}
我们又容易知道
\mathrm{P}(A) = \frac{1}{n \cdot \mathrm{C}_{n-1}^{k-1}}
直接就可以算出来了,结果就是
\mathrm{E}(x_1|A) = \frac{k}{n+1}
// Code by KSkun, 2018/3
#include <cstdio>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
const int MAXN = 13, MAXM = 50;
int n, m, ut, vt;
LL gra[MAXN], cnt[1 << MAXN], ecnt[1 << MAXN], f[1 << MAXN][MAXM], g[1 << MAXN][MAXM];
LL C[MAXM][MAXM];
inline void calc() {
for(int i = 0; i <= m; i++) {
C[i][0] = 1;
for(int j = 1; j <= i; j++) {
C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
}
}
int main() {
n = readint(); m = readint();
calc();
for(int i = 0; i < m; i++) {
ut = readint(); vt = readint();
gra[ut] |= (1 << (vt - 1));
gra[vt] |= (1 << (ut - 1));
}
for(int i = 0; i < 1 << n; i++) {
cnt[i] = cnt[i >> 1] + (i & 1);
}
for(int i = 0; i < 1 << n; i++) {
for(int j = 1; j <= n; j++) {
if(i & (1 << (j - 1))) {
ecnt[i] += cnt[gra[j] & i];
}
}
ecnt[i] >>= 1;
}
for(int i = 0; i < 1 << n; i++) {
if(cnt[i] == 1) {
g[i][0] = 1;
continue;
}
int t = i & (-i);
for(int j = (i - 1) & i; j; j = (j - 1) & i) {
if(j & t) {
for(int k1 = 0; k1 <= ecnt[j]; k1++) {
for(int k2 = 0; k2 <= ecnt[i ^ j]; k2++) {
f[i][k1 + k2] += g[j][k1] * C[ecnt[i ^ j]][k2];
}
}
}
}
for(int j = 0; j <= ecnt[i]; j++) {
g[i][j] = C[ecnt[i]][j] - f[i][j];
}
}
double ans = 0;
for(int i = 0; i <= m; i++) {
ans += double(f[(1 << n) - 1][i]) / C[m][i];
}
printf("%.6lf", ans / (m + 1));
return 0;
}