莫比乌斯反演学习笔记
前言
莫比乌斯反演是数论中的重要内容
对于一些函数\(f(n)\)，如果很难直接求出它的值
而容易求出其倍数和或约数和\(g(n)\)，那么可以通过莫比乌斯反演简化运算，求得\(f(n)\)的值
前置知识：整除分块
对于类似于 \(\sum\limits_{i=1}^{n}\lfloor\frac{n}{i}\rfloor\) 的式子
我们可以很容易地用 \(O(n)\) 的效率求出它的值
但是在一些题目中，我们需要更优秀的时间效率，这时候就要用到整除分块
主要的思想是对于任意一个 \(i(i \leq n)\) 找到一个最大的 \(j\)
使得 \(i \leq j \leq n\) 并且 \(\lfloor\frac{n}{i}\rfloor=\lfloor\frac{n}{j}\rfloor\)
此时 \(j=\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor\)
因为对于 $\lfloor\frac{n}{i}\rfloor $，当 \(i \leq \sqrt{n}\) 时只有 \(\sqrt{n}\) 种，当 \(i>\sqrt{n}\)时，\(\lfloor\frac{n}{i}\rfloor< \sqrt{n}\)，也只有 \(\sqrt{n}\) 种取值，共计 \(2\sqrt{n}\) 种
所以该过程的时间复杂度为 \(O(\sqrt{n})\)
关于正确性的证明
首先显然 \(j \leq n\)
又有 \(j=\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor \geq \lfloor \frac{n}{\frac{n}{i}}\rfloor =i\)
设 \(k=\lfloor \frac{n}{i} \rfloor\),那么我们要证明的就是当 \(\lfloor \frac{n}{j} \rfloor=k\) 时，\(j\) 的最大值为 \(\lfloor \frac{n}{k} \rfloor\)
因为 \(j\) 为整数，所以 \(max(j)=\lfloor \frac{n}{k} \rfloor\)
练习题
分析
\(ans=\sum\limits_{i=1}^n(k\bmod i)=\sum\limits_{i=1}^nk-i\left\lfloor\frac{k}{i}\right\rfloor\)
代码
#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
int n,k,mmax;
long long ans=0;
int main(){
scanf("%d%d",&n,&k);
mmax=std::min(n,k);
for(int l=1,r;l<=mmax;l=r+1){
r=k/(k/l);
r=std::min(r,mmax);
ans+=1LL*(k/l)*(l+r)*(r-l+1)/2;
}
ans=1LL*n*k-ans;
printf("%lld\n",ans);
return 0;
}
拓展：二维整除分块
求 $ \sum\limits_{i=1}^{\min (n,m)}\left\lfloor\frac{n}{i} \right\rfloor\left\lfloor\frac{m}{i} \right\rfloor $的值
此时我们需要将代码中 \(r\) 的取值改为 \(min(m/(m/l),n/(n/l))\)
莫比乌斯函数
\(μ(d)\) 为莫比乌斯函数的函数名，是一个由容斥系数所构成的函数。
其定义为：
\(1\)、当 \(d=1\) 时， \(\mu(d)=1\) 。
\(2\)、当 \(d = \prod_{i=1}^k\ p_i\)，并且所有的 \(p_i\) 为不同的素数的时候， \(\mu(d)=(-1)^k\) 。
\(3\)、当 \(d\) 含有的任一质因子的幂大于等于 \(2\) 次，那么 \(\mu(d)=0\) 。
求法（线性筛）
const int maxn=4e5+5;
bool not_pri[maxn];
int t,mu[maxn],pri[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
}
性质一
对于任意正整数 \(n\)，\(\sum\limits_{d|n}\mu(d)=[n=1]\)
其中 \([n=1]\) 表示只有当 \(n=1\) 成立的时候返回值为 \(1\)，否则返回值为 \(1\)
在反演的时候，这一个式子经常会将 \(gcd(i,j)=1\) 这个条件代换掉
证明
当 \(n=1\)时，返回值显然为 \(1\)
当 \(n\neq1\)时，根据唯一分解定理，\(n\) 一定能划分为若干质因数的乘积
我们设这样的质因数一共有 \(k\) 种
如果任意一种质因数选择了大于一个，那么组成的数的 \(\mu\) 值一定为 \(0\)
所以每种质因子最多选择 \(1\) 个
所以选择 \(i\) 种质因子的方案数为 \(C_k^i\)，返回值为 \((-1)^i\)
贡献的总和为 \(\sum_{i=1}^k C_k^i(-1)^i\)
如果再加上一项 \(C_k^0(-1)^0\)
根据二项式定理，恰好是 \((1-1)^k=1\)
减去加上的一项，结果就是 \(0\)
性质二
对于任意正整数\(n\)，\(\sum\limits_{d|n}\frac{\mu(d)}{d}=\frac{\varphi(n)}{n}\)
这个性质阐述了莫比乌斯函数和欧拉函数之间的联系
证明
\(\varphi(n)=\sum\limits_{i=1}^n[\gcd(n,i)=1]=\sum\limits_{i=1}^{n} \sum\limits_{d|gcd(n,i)} \mu(d)=\sum\limits_{d|n} \mu(d) \times \frac{n}{d}\)
转换一下就变成了\(\sum\limits_{d|n}\frac{\mu(d)}{d}=\frac{\varphi(n)}{n}\)
莫比乌斯反演
形式一
\(F(n)=\sum\limits_{d|n}f(d)=>f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})\)
证明
\(f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})=\sum\limits_{d|n}\mu(d)\sum\limits_{k|\frac{n}{d}}f(k)=\sum\limits_{k|n}f(k)\sum\limits_{d|\frac{n}{k}}\mu(d)\)
根据莫比乌斯函数性质一，只有当 \(n=1\) 时，\(\sum\limits_{d|n} \mu(d)=1\)
所以只有当 \(n=k\) 时，求出的值才不为 \(0\)
所以右边等于 \(f(n)\)
形式二
\(F(n)=\sum\limits_{n|d}f(d)=>f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d)\)
证明
令\(k=\frac{d}{n}\)，
\(f(n)=\sum^{+\infty}\limits_{k=1}\mu(k)F(nk)=\sum^{+\infty}\limits_{k=1}\mu(k)\sum\limits_{nk|t}f(t)=\sum\limits_{n|t}f(t)\sum\limits_{k|\frac{t}{n}}\mu(k)\)
当且仅当 \(t=n\) 时，\(\sum\limits_{k|\frac{t}{n}}\mu(k)=1\)，其它情况都为 \(0\)
所以右边等于 \(f(n)\)
例题一、ZAP-Queries
分析
要求的式子是
\(\sum\limits_{i=1}^a \sum\limits_{j=1}^b[gcd(i,j)=k]\)
默认 \(a>b\)
把 \(a\) 和 \(b\) 的上界同时除以 \(k\)，可以得到
\(\sum\limits_{i=1}^{a/k} \sum\limits_{j=1}^{b/k}[gcd(i,j)=1]\)
把后面的部分用莫比乌斯函数带入，可以得到
\(\sum\limits_{i=1}^{a/k} \sum\limits_{j=1}^{b/k}\sum\limits_{d|gcd(i,j)}\mu(d)\)
把 \(d\) 提到前面
\(\sum\limits_{d=1}^{b/k}\mu(d)\sum\limits_{i=1}^{a/kd} \sum\limits_{j=1}^{b/kd} 1\)
也就是
\(\sum\limits_{d=1}^{b/k}\mu(d) \frac{a}{kd} \frac{b}{kd}\)
我们用整除分块处理一下，就可以做到单次询问 \(\sqrt{b}\) 的复杂度
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5;
bool not_pri[maxn];
int t,n,m,d,mu[maxn],pri[maxn],sum[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
for(rg int i=1;i<maxn;i++){
sum[i]=sum[i-1]+mu[i];
}
}
int main(){
xxs();
t=read();
while(t--){
n=read(),m=read(),d=read();
if(n>m) std::swap(n,m);
n/=d,m/=d;
rg int mmax=n;
rg long long ans=0;
for(rg int l=1,r;l<=mmax;l=r+1){
r=std::min(n/(n/l),m/(m/l));
if(r>mmax) r=mmax;
ans+=1LL*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
例题二、YY的GCD
分析
要求的式子是
\(\sum\limits_{i=1}^n \sum\limits_{j=1}^m[gcd(i,j)=k](k \in prime)\)
默认 \(n>m\)
和上一道题一样，我们把 \(k\) 提出来
\(\sum\limits_{k=1}^m\sum\limits_{i=1}^{n/k} \sum\limits_{j=1}^{m/k}[gcd(i,j)=1](k \in prime)\)
把后面的部分用莫比乌斯函数带入，可以得到
\(\sum\limits_{k=1}^m\sum\limits_{i=1}^{n/k} \sum\limits_{j=1}^{m/k}\sum\limits_{d|gcd(i,j)}\mu(d)(k \in prime)\)
把 \(d\) 扔到前面
\(\sum\limits_{k=1}^m\sum\limits_{d=1}^m \mu(d)\sum\limits_{i=1}^{n/kd}\sum\limits_{j=1}^{m/kd}(k \in prime) 1\)
其实就是
\(\sum\limits_{k=1}^m\sum\limits_{d=1}^m \mu(d)\frac{n}{kd}\frac{m}{kd}(k \in prime)\)
设\(kd=T\)，改为枚举 \(T\)，则有
\(\sum\limits_{T=1}^m \sum\limits_{k=1}^m\mu(T/k)\frac{n}{T}\frac{m}{T}(k \in prime)\)
把求和符号换一下位置
\(\sum\limits_{T=1}^m \frac{n}{T}\frac{m}{T}\sum\limits_{k=1}^m\mu(T/k)(k \in prime)\)
前半部分可以整除分块计算，后半部分可以用 \(Dirichlet\) 前缀和预处理
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
typedef long long ll;
const int maxn=1e7+5;
bool not_pri[maxn];
int n,m,mu[maxn],pri[maxn];
ll sum[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
for(rg int i=1;i<=pri[0];i++){
for(rg int j=1;1LL*j*pri[i]<maxn;j++){
sum[j*pri[i]]+=mu[j];
}
}
for(rg int i=1;i<maxn;i++){
sum[i]+=sum[i-1];
}
}
void solve(int n,int m){
ll ans=0;
for(rg int l=1,r;l<=n;l=r+1){
r=std::min(n/(n/l),m/(m/l));
ans+=1LL*(sum[r]-sum[l-1])*(n/l)*(m/l);
}
printf("%lld\n",ans);
}
int main(){
int n,m,t;
scanf("%d",&t);
xxs();
while(t--){
scanf("%d%d",&n,&m);
if(n>m) std::swap(n,m);
solve(n,m);
}
return 0;
}
例题三、gcd
题目描述
有 \(n\) 个正整数 \(x_1 \sim x_n\)，初始时状态均为未选。
有 \(m\) 个操作，每个操作给定一个编号 \(i\)，将 \(x_i\) 的选取状态取反。
每次操作后，你需要求出选取的数中有多少个互质的无序数对。
输入格式
第一行两个整数 \(n,m\)。第二行 \(n\) 个整数 \(x_1 \sim x_n\)。接下来 \(m\) 行每行一个整数。
输出格式
\(m\) 行，每行一个整数表示答案。
样例
样例输入
4 5
1 2 3 4
1
2
3
4
1
样例输出
0
1
3
5
2
数据范围与提示
对于 \(20\%\) 的数据，\(n,m<=1000\)。
对于另外 \(30\%\) 的数据，\(x_i<=100\)。
对于 \(100\%\) 的数据，\(n,m<=200000,x_i<=500000,1<=i<=n\)。
分析
我们设 \(s[i]\) 为当前选取的数中有多少数是 \(i\) 的倍数
设 \(f[i]\) 为 \(gcd\) 是 \(i\) 的倍数的数对的个数
那么就有 \(f[i]=\frac{s[i](s[i]+1)}{2}\)
设 \(g[i]\) 为 \(gcd\) 是 \(i\) 的数对的个数
我们要求的就是 \(g[1]\)
根据莫比乌斯反演的形式二
\(f[i]=\sum\limits_{i|d} g[d]\)
则
\(g[i]=\sum\limits_{i|d} \mu(\frac{d}{i})f[d]\)
而 \(f[d]\) 可以在 \(O(\sqrt{n})\) 的时间内更新
所以我们只要开一个变量一直记录答案即可
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e6+5;
bool not_pri[maxn],vis[maxn];
int n,m,mu[maxn],pri[maxn],a[maxn],cnt[maxn];
long long g[maxn],ans;
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
}
int main(){
xxs();
n=read(),m=read();
for(rg int i=1;i<=n;i++){
a[i]=read();
}
rg int aa,bb,cc,now;
for(rg int i=1;i<=m;i++){
aa=read();
vis[aa]^=1;
cc=a[aa];
bb=sqrt(cc);
if(vis[aa]){
for(rg int j=1;j<=bb;j++){
if(cc%j==0){
now=j;
cnt[now]++;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
if(j*j!=cc){
now=cc/j;
cnt[now]++;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
}
}
}
} else {
for(rg int j=1;j<=bb;j++){
if(cc%j==0){
now=j;
cnt[now]--;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
if(j*j!=cc){
now=cc/j;
cnt[now]--;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
}
}
}
}
printf("%lld\n",ans);
}
return 0;
}
例题四、数表
分析
要求的式子是
而且还要满足 \(\sigma(gcd(i,j)) \leq a\)
我们先不考虑 \(a\) 的限制
令 \(n<m\)，开始推式子
令 \(T=dk\)
这个式子前半部分可以整除分块计算，关键是后半部分
设 \(f(t)=\sum\limits_{d|T}\sigma(d) \mu(\frac{T}{d})\)
要满足的条件是 \(\sigma(d) \leq a\)
我们会发现，随着 \(a\) 的增大，满足条件的 \(d\) 会越来越多
所以我们可以把询问离线下来，按照 \(a\) 从小到大排序
每次把新的能够做出贡献的 \(d\) 的答案更新至 \(f(t)\)
因为分块时还要询问前缀和，所以采用数状数组维护
更新的次数是 \(\sum\limits_{i=1}^{n}\frac{n}{i}=nlogn\)
总复杂度 \(O(Q\sqrt{n}logn+nlog^2n)\)
代码
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#include<cmath>
#include<iostream>
#define rg register
inline int read(){
rg int x=0,fh=1;
char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e6+5;
const long long mod=1LL<<31;
int pri[maxn],mu[maxn],q,mmax;
long long f[maxn];
bool not_pri[maxn];
struct asd{
int n,m,a,id;
}b[maxn];
bool cmp(asd aa,asd bb){
return aa.a<bb.a;
}
struct jie{
int id;
long long val;
}c[maxn];
bool cmp2(jie aa,jie bb){
return aa.val<bb.val;
}
void xxs(){
not_pri[0]=not_pri[1]=1;
f[1]=mu[1]=1;
for(rg int i=2;i<=mmax;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
f[i]=i+1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<=mmax;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0){
mu[i*pri[j]]=0;
f[i*pri[j]]=f[i]*f[pri[j]]-pri[j]*f[i/pri[j]];
break;
} else {
mu[i*pri[j]]=-mu[i];
f[i*pri[j]]=f[i]*f[pri[j]];
}
}
}
for(rg int i=1;i<=mmax;i++){
c[i].id=i;
c[i].val=f[i];
}
std::sort(c+1,c+1+mmax,cmp2);
}
int lb(rg int xx){
return xx&-xx;
}
long long tr[maxn];
void ad(rg int wz,rg long long val){
for(rg int i=wz;i<maxn;i+=lb(i)){
tr[i]+=val;
if(tr[i]>=mod) tr[i]-=mod;
}
}
long long cx(rg int wz){
rg long long nans=0;
for(rg int i=wz;i>0;i-=lb(i)){
nans+=tr[i];
if(nans>=mod) nans-=mod;
if(nans<0) nans+=mod;
}
return nans;
}
long long cxqj(rg int l,rg int r){
rg long long nans=cx(r)-cx(l-1);
if(nans<0) nans+=mod;
return nans;
}
long long ans[maxn];
void updat(rg int now){
rg long long nans=f[now]%mod;
for(rg int i=now;i<=mmax;i+=now){
ad(i,nans*mu[i/now]%mod);
}
}
int main(){
q=read();
for(rg int i=1;i<=q;i++){
b[i].n=read(),b[i].m=read(),b[i].a=read();
b[i].id=i;
mmax=std::max(mmax,std::max(b[i].n,b[i].m));
}
xxs();
std::sort(b+1,b+1+q,cmp);
rg int head=1;
for(rg int i=1;i<=q;i++){
while(c[head].val<=b[i].a && head<=mmax){
updat(c[head].id);
head++;
}
if(b[i].n>b[i].m) std::swap(b[i].n,b[i].m);
for(rg int l=1,r;l<=b[i].n;l=r+1){
r=std::min(b[i].n/(b[i].n/l),b[i].m/(b[i].m/l));
ans[b[i].id]+=(b[i].n/l)%mod*(b[i].m/l)%mod*cxqj(l,r)%mod;
if(ans[b[i].id]>=mod) ans[b[i].id]-=mod;
}
}
for(rg int i=1;i<=q;i++){
printf("%lld\n",ans[i]);
}
return 0;
}
例题五、数字表格
分析
求
设 \(n<m\)，开始推式子
设 \(T=dk\)，则
和上一道题一样，枚举每一个 \(d\) ，然后把贡献加到它的倍数里
要注意的是，\(f(d)\) 的 \(-1\) 次幂要提前处理出来，否则复杂度会多一个 \(log\)
代码
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#include<cmath>
#include<iostream>
#define rg register
inline int read(){
rg int x=0,fh=1;
char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e6+5,mod=1e9+7;
int pri[maxn],mu[maxn],t,mmax,n,m,f[maxn],sum[maxn],ny[maxn],ans[maxn];
bool not_pri[maxn];
int ksm(rg int ds,rg long long zs){
rg int nans=1;
while(zs){
if(zs&1LL) nans=1LL*nans*ds%mod;
ds=1LL*ds*ds%mod;
zs>>=1LL;
}
return nans;
}
void xxs(){
not_pri[0]=not_pri[1]=1;
f[1]=mu[1]=1;
for(rg int i=2;i<=mmax;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<=mmax;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0){
mu[i*pri[j]]=0;
break;
} else {
mu[i*pri[j]]=-mu[i];
}
}
}
for(rg int i=2;i<=mmax;i++){
f[i]=f[i-1]+f[i-2];
if(f[i]>=mod) f[i]-=mod;
}
for(rg int i=1;i<=mmax;i++){
ny[i]=ksm(f[i],mod-2);
}
for(rg int i=1;i<=mmax;i++){
ans[i]=1;
}
rg int now;
for(rg int i=1;i<=mmax;i++){
for(rg int j=i;j<=mmax;j+=i){
now=mu[j/i];
if(now==-1) ans[j]=1LL*ans[j]*ny[i]%mod;
else if(now==1) ans[j]=1LL*ans[j]*f[i]%mod;
}
}
sum[0]=1;
for(rg int i=1;i<=mmax;i++){
sum[i]=1LL*sum[i-1]*ans[i]%mod;
}
}
int getsum(rg int l,rg int r){
return 1LL*sum[r]*ksm(sum[l-1],mod-2)%mod;
}
int main(){
mmax=1000000;
xxs();
t=read();
while(t--){
n=read(),m=read();
if(n>m) std::swap(n,m);
rg int nans=1;
for(rg int l=1,r;l<=n;l=r+1){
r=std::min(n/(n/l),m/(m/l));
nans=1LL*nans*ksm(getsum(l,r),1LL*(n/l)*(m/l))%mod;
}
printf("%d\n",nans);
}
return 0;
}
例题六、Crash的数字表格 / JZPTAB
分析
求
设 \(n<m\)，开始推式子
设 \(T=dk\)，则
前面的可以整除分块进行计算
设 \(f(T)=\sum\limits_{d|T}d\mu(d)\)
考虑如何求出 \(f(T)\)
这个东西可以线性筛 \(O(n)\) 预处理
当\(T=1\)时，\(f(T)=1\)
当 \(T\) 为质数时，\(f(T)=1-T\)
当 \(gcd(a,b)=1\) 时 \(f(ab)=f(a)\times f(b)\)
当 \(b\) 为质数且 \(a\%b=0\) 时，\(f(ab)=f(a)\)
代码
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#include<cmath>
#include<iostream>
#define rg register
inline int read(){
rg int x=0,fh=1;
char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e7+5,mod=20101009;
int pri[maxn],mu[maxn],t,mmax,n,m,f[maxn],sum[maxn];
bool not_pri[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
f[1]=mu[1]=1;
for(rg int i=2;i<=mmax;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
f[i]=1-i;
f[i]+=mod;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<=mmax;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0){
mu[i*pri[j]]=0;
f[i*pri[j]]=f[i];
break;
} else {
mu[i*pri[j]]=-mu[i];
f[i*pri[j]]=1LL*f[i]*f[pri[j]]%mod;
}
}
}
for(rg int i=1;i<=mmax;i++){
f[i]=1LL*f[i]*i%mod;
}
for(rg int i=1;i<=mmax;i++){
sum[i]=sum[i-1]+f[i];
if(sum[i]>=mod) sum[i]-=mod;
}
}
int getsum(rg int now){
return 1LL*(now)*(now+1)/2LL%mod;
}
int main(){
mmax=10000000;
xxs();
n=read(),m=read();
if(n>m) std::swap(n,m);
rg int nans=0,cs;
for(rg int l=1,r;l<=n;l=r+1){
r=std::min(n/(n/l),m/(m/l));
cs=sum[r]-sum[l-1];
if(cs<0) cs+=mod;
nans+=1LL*cs%mod*getsum(n/l)%mod*getsum(m/l)%mod;
if(nans>=mod) nans-=mod;
}
printf("%d\n",nans);
return 0;
}
