博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
HDU 5321 Beautiful Set 漂亮集合
阅读量:5944 次
发布时间:2019-06-19

本文共 5707 字,大约阅读时间需要 19 分钟。

题意:给定一个集合,含有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 
#include
#include
#include
#include
#include
#include
#include
#include
#include
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 = 0; for(int i = 0; i < n; i++) { scanf("%d", &x); cnt_temp[x]++; max_value = max(max_value, x); }}void solve(){ for(int i = 1; i <= max_value; i++) //计算cnt数组 { for(int j = i; j <= max_value; j += i) //枚举i的倍数 cnt[i] += cnt_temp[j]; } P[0] = 0; //不要忘了特殊值。

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; }




你可能感兴趣的文章
double free or corruption (fasttop)
查看>>
Mac 常用快捷键
查看>>
阿里云CentOS7安装Oracle11GR2
查看>>
nginc+memcache
查看>>
从拼多多优惠券事件看到的一些反思
查看>>
mac下完全卸载postgresql的方法
查看>>
20个纯css3写的logo
查看>>
四周第四次课 6.1 压缩打包介绍 6.2 gzip压缩工具 6.3 bzip2压缩工具 6.4 x
查看>>
交换机自动学习vlan
查看>>
三层交换配置与原理
查看>>
Nginx ssl、rewrite配置
查看>>
thinkphp-查询数据-基本查询
查看>>
bootstrap-自适应导航
查看>>
SQL-4查找所有已经分配部门的员工的last_name和first_name(自然连接)
查看>>
查找最近修改的SP
查看>>
linux交换空间
查看>>
加入马帮,马到功成
查看>>
使用cpau.exe让不是管理员的用户也有权限运行哪些需要管理员权限的软件。
查看>>
编译安装mariadb-10.0.10
查看>>
UML类图
查看>>