hdu4675——GCD of Sequence(莫比乌斯反演+组合数取模+乘法逆元+快速幂)

来源：互联网发布：手机登录淘宝达人编辑：程序博客网时间：2024/06/15 17:36

GCD of Sequence

Time Limit: 6000/3000 MS (Java/Others) Memory Limit: 65535/65535 K (Java/Others)
Total Submission(s): 1595 Accepted Submission(s): 540

Problem Description

Alice is playing a game with Bob.
Alice shows N integers a₁, a₂, …, a_N, and M, K. She says each integers 1 ≤ a_i ≤ M.
And now Alice wants to ask for each d = 1 to M, how many different sequences b₁, b₂, …, b_N. which satisfies :
1. For each i = 1…N, 1 ≤ b[i] ≤ M
2. gcd(b₁, b₂, …, b_N) = d
3. There will be exactly K position i that ai != bi (1 ≤ i ≤ n)

Alice thinks that the answer will be too large. In order not to annoy Bob, she only wants to know the answer modulo 1000000007.Bob can not solve the problem. Now he asks you for HELP!
Notes: gcd(x₁, x₂, …, x_n) is the greatest common divisor of x₁, x₂, …, x_n

Input

The input contains several test cases, terminated by EOF.
The first line of each test contains three integers N, M, K. (1 ≤ N, M ≤ 300000, 1 ≤ K ≤ N)
The second line contains N integers: a₁, a₂, ..., a_n (1 ≤ a_i ≤ M) which is original sequence.

Output

For each test contains 1 lines :
The line contains M integer, the i-th integer is the answer shows above when d is the i-th number.

Sample Input

3 3 33 3 33 5 31 2 3

Sample Output

7 1 059 3 0 1 1Hint
In the first test case :when d = 1, {b} can be :(1, 1, 1)(1, 1, 2)(1, 2, 1)(1, 2, 2)(2, 1, 1)(2, 1, 2)(2, 2, 1)when d = 2, {b} can be :(2, 2, 2)And because {b} must have exactly K number(s) different from {a}, so {b} can't be (3, 3, 3), so Answer = 0

Source

2013 Multi-University Training Contest 7

碰上一道数论大综合，顿时学了一波知识点。

这题好难读，读懂之后很难下手，因为不会分析复杂度...感觉就是反演的基础上加循环。

题意：

给定一个长度为n的序列a，且 1<=a[i]<=m,求有多少个序列b,使得GCD(b[1],b[2],...b[n])=d (1<=d<=m),且

正好有k个b[i]！=a[i].

也就是说，对于每个d，要给出一个答案。

f(d)表示d=gcd(b1,b2,...bn) 的所有满足的组合的总数。即最大公约数为d的组数
F(d)表示d|gcd(b1,b2,...bn)b 的所有满足的组合的总数。即公约数为d的组数

题目要求，求出每一个f(x)值。也就是公约数问题转最大公约数问题

f(x)=∑x|dμ(dx)∗F(d)。

只要求出所有的F(d)即可。

F(d)=(tot(d)n−k)∗(Md−1)tot(d)+k−n∗(Md)n−tot(d)。那个组合数可千万要注意！是C(tot(d),n-k) tot(d)>n-k

在做的时候我们先求F[d]我们枚举gcd,在原来的序列中找到d=gcd倍数的个数，设为tot.因为题目要求恰好与k个不同而且，gcd = d，那么新的序列肯定都得是d的倍数，原序列中不是d的倍数的个数为 n-tot.这些数十肯定需要修改的，如果

1)n-tot>k那么肯定是不可能的了。

2)n-tot<=k 那么这n-tot个数要换成d的倍数，每个数有m/d种，所有的就是 (m/d)^(n-tot),剩下的tot个数中还要选出k-(n-tot)个数变成不等于他们自己本身的d的倍数，每个数有(m/d-1)种，那么这种情况有

C(tot,k-(n-tot))*(m/d-1)^(k-(n-tot))种。G[d]=(m/d)^(n-tot)*C(tot,k-(n-tot))*(m/d-1)^(k-(n-tot))%mod

在求组合数的时候我们要先预处理一下阶乘的逆元。因为数据量比较大我们需要用递推法来预处理逆元

inv[n] = inv[mod%n]*(mod-mod/n)%mod;具体的证明：传送门

这里逆元我整理了一些方法，

