【BZOJ3434】[Wc2014]时空穿梭 莫比乌斯反演

【BZOJ3434】[Wc2014]时空穿梭

Description

Input

第一行包含一个正整数T，表示有T组数据求解
每组数据包含两行，第一行包含两个正整数N,C(c>=2),分别表示空间的
维数和需要选择的暂停点个数
第二行包含N个正整数，依次表示M1,M2....Mn

Output

有T行，每行一个非负整数，依次对应每组数据的答案。

Sample Input

3 
2 3 
3 4 
3 3 
3 4 4 
4 4 
5 9 7 8

Sample Output

2 
4 
846

HINT

样例数据第一组共有两种可行方案：一种是选择（1,1），（ 2,2） ,（ 3,3） ，另一种是选择 （ 1,2） ,（ 2,3） ,（ 3,4） 。

T<=1000,N<=11,C<=20,Mi<=100000

题解：BZOJ3518的加强版，推导过程倒不难，就是巨麻烦~~~

我们设由1号点指向c号点所构成的向量为（p1,p2...pn），那么我们枚举p1,p2...pn的值，于是中间的点可以在这条线段经过的所有整点中取，设d=gcd(p1,p2...pn)，那么这条线段经过的整点数就是d-1(不算起点和终点)。并且这条线段还能在空间中不断平移，比如起点的第i位坐标可以在[1..mi-pi]中随便取。于是经过一番整理可以得到所求：

$ans=\sum\limits_{p_1=1}^{m_1}(m_1-p_1)\sum\limits_{p_2=1}^{m_2}(m_2-p_2)...\sum\limits_{p_n=1}^{m_n}(m_n-p_n)\times C_{gcd(p_1,p_2...p_n)-1}^{c-2}$

然后就是熟悉的莫比乌斯反演啦！我们设gcd=d，然后枚举d：

$ans=\sum\limits_{p_1=1}^{m_1}(m_1-p_1)\sum\limits_{p_2=1}^{m_2}(m_2-p_2)...\sum\limits_{p_n=1}^{m_n}(m_n-p_n)\times\sum\limits_{d=1}^{m_1}[d==gcd(p_1,p_2...p_n)] C_{d-1}^{c-2}\\=\sum\limits_{d=1}^{m_1}C_{d-1}^{c-2}\sum\limits_{p_1=1}^{\lfloor \frac {m_1} d \rfloor}(m_1-d\times p_1)\sum\limits_{p_2=1}^{\lfloor \frac {m_2} d \rfloor}(m_2-d\times p_2)...\sum\limits_{p_n=1}^{\lfloor \frac {m_n} d \rfloor}(m_n-d\times p_n)[gcd(p_1,p_2..p_n)==1]\\=\sum\limits_{d=1}^{m_1}C_{d-1}^{c-2}\sum\limits_{e=1}^{\lfloor \frac {m_1} d \rfloor}\mu(e)\sum\limits_{p_1=1}^{\lfloor \frac {m_1} {de} \rfloor}(m_1-de\times p_1)\sum\limits_{p_2=1}^{\lfloor \frac {m_2} {de} \rfloor}(m_2-de\times p_2)...\sum\limits_{p_n=1}^{\lfloor \frac {m_n} {de} \rfloor}(m_n-de\times p_n)$

令D=de：

$=\sum\limits_{D=1}^{m_1}\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)\prod\limits_{i=1}^n\sum\limits_{p_i=1}^{\lfloor \frac {m_i} {D} \rfloor}(m_i-D\times p_i)$

接着上等差数列求和公式：

$=\sum\limits_{D=1}^{m_1}\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)\prod\limits_{i=1}^n\frac {(2m_i-({\lfloor \frac {m_i} {D} \rfloor}+1)\times D){\lfloor \frac {m_i} {D} \rfloor}} 2$

