[WC2014]时空穿梭(莫比乌斯反演)-CSDN博客

https://www.cnblogs.com/CQzhangyu/p/7891363.html

不难推到$\sum\limits_{D=1}^{m_1}\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)\prod\limits_{i=1}^n\frac {(2m_i-({\lfloor \frac {m_i} {D} \rfloor}+1)\times D){\lfloor \frac {m_i} {D} \rfloor}}{2}$。

$O(Tnm)$，可以拿80甚至100。

我们发现，求和部分与累积部分都含有D，这使分块加速变得困难。

化一下式子发现，当将$\lfloor\frac{m}{D}\rfloor$看作常数时，右边可以化成一个n次多项式。

$O(cm\log m+nmc)$预处理出多项式系数，$O(n\sqrt{m})$整除分块，$O(n^2)$暴力求多项式系数即可。

常数过大BZ被卡。

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
 5 using namespace std;
 6 
 7 const int N=100010,M=100000,mod=10007,inv2=5004;
 8 bool b[N];
 9 int T,c,n,tot,mn,ans,m[12],p[N],mu[N],C[N][22],s[N][21][12],sj[N][12],f[12],g[N][22];
10 
11 int main(){
12     freopen("space.in","r",stdin);
13     freopen("space.out","w",stdout);
14     mu[1]=1;
15     rep(i,2,M){
16         if (!b[i]) p[++tot]=i,mu[i]=-1;
17         for (int j=1; j<=tot && i*p[j]<=M; j++){
18             b[i*p[j]]=1;
19             if (i%p[j]==0) break;
20             mu[i*p[j]]=-mu[i];
21         }
22     }
23     rep(i,0,M){
24         C[i][0]=1;
25         rep(j,1,min(i,20)) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
26     }
27     rep(i,1,M){
28         int tmp=1;
29         rep(k,0,11) sj[i][k]=tmp,tmp=1ll*tmp*i%mod;
30     }
31     rep(j,0,20) rep(i,1,M) if (mu[i])
32         for (int k=i; k<=M; k+=i) g[k][j]=(g[k][j]+1ll*mu[i]*C[k/i-1][j]%mod+mod)%mod;
33     rep(i,1,M) rep(j,0,20) rep(k,0,11) s[i][j][k]=(s[i-1][j][k]+1ll*g[i][j]*sj[i][k])%mod;
34     for (scanf("%d",&T); T--; ){
35         scanf("%d%d",&n,&c); mn=N; ans=0;
36         rep(i,1,n) scanf("%d",&m[i]),mn=min(mn,m[i]);
37         for (int i=1,lst; i<=mn; i=lst+1){
38             lst=mn; rep(j,1,n) lst=min(lst,m[j]/(m[j]/i));
39             int tmp=1; rep(j,1,n) tmp=1ll*tmp*(m[j]/i)%mod*inv2%mod;
40             memset(f,0,sizeof(f)); f[0]=1;
41             rep(j,1,n) for (int k=j; ~k; k--) f[k]=(2ll*f[k]*m[j]%mod-1ll*f[k-1]*(m[j]/i+1)%mod+mod)%mod;
42             rep(j,0,n) ans=(ans+1ll*tmp*f[j]%mod*(s[lst][c-2][j]-s[i-1][c-2][j]+mod)%mod)%mod;
43         }
44         printf("%d\n",ans);
45     }
46     return 0;
47 }