首页 > 代码库 > HDU 5321 Beautiful Set 漂亮集合
HDU 5321 Beautiful Set 漂亮集合
题意:给定一个集合,含有n个数。浙理工先生和杭电先生各自有计算这个集合漂亮值的方法。
浙理工先生的计算方法是:对于这个n个数的某个排列,此排列的漂亮值为这个排列全部的区间最大公约数之和。然后这个集合的漂亮值为n个数的全部排列的漂亮值之和。
杭电先生的计算方法是:在这个n个数中选出k(1 ≤ k ≤ n)个数。对于某种选取方案。这样的方案的漂亮值为k个数的最大公约数乘上k。然后这个集合的漂亮值为全部选数方案的漂亮值之和。
然后他们想比比谁得到的漂亮值更大。由于数非常大,所以他们仅仅比較各自的结果对258280327取余的大小。
做此题的过程真是艰难。。。
。。我来分享下。
。
。
。。
看起来全然没法做。
只是数论题么,总应该是能够怎么乱搞一下或者推公式的。
。。。。
推出公式来要非常半天。然后写程序10分钟基本上就这样。
。
。。
。
当然。依照题意去写公式非常好写,然而公式是(n^2)的复杂度。须要利用各种数论定理或者化简技巧去化简。
我们先看浙理工先生的计算方法。
对随意一个排列。计算全部区间GCD的和。所以直接枚举排列再枚举区间不就完了嘛。然而这显然是不可能的。于是我们转化一下方法。既然归根结底是求全部的区间GCD,我们最好还是看看有哪些的区间能够用来计算。
显然,在n个数中任取k个数,这k个数连续的放在一起。能够构成一个长度为k区间。由于k个数随意排列都不影响整个区间的GCD,所以我们把这整个区间看成一个物品,和剩下n - k个数一共组成n - k + 1个物品,又有(n - k + 1)!种排列方法(也就是捆绑法思想)。所以这个区间在全部的排列中。一共会被计算k!*(n - k + 1)!次。而每次计算都会加上这个区间的GCD值,最好还是设为x,于是GCD为x,且长度为k的这样一个区间对整个漂亮值的贡献即为x*k!*(n - k + 1)!。我们最好还是设n个数中最大的数为MAXN。我们定义一个函数G(k, x)表示从n个数中选k个数且这个k个数最大公约数为x的方案数。
于是我们得到了计算浙理工先生的计算公式:
ZSTU = sigma(1 ≤ k ≤ n。sigma(1 ≤ x ≤ MAXN,G(k,x)*x*k!*(n - k + 1)! ))。
然而。这样暴力枚举k和x是n^2复杂度。必定须要将公式化简。直观看来,两个sigma求和似乎都找不到什么好的公式。
于是我们对G(k,x)下手。设法计算G(k,x)的值,看能不能有什么突破。
考虑到G(k,x)是表示k个数最大公约数恰为x的方案数,恰为x不好计算,我们能够定义一个辅助函数——F(k,x),表示n个数中选k个数的最大公约数为x的倍数的方案数。为什么用这个辅助函数呢。由于这个函数非常好计算——既然gcd是x的倍数。那么就在全部x的倍数中选k个不就完了嘛。
于是我们定一个cnt[x]表示n个数中x的倍数的个数。就得到F(k,x) = C(cnt[x],k)(C是组合数,下同)。我们再来看看F函数和G函数的关系——非常easy,事实上就是F(k,x) = sigma(x | d,G(k,d))看到这样的形式的等式。不由得就想到了莫比乌斯反演。于是G函数的计算就被转化为F了。
然而。貌似好像前进了一大步。只是,剩下来的问题更加困难——怎么化简公式。打字太麻烦。
。。。。我就直接上图。。。。
。计算中的μ函数是莫比乌斯函数。
注意这里第二次变形的时候把x和d的地位互换了,为什么呢。
由于枚举约数x计算倍数d有个限制,那就是d不能超过MAXN。这样就比較麻烦,而转化为枚举倍数计算约数时。每一个数的约数是固定的。于是计算起来就比較方便。所以对上式最后一个sigma,我们能够对任一的d直接预处理。我们记最后一个sigma所表示的值为H(d)。这个H(d)是能够离线打表计算的。继续化简——
化简到这里,貌似还是个(n^2)公式。只是内层的sigma已经是一个很有规律的公式了——我们把它写出来看看:答案记为P(cnt[d])。然后我们来求这个P函数。为了简单点写。
。。。。就临时用x取代cnt[d]。
P(x)
= sigma(1 ≤ k ≤ x。(n - k + 1)! / (x - k)!
= n! / (x - 1)! + (n - 1)! / (x - 2)! + ...... + (n - x + 1)! / 0!
= n*(n - 1)*(n - 2)* ...... *x +
(n - 1)*(n - 2)* ...... *(x - 1) +
(n - 2)*(n - 3)* ...... *(x - 2) +
(n - x + 1)*(n - x)* ...... *1
能够发现此和式的每一项都是连续(n - x + 1)个数相乘的结果,且每一项的连续数的开头都等于上一项的开头 + 1。
对于这样的类型的和式,有一种很巧妙的求和办法——我举个简单样例:求 1*2*3 + 2*3*4 + ..... + n*(n + 1)*(n + 2)。
我们记答案为ans。
对于每一个i*(i + 1)*(i + 2)这一项,我们都让它乘上4(4的由来是每一项中连乘数的个数 + 1)然后。对于每一项中的4。都拆成( (i + 3) - (i - 1) )。为什么要这么拆呢?拆完带进去你就发现答案了——
4*ans
= 1*2*3*(4 - 0) + 2*3*4*(5 - 1) + ...... +(n*(n + 1)*(n + 2)*( (n + 3) - (n - 1) )
= 1*2*3*4 - 0*1*2*3 + 2*3*4*5 - 1*2*3*4 + 3*4*5*6 - 2*3*4*5 + ...... + n*(n + 1)*(n + 2)*(n + 3) - (n - 1)*n*(n + 1)*(n + 2)
中间项全都抵消了。于是答案为ans = n*(n + 1)*(n + 2)*(n + 3)/4。
这样,我们得到了P(x)的公式——P(x) = (n + 1)*n*(n - 1)* ...... *x / (n - x + 1) = (n + 1)! / (x - 1)! / (n - x + 1) —— 在预处理出逆元、阶乘取模、阶乘逆元后。这是个O(1)的公式,除法能够用逆元计算。至于逆元的求法不要用扩欧,由于会反复用到非常多个数的逆元。这里我们用O(n)递推打表的公式求每一个数的逆元inv[i]——公式为inv[i] = (MOD - MOD/i)*inv[MOD%i]%MOD,约定inv[0] = inv[1] = 1。这个公式我就不证了。。
。
。。这样,我们先预处理出莫比乌斯函数、H函数、逆元、阶乘取模和阶乘逆元等,最后的公式就全然成了个O(n)的公式。至此浙理工先生的答案已成功解出。
那么杭电先生的公式呢。比較简单。我就直接写了——当中内层G函数还是转化为F函数去求,方法和计算浙理工先生时一模一样。
HDU
= sigma(1 ≤ k ≤ n,sigma(1 ≤ x ≤ MAXN。G(k。x)*k))
= sigma(1 ≤ d ≤ n,H(d)*sigma(1 ≤ k ≤ cnt[d],k*C(cnt[d],k)))
由于k*C(cnt[d],k) = cnt[d]*C(cnt[d] - 1, k - 1)。所以内层sigma
= sigma(1 ≤ k ≤ cnt[d]。k*C(cnt[d]。k)))
= sigma(1 ≤ k ≤ cnt[d]。cnt[d]*C(cnt[d] - 1。k - 1))
= cnt[d]*(sigma(0 ≤ k ≤ cnt[d] - 1,C(cnt[d] - 1。k))
= cnt[d]*2^(cnt[d] - 1)。
带回去后得到HDU
= sigma(1 ≤ d ≤ n。H(d)*cnt[d]*2^(cnt[d] - 1)。
OK,两位先生的答案都算出来了,再比較一下输出就能够了!
#include <iostream> #include <cstdio> #include <cmath> #include <algorithm> #include <cstring> #include <string> #include <vector> #include <map> #include <stack> #include <queue> using namespace std; const int MAX = 1e5 + 5; const int MOD = 258280327; bool u[MAX]; //素数筛 int prime[MAX], num; //素数表 int mobius[MAX]; //莫比乌斯函数表 long long inv[MAX]; //逆元表 long long factorial[MAX]; //阶乘表 long long fac_inv[MAX]; //阶乘逆元表 long long pow2[MAX]; //2的次方表 long long sum[MAX]; //H函数表 int n, max_value; int cnt[MAX], cnt_temp[MAX]; //cnt_temp存储每一个数出现了多少次,方便cnt数组的计算 long long P[MAX]; void initial() //预处理各种表 { memset(u, true, sizeof(u)); num = 0; mobius[1] = 1; for(int i = 2; i < MAX; i++) { if(u[i]) { prime[num++] = i; mobius[i] = -1; } for(int j = 0; j < num; j++) { if(i*prime[j] >= MAX) break; u[i*prime[j]] = false; if(i%prime[j] == 0) { mobius[i*prime[j]] = 0; break; } else mobius[i*prime[j]] = -mobius[i]; } } inv[0] = inv[1] = 1; factorial[0] = factorial[1] = 1; fac_inv[0] = fac_inv[1] = 1; for(int i = 2; i < MAX; i++) { inv[i] = (MOD - MOD/i)*inv[MOD%i]%MOD; factorial[i] = factorial[i - 1]*i%MOD; fac_inv[i] = fac_inv[i - 1]*inv[i]%MOD; } memset(sum, 0, sizeof(sum)); sum[1] = 1; for(int i = 2; i < MAX; i++) { for(int j = 1; j*j <= i; j++) { if(i%j == 0) { sum[i] += j*mobius[i/j]; if(j*j != i) sum[i] += i/j*mobius[j]; } } } pow2[0] = 1; for(int i = 1; i < MAX; i++) pow2[i] = pow2[i - 1]*2%MOD; } void input() { int x; memset(cnt_temp, 0, sizeof(cnt_temp)); memset(cnt, 0, sizeof(cnt)); max_value = http://www.mamicode.com/0;"%d", &x); cnt_temp[x]++; max_value = http://www.mamicode.com/max(max_value, x);>for(int i = 1; i <= n; i++) P[i] = factorial[n + 1]*fac_inv[i - 1]%MOD*inv[n - i + 2]%MOD; long long ZSTU = 0, HDU = 0; for(int i = 1; i <= max_value; i++) { HDU = (HDU + sum[i]*cnt[i]%MOD*(pow2[cnt[i] - 1]))%MOD; ZSTU = (ZSTU + sum[i]*factorial[cnt[i]]%MOD*P[cnt[i]]%MOD)%MOD; } if(ZSTU > HDU) printf("Mr. Zstu %I64d\n", ZSTU); if(HDU > ZSTU) printf("Mr. Hdu %I64d\n", HDU); if(ZSTU == HDU) printf("Equal %I64d\n", HDU); } int main() { //for(int i = 2; i*i <= MOD; i++)if(MOD%i == 0){printf("NO\n");return 0;}printf("YES\n"); 好吧它是个质数 ^_^ initial(); while(scanf("%d", &n) != EOF) { input(); solve(); } return 0; }
??
HDU 5321 Beautiful Set 漂亮集合