HDU 6439 Congruence equation（莫比乌斯反演）

Description

给出一个长度为$k$的序列$A_1,...,A_k$，定义$f(m)$为满足以下条件的$C$序列的数量：

1.若$A_i=-1$则$C_i$可以取$[0,m)$中任意整数，否则$C_i=A_i\%m$

2.同余方程$\sum _{i=1}^k{C}_i{x}_i\equiv 1(modm)$有解

给出一整数$n$，求$\sum _{m=1}^{n}f(m)$

Input

第一行一整数$T$表示用例组数，每组用例首先输入两个整数$k,n$，之后输入一长度为$k$的序列$A_i$

$(1\le T\le 250,1\le k\le 50,1\le n\le 10^9,-1\le A_i\le 10^9)$

Output

输出答案，结果模$10^9+7$

Sample Input

2
5 10
-1 -1 8 -1 -1
3 20
-1 6 18

Sample Output

Case #1: 24354
Case #2: 140

Solution

方程有解当且仅当$gcd(C,m)=1$，假设$A$序列中已经确定值的$gcd$为$g$，尚未确定的数字有$x$个

1.若$g̸=0$，记$f(d)$为使得$gcd(C,m)=d$的序列$C$的个数，$F(d)$为使得$d|gcd(C,m)$的序列$C$的个数，则显然有
$$F(d)=(\frac{m}{d}{)}^x$$
且由$F(d)=\sum _{d|i}f(i)$，莫比乌斯反演有
$$f(1)=\sum _{d|gcd(g,m)}\mu (d)F(d)=\sum _{d|gcd(g,m)}\mu (d)(\frac{m}{d}{)}^x$$
进而有
$$ans=\sum _{m=1}^n\sum _{d|gcd(g,m)}\mu (d)(\frac{m}{d}{)}^x=\sum _{d|g}\mu (d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }{i}^x=\sum _{d|g}\mu (d)S(\lfloor \frac{n}{d}\rfloor ,x)$$
其中$S(n,k)=\sum _{i=1}{i}^k$，通过伯努利序列即可$O(k^2)$预处理$O(k)$得到$S(n,k)$，那么该种情况即可$O(x\sqrt{g})$解决

2.若$g=0$，将之前的式子稍作修改有
$$ans=\sum _{m=1}^n\sum _{d|m}\mu (d)(\frac{m}{d}{)}^x=\sum _{d=1}^n\mu (d)S(\lfloor \frac{n}{d}\rfloor ,x)$$
由于$\lfloor \frac{n}{d}\rfloor$之后$O(\sqrt{n})$种取值，只要对于每种取值对应的$d$所处区间$[l,r]$，知道该区间的$\mu$之和即可，也即求出莫比乌斯函数前缀和即可

令$A(n)=\sum _{i=1}^n\mu (i)$，由于$\sum _{i|d}\mu (i)=[d=1]$，故有
$$1=\sum _{d=1}^n\sum _{i|d}\mu (i)=\sum _{i=1}^n\mu (i)\lfloor \frac{n}{i}\rfloor =\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu (j)=\sum _{i=1}^{n}A(\lfloor \frac{n}{i}\rfloor )$$
进而有
$$A(n)=1-\sum _{i=2}^{n}A(\lfloor \frac{n}{i}\rfloor )$$
同样分块加速即可，小范围前缀和可以线性筛预处理，较大值每次递归记忆化

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=1000005;
#define mod 1000000007
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int mu[maxn],p[maxn],res=0,vis[maxn],sum[maxn];
void pre(int n=1e6)
{
	mu[1]=sum[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])p[res++]=i,mu[i]=-1;
		sum[i]=sum[i-1]+mu[i];
		if(sum[i]<0)sum[i]+=mod;
		for(int j=0;j<res&&i*p[j]<=n;j++)
		{
			vis[i*p[j]]=1;
			if(i%p[j]==0)
			{
				mu[i*p[j]]=0;
				break;
			}
			mu[i*p[j]]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++)
		if(mu[i]==-1)mu[i]=mod-1;
}
int inv[60],B[60],C[60][60];
void init(int n=55)
{
	C[0][0]=1;
	for(int i=1;i<=n;i++)
	{
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;j++)C[i][j]=add(C[i-1][j-1],C[i-1][j]);
	} 
	inv[1]=1;
	for(int i=2;i<=n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
	B[0]=1;
	for(int i=1;i<n;i++)
	{
		for(int j=0;j<i;j++)B[i]=add(B[i],mul(C[i+1][j],B[j]));
		B[i]=mod-mul(B[i],inv[i+1]);
	}
}
int S(int n,int k)
{
	if(k==0)return n%mod;
	int ans=0,temp=1;
	for(int i=1;i<=k+1;i++)
	{
		temp=mul(temp,n+1);
		ans=add(ans,mul(mul(C[k+1][i],B[k+1-i]),temp));
	}
	return mul(ans,inv[k+1]);
}
int gcd(int a,int b)
{
	return b?gcd(b,a%b):a;
}
int get_mu(int n)
{
	if(n<=1e6)return mu[n];
	int ans=1;
	for(int i=2;i*i<=n;i++)
		if(n%i==0)
		{
			int num=0;
			while(n%i==0)n/=i,num++;
			if(num>1)return 0;
			ans*=-1;
		}
	if(n>1)ans*=-1;
	return add(ans,mod);
}
map<int,int>M;
int A(int n)
{
	if(n<=1e6)return sum[n];
	if(M.find(n)!=M.end())return M[n];
	int ans=1;
	for(int i=2,nxt;i<=n;i=nxt+1)
	{
		nxt=n/(n/i);
		ans=add(ans,mod-mul(nxt-i+1,A(n/i)));
	}
	return M[n]=ans;
}
int Solve(int l,int r,int n,int k)
{
	if(l==r)return mul(get_mu(l),S(n,k));
	return mul(add(A(r),mod-A(l-1)),S(n,k));
}
int Case=1,T,k,n,a[maxn];
int main()
{
	pre();
	init();
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d",&k,&n);
		int g=0,m=0,ans=0;
		for(int i=1;i<=k;i++)
		{
			scanf("%d",&a[i]);
			if(a[i]!=-1)g=gcd(g,a[i]);
			else m++;
		}
		if(g)
		{
			for(int i=1;i*i<=g;i++)
				if(g%i==0)
				{
					ans=add(ans,Solve(i,i,n/i,m));
					if(i*i!=g)ans=add(ans,Solve(g/i,g/i,n/(g/i),m));
				}
		}
		else
		{
			for(int i=1,nxt;i<=n;i=nxt+1)
			{
				nxt=n/(n/i);
				ans=add(ans,mul(add(A(nxt),mod-A(i-1)),S(n/i,m)));
			}
		}
		printf("Case #%d: %d\n",Case++,ans);
	}
	return 0;
}