【XSY2721】求和 杜教篩
題目描述
設\(n=\prod a_i^{p_i}\)，那麽定義\(f_d(n)=\prod{(-1)^{p_i}[p_i\leq d]}\)。特別的，\(f_1(n)=\mu(n)\)。
給你\(n,k\)，求
\[
\sum_{i=1}^n\sum_{j=1}^n\sum_{d=1}^kf_d(\gcd(i,j))
\]
\(n\leq {10}^{10},k\leq 40\)
題解
先做一些簡單的處理
\[
\begin{align}
ans&=\sum_{i=1}^n\sum_{j=1}^n\sum_{d=1}^kf_d(\gcd(i,j))\&=\sum_{i=1}^n\sum_{d=1}^kf_d(i)(2(\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\varphi(j))-1)
\end{align}
\]
後面\(\varphi\)用杜教篩可以在\(O(n^\frac{2}{3})\)內搞出來。
設\(\lambda(n)=f_\infty(n)=\prod{(-1)}^{p_i}\)
考慮容斥，有：
\[
f_d(n)=\lambda(n)\sum_{i^d|n}\mu(i)
\]
\[ \begin{align} F_d(n)&=\sum_{i=1}^nf_d(i)\&=\sum_{i=1}^n\lambda(i)\sum_{j^d|i}\mu(j)\&=\sum_{i=1}^n\mu(i)\sum_{j=1}^{\lfloor\frac{n}{i^{d+1}}\rfloor}\lambda(i^{d+1}j)\&=\sum_{i=1}^{\lfloor\sqrt[d+1]{n}\rfloor}\lambda^{d+1}(i)\mu(i)\Lambda(\lfloor\frac{n}{i^{d+1}}\rfloor) \end{align} \]
\(n\leq {10}^7\)的部分預處理，其他的每次枚舉。這部分每次枚舉是\(\sqrt{n}\)的。總的是\(O(n^\frac{2}{3})\)的。（和杜教篩的分析方法一樣。）
\[ \begin{align} \sum_{j|i}\lambda(j)&=[i是完全平方數]\\sum_{i=1}^n\sum_{j|i}\lambda(j)&=\sqrt{n}\\sqrt{n}=\sum_{i=1}^n\sum_{j}[j|i]\lambda(i)&=\sum_{\frac{i}{j}=1}^n\sum_{j=1}^{\lfloor\frac{n}{\frac{i}{j}}\rfloor}\lambda(j) =\sum_{i=1}^n\Lambda(\lfloor\frac{n}{i}\rfloor)\\Lambda(n)&=\sqrt {n}-\sum_{i=2}^n\Lambda(\lfloor\frac{n}{i}\rfloor) \end{align} \]
後面\(\Lambda(n)\)用杜教篩可以在\(O(n^\frac{2}{3})\)內搞出來
反正總的是\(O(n^\frac{2}{3})\)的就對了。。。
時間復雜度：\(O(n^\frac{2}{3})\)
代碼
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int getsqrt(ll n)
{
int l=1,r=1000000;
while(l<r)
{
int mid=(l+r+1)>>1;
if((ll)mid*mid>n)
r=mid-1;
else
l=mid;
}
return l;
}
ll n;
ll _sqrt;
namespace hashset
{
int getnum(ll x)
{
return n/x;
}
}
using hashset::getnum;
int miu[10000010];
int phi[10000010];
int c[10000010];
int cs[10000010];
const int maxn=10000000;
int b[10000010];
int pri[1000010];
int cnt;
int d[10000010];
int e[10000010];
int f[10000010];
int k;
int vis[10000010];
int qp[10000010];
int qc[10000010];
void init()
{
c[1]=phi[1]=miu[1]=f[1]=e[1]=1;
d[1]=f[1]=0;
int i,j;
for(i=2;i<=maxn;i++)
{
if(!b[i])
{
miu[i]=-1;
phi[i]=i-1;
c[i]=-1;
pri[++cnt]=i;
d[i]=e[i]=1;
f[i]=1;
}
for(j=1;j<=cnt&&i*pri[j]<=maxn;j++)
{
b[i*pri[j]]=1;
c[i*pri[j]]=-c[i];
f[i*pri[j]]=f[i]+1;
if(i%pri[j]==0)
{
miu[i*pri[j]]=0;
phi[i*pri[j]]=phi[i]*pri[j];
d[i*pri[j]]=d[i]+1;
e[i*pri[j]]=max(d[i*pri[j]],e[i]);
break;
}
d[i*pri[j]]=1;
e[i*pri[j]]=e[i];
miu[i*pri[j]]=-miu[i];
phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for(i=1;i<=maxn;i++)
{
if(e[i]>k)
f[i]=0;
else
f[i]=(f[i]&1?-1:1)*(k-e[i]+1);
f[i]+=f[i-1];
// miu[i]+=miu[i-1];
phi[i]+=phi[i-1];
cs[i]=cs[i-1]+c[i];
}
}
int getphi(ll n)
{
if(n<=maxn)
return phi[n];
int x=getnum(n);
if(vis[x]&1)
return qp[x];
vis[x]|=1;
ll i,j;
int s=n*(n+1)>>1;
for(i=2;i<=n;i=j+1)
{
j=n/(n/i);
s-=(j-i+1)*getphi(n/i);
}
qp[x]=s;
return s;
}
int getc(ll n)
{
if(n<=maxn)
return cs[n];
int x=getnum(n);
if(vis[x]&2)
return qc[x];
vis[x]|=2;
int s=getsqrt(n);
ll i,j;
for(i=2;i<=n;i=j+1)
{
j=n/(n/i);
s-=(j-i+1)*getc(n/i);
}
qc[x]=s;
return s;
}
ll pw[1000010];
int pw2[1000010];
int pw3[1000010];
int getfd(ll n)
{
if(n<=maxn)
return f[n];
int i,j;
for(i=1;(ll)i*i<=n;i++)
{
pw[i]=i;
pw2[i]=pw3[i]=c[i];
}
int m=i-1;
int s=0;
for(j=1;j<=k;j++)
{
for(i=1;i<=m;i++)
{
pw[i]*=i;
if(pw[i]>n)
break;
pw2[i]*=pw3[i];
s+=miu[i]*pw2[i]*getc(n/pw[i]);
}
}
return s;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("b.in","r",stdin);
freopen("b.out","w",stdout);
#endif
scanf("%lld%d",&n,&k);
_sqrt=getsqrt(n);
init();
int s=0;
ll i,j;
int now,last=0;
int ans=0;
for(i=1;i<=n;i=j+1)
{
j=n/(n/i);
now=getfd(j);
ans+=(now-last)*(2*getphi(n/i)-1);
last=now;
}
ans&=(1<<30)-1;
printf("%d\n",ans);
return 0;
}【XSY2721】求和 杜教篩