【题解】TopCoder 11351 TheCowDivOne (单位根反演)

【题解】TopCoder 11351 TheCowDivOne (单位根反演)

传送门

题目大意:

你有\(n\)头牛，编号\(0\sim n-1\)。问你从中选出\(k\)头牛，且选出牛的编号的和为\(n\)的倍数，问有多少方案，答案对1e9+7取模。\(k\le 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}{1\over n}f(\omega_n^i) \]

然而现在还是不好搞，我们把\(f(x)\)展开

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

根据单位根的性质，可以发现\((1+\omega_n^{ij})\)是有循环节的，这个的循环节的情况等价于$\text{for n>j>=0 , } i\times j \mod n $的情况。

一个好结论:^[1]

设\(a_j=ij\mod n,b_j=dj\mod n\)，其中\(d=\gcd(i,n),i=kd\)那么有\(k\perp n,\therefore \exist k^{-1} (\mod n)\)。此外\(a_j=a_{j+{n\over 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_{n\over d}^j) \]

因此还按循环节长度\(d\)枚举，并且乘上贡献系数(有多少个\(i\)满足循环节长度是\(d\)，也就是\((i,n)={n\over d},i<d\)的个数，这个数量=\(\phi(d)\))。

\[=[y^k]{1\over n} \sum_{d|n} \phi(d)\prod_{j=0}^{d-1} (1+\omega_d^{j}y)^{n\over d} \]

二项式定理化开幂

\[=[y^k]{1\over n} \sum_{d|n} \phi(d)\prod_{j=0}^{d-1} \sum_i^{n\over d} {{n\over d }\choose i} \omega_d^{ij} y^i \]

\(j\)被孤立，交换prod

\[=[y^k]{1\over n} \sum_{d|n} \phi(d) \sum_i^{n\over d} {{n\over d }\choose i} y^{di} \prod_{j=0}^{d-1}\omega_d^{ij} \]

化简得

\[=[y^k]{1\over n} \sum_{d|n} \phi(d) \sum_i^{n\over d} {{n\over d }\choose i} y^{di} \omega_d^{{d(d-1)\over 2}\times i} \]

又因为

\[\omega_d^{d(d-1)\over 2}=e^{2\pi i \times {d(d-1)\over 2}\over d}=e^{\pi i(d-1)}=\cos(\pi(d-1))+i\sin(\pi(d-1))=(-1)^{d+1} \]

那么

\[={1\over n} \sum_{d|n} \phi(d) {{n\over d}\choose {k\over d}} (-1)^{(d+1){k\over d}}[d|k] \]

用那个根号算\(\phi(x)\)，复杂度\(O(k\sqrt 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\n",ans);
	}
	return 0;
}

---

1. 试下这个东西 ↩︎