【题解】TopCoder 11351 TheCowDivOne (单位根反演)
题目大意:
你有(n)头牛，编号(0sim n-1)。问你从中选出(k)头牛，且选出牛的编号的和为(n)的倍数，问有多少方案，答案对1e9+7取模。(kle n le 10^9)
这个可以用一个二元生成函数表示
[ans=sum_i^{n(n-1)over 2} [n|i][x^i][y^k]prod_{j=0}^{n-1} (1+x^jy)
]
令(f(x)=[y^k]prod_{j=0}^{n-1} (1+x^jy))，直接考虑单位根反演，原式等于
[=sum_{i=0}{1over n}f(omega_n^i)
]
然而现在还是不好搞，我们把(f(x))展开
[=[y^k]{1over n} sum_{i=0} prod_j^{n-1} (1+omega_n^{ij}y)
]
根据单位根的性质，可以发现((1+omega_n^{ij}))是有循环节的，这个的循环节的情况等价于$ ext{for n>j>=0 , } i imes j mod n $的情况。
一个好结论:[1]
设(a_j=ijmod n,b_j=djmod n)，其中(d=gcd(i,n),i=kd)那么有(kperp n, herefore exist k^{-1} (mod n))。此外(a_j=a_{j+{nover d}}) ，
因此(a_j=b_{jk^{-1} mod n},b_j=a_{jk})，构成一一对应。由此，有
[prod (1+omega_n^{ij})=prod (1+omega_n^{dj})=prod (1+omega_{nover d}^j) ]
因此还按循环节长度(d)枚举，并且乘上贡献系数(有多少个(i)满足循环节长度是(d)，也就是((i,n)={nover d},i<d)的个数，这个数量=(phi(d)))。
[=[y^k]{1over n} sum_{d|n} phi(d)prod_{j=0}^{d-1} (1+omega_d^{j}y)^{nover d}
]
二项式定理化开幂
[=[y^k]{1over n} sum_{d|n} phi(d)prod_{j=0}^{d-1} sum_i^{nover d} {{nover d }choose i} omega_d^{ij} y^i
]
(j)被孤立，交换prod
[=[y^k]{1over n} sum_{d|n} phi(d) sum_i^{nover d} {{nover d }choose i} y^{di} prod_{j=0}^{d-1}omega_d^{ij}
]
化简得
[=[y^k]{1over n} sum_{d|n} phi(d) sum_i^{nover d} {{nover d }choose i} y^{di} omega_d^{{d(d-1)over 2} imes i}
]
又因为
[omega_d^{d(d-1)over 2}=e^{2pi i imes {d(d-1)over 2}over d}=e^{pi i(d-1)}=cos(pi(d-1))+isin(pi(d-1))=(-1)^{d+1}
]
那么
[={1over n} sum_{d|n} phi(d) {{nover d}choose {kover d}} (-1)^{(d+1){kover d}}[d|k]
]
用那个根号算(phi(x))，复杂度(O(ksqrt n+n^{0.75}))
本题主要利用了那个结论，也就是观察到这个式子只和(gcd(i,n))有关，把状态数直接将至(sqrt n)数量级
//@winlere
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std; typedef long long ll;
inline int qr(){
int ret=0,f=0,c=getchar();
while(!isdigit(c))f|=c==45,c=getchar();
while(isdigit(c)) ret=ret*10+c-48,c=getchar();
return f?-ret:ret;
}
const int mod=1e9+7;
const int maxn=1005;
int inv[maxn];
int phi(int x){
if(x==1) return 1;
int ret=x;
for(int t=2;1ll*t*t<=x;++t){
if(x%t==0){
while(x%t==0) x/=t;
ret-=ret/t;
}
}
if(x>1) ret-=ret/x;
return ret;
}
int c(int n,int m){
int ret=1;
for(int t=1;t<=m;++t)
ret=1ll*ret*(n-t+1)%mod;
return 1ll*ret*inv[m]%mod;
}
void pre(const int&n){
inv[1]=1; inv[0]=1;
for(int t=2;t<=n;++t) inv[t]=1ll*(mod-mod/t)*inv[mod%t]%mod;
for(int t=2;t<=n;++t) inv[t]=1ll*inv[t]*inv[t-1]%mod;
}
int ksm(const int&ba,const int&p){
int ret=1;
for(int t=p,b=ba;t;t>>=1,b=1ll*b*b%mod)
if(t&1) ret=1ll*ret*b%mod;
return ret;
}
int getAns(int n,int d,int k){
int ret=1ll*phi(d)*c(n/d,k)%mod;
return (1ll*(d+1)*k)%2?mod-ret:ret;
}
int main(){
pre(1e3);
int n,k;
while(~scanf("%d%d",&n,&k)){
int ans=0;
for(int t=1;1ll*t*t<=n;++t){
if(n%t==0){
if(k%t==0) ans=(ans+getAns(n,t,k/t))%mod;
if(t*t!=n&&k%(n/t)==0) ans=(ans+getAns(n,n/t,k/(n/t)))%mod;
}
}
ans=1ll*ans*ksm(n,mod-2)%mod;
printf("%d
",ans);
}
return 0;
}
试下这个东西 ↩︎