https://www.cnblogs.com/CQzhangyu/p/7891363.html
不难推到$sumlimits_{D=1}^{m_1}sumlimits_{d|D}C_{d-1}^{c-2}mu(frac D d)prodlimits_{i=1}^nfrac {(2m_i-({lfloor frac {m_i} {D} floor}+1) imes D){lfloor frac {m_i} {D} floor}}{2}$。
$O(Tnm)$，可以拿80甚至100。
我们发现，求和部分与累积部分都含有D，这使分块加速变得困难。
化一下式子发现，当将$lfloorfrac{m}{D} floor$看作常数时，右边可以化成一个n次多项式。
$O(cmlog m+nmc)$预处理出多项式系数，$O(nsqrt{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 ",ans); 45 } 46 return 0; 47 }