BZOJ 3434 [WC2014]时空穿梭 (莫比乌斯反演)
好难啊..反演的终极题目
首先，本题的突破口在于直线的性质。不论是几维的空间，两点一定能确定一条直线
选取两个点作为最左下和最右上的点！
假设现在是二维空间，选取了$(x1,y1)$和$(x2,y2)$两个点，那么它们连线上经过的点数就是$gcd(x2-x1,y2-y1)-1$
选取的方案数为$\left ( _{gcd-1}^{c-2} \right )$
发现这个方案数只和坐标的差值有关，我们直接枚举差值就行
接下来就是反演了，为了方便叙述，以$n=2$为例
$\sum\limits_{x=1}^{m_{1}}\sum\limits_{y=1}^{m_{2}} \left ( _{c-2}^{gcd(x,y)-1} \right ) (m_{1}-x)(m_{2}-y)$
$\sum\limits_{k=1}^{m} \left ( _{c-2}^{k-1} \right ) \sum\limits_{x=1}^{m_{1}} \sum\limits_{y=1}^{m_{2}} [gcd(x,y)==k] (m_{1}-x)(m_{2}-y)$
$\sum\limits_{k=1}^{m} \left ( _{c-2}^{k-1} \right ) \sum\limits_{x=1}^{\left \lfloor \frac{m_{1}}{k} \right \rfloor} \sum\limits_{y=1}^{\left \lfloor \frac{m_{2}}{k} \right \rfloor} [gcd(x,y)==1] (m_{1}-xk)(m_{2}-yk)$
$\sum\limits_{k=1}^{m} \left ( _{c-2}^{k-1} \right ) \sum\limits_{x=1}^{\left \lfloor \frac{m_{1}}{k} \right \rfloor} \sum\limits_{y=1}^{\left \lfloor \frac{m_{2}}{k} \right \rfloor} \sum\limits_{d|k}\mu(d)(m_{1}-xk)(m_{2}-yk)$
$\sum\limits_{k=1}^{m} \left ( _{c-2}^{k-1} \right ) \sum\limits_{d=1}^{\left \lfloor \frac{m_{1}}{k} \right \rfloor} \mu(d) \sum\limits_{x=1}^{\left \lfloor \frac{m_{1}}{kd} \right \rfloor} (m_{1}-xkd) \sum\limits_{y=1}^{ \left \lfloor \frac{m_{2}}{kd} \right \rfloor} (m_{2}-ykd)$
令$Q=kd$
$\sum\limits_{Q=1}^{m} \left ( \sum\limits_{k|Q} \left ( _{c-2}^{k-1} \right ) \mu(\frac{Q}{k}) \right ) \left ( \sum\limits_{x=1}^{\left \lfloor \frac{m_{1}}{Q} \right \rfloor} (m_{1}-xQ) \right ) \left ( \sum\limits_{y=1}^{ \left \lfloor \frac{m_{2}}{Q} \right \rfloor} (m_{2}-yQ) \right ) $
利用等差数列公式 $\left ( \sum\limits_{x=1}^{\left \lfloor \frac{m}{Q} \right \rfloor} (m-xQ) \right ) = \left ( \left \lfloor \frac{m}{Q} \right \rfloor m- \frac{ (1+\left \lfloor \frac{m}{Q} \right \rfloor )\left \lfloor \frac{m}{Q} \right \rfloor Q}{2} \right ) $
$\sum\limits_{Q=1}^{m} \left ( \sum\limits_{k|Q} \left ( _{c-2}^{k-1} \right ) \mu(\frac{Q}{k}) \right ) \left ( \left \lfloor \frac{m_{1}}{Q} \right \rfloor m_{1}- \frac{ (1+\left \lfloor \frac{m_{1}}{Q} \right \rfloor )\left \lfloor \frac{m_{1}}{Q} \right \rfloor Q}{2} \right ) \left ( \left \lfloor \frac{m_{2}}{Q} \right \rfloor m_{2}- \frac{ (1+\left \lfloor \frac{m_{2}}{Q} \right \rfloor )\left \lfloor \frac{m_{2}}{Q} \right \rfloor Q}{2} \right ) $
推广到更高维上
$\sum\limits_{Q=1}^{m} \left ( \sum\limits_{k|Q} \left ( _{c-2}^{k-1} \right ) \mu(\frac{Q}{k}) \right ) \prod_{t=1}^{n} \left ( \left \lfloor \frac{m_{t}}{Q} \right \rfloor m_{t}- \frac{ (1+\left \lfloor \frac{m_{t}}{Q} \right \rfloor )\left \lfloor \frac{m_{t}}{Q} \right \rfloor Q}{2} \right ) $
我们不能把$Q$这一项带入到整除分块里去
所以每次$O(n^{2})$暴力计算出$Q^{k}(k=0,1,2...n)$的系数
预处理出$f(Q,c)=\sum\limits_{k|Q} \left ( _{c-2}^{k-1} \right ) Q^{u}$
用整除分块即可
时间$O(Tn^{3}\sqrt {m})$
人傻常数大，BZOJ过不去，洛谷上开O2过了，懒得优化了
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #define ll long long 5 #define N1 100100 6 using namespace std; 7 const int inf=0x3f3f3f3f; 8 9 int gint() 10 { 11 int ret=0,fh=1;char c=getchar(); 12 while(c<'0'||c>'9'){if(c=='-')fh=-1;c=getchar();} 13 while(c>='0'&&c<='9'){ret=ret*10+c-'0';c=getchar();} 14 return ret*fh; 15 } 16 const int zwz=10007; 17 int qpow(int x,int y) 18 { 19 int ans=1; 20 for(;y;x=x*x%zwz,y>>=1) if(y&1) ans=ans*x%zwz; 21 return ans; 22 } 23 24 int T,n,c,cnt; 25 int mu[N1],pr[N1],use[N1], f[21][N1][21],m[N1],C[N1][21],mul[N1],_mul[N1]; 26 int maxn=100000; 27 int Lucas(int a,int b) 28 { 29 if(b>a) return 0; 30 if(a<zwz&&b<zwz) return mul[a]*_mul[b]%zwz*_mul[a-b]%zwz; 31 else return Lucas(a%zwz,b%zwz)*Lucas(a/zwz,b/zwz)%zwz; 32 } 33 void Pre() 34 { 35 int i,j,k,l; 36 mu[1]=1; 37 for(i=2;i<=maxn;i++) 38 { 39 if(!use[i]) pr[++cnt]=i,mu[i]=-1; 40 for(j=1;(j<=cnt)&&(i*pr[j]<=maxn);j++) 41 { 42 use[i*pr[j]]=1; 43 if(i%pr[j]) mu[i*pr[j]]=-mu[i]; 44 else mu[i*pr[j]]=0; 45 } 46 } 47 mul[0]=mul[1]=_mul[0]=_mul[1]=1; 48 for(i=2;i<zwz;i++) 49 { 50 mul[i]=mul[i-1]*i%zwz; 51 _mul[i]=qpow(mul[i],zwz-2); 52 } 53 for(i=0;i<=maxn;i++) for(k=0;k<=min(i,18);k++) 54 C[i][k]=Lucas(i,k); 55 for(i=1;i<=maxn;i++) for(j=i;j<=maxn;j+=i) 56 { 57 for(k=2;k<=20;k++) 58 (f[0][j][k]+=C[i-1][k-2]*mu[j/i]%zwz+zwz)%=zwz; 59 } 60 for(l=1;l<=20;l++) for(i=1;i<=maxn;i++) for(k=2;k<=20;k++) f[l][i][k]=i*f[l-1][i][k]%zwz; 61 for(l=0;l<=20;l++) for(i=1;i<=maxn;i++) for(k=2;k<=20;k++) (f[l][i][k]+=f[l][i-1][k])%=zwz; 62 } 63 int p[2][22]; 64 65 int main() 66 { 67 scanf("%d",&T); 68 int i,la,j,x,y,k,t,mi,ans,tmp,inv2,now,pst;// ll ans=0; 69 Pre(); 70 for(t=1;t<=T;t++){ 71 72 scanf("%d%d",&n,&c); ans=0; inv2=qpow(2,zwz-2); 73 for(i=1,mi=inf;i<=n;i++) m[i]=gint(), mi=min(mi,m[i]); 74 for(i=1;i<=mi;i=la+1) 75 { 76 for(j=1,la=inf;j<=n;j++) la=min(la,m[j]/(m[j]/i)); 77 memset(p,0,sizeof(p)); now=1; pst=0; p[pst][0]=1; 78 for(j=1;j<=n;j++) 79 { 80 memset(p[now],0,sizeof(p[now])); 81 for(k=0;k<n;k++) 82 p[now][k+1]=(p[now][k+1]-p[pst][k]*(1+m[j]/i)%zwz*(m[j]/i)%zwz*inv2%zwz+zwz)%zwz, 83 p[now][k]=(p[now][k]+p[pst][k]*(m[j]/i)%zwz*m[j])%zwz; 84 swap(now,pst); 85 } 86 for(k=0;k<=n;k++) 87 ans=(ans+p[pst][k]*(f[k][la][c]-f[k][i-1][c]+zwz)%zwz)%zwz; 88 } 89 printf("%d\n",ans); 90 91 } 92 return 0; 93 }
