[WC2014]时空穿梭-题解

【题目地址】[IN]

题目意思见原题面。

---

我们首先来看一些前置知识：

• 对于一个线段，我们知道其两个端点坐标为$(x_1,y_1),(x_2,y_2)$，那么它上面的整点个数为$gcd(|x_1-x_2|,|y_1-y_2|)+1$（+1是开始的端点）。

我们很容易的将其推广到多维的情况上面。

那么我们考虑枚举坐标每维的最大差（也就是端点坐标差），那么一条直线的上的点的个数为$gcd(\Delta x_i)+1$，除去端点上，其余可选点数为$gcd(\Delta x_i)-1$，那么我们还需选择$c-2$个点，所以对于一条直线的所有方案数就为$C_{gcd(\Delta x_i)-1}^{c-2}$，然后我们考虑对于一种坐标差的情况，直线有多少条，显然对于每维坐标选择的方案数就是$m_i-\Delta x_i$（不加1是因为坐标大于等于1）,所以总的直线选择方案数为

$$\prod _i^n(m_i-\Delta x_i)$$

所以题目就变成了我们要求这个式子：

$$\sum _{\Delta x_1=1}^{m_1}\sum _{\Delta x_2=1}^{m_2}\cdots \sum _{\Delta x_n=1}^{m_n}{C}_{gcd(\Delta x_i)-1}^{c-2}\prod _i^n(m_i-\Delta x_i)$$

暴力求的话，复杂度非常之高，所以我们要考虑优化。

看到有gcd，一般可以想到用莫比乌斯反演的套路，我们来枚举gcd，则原式变成了：（令$m=\min (m_i)$）
$$\sum _{g=1}^m{C}_{g-1}^{c-2}\sum _{\Delta x_1=1}^{\frac{m_1}{g}}\sum _{\Delta x_2=1}^{\frac{m_2}{g}}\cdots \sum _{\Delta x_n=1}^{\frac{m_n}{g}}\prod _i^n(m_i-\Delta x_i)[gcd(\Delta x_i)=1]$$

我们用莫比乌斯函数来代替$[gcd(\Delta x_i)=1]$则可以得到：

$$\sum _{g=1}^m{C}_{g-1}^{c-2}\sum _{\Delta x_1=1}^{\frac{m_1}{g}}\sum _{\Delta x_2=1}^{\frac{m_2}{g}}\cdots \sum _{\Delta x_n=1}^{\frac{m_n}{g}}\prod _i^n(m_i-\Delta x_{i}g)\sum _{t|gcd(\Delta x_i)}\mu (t)$$

我们交换求和顺序，将$\mu$放在前面，则可以得到：

$$\sum _{g=1}^m{C}_{g-1}^{c-2}\sum _{t=1}^{\frac{m}{g}}\mu (t)\sum _{\Delta x_i=1}^{\frac{m_i}{gt}}\prod _{j=1}^n(m_j-\Delta x_{j}gt)$$

交换一下后面的求和顺序：

$$\sum _{g=1}^m{C}_{g-1}^{c-2}\sum _{t=1}^{\frac{m}{g}}\mu (t)\prod _{i=1}^n\sum _{x_i=1}^{\frac{m_i}{gt}}(m_i-\Delta x_{i}gt)$$

我们发现 ：
$$\sum _{x_i=1}^{\frac{m_i}{gt}}(m_i-\Delta x_{i}gt)$$
这个可以直接用等差数列求和公式$O(1)$求出，我们令：
$$f(m_i,gt)=\lfloor \frac{m_i}{gt}\rfloor m_i-gt\frac{\lfloor \frac{m_i}{gt}\rfloor \left(\lfloor \frac{m_i}{gt}\rfloor +1\right)}{2}$$

原式就变成：

$$\sum _{g=1}^m{C}_{g-1}^{c-2}\sum _{t=1}^{\frac{m}{g}}\mu (t)\prod _{i=1}^{n}f(m_i,gt)$$

但是此时暴力求的话复杂度还是$O(n{m}^2)$，所以我们考虑进一步优化，根据套路，我们令$h=gt$然后来枚举$h$可以得：

$$\sum _{h=1}^m\sum _{g|h}{C}_{g-1}^{c-2}\mu \left(\frac{h}{g}\right)\prod _{i=1}^{n}f(m_i,h)$$

我们显然可以枚举倍数，在$O(mlogm)$的时间里预处理出${\sum }_{g|h}{C}_{g-1}^{c-2}\mu \left(\frac{h}{g}\right)$，然后复杂度就优化成了$O(nm)$，大概可以过$50\sim 60pts$，所以我们还要考虑优化。

我们看$f(m,h)$中的$\lfloor \frac{m}{h}\rfloor$在所有的$\left(\frac{m_i}{h}\right)$中根据除法分块的理论最多只有$n\sqrt{m}$种取值，所以我们可以将其分成$n\sqrt{m}$段来求取。

