您的位置:首页 > 其它

CodeForces 329 E.Devu and Birthday Celebration(莫比乌斯反演+组合数学)

2018-01-08 21:57 417 查看
Description

问满足gcd(a1,a2,...,af)=1,∑i=1fai=n的序列a的个数

Input

第一行一整数T表示用例组数,每组用例输入两个整数n,f表示一组查询(1≤T≤105,1≤f≤n≤105)

Output

对于每组查询,输出满足条件的a序列个数

Sample Input

5

6 2

7 2

6 3

6 4

7 4

Sample Output

2

6

9

10

20

Solution

令f(d)为满足gcd(a1,a2,...,af)=d,∑i=1fai=n的方案数,F(d)为满足d|gcd(a1,a2,...,af),∑i=1fai=n的方案数,则F(d)=∑d|kf(k),由莫比乌斯反演,f(d)=∑d|kμ(kd)F(k)

现在求F(k),显然k|n,问题转化为将nk拆成f个正整数bi,由插板法得F(k)=Cf−1nk−1

进而得知答案为f(1)=∑d|nμ(d)Cf−1nd−1,时间复杂度O(Tn√)

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=100005;
#define mod 1000000007
int prime[maxn],tot,mu[maxn],check[maxn];
void Moblus(int n=1e5)
{
mu[1]=1;
tot=0;
for(int i=2;i<=n;i++)
{
if(!check[i])
{
prime[tot++]=i;
mu[i]=-1;
}
for(int j=0;j<tot&&i*prime[j]<=n;j++)
{
check[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
}
int fact[maxn],inv[maxn];
void init(int n=1e5)
{
fact[0]=1;
for(int i=1;i<=n;i++)fact[i]=(ll)i*fact[i-1]%mod;
inv[1]=1;
for(int i=2;i<=n;i++)inv[i]=mod-(ll)(mod/i)*inv[mod%i]%mod;
inv[0]=1;
for(int i=1;i<=n;i++)inv[i]=(ll)inv[i-1]*inv[i]%mod;
}
int C(int n,int m)
{
return (ll)fact
*inv[m]%mod*inv[n-m]%mod;
}
void add(int &x,int y)
{
x=x+y>=mod?x+y-mod:x+y;
}
int T,n,f;
int Solve(int i)
{
if(mu[i]==1&&n/i>=f)return C(n/i-1,f-1);
if(mu[i]==-1&&n/i>=f)return mod-C(n/i-1,f-1);
return 0;
}
int main()
{
Moblus();
init();
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&f);
int ans=0;
for(int i=1;i*i<=n;i++)
if(n%i==0)
{
add(ans,Solve(i));
if(i*i!=n)add(ans,Solve(n/i));
}
printf("%d\n",ans);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: