「学习笔记」莫比乌斯反演
前置芝士：数论分块。（之后再说。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}\)
杜教筛
杜教筛代码真是好短啊。
一种时间复杂度为 \(O(n^{\frac{2}{3}})\) 的筛，可以求出一些积性函数的前缀和。
具体来说，求 \(S(n)=\sum\limits_{i=1}^{n} f(i)\)，\(f\) 是积性函数。
构造函数 \(h,g\) 满足 \(h=f*g\)。
那么：
\[\begin{aligned}\sum\limits_{i=1}^nh(i)&=\sum\limits_{i=1}^n\sum\limits_{d|i}f(i)g(\frac{n}{i})\\&=\sum\limits_{d=1}^ ng(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f(i)\\&=\sum\limits_{d=1}^ ng(d)S(\left\lfloor\frac{n}{d}\right\rfloor)\end{aligned}
\]
那么：
\[g(1)S(n)=\sum\limits_{i=1}^nh(i)-\sum\limits_{d=2}^ ng(d)S(\left\lfloor\frac{n}{d}\right\rfloor)
\]
（这里就是把 \(d=1\) 的给提出来了）
\[S(n)=\frac{\sum\limits_{i=1}^nh(i)-\sum\limits_{d=2}^ ng(d)S(\left\lfloor\frac{n}{d}\right\rfloor)}{g(1)}
\]
对于 \(\mu\) 函数，知道 \(\mu * I = \epsilon\)，而 \(\epsilon\) 和 \(\mu\) 的前缀和都很好求，所以可以用杜教筛。
而对于 \(\varphi\) 函数，知道 \(\varphi * I = id\)，也可以用杜教筛。
P4213 【模板】杜教筛（Sum）
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
int T,n;
int phi[MAXN],mu[MAXN],sumphi[MAXN],summu[MAXN];
bool vis[MAXN];
vector<int> prm;
map<int,int> mphi,mmu;
void sieve(){
mu[1]=phi[1]=sumphi[1]=summu[1]=1;
for(int i=2;i<=1e7;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;mu[i]=-1;
}
sumphi[i]=sumphi[i-1]+phi[i];
summu[i]=summu[i-1]+mu[i];
for(auto j:prm){
if(i*j>1e7) break;
vis[i*j]=true;
if(i%j==0){
phi[i*j]=phi[i]*j;
break;
}
phi[i*j]=phi[i]*phi[j];
mu[i*j]=-mu[i];
}
}
}
int getmu(int x){
if(x<=1e7) return summu[x];
if(mmu[x]) return mmu[x];
int nem=1,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
nem-=(r-l+1)*getmu(x/l);
}
return mmu[x]=nem;
}
int getphi(int x){
if(x<=1e7) return sumphi[x];
if(mphi[x]) return mphi[x];
int nem=x*(x+1)/2,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
nem-=(r-l+1)*getphi(x/l);
}
return mphi[x]=nem;
}
signed main(){
cin>>T;
sieve();
while(T--){
cin>>n;
cout<<getphi(n)<<" "<<getmu(n)<<"\n";
}
return 0;
}
莫比乌斯反演
P2398 GCD SUM
（GCDMAT - GCD OF MATRIX 双倍经验）
求：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n \gcd(i,j)
\]
运用欧拉反演，可得：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n \sum\limits_{d|\gcd(i,j)} \varphi(d)
\]
转换，得：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n \sum\limits_{d=1} ^n\varphi(d)[d|i][d|j]
\]
\[\sum\limits_{d=1} ^n\varphi(d)\sum\limits_{i=1}^n[d|i]\sum\limits_{j=1}^n [d|j]
\]
\[\sum\limits_{d=1}^n\varphi(d)\lfloor\frac{n}{d}\rfloor^2
\]
线性欧拉筛即可。（可以用数论分块，也可以不用，因为 \(n\) 范围太小了）
不用数论分块：
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;
}
for(int j=0;j<prm.size();j++){
if(i*prm[j]<=n){
int nem=i*prm[j];
vis[nem]=true;
if(i%prm[j]==0){
phi[nem]=phi[i]*prm[j];
break;
}else{
phi[nem]=phi[i]*phi[prm[j]];
}
}else break;
}
}
}
signed main(){
cin>>n;
prime();
for(int i=1;i<=n;i++){
ans+=phi[i]*(n/i)*(n/i);
}
cout<<ans;
return 0;
}
用数论分块：
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 100010
#define int long long
using namespace std;
int ans;
int n,sum[MAXN],phi[MAXN];
vector<int> prm;
bool vis[MAXN];
void prime(){
sum[1]=phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;
}
for(int j=0;j<prm.size();j++){
if(i*prm[j]<=n){
int nem=i*prm[j];
vis[nem]=true;
if(i%prm[j]==0){
phi[nem]=phi[i]*prm[j];
break;
}else{
phi[nem]=phi[i]*phi[prm[j]];
}
}else break;
}
sum[i]=sum[i-1]+phi[i];
}
}
signed main(){
cin>>n;
prime();
int r;
for(int l=1;l<=n;l=r+1){
int nem=n/l;
r=n/nem;
ans+=(sum[r]-sum[l-1])*(n/l)*(n/l);
}
cout<<ans;
return 0;
}
P2303 [SDOI2012] Longge 的问题
\[\sum\limits_{i=1}^n \gcd(i, n)
\]
\[\sum\limits_{k|n}k\sum\limits_{i=1}^n [\gcd(i, n)=k]
\]
\[\sum\limits_{k|n}k\sum\limits_{i=1}^n [\gcd(\frac{i}{k}, \frac{n}{k})=1]
\]
\[\sum\limits_{k|n}k\sum\limits_{i=1}^{\frac{n}{k}} [\gcd(i, \frac{n}{k})=1]
\]
\[\sum\limits_{k|n}k \times \varphi(\frac{n}{k})
\]
单次求欧拉函数即可。
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define int long long
using namespace std;
int n,ans;
int phi(int x){
int num=x;
for(int i=2;i*i<=x;i++){
if(x%i==0) num=num/i*(i-1);
while(x%i==0) x/=i;
}
if(x>1) num=num/x*(x-1);
return num;
}
signed main(){
cin>>n;
for(int i=1;i*i<=n;i++){
if(i*i==n) ans+=i*phi(n/i);
else if(n%i==0) ans+=i*phi(n/i)+(n/i)*phi(i);
}
cout<<ans;
return 0;
}
P1390 公约数的和
UVA11417 GCD / UVA11424 GCD - Extreme (I) /P1390 公约数的和/ SP3871 GCDEX - GCD Extreme /UVA11426 拿行李（极限版） GCD - Extreme (II) （5 倍经验，按数据范围从小到大排序）
求：
\[\sum\limits_{i=1}^n\sum\limits_{j=i+1}^n\gcd(i,j)
\]
转化为：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^{i-1}\gcd(i,j)
\]
后面的直接套用上一道题的式子，可以提前线性筛欧拉函数，然后提前把答案都求出来，这样就可以把前两倍经验给切了。QWQ
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 200010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
vector<int> prm;
void sieve(){
phi[1]=1;
for(int i=2;i<=2e5+1;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;
}
for(int j=0;j<prm.size() and i*prm[j]<=2e5+1;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0){
phi[i*prm[j]]=phi[i]*prm[j];
break;
}
phi[i*prm[j]]=phi[i]*phi[prm[j]];
}
}
}
long long ans[MAXN],num;
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
sieve();
for(int i=1;i<=2e5+1;i++){
for(int j=1;j*j<=i;j++){
if(j*j==i) num+=j*phi[j];
else if(i%j==0) num+=1ll*j*phi[i/j]+(i/j)*phi[j];
}
num-=i;
ans[i]=num;
}
while(true){
cin>>n;if(n==0) break;
cout<<ans[n]<<"\n";
}
return 0;
}
时间复杂度预处理 \(O(n\sqrt{n})\)，查询 \(O(1)\)。
考虑优化。
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^{i-1}\gcd(i,j)
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^{i}\gcd(i,j)-i
\]
\[\sum\limits_{i=1}^n\sum\limits_{k|i}k \times \varphi(\frac{i}{k})-\frac{(n+1)\times n}{2}
\]
看左面的式子。
\[\sum\limits_{i=1}^n\sum\limits_{k=1}^nk \times \varphi(\frac{i}{k})[k|i]
\]
\[\sum\limits_{k=1}^nk\sum\limits_{i=1}^n\varphi(\frac{i}{k})[k|i]
\]
\[\sum\limits_{k=1}^nk\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\varphi(i)
\]
还是预处理线性筛欧拉函数，查询直接数论分块即可。
在代码里，为防止爆 long long，我把要减的 \(\frac{(n+1)\times n}{2}\) 直接在数论分块中算了。QWQ
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 4000010
using namespace std;
int n;
bool vis[MAXN];
int phi[MAXN];
long long num[MAXN];
vector<int> prm;
void sieve(){
phi[1]=1;num[1]=1;
for(int i=2;i<=4e6;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;
}
num[i]=num[i-1]+phi[i];
for(int j=0;j<prm.size() and i*prm[j]<=4e6;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0){
phi[i*prm[j]]=phi[i]*prm[j];
break;
}
phi[i*prm[j]]=phi[i]*phi[prm[j]];
}
}
}
long long ans;
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
sieve();
while(true){
ans=0;
cin>>n;
if(n==0) break;
int r;
for(int l=1;l<=n;l=r+1){
r=n/(n/l);
ans+=(1ll*(r+1)*r/2-1ll*(l-1)*l/2)*(1ll*num[n/l]-1);
}
cout<<ans<<"\n";
}
return 0;
}
说句闲话，这篇题解的式子（是我上面推的式子的方法）有点问题。QWQ
膜运算
这是一个裸题，因为显然可以用欧几里得定理搞成原题。
而我才不会告诉你我想的有亿点复杂。
推式子。
P3455 [POI2007] ZAP-Queries
（P4450 双亲数 双倍经验）
求：
\[\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;
}
P2522 [HAOI2011] Problem b
就是用容斥原理搞一下，之后就和上面的一样了。
点击查看代码
#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
下面正式推式子。
求：
\[\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;
}
P3327 [SDOI2015] 约数个数和
求：
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m d(ij)
\]
首先，要知道 \(d\) 函数的一个性质：
\[d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1]
\]
所以原式就变成了：
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1]
\]
\[\sum\limits_{x=1}^n\left\lfloor\frac{n}{x}\right\rfloor \sum\limits_{y=1}^m \left\lfloor\frac{m}{y}\right\rfloor[\gcd(x,y)=1]
\]
\[\sum\limits_{x=1}^n\left\lfloor\frac{n}{x}\right\rfloor \sum\limits_{y=1}^m \left\lfloor\frac{m}{y}\right\rfloor\sum\limits_{d=1}^{\min(n,m)}\mu(d)[d|x][d|y]
\]
\[\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{x=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{xd}\right\rfloor \sum\limits_{y=1}^{\left\lfloor\frac{m}{d}\right\rfloor} \left\lfloor\frac{m}{yd}\right\rfloor
\]
设 \(p=\left\lfloor\frac{n}{d}\right\rfloor\)，\(q=\left\lfloor\frac{m}{d}\right\rfloor\)，则:
\[\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{x=1}^p{\left\lfloor\frac{p}{x}\right\rfloor} \sum\limits_{y=1}^q \left\lfloor\frac{q}{y}\right\rfloor
\]
可以发现 \(\sum\limits_{x=1}^p{\left\lfloor\frac{p}{x}\right\rfloor}\) 和 \(\sum\limits_{y=1}^q \left\lfloor\frac{q}{y}\right\rfloor\) 可以预处理时存在数组里，我们设 \(f(x)=\sum\limits_{i=1}^x{\left\lfloor\frac{x}{i}\right\rfloor}\)。
\[\sum\limits_{d=1}^{\min(n,m)}\mu(d)f(p) f(q)
\]
\[\sum\limits_{d=1}^{\min(n,m)}\mu(d) f(\left\lfloor\frac{n}{d}\right\rfloor) f(\left\lfloor\frac{m}{d}\right\rfloor)
\]
于是查询时数论分块即可。（相当于数论分块套数论分块 QWQ）
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
#define int long long
using namespace std;
int T,n,m,ans;
bool vis[MAXN];
int mu[MAXN],sum[MAXN],num[MAXN];
vector<int> prm;
void sieve(){
mu[1]=1;sum[1]=1;
for(int i=2;i<=5e4;i++){
if(!vis[i]){
prm.push_back(i);
mu[i]=-1;
}
for(int j=0;j<prm.size() and prm[j]*i<=5e4;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0) break;
mu[i*prm[j]]=-mu[i];
}
sum[i]=sum[i-1]+mu[i];
}
for(int i=1;i<=5e4;i++){
int r,nem=0;
for(int l=1;l<=i;l=r+1){
r=i/(i/l);
nem+=(r-l+1)*(i/l);
}
num[i]=nem;
}
}
void calc(){
int r;
for(int l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(sum[r]-sum[l-1])*num[n/l]*num[m/l];
}
}
signed main(){
cin>>T;
sieve();
while(T--){
ans=0;
cin>>n>>m;
if(n>m) swap(n,m);
calc();
cout<<ans<<"\n";
}
return 0;
}
P3768 简单的数学题
求：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}ij\gcd(i,j)
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}ij\sum\limits_{d=1}^n\varphi(d)[d|i][d|j]
\]
\[\sum\limits_{d=1}^n\varphi(d)\sum\limits_{i=1}^ni[d|i]\sum\limits_{j=1}^{n}j[d|j]
\]
\[\sum\limits_{d=1}^n\varphi(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}i\times d\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}j\times d
\]
为方便，设 \(m=\left\lfloor\frac{n}{d}\right\rfloor\)。
\[\sum\limits_{d=1}^n\varphi(d)d^2\times\frac{m^2(m+1)^2}{4}
\]
\(\varphi(n)n^2\) 这个东西要用杜教筛。
设 \(f=\varphi\times id \times id\)，\(g=id \times id\)，则：
\[\begin{aligned}(h*g)(n)&=\sum\limits_{d|n}f(d)g(\frac{n}{d})\\&=\sum\limits_{d|n}\varphi(d)\times d^2\times (\frac{n}{d})^2\\&=n^2\times \sum\limits_{d|n}\varphi(d)\\&=n^3\end{aligned}
\]
所以用杜教筛后数论分块即可。
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010
#define int long long
using namespace std;
int p,n,ans;
int inv4,inv6;
int phi[MAXN],sumphi[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
phi[1]=sumphi[1]=1;
for(int i=2;i<=5e6;i++){
if(!vis[i]){
prm.push_back(i);
phi[i]=i-1;
}
sumphi[i]=(sumphi[i-1]+phi[i]*i%p*i%p)%p;
for(auto j:prm){
if(i*j>5e6) break;
vis[i*j]=true;
if(i%j==0){
phi[i*j]=phi[i]*j;
break;
}
phi[i*j]=phi[i]*phi[j];
}
}
}
int power(int x,int y){
int nem=1;
while(x){
if(x&1) nem=(nem*y)%p;
y=(y*y)%p;x>>=1;
}
return nem;
}
int pow2(int x){
x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;
}
int pow3(int x){
x%=p;return x*x%p*(x+1)%p*(x+1)%p*inv4%p;
}
map<int,int> mphi;
int getphi(int x){
if(x<=5e6) return sumphi[x];
if(mphi[x]) return mphi[x];
int nem=pow3(x),r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
nem=(nem-((pow2(r)-pow2(l-1)+p)%p*getphi(x/l)%p)%p+p)%p;
}
return mphi[x]=(nem%p+p)%p;
}
signed main(){
cin>>p>>n;
sieve();
inv4=power(p-2,4);
inv6=power(p-2,6);
int r;
for(int l=1;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+((getphi(r)-getphi(l-1))%p+p)%p*pow3(n/l))%p;
}
cout<<(ans%p+p)%p;
return 0;
}
P1891 疯狂 LCM 和 LCMSUM - LCM Sum
\[\sum\limits_{i=1}^n \operatorname{lcm}(n,i)
\]
\[\sum\limits_{i=1}^n\sum\limits_{d|n}\frac{n\times i}{d}[\gcd(n,i)=d]
\]
\[n\sum\limits_{d|n}\sum\limits_{i=1}^{\frac{n}{d}}i[\gcd(i,\frac{n}{d})=1]
\]
有一个式子： \(\sum\limits_{i=1}^ni[\gcd(n,i)=1]=\frac{\varphi(n)n}{2}\)。
\[n\sum\limits_{d|n}\frac{\varphi(d)d}{2}
\]
提前求出 \(\sum\limits_{d|n}\frac{\varphi(d)d}{2}\) 的值即可。（时间复杂度 \(O(n\ln n)\)）
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 1000010
#define ll long long
using namespace std;
int t,n;
int phi[MAXN];
ll ans[MAXN];
vector<int> prm;bool vis[MAXN];
void sieve(){
phi[1]=1;
for(int i=2;i<=1e6;i++){
if(!vis[i]) prm.push_back(i),phi[i]=i-1;
for(auto j:prm){
if(i*j>1e6) break;
vis[i*j]=true;
if(i%j==0){
phi[i*j]=phi[i]*j;
break;
}
phi[i*j]=phi[i]*phi[j];
}
}
for(int d=1;d<=1e6;d++){
for(int i=d;i<=1e6;i+=d){
if(d==1){ans[i]++;continue;}
ans[i]+=1ll*phi[d]*d/2;
}
}
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
sieve();
cin>>t;
while(t--){
cin>>n;
cout<<1ll*n*ans[n]<<"\n";
}
return 0;
}
P1829 [国家集训队] Crash的数字表格 / JZPTAB
设 \(n<m\)。求：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\operatorname{lcm}(i,j)
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d=1}^n\frac{i\times j}{d}[\gcd(i,j)=d]
\]
为方便，设 \(x=\left\lfloor\frac{n}{d}\right\rfloor\)，\(y=\left\lfloor\frac{m}{d}\right\rfloor\)。
\[\sum\limits_{d=1}^n\sum\limits_{i=1}^x\sum\limits_{j=1}^y d\times i\times j[\gcd(i,j)=1]
\]
\[\sum\limits_{d=1}^nd\sum\limits_{i=1}^x\sum\limits_{j=1}^y i\times j\sum\limits_{k|i,k|j}\mu(k)
\]
\[\sum\limits_{d=1}^nd\sum\limits_{k=1}^x\mu(k)\times k\times k\sum\limits_{i=1}^{\left\lfloor\frac{x}{d}\right\rfloor}i\sum\limits_{j=1}^{\left\lfloor\frac{y}{d}\right\rfloor}j
\]
设 \(\operatorname{sum}(x)=\sum\limits_{i=1}^x i\)，\(f(x,y)=\sum\limits_{k=1}^x\mu(k)\times k\times k\times \operatorname{sum}(\left\lfloor\frac{x}{d}\right\rfloor)\times \operatorname{sum}(\left\lfloor\frac{y}{d}\right\rfloor)\)。
\[\sum\limits_{d=1}^nd\times f(x,y)
\]
数论分块套数论分块即可。
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
const int mod=20101009;
int t,n,m,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<=1e7;i++){
if(!vis[i]){
prm.push_back(i);
mu[i]=-1;
}
num[i]=(num[i-1]+mu[i]*i%mod*i%mod)%mod;
for(int j=0;j<prm.size() and i*prm[j]<=1e7;j++){
vis[i*prm[j]]=true;
if(i%prm[j]==0) break;
mu[i*prm[j]]=-mu[i];
}
}
}
int getnum(int l,int r){return ((r+1)*r/2%mod-(l-1)*l/2%mod+mod)%mod;}
int f(int x,int y){
int r,nem=0;
for(int l=1;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
nem=(nem+(num[r]-num[l-1])%mod*getnum(1,x/l)%mod*getnum(1,y/l))%mod;
}
return nem%mod;
}
void calc(){
int r;
for(int l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+getnum(l,r)*f(n/l,m/l))%mod;
}
}
signed main(){
sieve();
cin>>n>>m;
if(n>m) swap(n,m);
calc();
cout<<(ans%mod+mod)%mod;
return 0;
}
P4917 天守阁的地板 和 P5221 Product
Product 凉心出题人又卡时间又卡空间
P3172 [CQOI2015] 选数
\[\sum\limits_{a_1=L}^R\sum\limits_{a_2=L}^R...\sum\limits_{a_n=L}^R[\gcd\limits_{i=1}^na_i=k]
\]
设 \(l=\left\lceil\frac{L}{k}\right\rceil\)，\(r=\left\lfloor\frac{R}{k}\right\rfloor\)。
\[\sum\limits_{a_1=l}^r\sum\limits_{a_2=l}^r...\sum\limits_{a_n=l}^r[\gcd\limits_{i=1}^na_i=1]
\]
\[\sum\limits_{d=1}^r\mu(d)\sum\limits_{a_1=l}^r[d|a_1]\sum\limits_{a_2=l}^r[d|a_2]...\sum\limits_{a_n=l}^r[d|a_n]
\]
\[\sum\limits_{d=1}^r\mu(d)(\left\lfloor\frac{r}{d}\right\rfloor-\left\lfloor\frac{l-1}{d}\right\rfloor)^n
\]
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
const int mod=1e9+7;
int n,k,L,R,ans;
int mu[MAXN],musum[MAXN];bool vis[MAXN];
vector<int> prm;
map<int,int> mmu;
void sieve(){
mu[1]=musum[1]=1;
for(int i=2;i<=5e6;i++){
if(!vis[i]) prm.push_back(i),mu[i]=-1;
musum[i]=musum[i-1]+mu[i];
for(auto j:prm){
if(i*j>5e6) break;
vis[i*j]=true;
if(i%j==0) break;
mu[i*j]=-mu[i];
}
}
}
int getmu(int x){
if(x<=5e6) return musum[x];
if(mmu.count(x)) return mmu[x];
int nem=1,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
nem-=getmu(x/l)*(r-l+1);
}
return mmu[x]=nem;
}
int power(int x,int y){
int nem=1;
while(y){
if(y&1) nem=nem*x%mod;
x=x*x%mod;y>>=1;
}
return nem;
}
void calc(){
int r;L--;
for(int l=1;l<=R;l=r+1){
if(!(L/l)) r=R/(R/l);
else r=min(R/(R/l),L/(L/l));
ans=(ans+power(R/l-L/l,n)*(getmu(r)-getmu(l-1))%mod)%mod;
}
}
signed main(){
sieve();
cin>>n>>k>>L>>R;
L=L/k+(L%k?1:0);R=R/k;
calc();
cout<<(ans%mod+mod)%mod;
return 0;
}
P6055 [RC-02] GCD
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{j}\right\rfloor}\sum\limits_{q=1}^{\left\lfloor\frac{n}{j}\right\rfloor} [\gcd(i,j)=1][\gcd(p,q)=1]
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{p=1}^j\sum\limits_{q=1}^j [\gcd(i,j)=1][\gcd(p,q)=j]
\]
\[\sum\limits_{i=1}^n\sum\limits_{p=1}^n\sum\limits_{q=1}^n [\gcd(i,p,q)=1]
\]
\[\sum\limits_{d=1}^n\mu(d)\left\lfloor\frac{n}{d}\right\rfloor^3
\]
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010
#define int long long
using namespace std;
const int mod=998244353;
int n;
int mu[MAXN],musum[MAXN];bool vis[MAXN];
vector<int> prm;
map<int,int> mmu;
void sieve(){
mu[1]=musum[1]=1;
for(int i=2;i<=5e6;i++){
if(!vis[i]) prm.push_back(i),mu[i]=-1;
musum[i]=musum[i-1]+mu[i];
for(auto j:prm){
if(i*j>5e6) break;
vis[i*j]=true;
if(i%j==0) break;
mu[i*j]=-mu[i];
}
}
}
int getmu(int x){
if(x<=5e6) return musum[x];
if(mmu.count(x)) return mmu[x];
int nem=1,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
nem-=(r-l+1)*getmu(x/l);
}
return mmu[x]=nem;
}
int ans;
signed main(){
sieve();
cin>>n;
int r;
for(int l=1;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+(getmu(r)-getmu(l-1))%mod*(n/l)%mod*(n/l)%mod*(n/l)%mod)%mod;
}
cout<<(ans%mod+mod)%mod;
return 0;
}
P4449 于神之怒加强版
假设 \(n<m\)。
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j)^k
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d=1}^nd^k[\gcd(i,j)=d]
\]
\[\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]
\]
\[\sum\limits_{d=1}^nd^k\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[p|i]\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[p|j]
\]
\[\sum\limits_{d=1}^nd^k\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\left\lfloor\frac{n}{dp}\right\rfloor\left\lfloor\frac{m}{dp}\right\rfloor
\]
设 \(T=dp\)。
\[\sum\limits_{T=1}^n\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{d|T}d^k\mu(\frac{T}{d})
\]
不难发现 \(\sum\limits_{d|T}d^k\mu(\frac{T}{d})\) 是个积性函数，所以可以像求欧拉函数一样用线性筛求出。
具体来说，设 \(f(T)=\sum\limits_{d|T}d^k\mu(\frac{T}{d})\)，则：
\[f(ab)=\begin{cases}f(a)\times f(b)\ \gcd(a,b)=1\\ f(a)\times b^k\ \gcd(a,b)\ne 1,b \in \mathbb{P}\end{cases}
\]
然后知道 \(f(1)=1\)，\(f(p)=p^k-1\)（\(p\in \mathbb{P}\)），所以就可以提前线性筛出来了。
点击查看代码
#include<bits/stdc++.h>
#define XD 114514
#define MAXN 5000010
#define int long long
using namespace std;
const int mod=1e9+7;
int t,n,m,k,ans;
vector<int> prm;bool vis[MAXN];
int sum[MAXN],num[MAXN];
int power(int x,int y){
int nem=1;
while(y){
if(y&1) nem=nem*x%mod;
x=x*x%mod;y>>=1;
}
return nem;
}
void sieve(){
num[1]=sum[1]=1;
for(int i=2;i<=5e6;i++){
if(!vis[i]){
prm.push_back(i);
num[i]=power(i,k)-1;
}
sum[i]=(sum[i-1]+num[i])%mod;
for(auto j:prm){
if(i*j>5e6) break;
vis[i*j]=true;
if(i%j==0){
num[i*j]=num[i]*(num[j]+1)%mod;
break;
}
num[i*j]=num[i]*num[j]%mod;
}
}
}
void calc(){
int r;ans=0;
for(int l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+(sum[r]-sum[l-1])*(n/l)%mod*(m/l)%mod+mod)%mod;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>t>>k;
sieve();
while(t--){
cin>>n>>m;
if(n>m) swap(n,m);
calc();
cout<<(ans+mod)%mod<<"\n";
}
return 0;
}
学习建议
建议看以下的博客，讲的非常好，反正我的也看不懂。QWQ
（我本人最推荐这两个。QWQ）
浙公网安备 33010602011771号