[WC2014]时空穿梭(莫比乌斯反演)

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被卡。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=100010,M=100000,mod=10007,inv2=5004;
bool b[N];
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];

int main(){
    freopen("space.in","r",stdin);
    freopen("space.out","w",stdout);
    mu[1]=1;
    rep(i,2,M){
        if (!b[i]) p[++tot]=i,mu[i]=-1;
        for (int j=1; j<=tot && i*p[j]<=M; j++){
            b[i*p[j]]=1;
            if (i%p[j]==0) break;
            mu[i*p[j]]=-mu[i];
        }
    }
    rep(i,0,M){
        C[i][0]=1;
        rep(j,1,min(i,20)) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
    rep(i,1,M){
        int tmp=1;
        rep(k,0,11) sj[i][k]=tmp,tmp=1ll*tmp*i%mod;
    }
    rep(j,0,20) rep(i,1,M) if (mu[i])
        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;
    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;
    for (scanf("%d",&T); T--; ){
        scanf("%d%d",&n,&c); mn=N; ans=0;
        rep(i,1,n) scanf("%d",&m[i]),mn=min(mn,m[i]);
        for (int i=1,lst; i<=mn; i=lst+1){
            lst=mn; rep(j,1,n) lst=min(lst,m[j]/(m[j]/i));
            int tmp=1; rep(j,1,n) tmp=1ll*tmp*(m[j]/i)%mod*inv2%mod;
            memset(f,0,sizeof(f)); f[0]=1;
            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;
            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;
        }
        printf("%d
",ans);
    }
    return 0;
}