而对于每一段，我们考虑这样的式子：（令$g(h)={\sum }_{g|h}{C}_{g-1}^{c-2}\mu \left(\frac{h}{g}\right),l(h)={\prod }_{i=1}^{n}f(m_i,h)$）
$$\sum _{h=l}^{r}g(h)l(h)$$

如果直接计算，每次计算需要$O(nm)$的时间，而总共有$n\sqrt{m}$段，复杂度还是$O(n^{2}m\sqrt{m})$的，多组询问下是过不去的。

所以我们考虑进一步优化，我们可以将$l$看成是一个关于$h$的函数，那么将$h$看成未知数，我们可以得到$l(x)={\sum }_{i=0}^n{a}_i{x}^i$，那么带回原式则可以得到：

$$\sum _{h=l}^{r}g(h)\sum _{i=0}^n{a}_i{h}^i$$

交换一下求和顺序可以得到：

$$\sum _{i=0}^n{a}_i\sum _{h=l}^{r}g(h)h^i$$

我们同样可以预处理出$g(h)h^i$前缀和，那么每次计算的复杂度就是$O(n)$的了，求系数$a_i$我们可以考虑将$f$暴力乘起来，化简便可以得到，复杂度为$O(n^2)$（你要用$FFT$优化也可以，不过太麻烦了）。所以此时每次询问的复杂度就变成了$O(n^3\sqrt{m})$了，就可以过了。

总的复杂度为$O(nmc+cmlogm+T{n}^3\sqrt{m})$

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int inf=2e9;
const int M=25,Mod=10007,U=1110,W=16;
int n[U],c[U],T,m[U][W],Maxv,Maxc;
const int N=1e5+20;
bool vis[N];
int prime[N],mu[N],cnt;
int g[M][N],f[M][W][N],C[N][M];
void init(){
    C[0][0]=1;
    for(int i=1;i<=Maxv;i++){
        C[i][0]=1;if(i<=Maxc-2)C[i][i]=1;
        for(int j=1;j<min(Maxc-1,i);j++){
            C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
        }
    }	mu[1]=1;
    for(int i=2;i<=Maxv;i++){
        if(!vis[i]){
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=Maxv;j++){
            vis[i*prime[j]]=1;
            if(!(i%prime[j])){
                break;
            }	mu[i*prime[j]]=-mu[i];
        }
    }
    for(int cv=2;cv<=Maxc;cv++){
        for(int i=1;i<=Maxv;i++){
            for(int j=i,k=1;j<=Maxv;j+=i,++k){
                (g[cv][j]+=C[i-1][cv-2]*mu[k]%Mod)%=Mod;
                if(g[cv][j]<0)g[cv][j]+=Mod;
            }
        }
    }
    for(int hv=1;hv<=Maxv;hv++){
        for(int cv=1;cv<=Maxc;cv++){
            for(int j=0,x=1;j<=11;j++){
                f[cv][j][hv]=(g[cv][hv]*x%Mod+f[cv][j][hv-1])%Mod;
                x=(x*hv)%Mod;
                if(f[cv][j][hv]<0)f[cv][j][hv]+=Mod;
            }
        }
    }
} 
int A[W];
void Pre(int h,int *mv,int nn){
    int x1,x0,t;memset(A,0,sizeof(A));A[0]=1;
    for(int i=1;i<=nn;i++){
        t=(mv[i]/h);
        x1=(-((1ll*t*(t+1))>>1))%Mod;
        x0=(1ll*mv[i]*t)%Mod;
        for(int j=nn;j>=1;j--){
            A[j]=(A[j]*x0%Mod+A[j-1]*x1%Mod)%Mod;
            if(A[j]<0)A[j]+=Mod;
        }	A[0]=(A[0]*x0)%Mod;if(A[0]<0)A[0]+=Mod;
    }
}
int Minv,ans;
int main(){
    scanf("%d",&T);
    for(int i=1;i<=T;i++){
        scanf("%d%d",&n[i],&c[i]);if(Maxc<c[i])Maxc=c[i];
        for(int j=1;j<=n[i];j++){
            scanf("%d",&m[i][j]);if(Maxv<m[i][j])Maxv=m[i][j];
        }
    }
    init();
    for(int w=1;w<=T;w++){
        Minv=inf;ans=0;
        for(int i=1;i<=n[w];i++)if(m[w][i]<Minv)Minv=m[w][i];
        for(int i=1,j;i<=Minv;i=j+1){
            j=Minv;for(int k=1;k<=n[w];k++)j=min(j,m[w][k]/(m[w][k]/i));
            Pre(i,m[w],n[w]);for(int k=0;k<=n[w];k++) ans=(ans+A[k]*(f[c[w]][k][j]-f[c[w]][k][i-1])%Mod)%Mod;
            if(ans<0)ans+=Mod;
        }	printf("%d\n",(ans%Mod+Mod)%Mod);
    }
    return 0;
}