[莫比乌斯反演][min_25筛]任凭风浪起，稳坐钓鱼台（续）
Description
\[\sum_{j=1}^n\sum_{k=1}^n\mu(j k)\pmod{2^{32}}
\]
其中\(n\le10^9.\)
Solution
首先是莫比乌斯反演经典套路。
\[\begin{align*}
\sum_{j=1}^n\sum_{k=1}^n\mu(j k)&=\sum_{j=1}^n\sum_{k=1}^n\mu(j)\mu(k)[\gcd(j,k)=1]\\
&=\sum_{j=1}^n\sum_{k=1}^n\mu(j)\mu(k)\sum_{d|j,d|k}\mu(d)\\
&=\sum_{d=1}^n\mu(d)\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(j d)\sum_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(k d)\\
&=\sum_{d=1}^n\mu(d)\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(j)[\gcd(j,d)=1]\sum_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(k)[\gcd(k,d)=1]
\end{align*}
\]
我们记
\[f(n,a)=\sum_{k=1}^n\mu(k)[\gcd(k,a)=1]
\]
那么就有
\[\sum_{j=1}^n\sum_{k=1}^n\mu(j k)=\sum_{d=1}^n\mu(d)f\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)^2
\]
哪怕不考虑求\(f\)的开销，这也是\(O(n)\)，必须进行优化。
观察到当\(d\)很大时，\(\left\lfloor\frac{n}{d}\right\rfloor\)很小，不妨设一个阈值\(B\)，和端点\(L=\left\lfloor\frac{n}{B+1}\right\rfloor\)
\[\begin{align*}
\sum_{j=1}^n\sum_{k=1}^n\mu(j k)&=\sum_{d=1}^B\mu(d)f\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)^2+\sum_{B+1}^n\mu(d)\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(j)[\gcd(j,d)=1]\sum_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(k)[\gcd(k,d)=1]\\
&=\sum_{d=1}^B\mu(d)f\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)^2+\sum_{j=1}^L\sum_{k=1}^L\mu(j)\mu(k)\sum_{i=B+1}^{\min{\{n/j,n/k\}}}\mu(i)[\gcd(i,j)=1][\gcd(i,k)=1]\\
&=\sum_{d=1}^B\mu(d)f\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)^2+2\sum_{j=1}^L\sum_{k=1}^{j-1}\mu(j)\mu(k)\sum_{i=B+1}^{\left\lfloor\frac{n}{j}\right\rfloor}\mu(i)[\gcd(i,jk)=1]\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j=1}^L\mu(j)^2\sum_{i=B+1}^{\left\lfloor\frac{n}{j}\right\rfloor}\mu(i)[\gcd(i,j)=1] \\
&=\sum_{d=1}^B\mu(d)f\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)^2+2\sum_{j=1}^L\sum_{k=1}^{j-1}\mu(j)\mu(k)\left(f\left(\left\lfloor\frac{n}{d}\right\rfloor,jk\right)-f(B,jk)\right)\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j=1}^L\mu(j)^2\left(f\left(\left\lfloor\frac{n}{d}\right\rfloor,j\right)-f(B,j)\right)
\end{align*}
\]
总时间复杂度(不考虑\(f\))为\(O(B+\left(\frac{n}{B}\right)^2)\)，显然当\(B=n^\frac{2}{3}\)时取到最小值，为\(O(n^\frac{2}{3}).\)
那么问题就变成如何快速地求出\(f(n,a).\)
\[\begin{align*}
f(n,a)&=\sum_{k=1}^n\mu(k)[\gcd(k,a)=1]\\
&=\sum_{k=1}^n\mu(k)\sum_{t|k,t|a}\mu(t)\\
&=\sum_{t|a}\mu(t)\sum_{k=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\mu(k t)\\
&=\sum_{t|a}\mu(t)^2\sum_{k=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\mu(k)[\gcd(k,t)=1]\\
&=\sum_{t|a}\mu(t)^2f\left(\left\lfloor\frac{n}{t}\right\rfloor,t\right)
\end{align*}
\]
用这个递推式就可以比较快速地计算\(f\)了。
对于边界：\(f(n,1)\)，就是\(\mu\)的前缀和，用min_25筛预处理出整除分块得出的所有\(2\sqrt{n}\)个位置的前缀和。
（原因：\(\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{a b}\right\rfloor\)）
同时，我们可以在\(O(A \log A)\)的时间内预处理出所有\(n a\le A\)的\(f(n,a)\)，以卡常。
至此，问题解决。总复杂度是亚线性的。
Code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e9+1,M=1e6+1;
int n,b,l,cnt,pri[M],frm[M];
bool cm[M];
uint ans,mu[M],mt[1<<8];
vector<uint> pre[M];
inline uint sqr(uint x){return x*x;}
void sieve(){
mu[1]=1;
for(int i=2;i<M;i++){
if(!cm[i]){pri[++cnt]=i;frm[i]=i;mu[i]=-1;}
for(int j=1;j<=cnt&&i*pri[j]<M;j++){
cm[i*pri[j]]=true;
frm[i*pri[j]]=pri[j];
if(i%pri[j])mu[i*pri[j]]=-mu[i];
else break;
}
}
}
namespace min_25{
int m,w[M],id1[M],id2[M];
uint g[M],s[M];
inline int id(int x){
if(x<M)return id1[x];
else return id2[n/x];
}
void work(){
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
g[m]=-w[m]+1;
if(w[m]<M)id1[w[m]]=m;
else id2[r]=m;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&(ll)pri[j]*pri[j]<=w[i];i++){
int k=id(w[i]/pri[j]);
g[i]-=g[k]+j-1;
}
}
for(int i=1;i<=m;i++)s[i]=g[i];
for(int j=cnt;j>=1;j--)for(int i=1;i<=m&&(ll)pri[j]*pri[j]<=w[i];i++)s[i]-=s[id(w[i]/pri[j])]+j;
}
inline uint smu(int x){return s[id(x)]+1;}
}
void init(){
for(int i=1;i<M;i++)pre[i].resize((M-1)/i+1);
pre[1][0]=1;
for(int i=1;i<M;i++)for(int j=1;j<=(M-1)/i;j++)pre[i][j]=(i<=j?pre[i][j-i]:pre[j][i]);
for(int i=1;i<M;i++)for(int j=0;j<=(M-1)/i;j++)pre[i][j]=pre[i][j-1]+mu[j]*pre[i][j];
}
uint f(int n,int a){
if((ll)a*n<=(M-1))return pre[a][n];
if(n==1||a==1)return min_25::smu(n);
vector<int>sum;
while(a!=1){
int p=frm[a];
sum.push_back(p);
while(a%p==0)a/=p;
}
int sz=sum.size();
mt[0]=1;
uint ans=min_25::smu(n);
for(int i=1;i<(1<<sz);i++){
int x=i&(-i);
mt[i]=mt[i^x]*sum[__builtin_ctz(i)];
ans+=f(n/mt[i],mt[i]);
}
return ans;
}
int main(){
sieve();
scanf("%d",&n);
min_25::work();
init();
b=n/(n/min(n,1000000)+1);
l=n/(b+1);
for(int i=1;i<=b;i++)if(mu[i])ans+=mu[i]*sqr(f(n/i,i));
for(int i=1;i<=l;i++)for(int j=1;j<i;j++)if(mu[i]&&mu[j])ans+=2*mu[i]*mu[j]*(f(n/i,i*j)-f(b,i*j));
for(int i=1;i<=l;i++)if(mu[i])ans+=sqr(mu[i])*(f(n/i,i)-f(b,i));
printf("%u\n",ans);
}
