BZOJ 3434 [Wc2014]时空穿梭

题目链接

https://lydsy.com/JudgeOnline/problem.php?id=3434

题解

枚举选取的第一个点和最后一个点的坐标差
$$\sum _{\Delta P=(\Delta x_1,\Delta x_2,\cdots ,\Delta x_n)}\left(\genfrac{}{}{0ex}{}{\gcd (\Delta P)-1}{c-2}\right)\prod _{k=1}^n(m_i-\Delta x_i+1)$$
反演得到
$$\sum _{T=1}^{{\min }_{i=1}^n{m}_i}g(c,T)\prod _{i=1}^{n}f(m_i,T)$$
其中
$$\begin{gathered} f(m,T)=-\frac{\lfloor \frac{m}{T}\rfloor (\lfloor \frac{m}{T}\rfloor +1)}{2}T+m\lfloor \frac{m}{T}\rfloor \\ g(c,T)=\sum _{d|T}\left(\genfrac{}{}{0ex}{}{d-1}{c-2}\right)\mu (\frac{T}{d}) \end{gathered}$$
观察发现
$$\prod _{i=1}^{n}f(m_i,T)$$
可以看作一个关于$T$的$n$次多项式，每一项的系数最多只有${\sum }_{i=1}^n\sqrt{m_i}$种，因此可以整除分块分别求系数。假设求出的系数分别为$a_0,a_1,\cdots ,a_n$，那么答案就是
$$\sum _{i=0}^n{a}_{i}H(c,i,\underset{i=1}{\overset{n}{\min }}{m}_i)$$
其中
$$H(c,i,k)=\sum _{T=1}^{k}g(c,T)T^i$$
可以预处理出$H$。

时间复杂度$O(ncm+\sum n^2{\sum }_{i=1}^n\sqrt{m_i})$。

代码

#include <cstdio>
#include <algorithm>

int read()
{
  int x=0,f=1;
  char ch=getchar();
  while((ch<'0')||(ch>'9'))
    {
      if(ch=='-')
        {
          f=-f;
        }
      ch=getchar();
    }
  while((ch>='0')&&(ch<='9'))
    {
      x=x*10+ch-'0';
      ch=getchar();
    }
  return x*f;
}

const int maxn=11;
const int maxc=20;
const int maxm=100000;
const int mod=10007;
const int inf=0x3f3f3f3f;

int p[maxm+10],prime[maxm+10],cnt,mu[maxm+10],C[maxm+10][maxc+2],g[maxc+2][maxm+10],h[maxc+2][maxn+2][maxm+10],power[maxm+10][maxn+2];

int getprime()
{
  p[1]=mu[1]=1;
  for(int i=2; i<=maxm; ++i)
    {
      if(!p[i])
        {
          prime[++cnt]=i;
          mu[i]=mod-1;
        }
      for(int j=1; (j<=cnt)&&(i*prime[j]<=maxm); ++j)
        {
          int x=i*prime[j];
          p[x]=1;
          if(i%prime[j]==0)
            {
              mu[x]=0;
              break;
            }
          mu[x]=mod-mu[i];
        }
    }
  for(int i=0; i<=maxm; ++i)
    {
      C[i][0]=1;
      for(int j=1; (j<=maxc)&&(j<=i); ++j)
        {
          C[i][j]=C[i-1][j]+C[i-1][j-1];
          if(C[i][j]>=mod)
            {
              C[i][j]-=mod;
            }
        }
    }
  for(int i=1; i<=maxm; ++i)
    {
      power[i][0]=1;
      for(int j=1; j<=maxn; ++j)
        {
          power[i][j]=power[i][j-1]*i%mod;
        }
    }
  for(int c=2; c<=maxc; ++c)
    {
      for(int d=1; d<=maxm; ++d)
        {
          for(int T=d; T<=maxm; T+=d)
            {
              g[c][T]=(g[c][T]+C[d-1][c-2]*mu[T/d])%mod;
            }
        }
    }
  for(int c=2; c<=maxc; ++c)
    {
      for(int i=0; i<=maxn; ++i)
        {
          for(int k=1; k<=maxm; ++k)
            {
              h[c][i][k]=(h[c][i][k-1]+g[c][k]*power[k][i])%mod;
            }
        }
    }
  return 0;
}

int T,n,c,m[maxn+2],minm,a[maxn+2];

int main()
{
  getprime();
  T=read();
  while(T--)
    {
      n=read();
      c=read();
      minm=inf;
      for(int i=1; i<=n; ++i)
        {
          m[i]=read();
          minm=std::min(minm,m[i]);
        }
      int ans=0;
      for(int l=1,r; l<=minm; l=r+1)
        {
          r=inf;
          for(int i=1; i<=n; ++i)
            {
              r=std::min(r,m[i]/(m[i]/l));
            }
          for(int i=0; i<=n; ++i)
            {
              a[i]=(i==0);
            }
          for(int i=1; i<=n; ++i)
            {
              int k=m[i]/l,x=1ll*m[i]*k%mod,y=(1ll*k*(k+1)/2)%mod*(mod-1)%mod;
              for(int j=n; j; --j)
                {
                  a[j]=(x*a[j]+y*a[j-1])%mod;
                }
              a[0]=a[0]*m[i]%mod*k%mod;
            }
          for(int i=0; i<=n; ++i)
            {
              ans=(ans+a[i]*(h[c][i][r]-h[c][i][l-1]+mod))%mod;
            }
        }
      printf("%d\n",ans);
    }
  return 0;
}