刷题笔记：莫比乌斯反演
入门紫题：yy的GCD
一句话题意：
求：
\[\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j) \in Prime]
\]
首先强制令\(n \le m\)。
考虑枚举\(gcd\)：
\[\sum_{k=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]*[k \in Prime]
\]
\[=
\]
\[\sum_{k=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor} [gcd(i,j)=1]*[k \in Prime]
\]
可以知道莫比乌斯函数有一个性质：
\[\sum_{d|n} \mu(d)=[n=1]
\]
所以原式等价于：
\[\sum_{k=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor}\sum_{d|gcd(i,j)}\mu(d) *[k \in Prime]
\]
\(d\)|\(gcd(i,j)\),那么，\(i,j\)都是\(d\)的倍数，考虑枚举\(d\)，则原式等价于：
\[\sum_{k=1}^{n}\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)*{\lfloor \frac{n}{kd} \rfloor}*{\lfloor \frac{m}{kd} \rfloor}*[k \in Prime]
\]
这里是考虑\(d\)会给范围内所有的倍数贡献，所以计算出倍数的个数即可将循环省去。
这里的复杂度会炸掉。
有一个套路，就是考虑变换枚举的对象。
另\(t=kd\)，枚举\(t\)，易知\(t \in [1,n]\)。
那么原式等价于：
\[\sum_{t=1}^{n} {\lfloor \frac{n}{t} \rfloor}*{\lfloor \frac{m}{t} \rfloor} \sum_{x|t,x \in Prime} \mu(\frac{t}{x})
\]
后者可以处理前缀和，前者可以除法分块。
总复杂度：\(O(q\sqrt{n})\).
Code
#include<bits/stdc++.h>
using namespace std;
namespace STD
{
#define rr register
using ll=long long ;
const int N=1e7+4;
int t,n,m;
inline void swap(int & x,int & y){x^=y,y^=x,x^=y;}
int read()
{
rr int x_read=0,y_read=1;
rr char c_read=getchar();
while(c_read<'0'||c_read>'9')
{
if(c_read=='-') y_read=-1;
c_read=getchar();
}
while(c_read<='9'&&c_read>='0')
{
x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
c_read=getchar();
}
return x_read*y_read;
}
int cnt;
int mu[N],f[N];
ll sum[N];
int prime[N/10];
bool is_not_prime[N];
void pre()
{
mu[1]=1;
for(int i=2;i<=N-4;i++)
{
if(!is_not_prime[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*prime[j]<=N-4;j++)
{
is_not_prime[i*prime[j]]=1;
if(i%prime[j]) mu[i*prime[j]]=-mu[i];
else
{
mu[i*prime[j]]=0;
break;
}
}
}
for(int i=1;i<=cnt;i++)
for(int j=1;j*prime[i]<=N-4;j++)
f[j*prime[i]]+=mu[j];
for(int i=1;i<=N-4;i++)
sum[i]=sum[i-1]+f[i];
}
};
using namespace STD;
int main()
{
pre();
t=read();
while(t--)
{
n=read(),m=read();
if(n>m) swap(n,m);
ll ans=0;
int l=1,r=0;
while(r<n)
{
r=min(n/(n/l),m/(m/l));
ans+=((ll)n/l)*((ll)m/l)*(ll)(sum[r]-sum[l-1]);
l=r+1;
}
printf("%lld\n",ans);
}
}
SDOI2014： 数表
离线询问+莫比乌斯反演+\(BIT\);
一句话题意：
求：
\[\sum_{i=1}^{n}\sum_{j=1}^{m} \sum_{d|gcd(i,j)} d
\]
强制令\(n \le m\)
首先记\(f(x)\)表示\(x\)所有的因数的和。
那么原式就是：
\[\sum_{i=1}^{n}\sum_{j=1}^{m} f(gcd(i,j))
\]
\[=
\]
\[\sum_{k=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor } [gcd(i,j)=1]*f(k)*[f(k) \le a]
\]
\[=
\]
\[\sum_{k=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor } \sum_{d|gcd(i,j)} \mu(d) * f(k) * [f(k) \le a]
\]
和上提一样的套路可以化为：
\[\sum_{t=1}^{n} {\lfloor \frac{n}{t} \rfloor}*{\lfloor \frac{m}{t} \rfloor} \sum_{d|t} {\frac{t}{d}} f(d) * [f(d) \le a]
\]
可以考虑一件事就是更小的\(a\)的值是可以给更大的\(a\)用的，将询问离线，然后按\(a\)升序排序。
每次将新加入的值加给它所有比\(n_{max}\)小的倍数，用树状数组维护前缀和，然后除法分块即可。
由于每个数只会被加一次，每次枚举倍数，所以这一步复杂度为\(N*ln N*log N\)
每次回答的复杂度为\(\sqrt{N}*logN\)。
所以总复杂度为\(O(N*ln N*log N)\)
Code
#include<bits/stdc++.h>
using namespace std;
namespace STD
{
#define rr register
using ll=long long;
const ll mod=1ll<<31;
const int N=1e5;
const int Q=2e4+4;
int q,last;
ll ans[Q];
void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
struct req
{
int n,m,a,id;
req():n(0),m(0),a(0),id(0) {}
bool operator<(const req x) const {return a<x.a;}
}a[Q];
template<typename type>
inline type cmin(type x,type y){return x<y?x:y;}
int read()
{
rr int x_read=0,y_read=1;
rr char c_read=getchar();
while(c_read<'0'||c_read>'9')
{
if(c_read=='-') y_read=-1;
c_read=getchar();
}
while(c_read<='9'&&c_read>='0')
{
x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
c_read=getchar();
}
return x_read*y_read;
}
class BIT
{
private:
ll t[N+4];
inline int lowbit(int x){return x&(-x);}
public:
void insert(int,int);
ll query(int,int);
}t;
void BIT::insert(int pos,int val)
{
while(pos<=N)
{
t[pos]=t[pos]+val;
if(t[pos]>=mod) t[pos]-=mod;
pos+=lowbit(pos);
}
}
ll BIT::query(int l,int r)
{
ll sum=0;
while(r)
{
sum=sum+t[r];
if(sum>=mod) sum-=mod;
r-=lowbit(r);
}
while(l)
{
sum=sum-t[l];
if(sum<0) sum+=mod;
l-=lowbit(l);
}
return sum;
}
struct pack
{
ll f;
int num;
pack():f(0ll),num(0) {}
bool operator<(const pack x) const {return f<x.f;}
}b[N+4];
int cnt;
int prime[(N+4)/10],mu[N+4];
bool tag[N+4];
void pre()
{
mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
{
tag[i*prime[j]]=1;
if(i%prime[j]) mu[i*prime[j]]=-mu[i];
else
{
mu[i*prime[j]]=0;
break;
}
}
}
for(int i=1;i<=N;i++)
{
b[i].num=i;
for(int j=1;(ll)(j*i)<=(ll)N;j++)
{
b[i*j].f+=i;
if(b[i*j].f>=mod)
b[i*j].f-=mod;
}
}
sort(b+1,b+1+N);
}
int find(ll x)
{
int l=1,r=N;
while(l<r)
{
int mid=(l+r+1)>>1;
if(b[mid].f<=x) l=mid;
else r=mid-1;
}
return l;
}
};
using namespace STD;
int main()
{
q=read();
pre();
for(int i=1;i<=q;i++)
{
a[i].n=read(),a[i].m=read();
a[i].a=read();
a[i].id=i;
if(a[i].n>a[i].m) my_swap(a[i].n,a[i].m);
}
sort(a+1,a+1+q);
for(int i=1;i<=q;i++)
{
if(a[i].a<=0) continue;
int pos=find(a[i].a);
for(int j=last+1;j<=pos;j++)
for(int k=1;(ll)(k*b[j].num)<=(ll)N;k++)
t.insert(k*b[j].num,mu[k]*b[j].f);
int l=1,r=0;
while(r<a[i].n)
{
r=cmin(a[i].n/(a[i].n/l),a[i].m/(a[i].m/l));
ans[a[i].id]=(ans[a[i].id]+((a[i].n/l)*(a[i].m/l)%mod)*t.query(l-1,r)%mod)%mod;
l=r+1;
}
last=pos;
}
for(int i=1;i<=q;i++) printf("%lld\n",ans[i]);
}
DZY Loves Math
线性筛的拓展使用+莫比乌斯反演。
首先考虑如何快速地求出题中给出的\(f\)。
传统做法比较无脑，就是直接\(O(n \sqrt n)\)枚举所有约数，边除边加加
有一个\(O(n)\)做法，就是利用线性筛，记录两个数组\(f\)与\(f'\)。\(f\)就是要求的数组，\(f'\)就是自己最小的质因子的幂数。
转移就有：
f[iprime[j]]=max(f'[i]+1,f[i]),f'[iprime[j]]=f'[i]+1;
其中\(prime[j]\)是线性筛枚举的素数。
最后的公式就是：
\[\sum_{t=1}^{n}\sum_{d|t} \mu(d)*f[\frac{t}{d}]*{\lfloor \frac{n}{t} \rfloor}{\lfloor \frac{m}{t} \rfloor }
\]
\(m,n\)老规矩。
推导略，因为太水了。
Code
#include<bits/stdc++.h>
using namespace std;
namespace STD
{
#define rr register
using ll=long long;
const int N=1e7;
int t;
int f[N+4],f_[N+4];
ll g[N+4];
template<typename type>
inline type cmax(type x,type y){return x>y?x:y;}
template<typename type>
inline type cmin(type x,type y){return x<y?x:y;}
inline void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
inline int read()
{
rr int x_read=0,y_read=1;
rr char c_read=getchar();
while(c_read<'0'||c_read>'9')
{
if(c_read=='-') y_read=-1;
c_read=getchar();
}
while(c_read<='9'&&c_read>='0')
{
x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
c_read=getchar();
}
return x_read*y_read;
}
int cnt;
int prime[N+4],mu[N+4];
bool tag[N+4];
void pre()
{
mu[1]=1;
int i,j,en;
for(i=2;i<=N;i++)
{
if(!tag[i]) prime[++cnt]=i,mu[i]=-1,f[i]=f_[i]=1;
for(j=1;j<=cnt&&prime[j]*i<=N;j++)
{
tag[i*prime[j]]=1;
if(i%prime[j]) mu[i*prime[j]]=-mu[i],f[i*prime[j]]=f[i],f_[i*prime[j]]=1;
else
{
f[i*prime[j]]=cmax(f[i],f_[i]+1);
f_[i*prime[j]]=f_[i]+1;
break;
}
}
}
for(i=2;i<=N;i++)
{
en=cmin(N/i,i);
for(j=1;j<=en;j++)
g[i*j]+=mu[j]*f[i]+mu[i]*f[j];
if(i==en&&i*i<=N)
g[i*i]-=mu[i]*f[i];
g[i]+=g[i-1];
}
}
};
using namespace STD;
int main()
{
t=read();
pre();
int a,b;
while(t--)
{
a=read(),b=read();
if(a>b) my_swap(a,b);
int l=1,r=0;
ll ans=0ll;
while(r<a)
{
r=cmin(a/(a/l),b/(b/l));
ans+=1ll*(a/l)*(b/l)*(g[r]-g[l-1]);
l=r+1;
}
printf("%lld\n",ans);
}
}
约数个数和
首先要知道\(d(ij)\)的性质：
\[d(ij)=\sum_{x|i} \sum_{y|j} [gcd(x,y)=1]
\]
证明见这里。
Code
#include<bits/stdc++.h>
using namespace std;
namespace STD
{
#define rr register
typedef long long ll;
const int N=50000;
int t;
int read()
{
rr int x_read=0,y_read=1;
rr char c_read=getchar();
while(c_read<'0'||c_read>'9')
{
if(c_read=='-') y_read=-1;
c_read=getchar();
}
while(c_read<='9'&&c_read>='0')
{
x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
c_read=getchar();
}
return x_read*y_read;
}
int cnt;
int mu[N+4],prime[N+4];
bool tag[N+1];ll f[N+1];
void pre()
{
mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
{
tag[i*prime[j]]=1;
if(i%prime[j]) mu[prime[j]*i]=-mu[i];
else
{
mu[prime[j]*i]=0;
break;
}
}
mu[i]+=mu[i-1];
}
for(int i=1;i<=N;i++)
{
int l=1,r=0;
while(r<i)
{
r=i/(i/l);
f[i]+=(r-l+1)*(i/l);
l=r+1;
}
}
}
inline void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
template<typename type>
inline type cmin(type x,type y){return x<y?x:y;}
};
using namespace STD;
int main()
{
pre();
t=read();
while(t--)
{
int n=read(),m=read();
ll ans=0;
if(n>m) my_swap(m,n);
int l=1,r=0;
while(r<n)
{
r=cmin(m/(m/l),n/(n/l));
ans+=(mu[r]-mu[l-1])*f[m/l]*f[n/l];
l=r+1;
}
printf("%lld\n",ans);
}
}
数字表格
人生第一道黑题，虽然说是黑里发紫。
首先\(O(n)\)搞出\(fib\)。
然后颓式子：
\[\prod_{i=1}^{n} \prod_{j=1}^{m} f[gcd(i,j)
\]
\[=
\]
\[\prod_{k=1}^{n} \prod_{i=1}^{n} \prod_{j=1}^{m} f[k]^{[gcd(i,j)=k]}
\]
\[=
\]
\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} f[k]^{[gcd(i,j)=1]}
\]
\[=
\]
\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} f[k]^{\sum_{d|gcd(i,j)} \mu(d)}
\]
\[=
\]
\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} \prod_{d|gcd(i,j)} f[k]^{\mu(d)}
\]
\[=
\]
\[\prod_{k=1}^{n} \prod_{d=1}^{\lfloor \frac{n}{k} \rfloor} (f[k]^{\mu(d)})^{{\lfloor \frac{n}{kd} \rfloor}{\lfloor \frac{m}{kd} \rfloor}}
\]
\[=
\]
\[\prod_{t=1}^{n} \prod_{d|t} (f[d]^{\mu(\frac{t}{d})})^{{\lfloor \frac{n}{t} \rfloor }{\lfloor \frac{m}{t} \rfloor}}
\]
处理一下前缀积即可。
Code
#include<bits/stdc++.h>
using namespace std;
namespace STD
{
#define rr register
using ll=long long;
const int mod=1e9+7;
const int N=1e6;
int t;
int read()
{
rr int x_read=0,y_read=1;
rr char c_read=getchar();
while(c_read<'0'||c_read>'9')
{
if(c_read=='-') y_read=-1;
c_read=getchar();
}
while(c_read<='9'&&c_read>='0')
{
x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
c_read=getchar();
}
return x_read*y_read;
}
ll qpow(ll base,ll exp)
{
if(exp<0) return qpow(base,mod-2);
ll ret=1ll;
while(exp)
{
if(exp&1) ret=ret*base%mod;
base=base*base%mod;
exp>>=1;
}
return ret;
}
int cnt;
int prime[N+1],mu[N+1];
ll f[N+1],g[N+1];
bool tag[N+1];
void pre()
{
mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
{
tag[i*prime[j]]=1;
if(!(i%prime[j])) break;
mu[i*prime[j]]=-mu[i];
}
}
f[1]=f[2]=1;
for(int i=3;i<=N;i++)
f[i]=(f[i-1]+f[i-2])%mod;
for(int i=0;i<=N;i++) g[i]=1ll;
for(int i=1;i<=N;i++)
{
for(int j=1;j*i<=N;j++)
g[i*j]=g[i*j]*qpow(f[j],mu[i])%mod;
g[i]=g[i]*g[i-1]%mod;
}
}
inline void my_swap(ll & x,ll & y){x^=y,y^=x,x^=y;}
template<typename type>
inline type cmin(type x,type y){return x<y?x:y;}
};
using namespace STD;
int main()
{
pre();
t=read();
while(t--)
{
ll n=read(),m=read();
if(n>m) my_swap(m,n);
int l=1,r=0;
ll ans=1ll,temp;
while(r<n)
{
r=cmin(n/(n/l),m/(m/l));
temp=g[r]*qpow(g[l-1],mod-2)%mod;
ans=(ans*qpow(temp,(n/l)*(m/l)))%mod;
l=r+1;
}
printf("%lld\n",ans);
}
}