到这里，我们就已经得到了一个$O(Tnm)$的做法，不过为了用前缀和来维护D，我们需要将右侧拆成关于D的多项式。因为我们过一会要对${\lfloor \frac {m_i} {D} \rfloor}$进行分块，所以我们可以将此时的${\lfloor \frac {m_i} {D} \rfloor}$看成常数，所以，右侧实质上就是一个n次的多项式，我们只需$O(n^2)$求出每一项的系数，然后套上D的j次方的前缀和即可。具体地，我们设第j项的系数为f(j)：

$=\sum\limits_{D=1}^{m_1}\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)\prod\limits_{i=1}^n\frac {{\lfloor \frac {m_i} {D} \rfloor}} 2\sum\limits_{j=0}^nD^jf(j)\\=\sum\limits_{D=1}^{m_1}\prod\limits_{i=1}^n\frac {{\lfloor \frac {m_i} {D} \rfloor}} 2\sum\limits_{j=0}^nf(j)D^j\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)$

所以呢，我们要维护的就是$D^j\sum\limits_{d|D}C_{d-1}^{c-2}\mu(\frac D d)$的前缀和，设其为s[D][c-2][j]，再套上分块就做完啦！

其实思路还是挺乱的，我们回顾一下全过程。

1.先预处理s数组
2.然后分块，对于每块，求出f数组
3.代入到多项式中计算答案

时间复杂度呢？

预处理s数组的复杂度=$O(cm\times ln_m+nmc)$
由于有n个m，所以我们一共分成了$n\sqrt m $块，每块我们都要用$n^2$的时间处理出f数组并用n的时间求和，所以总复杂度=$O(cm\times ln_m+nmc+Tn^3\sqrt m)$（这也能过？）

然而BZ上TLE了，不过UOJ上是AC的，辣鸡BZ。

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int P=10007;
const int ine=5004;
const int N=100000;
int num;
bool np[N+5];
int pri[N],mu[N+5];
int c[N+5][22],s[N+5][21][12],sj[N+5][12],f[12],g[N+5][22];
int T,C;
int n,m[12],mn,ans;
int main()
{
	register int i,j,k,last,tmp;
	mu[1]=1;
	for(i=2;i<=N;i++)
	{
		if(!np[i])	pri[++num]=i,mu[i]=-1;
		for(j=1;j<=num&&i*pri[j]<=N;j++)
		{
			np[i*pri[j]]=1;
			if(i%pri[j]==0)	break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(i=0;i<=N;i++)	for(c[i][0]=j=1;j<=i&&j<=20;j++)	c[i][j]=(c[i-1][j-1]+c[i-1][j])%P;
	for(i=1;i<=N;i++)	for(k=0,tmp=1;k<=11;k++,tmp=tmp*i%P)	sj[i][k]=tmp;
	for(j=0;j<=20;j++)	for(i=1;i<=N;i++)	if(mu[i])	for(k=i;k<=N;k+=i)	g[k][j]=g[k][j]+mu[i]*c[k/i-1][j];
	for(i=1;i<=N;i++)	for(j=0;j<=20;j++)	for(g[i][j]%=P,k=0;k<=11;k++)
		s[i][j][k]=(s[i-1][j][k]+g[i][j]*sj[i][k])%P;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d",&n,&C),mn=N,ans=0;
		for(i=1;i<=n;i++)	scanf("%d",&m[i]),mn=min(mn,m[i]);
		for(i=1;i<=mn;i=last+1)
		{
			for(last=mn,j=1;j<=n;j++)	last=min(last,m[j]/(m[j]/i));
			for(tmp=j=1;j<=n;j++)	tmp=tmp*(m[j]/i)%P*ine%P;
			memset(f,0,sizeof(f)),f[0]=1;
			for(j=1;j<=n;j++)	for(k=j;k>=0;k--)	f[k]=(f[k]*m[j]*2%P-f[k-1]*(m[j]/i+1))%P;
			for(j=0;j<=n;j++)	ans=(ans+tmp*f[j]%P*(s[last][C-2][j]-s[i-1][C-2][j]))%P;
		}
		printf("%d\n",(ans+P)%P);
	}
	return 0;
}//1 2 3 99997 99953