TopCoder 11351 TheCowDivOne

vj

考虑生成函数,答案可以写成\([x^0y^k]\prod_{i=0}^{n-1}(1+x^iy)\pmod {x^n-1}\)

\(=\sum_i[x^{in}y^k]\prod_{j=0}^{n-1}(1+x^jy)\)

\(=\sum_i[x^iy^k]\prod_{j=0}^{n-1}(1+x^jy)[n|i]\)

这里单位根反演,原式

\(=\frac{1}{n}[y^k]\sum_{i=0}^{n-1}\prod_{j=0}^{n-1}(1+\omega_n^{ij}y)\)

然后对于每个\(i\),后面那坨多项式的形式只和\(\gcd(i,n)\)有关,这里改为枚举这个\(\gcd\)的值,原式

\(=\frac{1}{n}\sum_{d|n}[y^k](\prod_{j=0}^{n/d-1}(1+\omega_n^{jd}y))^d\sum_{j=1}^{n}[\gcd(j,n)=d]\)

\(=\frac{1}{n}\sum_{d|n}[y^k](\prod_{j=0}^{d-1}(1+\omega_d^jy))^{n/d}\sum_{j=1}^{n}[\gcd(j,n)=\frac{n}{d}]\)

\(=\frac{1}{n}\sum_{d|n}\varphi(d)[y^k](1+(-1)^{d+1}y^d)^{n/d}\)

然后二项式定理,求出对应组合数即可知道对应的\(y^k\)的系数.这里的\(d\)也要满足\(d|k\)

btw,下面这个代码不是忘topcoder上交的(

#include<bits/stdc++.h>
#define LL long long
#define db double

using namespace std;
const int N=2000+10,mod=1e9+7;
int rd()
{
    int x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
    return x*w;
}
void ad(int &x,int y){x+=y,x-=x>=mod?mod:0;}
int fpow(int a,int b){int an=1;while(b){if(b&1) an=1ll*an*a%mod;a=1ll*a*a%mod,b>>=1;}return an;}
int ginv(int a){return fpow(a,mod-2);}
int n,k,inv[N],d[N],td,phi[N],p[N][2],tt;
void dfs(int o,int s,int ph)
{
    if(o>tt){d[++td]=s,phi[td]=ph;return;}
    int ft=1;
    for(int i=0;i<=p[o][1];++i)
	dfs(o+1,s,ph),s*=p[o][0],ph*=p[o][0]-ft,ft=0;
}
int C(int a,int b)
{
    if(b<0||a<b) return 0;
    int an=1;
    for(int i=1;i<=b;++i) an=1ll*an*(a-i+1)%mod*inv[i]%mod;
    return an;
}

int main()
{
    inv[0]=inv[1]=1;
    for(int i=2;i<=N-5;++i) inv[i]=(mod-1ll*mod/i*inv[mod%i]%mod)%mod;
    while(~scanf("%d%d",&n,&k))
    {
    	tt=0;
    	int x=n,sqt=sqrt(n);
    	for(int i=2;i<=sqt;++i)
    	    if(x%i==0)
    	    {
    	    	p[++tt][0]=i,p[tt][1]=0;
    	    	while(x%i==0) x/=i,++p[tt][1];
    	    }
    	if(x>1) p[++tt][0]=x,p[tt][1]=1;
    	td=0,dfs(1,1,1);
    	int ans=0;
    	for(int i=1;i<=td;++i)
    	{
    	    int x=d[i];
    	    if(k%x) continue;
    	    int dx=1ll*phi[i]*C(n/x,k/x)%mod;
    	    ad(ans,((x^1)&(k/x)&1)?mod-dx:dx);
    	}
    	ans=1ll*ans*ginv(n)%mod;
    	printf("%d\n",ans);
    }
    return 0;
}