开新坑了。QWQ
前置芝士：数论分块。（之后再说。QWQ）
积性函数
定义
一个数论函数 \(f(n)\) 满足 \(f(xy)=f(x) \times f(y)\)（\(\gcd(x,y)=1\)），则称 \(f(n)\) 是积性函数。
莫比乌斯函数：
\(\mu(n) = \begin{cases}1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\)
求：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k]
\]
经典把 \(\gcd(i,j)=k\) 转化为 \(\gcd(i,j)=1\)。
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(\frac{i}{k},\frac{j}{k})=1][i|k][j|k]
\]
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j)=1]
\]
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum\limits_{d|\gcd(i,j)}\mu(d)
\]
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum\limits_{d|i,d|j}\mu(d)
\]
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum\limits_{d=1}^{\min(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor)}\mu(d)[d|i][d|j]
\]
\[\sum\limits_{d=1}^{\min(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}[d|i]\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}[d|j]
\]
\[\sum\limits_{d=1}^{\min(\left\lfloor\frac{n}{k}\right\rfloor,\left\lfloor\frac{m}{k}\right\rfloor)}\mu(d)\left\lfloor\frac{\left\lfloor\frac{n}{k}\right\rfloor}{d}\right\rfloor\left\lfloor{\frac{\left\lfloor\frac{m}{k}\right\rfloor}{d}}\right\rfloor
\]
线性筛一遍莫比乌斯函数再用数论分块即可。
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,k,ans;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
mu[1]=1;num[1]=1;
for(int i=2;i<=5e4;i++){
if(!vis[i]){
prm.push_back(i);
mu[i]=-1;
}
num[i]=num[i-1]+mu[i];
for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0) break;
mu[i*prm[j]]=-mu[i];
}
}
}
void calc(){
int r;
for(int l=1;l<=a;l=r+1){
r=min(a/(a/l),b/(b/l));
ans+=(num[r]-num[l-1])*(a/l)*(b/l);
}
}
int main(){
cin>>t;
sieve();
while(t--){
ans=0;
cin>>n>>m>>k;
a=n/k;b=m/k;
if(a>b) swap(a,b);
calc();
cout<<ans<<"\n";
}
return 0;
}
就是用容斥原理搞一下，之后就和上面的一样了。
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,c,d,k;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
mu[1]=1;num[1]=1;
for(int i=2;i<=5e4;i++){
if(!vis[i]){
prm.push_back(i);
mu[i]=-1;
}
num[i]=num[i-1]+mu[i];
for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0) break;
mu[i*prm[j]]=-mu[i];
}
}
}
int calc(int x,int y){
n=x/k;m=y/k;
if(n>m) swap(n,m);
int r,ans=0;
for(int l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(num[r]-num[l-1])*(n/l)*(m/l);
}
return ans;
}
int main(){
sieve();
cin>>t;
while(t--){
cin>>a>>b>>c>>d>>k;
cout<<calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1)<<"\n";
}
return 0;
}
P2257 YY的GCD 和 PGCD - Primes in GCD Table
说句题外话，我是先写的这道题，再写的上面那两道题。QWQ
用图画手推的式子，很抽象，就当图个乐就行了。
下面正式推式子。
求：
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m [\gcd(i,j) \in \operatorname{prime}]
\]
经典把 \(\gcd(i,j)=k\) 转化为 \(\gcd(i,j)=1\)。
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{p \in \operatorname{prime},p|i,p|j}[\gcd(\frac{i}{p},\frac{j}{p})=1]
\]
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{p \in \operatorname{prime}}[\gcd(\frac{i}{p},\frac{j}{p})=1][p|i][p|j]
\]
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} [\gcd(i,j)=1]
\]
然后用莫比乌斯反演。
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} \sum\limits_{d|\gcd(i,j)}\mu(d)
\]
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} \sum\limits_{d|i,d|j}\mu(d)
\]
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} \sum\limits_{d=1}^{\min(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor)}\mu(d)[d|i][d|j]
\]
这里就把 \(n\) 当做 \(\min(n,m)\)。
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} \sum\limits_{d=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\mu(d)[d|i][d|j]
\]
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{d=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\mu(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor} [d|i]\sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor} [d|j]
\]
\[\sum\limits_{p \in \operatorname{prime}}\sum\limits_{d=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\mu(d)\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor
\]
我们设 \(T=pd\)，并将式子转换成枚举 \(T\)。
\[\sum\limits_{T=1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{p \in \operatorname{prime},p|T}\mu(\frac{T}{p})
\]
可以发现右面的式子可以提前求出。左面的式子数论分块即可。
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
using namespace std;
int t,n,m;
int mu[MAXN],f[MAXN],sum[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
mu[1]=1;
for(int i=2;i<=1e7;i++){
if(!vis[i]){
prm.push_back(i);
mu[i]=-1;
}
for(int j=0;j<prm.size();j++){
if(i*prm[j]>1e7) break;
vis[i*prm[j]]=true;
if(i%prm[j]==0) break;
mu[i*prm[j]]=-mu[i];
}
}
for(int i=0;i<prm.size();i++){
for(int j=1;j<=1e7;j++){
if(prm[i]*j>1e7) break;
f[prm[i]*j]+=mu[j];
}
}
for(int i=1;i<=1e7;i++) sum[i]=sum[i-1]+f[i];
}
long long ans;
void solve(){
int r=0;if(n>m) swap(n,m);
for(int l=1;l<=n;l=r+1){
int nem1=n/l,nem2=m/l;
r=min(n/nem1,m/nem2);
ans+=1ll*(sum[r]-sum[l-1])*(n/l)*(m/l);
}
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>t;
sieve();
while(t--){
ans=0;
cin>>n>>m;
solve();
cout<<ans<<"\n";
}
return 0;
}
学习建议
建议看以下的博客，讲的非常好，反正我的也看不懂。QWQ
（我本人最推荐这两个。QWQ）