莫比乌斯反演
经典柿子
\[[gcd(i,j)=1]=\sum_{w\mid gcd(i,j)}\mu(w)=\sum_{w\mid i,w\mid j}\mu(w)
\]
POJ3904 Sky Code
description
给定\(n\) 个不超过\(10000\) 的正整数\(a_1,a_2,\cdots,a_n\) ，求
\[\sum_{1\le i<j<k<l\le n}[\gcd(a_i,a_j,a_k,a_l)=1]
\]
\(n\le 10000,T\le10000\)
solution
大力推柿子
\[\begin{aligned}
ans&=\sum_{1\le i<j<k<l\le n}[\gcd(a_i,a_j,a_k,a_l)=1]\\
&=\sum_{1\le i<j<k<l\le n}\sum_{w\mid a_i,w\mid a_j,w\mid a_k,w\mid a_l}\mu(w)\\
&=\sum_{w=1}^{10000}\mu(w)\sum_{1\le i<j<k<l\le n}[w\mid a_i][w\mid a_j][w\mid a_k][w\mid a_l]
\end{aligned}
\]
我们不妨假设有\(c(w)\) 个数是\(w\) 的倍数，那么后面那些东西就等于\(\dbinom{c(w)}{4}\) 。
复杂度\(\mathcal O(nT)\) 。
code
#include<cstdio>
#include<vector>
using namespace std;
const int N=1e4+5;
typedef long long ll;ll C[N];
bool flag[N];int mu[N];
vector<int>pr;
inline void pre()
{
int n=N-5;
flag[1]=1;mu[1]=1;
for(int i=2;i<=n;++i)
{
if(!flag[i]){pr.push_back(i);mu[i]=-1;}
for(int j=0;j<pr.size();++j)
{
int num=i*pr[j];if(num>n)break;
flag[num]=1;
if(i%pr[j])mu[num]=-mu[i];
else{mu[num]=0;break;}
}
}
for(int i=4;i<=n;++i)
C[i]=1ll*i*(i-1)*(i-2)*(i-3)/24;
}
int n,cnt[N];
int main()
{
pre();
while(scanf("%d",&n)==1)
{
int mx=0;
for(int i=1,d;i<=n;++i)
{
scanf("%d",&d);
for(int j=1;j*j<=d;++j)
if(d%j==0)
{
int w1=d/j,w2=j;
++cnt[w1];
if(w1^w2)++cnt[w2];
}
mx=max(mx,d);
}ll ans=0;
for(int i=1;i<=mx;++i)
{
if(!mu[i])continue;
ans+=mu[i]*C[cnt[i]];
}
printf("%lld\n",ans);
for(int i=1;i<=mx;++i)cnt[i]=0;
}
return 0;
}
可见格点
description
三维直角坐标系中三个坐标均非负且不超过\(n\) 的整点中有多少点可以从原点看到。
\(n\le 10^6\)
solution
原题相当于是求
\[\sum_{x\ge 0,y\ge0,z\ge0}[\gcd(x,y,z)=1]
\]
这是经典的莫反柿子。不过需要考虑\(x,y,z\) 中部分为零的情况，于是就有
\[\begin{aligned}
&=3+3\sum_{x>0,y>0}[\gcd(x,y)=1]+\sum_{x>0,y>0,z>0}[\gcd(x,y,z)=1]\\
&=3+3\sum_{x>0,y>0}\sum_{w\mid x,w\mid y}\mu(w)+\sum_{x>0,y>0,z>0}\sum_{w\mid x,w\mid y,w\mid z}\mu(w)\\
&=3+3\sum_{w=1}^n\mu(w)\lfloor\frac nw\rfloor^2+\sum_{w=1}^n\mu(w)\lfloor\frac nw\rfloor^3\\
&=3+\sum_{w=1}^n\mu(w)\lfloor\frac nw\rfloor^2(\lfloor\frac nw\rfloor+3)
\end{aligned}
\]
复杂度\(\mathcal O(nT)\) 。
code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1E+6+5,V=1e+6;
bool flag[N];vector<int>pr;
int mu[N];
inline void pre()
{
mu[1]=1,flag[1]=1;
for(int i=2;i<=V;++i)
{
if(!flag[i])pr.push_back(i),mu[i]=-1;
for(int c:pr)
{
int num=c*i;if(num>V)break;
flag[num]=1;
if(i%c)mu[num]=-mu[i];
else{mu[num]=0;break;}
}
}
}
int main()
{
pre();
int T;scanf("%d",&T);
while(T-->0)
{
int n;scanf("%d",&n);
ll ans=3;
for(int i=1;i<=n;++i)
{
if(!mu[i])continue;
int k=n/i;
ans+=1ll*mu[i]*k*k*(k+3);
}
printf("%lld\n",ans);
}
return 0;
}
数一数a x b
description
令
\[f(n)=\sum_{i=1}^n\sum_{j=1}^n[n\nmid ij]
\]
求
\[g(n)=\sum_{d\mid n}f(d)
\]
\(n\le 10^9,T\le 20000\)
solution
不整除显然是难以计算的，我们考虑简单容斥，即
\[f(n)=n^2-h(n)=n^2-\sum_{i=1}^n\sum_{j=1}^n[n\mid ij]
\]
因此就有
\[g(n)=\sum_{d\mid n}d^2-\sum_{d\mid n}\sum_{i=1}^d\sum_{j=1}^d[d\mid ij]
\]
分开考虑。我们假设\(n\) 的质因数分解为\(n=\prod_\limits{i=1}^kp_i^{a_i}\) ，那么就有
\[\sum_{d\mid n}d^2=\prod_{i=1}^k\sum_{j=0}^{a_i}(p_i^j)^2
\]
因为次数之和大概是\(\mathcal O(\log n)\) 级别的，因此这部分复杂度为\(\mathcal O(\log n)\) 。
对于后面部分，我们枚举\(w=\gcd(d,i)\) ，则有
\[\begin{aligned}
&=\sum_{d\mid n}\sum_{i=1}^d\sum_{j=1}^d[d\mid ij]\\
&=\sum_{w=1}^n\sum_{d\mid n,w\mid d}\sum_{w\mid i}^d\sum_{j=1}^d[\frac dw\mid \frac iw\cdot j][\gcd(\frac dw,\frac iw)=1]\\
&=\sum_{w\mid n}\sum_{d\mid \frac nw}\sum_{j=1}^{wd}[d\mid j]\sum_{i=1}^{d}[\gcd(d,i)=1]\\
&=\sum_{w\mid n}\sum_{d\mid \frac nw}\sum_{j=1}^{wd}[d\mid j]\varphi(d)\\
&=\sum_{w\mid n}\sum_{d\mid \frac nw}\sum_{j=1}^{w}\varphi(d)\\
&=\sum_{w\mid n}w\sum_{d\mid \frac nw}\varphi(d)\\
&=\sum_{w\mid n}w\frac nw\\
&=n\cdot \sigma_0(n)
\end{aligned}
\]
可以通过预处理质数来优化后续复杂度。至此，本题完结。
总复杂度为\(\mathcal O(\sqrt n+T\log n)\) 。
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
typedef unsigned long long ll;
bool flag[N];vector<int>pr;
inline void pre(int n)
{
flag[1]=1;
for(int i=2;i<=n;++i)
{
if(!flag[i])pr.push_back(i);
for(int j=0;j<pr.size();++j)
{
int num=i*pr[j];if(num>n)break;
flag[num]=1;
if(i%pr[j]==0)break;
}
}
}
int main()
{
pre(4e4);int T;scanf("%d",&T);
while(T-->0)
{
int n;scanf("%d",&n);ll g=1,h=n;
for(int i=0;i<pr.size();++i)
{
int p=pr[i];ll ret=0;
if(n<p)break;if(n%p)continue;
int k=0;while(n%p==0)n/=p,++k;
for(int j=0,cur=1;j<=k;++j,cur*=p)ret+=1ull*cur*cur;
g*=ret;h*=k+1;
}
if(n>1)h<<=1,g*=(1+1ull*n*n);
printf("%llu\n",g-h);
}
return 0;
}
墨菲斯
description
给定\(n,m,p\) ，求
\[\sum_{i=1}^n\sum_{j=1}^m[h(\gcd(i,j))\le p]
\]
其中\(h(x)\) 表示\(x\) 的唯一分解中质数个数。
solution
首先套路地推柿子，枚举\(\gcd\)
\[\begin{aligned}
ans&=\sum_{d=1}^n[h(d)\le p]\sum_{d\mid i}^n\sum_{d\mid j}^m[\gcd(\frac id,\frac jd)=1]\\
&=\sum_{d=1}^n[h(d)\le p]\sum_i\sum_j\sum_{w\mid i,w\mid j}\mu(w)\\
&=\sum_{d=1}^n[h(d)\le p]\sum_{w=1}^{\lfloor \frac nd\rfloor}\mu(w)\lfloor \frac n{dw}\rfloor\lfloor \frac m{dw}\rfloor\\
&=\sum_{T=1}^n\lfloor \frac nT\rfloor\lfloor \frac mT\rfloor\sum_{d\mid T}[h(d)\le p]\mu(\frac Td)
\end{aligned}
\]
似乎遇到了障碍。不过我们注意到\(2^{19}>5\times 10^5\) ，因此当\(p\ge 18\) 时答案就是\(nm\) 。我们只用考虑\(p<18\) 的情况。
令
\[f(p,d)=\sum_{d\mid T}[h(d)\le p]\mu(\frac Td)
\]
我们可以在\(\mathcal O(18n\log n)\) 的时间内求出所有的\(f(p,d)\) ，然后对其求前缀和以便在整除分块时可以快速计算区间和。
复杂度\(\mathcal O(18n\log n+q\sqrt n)\) 。
code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e5+5,MX=18;
bool flag[N];int mu[N],h[N];ll f[MX+5][N];
vector<int>pr;
inline void pre()
{
int n=N-5;
flag[1]=1;mu[1]=1;h[1]=0;
for(int i=2;i<=n;++i)
{
if(!flag[i])pr.push_back(i),mu[i]=-1,h[i]=1;
for(int j=0;j<pr.size();++j)
{
int num=i*pr[j];if(num>n)break;
flag[num]=1;h[num]=h[i]+1;
if(i%pr[j])mu[num]=-mu[i];
else{mu[num]=0;break;}
}
}
for(int p=0;p<=MX;++p)
{
for(int d=1;d<=n;++d)
{
if(h[d]>p)continue;
for(int t=d,c=1;t<=n;t+=d,++c)f[p][t]+=mu[c];
}
for(int i=1;i<=n;++i)f[p][i]+=f[p][i-1];
}
}
int main()
{
pre();int T;scanf("%d",&T);
while(T-->0)
{
int n,m,p;scanf("%d%d%d",&n,&m,&p);
if(p>MX){printf("%lld\n",1ll*n*m);continue;}
if(n>m)swap(n,m);ll ans=0;
for(int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*(f[p][r]-f[p][l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
NO PAIN NO GAIN
