BZOJ-4196 Lucas的數論(莫比烏斯反演+杜教篩)
阿新 • • 發佈:2020-11-23
題目描述
計算:
\[\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_0(ij) \]其中 \(\sigma_0(x)\) 表示 \(x\) 的約數個數,\(n\leq 10^9\),答案對 \(10^9+7\) 取模。
分析
首先有一個結論:
\[\sigma_0(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1] \]證明:
在等號左邊,對於 \(ij\) 的每個質因子 \(p\),設 \(i=i'\times p^{k_1},j=j'\times p^{k_2}\)(其中 \(k_1,k_2\) 可以為 \(0\)),則 \(ij=i'\times j'\times p^{k_1+k_2}\)
因此:
\[\begin{aligned}&\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma_0(ij)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]\\ =&\sum_{x=1}^{n}\sum_{y=1}^{m}\Big\lfloor\frac{n}{x}\Big\rfloor\Big\lfloor\frac{m}{y}\Big\rfloor[\gcd(x,y)=1]\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor[\gcd(i,j)=1]\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor\sum_{d\mid \gcd(i,j)}\mu(d)\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{m}[d\mid \gcd(i,j)]\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\Big\lfloor\frac{n}{di}\Big\rfloor\Big\lfloor\frac{m}{dj}\Big\rfloor\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\Big(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\Big\lfloor\frac{n}{di}\Big\rfloor\Big)\Big(\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\Big\lfloor\frac{m}{dj}\Big\rfloor\Big) \end{aligned} \]\(\mu\) 的字首和用杜教篩解決,後面兩部分數論分塊處理,然後分塊套分塊。
程式碼
#include<bits/stdc++.h> using namespace std; const int N=2e6+10;//n^(2/3) const int mod=1e9+7; long long prime[N+10],vis[N+10],mu[N+10],sum[N+10],cnt; void init() { mu[1]=1; for(int i=2;i<=N;i++) { if(!vis[i]) { prime[++cnt]=i; mu[i]=-1; } for(int j=1;j<=cnt&&i*prime[j]<=N;j++) { vis[i*prime[j]]=1; if(i%prime[j]==0) { mu[i*prime[j]]=0; break; } else mu[i*prime[j]]=-mu[i]; } } for(int i=1;i<=N;i++) sum[i]=sum[i-1]+mu[i]; } map<long long,long long> mp; long long S_mu(long long n) { if(n<N) return sum[n]; if(mp[n]) return mp[n]; long long ans=1; for(long long l=2,r;l<=n;l=r+1) { r=n/(n/l); ans=((ans-(r-l+1)*S_mu(n/l)%mod)%mod+mod)%mod; } mp[n]=(ans+mod)%mod; return mp[n]; } long long solve(long long n) { long long ans=0; for(long long l=1,r;l<=n;l=r+1) { r=n/(n/l); ans=(ans+(r-l+1)*(n/l)%mod)%mod; } return (ans+mod)%mod; } int main() { init(); long long n; cin>>n; long long ans=0; for(long long l=1,r;l<=n;l=r+1) { r=n/(n/l); ans=(ans+((S_mu(r)-S_mu(l-1))*solve(n/l)%mod*solve(n/l)%mod)%mod)%mod; } printf("%lld\n",(ans+mod)%mod); return 0; }