好难啊..反演的终极题目
首先，本题的突破口在于直线的性质。不论是几维的空间，两点一定能确定一条直线
选取两个点作为最左下和最右上的点！
假设现在是二维空间，选取了$(x1,y1)$和$(x2,y2)$两个点，那么它们连线上经过的点数就是$gcd(x2-x1,y2-y1)-1$
选取的方案数为$left ( _{gcd-1}^{c-2} ight )$
发现这个方案数只和坐标的差值有关，我们直接枚举差值就行
接下来就是反演了，为了方便叙述，以$n=2$为例
$sumlimits_{x=1}^{m_{1}}sumlimits_{y=1}^{m_{2}} left ( _{c-2}^{gcd(x,y)-1} ight ) (m_{1}-x)(m_{2}-y)$
$sumlimits_{k=1}^{m} left ( _{c-2}^{k-1} ight ) sumlimits_{x=1}^{m_{1}} sumlimits_{y=1}^{m_{2}} [gcd(x,y)==k] (m_{1}-x)(m_{2}-y)$
$sumlimits_{k=1}^{m} left ( _{c-2}^{k-1} ight ) sumlimits_{x=1}^{left lfloor frac{m_{1}}{k} ight floor} sumlimits_{y=1}^{left lfloor frac{m_{2}}{k} ight floor} [gcd(x,y)==1] (m_{1}-xk)(m_{2}-yk)$
$sumlimits_{k=1}^{m} left ( _{c-2}^{k-1} ight ) sumlimits_{x=1}^{left lfloor frac{m_{1}}{k} ight floor} sumlimits_{y=1}^{left lfloor frac{m_{2}}{k} ight floor} sumlimits_{d|k}mu(d)(m_{1}-xk)(m_{2}-yk)$
$sumlimits_{k=1}^{m} left ( _{c-2}^{k-1} ight ) sumlimits_{d=1}^{left lfloor frac{m_{1}}{k} ight floor} mu(d) sumlimits_{x=1}^{left lfloor frac{m_{1}}{kd} ight floor} (m_{1}-xkd) sumlimits_{y=1}^{ left lfloor frac{m_{2}}{kd} ight floor} (m_{2}-ykd)$
令$Q=kd$
$sumlimits_{Q=1}^{m} left ( sumlimits_{k|Q} left ( _{c-2}^{k-1} ight ) mu(frac{Q}{k}) ight ) left ( sumlimits_{x=1}^{left lfloor frac{m_{1}}{Q} ight floor} (m_{1}-xQ) ight ) left ( sumlimits_{y=1}^{ left lfloor frac{m_{2}}{Q} ight floor} (m_{2}-yQ) ight ) $
利用等差数列公式 $left ( sumlimits_{x=1}^{left lfloor frac{m}{Q} ight floor} (m-xQ) ight ) = left ( left lfloor frac{m}{Q} ight floor m- frac{ (1+left lfloor frac{m}{Q} ight floor )left lfloor frac{m}{Q} ight floor Q}{2} ight ) $
$sumlimits_{Q=1}^{m} left ( sumlimits_{k|Q} left ( _{c-2}^{k-1} ight ) mu(frac{Q}{k}) ight ) left ( left lfloor frac{m_{1}}{Q} ight floor m_{1}- frac{ (1+left lfloor frac{m_{1}}{Q} ight floor )left lfloor frac{m_{1}}{Q} ight floor Q}{2} ight ) left ( left lfloor frac{m_{2}}{Q} ight floor m_{2}- frac{ (1+left lfloor frac{m_{2}}{Q} ight floor )left lfloor frac{m_{2}}{Q} ight floor Q}{2} ight ) $
推广到更高维上
$sumlimits_{Q=1}^{m} left ( sumlimits_{k|Q} left ( _{c-2}^{k-1} ight ) mu(frac{Q}{k}) ight ) prod_{t=1}^{n} left ( left lfloor frac{m_{t}}{Q} ight floor m_{t}- frac{ (1+left lfloor frac{m_{t}}{Q} ight floor )left lfloor frac{m_{t}}{Q} ight floor Q}{2} ight ) $
我们不能把$Q$这一项带入到整除分块里去
所以每次$O(n^{2})$暴力计算出$Q^{k}(k=0,1,2...n)$的系数
预处理出$f(Q,c)=sumlimits_{k|Q} left ( _{c-2}^{k-1} ight ) 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 ",ans); 90 91 } 92 return 0; 93 }