HDU 6439 2018CCPC網路賽 Congruence equationI（杜教篩 + 莫比烏斯反演 + 伯努利數）
大致題意：給你一個長度為k的序列a。對於序列c，當 時，
；當
時，
取[0,m)中任意一個數字。令
表示滿足
的序列c的方案數。現在讓你求
。
首先，根據裴蜀定理，滿足的條件是
，那麼我們不妨分為兩種情況處理。對於
的數字，假設他們的gcd為g，那麼剩下的數字與g的gcd就要是1。設
的項有k個，加上這個m，設這k+1個數字的gcd為d，那麼gcd(d,g)要等於1。由於這k+1個數字裡面有一個定值m，所以這個d一定是m的因子。我們令f(d)表示這k+1個數字的gcd為d的方案數。那麼開始第一次莫比烏斯反演，有：
令，那麼最後的答案就是:
令，於是我們開始第二次莫比烏斯反演：
這樣，對於g，我們只需要用到gcd為1的，所以我們不妨把第二個引數去掉。整理一下，最後的答案就是：
對於這個g(x)，根據莫比烏斯函式的性質，有效的i肯定不會是任意一個質因子的2次冪及以上，所以i一定是x的質因子的線性組合，因此我們可以把這個預處理出來。然後
我們發現，如果令，那麼h(d)就是H(d)和
的迪利克雷卷積，而H(d)和
是顯然具有積性的。所以說我們可以用積性函式的性質，構造線性篩來求解h(d)的字首和。但是注意到，本題的資料範圍是
，而且還有多組資料，即使是線性的篩法也無法滿足條件。所以我們這裡考慮用杜教篩。
簡單來說就是，利用迪利克雷卷積的恆等變換，使得一個原本不容易求的積性數論函式的字首和，變成兩個容易求的積性數論函式運算，最後轉化成的形式，其中A(n)表示輔助的容易求和的積性數論函式的字首和。之後對於前
對於本題，令，我們可以這麼推導：
本題的話選用伯努利數的方法來計算冪和，即 。經過曲折之後，我們就可以在O(
)的時間複雜度內求出這個h(d)的字首和S(d)。最後的答案
，分成兩部分分塊求和即可。最後總的時間複雜度是
。具體見程式碼：
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
#define pb push_back
#define lb lower_bound
#define ub upper_bound
#define INF 0x3f3f3f3f
#define sf(x) scanf("%d",&x)
#include <ext/pb_ds/hash_policy.hpp>
#include <ext/pb_ds/assoc_container.hpp>
#define sc(x,y,z) scanf("%d%d%d",&x,&y,&z)
#define clr(x,n) memset(x,0,sizeof(x[0])*(n+5))
#define file(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
using namespace __gnu_pbds;
using namespace std;
const int N = 1e6 + 10;
const int K = 1e2 + 10;
int n,m,k,tot,S[N],C[K][K],p[N];
gp_hash_table<int,int> mp;
std::vector<int> pri_fac;
typedef pair<int,int> P;
LL B[K],invk,pw[N];
std::vector<P> fac;
bool isp[N];
LL qpow(LL x,int n)
{
LL res=1;
while(n)
{
if(n&1) res=res*x%mod;
x=x*x%mod; n>>=1;
}
return res;
}
void init()
{
C[0][0]=1;
for(int i=1;i<K;i++) C[i][0]=1;
for(int i=1;i<K;i++)
for(int j=1;j<K;j++)
C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
B[0]=1;
for(int i=1;i<K;i++)
{
for(int j=0;j<i;j++)
B[i]+=C[i+1][j]*B[j]%mod;
B[i]%=mod; B[i]=(mod-B[i])*qpow(i+1,mod-2)%mod;
}
}
void sieve(int n)
{
int sz=0; S[1]=1;
for(int i=2;i<n;i++)
{
if (!isp[i])
{
p[++sz]=i;
pw[i]=qpow(i,k);
S[i]=pw[i]-1;
}
for(int j=1;j<=sz&&p[j]*i<n;j++)
{
int x=i*p[j]; isp[x]=1;
if (i%p[j]==0)
{
S[x]=S[i]*pw[p[j]]%mod;
break;
} else S[x]=S[i]*(pw[p[j]]-1)%mod;
}
}
for(int i=2;i<n;i++) S[i]=(S[i-1]+S[i])%mod;
}
int powsum(int x)
{
if (k==0) return x;
LL res=0,pw=1;
for(int i=1;i<=k+1;i++)
{
pw=pw*(x+1)%mod;
res+=C[k+1][i]*B[k+1-i]%mod*pw%mod;
}
res%=mod;
return res*invk%mod;
}
int s(int x)
{
if (x<tot) return S[x];
if (mp[x]) return mp[x];
LL res=powsum(x);
for(int l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
res-=(r-l+1)*(LL)s(x/l)%mod;
}
res=res%mod+mod;
return mp[x]=res%mod;
}
LL cal(int x)
{
LL res=0;
for(auto i:fac)
{
if (i.first>x) break;
res+=(x/i.first)*(i.second&1?-1:1);
}
return res%mod;
}
int main()
{
int T; sf(T);
init(); int AC=1;
while(T--)
{
sf(n); sf(m);
int g=0; k=0; mp.clear();
for(int i=1;i<=n;i++)
{
int x; sf(x);
if(x==-1)k++;else g=__gcd(g,x);
}
sieve(tot=ceil(pow(m,2.0/3)));
invk=qpow(k+1,mod-2);
if (g==0)
{
printf("Case #%d: %d\n",AC++,s(m));
continue;
}
LL ans=0;
fac.clear();
pri_fac.clear();
for(int i=2;i*i<=g;i++)
{
if (g%i) continue;
pri_fac.pb(i);while(g%i==0) g/=i;
}
if (g>1) pri_fac.pb(g);
int up=1<<pri_fac.size();
for(int i=0;i<up;i++)
{
int cnt=0,d=1;
for(int j=0;j<pri_fac.size();j++)
if (i&(1<<j)) cnt++,d*=pri_fac[j];
fac.pb(P(d,cnt));
}
sort(fac.begin(),fac.end());
int pre=0,cur;
for(int l=1,r;l<=m;l=r+1)
{
r=m/(m/l); cur=s(r);
ans+=(cur-pre)*cal(m/l)%mod;
pre=cur;
}
ans%=mod;
printf("Case #%d: %lld\n",AC++,(ans+mod)%mod);
}
return 0;
}