[WC2014]时空穿梭
这才叫莫比乌斯反演题。
一、题目
二、解法
也没有什么好的思路,我们不妨把暴力柿子写出来,我们想枚举直线,但是这道题不能枚举直线的斜率,所以就要用整数来表示直线,我们不妨枚举出发点和终止点的向量差 \((x_1,x_2...x_n)\) ,那么起始点的方案数就是 \(\prod m-x_i\) ,剩余点的选取是在起始点的基础上加上这个向量 \((x_1/d,x_2/d...x_n/d)\) ,不难发现 \(d\) 必然满足是所有 \(x\) 的因数,也就是 \(\gcd(x)\) 的因数( \(d=0\) 也可以),所以直线上的点都是可以选的,由于已经确定了起终点,所以选剩下点的方案数是 \({\gcd(x)-1\choose c-2}\) ,写成柿子:
看到 \(\gcd\) 的柿子就一定要想莫比乌斯反演,我觉得网上的推导都不是很好,这种题是有套路的,按照我给的方法去推一定推得出来。第一步,我们在最前面枚举 \(\gcd(x)\) ,那么柿子变成这样:
第二步,把 \([\gcd(x)=1]\) 反演了,很多题都只需要反演这种结构:
第三步,由于 \(\mu\) 的求和结构放在里面很不舒服,我们把它提前,本题我们放在最外层求和的后边:
第四步,你发现出现 \(dj\) 这个整体结构,于是令 \(d=dj\) 结构一定会更简单(注意 \(d\) 的含义不同):
设 \(g(n)=\sum_{d|n}C(d-1,c-2)\mu(\frac{n}{d})\) 那么预处理 \(g\) 是需要 \(O(nc\log m)\) 的时间的(对于每一种可能的 \(c\) 我们都预处理),后面那一大块可以用乘法原理化简,那么柿子变成了这个样子:
暴力算他是 \(O(Tnm)\) 的,\(1e8\) 还是过不了,注意到出现了 \(m_i-x_id\) 这种结构,\(d\) 又是在前面枚举的,你难道对这个结构没有感觉吗?noip2020微信步数考过的啊!我们把 \(d\) 看成未知数,那么后面的那一块就可以变成一个关于 \(d\) 的多项式,对于 \(m_i/d\) 相同的那些 \(d\) ,他们的多项式是相同的,所以可以把前面改成一个分块,我们把它的系数给整出来,设系数为 \(dp_i\) ,\(d\in[l,r]\) 则柿子变为:
那么把 \(dp_i\) 放在前面去枚举:
后面的那部分相当于是一个前缀和,显然是可以预处理,预处理的数据要加上 \(n\) 的那维所以时间复杂度 \(O(cnm)\) ,现在来想一下算系数的复杂度,一共有 \(O(n\sqrt m)\) 个块(要保证 \(m_i/d\) 的值相同嘛),然后每个块的系数要 \(O(n^2)\) 去算,那么复杂度是 \(O(n^3\sqrt m)\) 的,\(n\) 特别小所以没关系。补充一种大佬的做法,由于系数是很多个 \(1\) 次多项式相乘来算的,根据分块的原理每次在端点处做退背包,这样做就可以优化到 \(O(n^2\sqrt m)\) (但是我太懒了没写)
然后计算一下分块的复杂度是 \(O(n^2\sqrt m)\) 的,所以我们就过掉了这道题!但是由于模数是 \(10007\) 所以算组合数要用杨辉三角哦,思路理清楚了写起来就不麻烦了!
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int N = 21;
const int M = 100005;
const int MOD = 10007;
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int T,n,m[M],c,up,ans,dp[N];
int cnt,mu[M],vis[M],p[M],C[M][N],g[M][N],sum[M][N][12];
void init(int n)
{
//算组合数
C[0][0]=1;
for(int i=1;i<=n;i++)
{
C[i][0]=1;
for(int j=1;j<20;j++)
C[i][j]=(C[i-1][j]+C[i-1][j-1])%MOD;
}
//欧拉筛
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt && i*p[j]<=n;j++)
{
vis[i*p[j]]=1;
if(i%p[j]==0) break;
mu[i*p[j]]=-mu[i];
}
}
//算g
for(int c=2;c<=20;c++)
{
for(int i=1;i<=n;i++)
{
for(int j=i;j<=n;j+=i)
g[j][c]=(g[j][c]+C[i-1][c-2]*mu[j/i])%MOD;
}
}
//算sum,也就是g(n,c)*n^k
for(int i=1;i<=n;i++)
for(int c=2;c<=20;c++)
{
int pw=1;
for(int k=0;k<=11;k++)
{
sum[i][c][k]=(sum[i-1][c][k]+g[i][c]*pw)%MOD;
pw=pw*i%MOD;
}
}
}
void work(int d)//算系数 -x_i*d+m_i
{
memset(dp,0,sizeof dp);dp[0]=1;
for(int i=1;i<=n;i++)
{//1-mi/d
int t=m[i]/d,x=-1ll*(t+1)*t/2%MOD,y=1ll*m[i]*t%MOD;
for(int j=i;j>=1;j--)
dp[j]=(dp[j]*y+dp[j-1]*x)%MOD;
dp[0]=dp[0]*y%MOD;//surprise!
}
}
void solve()
{
n=read();c=read();up=1e5;ans=0;
for(int i=1;i<=n;i++)
{
m[i]=read();
up=min(up,m[i]);
}
for(int l=1,r;l<=up;l=r+1)
{
r=1e5;
for(int i=1;i<=n;i++)
r=min(r,m[i]/(m[i]/l));
work(l);
for(int i=0;i<=n;i++)
ans=(ans+dp[i]*(sum[r][c][i]-sum[l-1][c][i]))%MOD;
}
printf("%d\n",(ans+MOD)%MOD);
}
signed main()
{
init(100000);
T=read();
while(T--)
{
solve();
}
}