F(d)=(cnt(d)n−k)∗(Md−1)cnt(d)+k−n∗(Md)n−cnt(d)。

算一下复杂度，预处理出所有x!即，x!的逆元，那么求(nx)的复杂度为o(1)，快速幂的复杂度为o(logn)，那么求F(d)的复杂度为o(nlogn)。

求一次f(x)的复杂度为o(nx)求所有f(x)的复杂度为o(∑nx=1nx)=o(nlogn)。

所以总的复杂度为o(nlogn)。

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
typedef long long ll;
const int maxn = 300005;
const int mod =1e9+7;
int n,m;

bool vis[maxn];
int prime[maxn];
int mu[maxn];
int tot;
int a[maxn];
int sum[300005];
ll fac[300005],inv[300005];//阶乘和逆元
ll F[maxn];
void Mobius()
{
memset(vis, 0, sizeof(vis));
mu[1] = 1;
tot = 0;
for(int i = 2; i < maxn; i ++)
{
if(!vis[i])
{
prime[tot ++] = i;
mu[i] = -1;
}
for(int j = 0; j < tot; j ++)
{
if(i * prime[j] >= maxn) break;
vis[i * prime[j]] = true;
if(i % prime[j] == 0)
{
mu[i * prime[j]] = 0;
break;
}
else
{
mu[i * prime[j]] = -mu[i];
}
}
}
}
ll qmod(ll a,ll b)
{
ll ans=1;
while(b)
{
if(b&1)ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans%mod;
}

/******递推求逆元一条龙********/
ll get_inv(ll n){
if(n==1) return 1;
return get_inv(mod%n)*(mod-mod/n)%mod;
}

ll Com(ll n,ll m){
return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
void init(){
inv[0]=fac[0]=1;
for(int i=1;i<maxn;i++){
fac[i]=fac[i-1]*i%mod;
inv[i]=get_inv(fac[i]);
}
}

/**************************/
int main()
{
int i,j,k,t;
init();
Mobius();
while(scanf("%d%d%d",&n,&m,&k)!=EOF)
{
memset(sum,0,sizeof(sum));
for(i=0; i<n; i++){scanf("%d",&a[i]);sum[a[i]]++;}
for(i=1;i<=m;i++){
for(j=i+i;j<=m;j+=i)//十分巧妙地统计倍数，用这一块代码，这样再从1遍历到m的时候，速度增加非常快，然后就不会超时。
sum[i]+=sum[j];
}

for(int d=1;d<=m;d++){
if(n-sum[d]>k){F[d]=0;continue;}
ll tmp=(inv[n-k]*inv[k-(n-sum[d])])%mod;//求逆元
tmp=(fac[sum[d]]*tmp)%mod; //求组合数
F[d]=tmp*(qmod(m/d-1,sum[d]-(n-k))*qmod(m/d,n-sum[d])%mod)%mod;//莫尼乌斯函数
}
ll ans;
for(int x=1;x<=m;x++){
ans=0;
for(int d=x;d<=m;d+=x){
ans=(ans+mu[d/x]*F[d]+mod)%mod; //莫比乌斯反演
}
if(x==m)printf("%I64d\n",ans);
else printf("%I64d ",ans);
}
}
return 0;
}

扩展欧几里得求逆元：

void Gcd(LL a,LL b,LL& d,LL& x,LL& y){
if(!b) {d = a;x = 1;y = 0;}
else{Gcd(b,a%b,d,y,x);y -= x*(a/b);}
}
LL Inv(LL a,LL n){
LL d,x,y;
Gcd(a,n,d,x,y);
return d == 1 ? (x+n)%n : -1;
}
void init(){
fac[0] = 1;
inv[0] = 1;
FOR(i,1,maxn){
fac[i] = fac[i-1]*i;
fac[i] %= Mod;
inv[i] = Inv(fac[i],Mod);
}
}

其他思路：F(d)=(cnt(d)n−k)∗(Md−1)cnt(d)+k−n∗(Md)n−cnt(d)。
http://www.acmerblog.com/hdu-4675-gcd-of-sequence-7720.html

F(d)=(cnt(d)n−k)∗(Md−1)cnt(d)+k−n∗(Md)n−cnt(d)。
http://blog.csdn.net/bigbigship/article/details/47679713

F(d)=(cnt(d)n−k)∗(Md−1)cnt(d)+k−n∗(Md)n−cnt(d)。

1 0