数学(2)——莫比乌斯反演做题记录
Luogu P1447 [NOI2010] 能量采集
求
可以转换为
考虑对蓝色部分莫比乌斯反演
设
这个 \(g(\mathrm d)\) 可以帮我们做数论分块。
首先可以看出前面的部分
我们可以直接枚举 \(i\) 而不是 \(i|d\)
所以有
答案就是
代码通过数论分块和线性筛实现可以就通过这题。并且首先保证 \(\color{red}n\leq m\) 。
#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
inline int read()
{
int x=0,w=1;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')ch=getchar(),w=-1;
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();
return x*w;
}
int n,m,tot;
const int N=100005;
ll phi[N],pri[N],sum[N];
bool vis[N];
void Phi()
{
phi[1]=1;
for(int i=2;i<=100000;i++)
{
if(!vis[i])pri[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&i*pri[j]<=100000;j++)
{
if(i*pri[j]>100000)continue;
vis[i*pri[j]]=1;
if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]];
else {phi[i*pri[j]]=phi[i]*pri[j];continue;}
}
}
}
int main()
{
n=read();m=read();
Phi();
for(int i=1;i<=100000;i++)sum[i]=sum[i-1]+phi[i];
ll ans=0;
ll l=1,r=0;
while(l<=min(n,m))
{
r=min(n/(n/l),m/(m/l));
ans+=(ll)(sum[r]-sum[l-1])*(n/l)*(m/l);
l=r+1;
}
printf("%lld\n",2ll*ans-(ll)n*m);
return 0;
}
Luogu P2257 YY的GCD
这个题有个神必的巧妙做法，参考这篇题解。
这个题是
(下文质数集合用 \(\mathbb P\) 表示且首先同样保证 \(\color{red}n\leq m\))
我们考虑直接枚举质数：
我们考虑莫比乌斯反演技巧，设：
那我们发现答案就是
再对答案进行莫比乌斯反演
消去 \(p\times \mathrm d\) 可以考虑换元为 \(k\)
这样数论分块，筛法，前缀和就可以实现了。
#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
inline int read()
{
int x=0;int w=1;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')ch=getchar(),w=-1;
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();
return x*w;
}
int n,m;
const int N=10000000;
const int M=10000005;
ll pri[M/10],phi[M],s[M],g[M];
bool vis[M];
ll mn;
ll ans,cnt;
void Phi()
{
phi[1]=1;
for(int i=2;i<=mn;i++)
{
if(!vis[i])phi[i]=-1,pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<=mn;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0)break;
else phi[i*pri[j]]=-phi[i];
}
}
for(int i=1;i<=cnt;i++)
{
for(int j=1;j*pri[i]<=mn;j++)
g[j*pri[i]]+=phi[j];
}
for(int i=1;i<=mn;i++)
{
s[i]=s[i-1]+g[i];
}
}
int main()
{
int T;
T=read();
mn=10000000;
Phi();ll ans=0;ll r;
while(T--)
{
n=read();m=read();
mn=(ll)min(n,m);ans=0;
for(int i=1;i<=mn;i=r+1)
{
r=min(n/(n/i),m/(m/i));
ans+=(ll)(s[r]-s[i-1])*(ll)(n/i)*(ll)(m/i);
}
printf("%lld\n",ans);
}
return 0;
}
Luogu P2261 [CQOI2007]余数求和
这个题好像很水，因为好像只是一个入门的除法分块练习题。
这样就可以做除法分块
就是对 \(\displaystyle {\color{blue}{\sum_{i=1}^n(i\times\lfloor\frac{k}{i}\rfloor)}}\) 做除法分块。
代码如下
for(ll i=1;i<=n;i=r+1)
{
if(k/i)
r=min((ll)n,k/(k/i));
else
r=n;
sum+=(ll)(k/i)*(ll)(r-i+1)*(i+r)/2ll;
}
Luogu P1390 公约数的和
这个题是求
我们先考虑求出
采用一些神必技巧，枚举一个数 \(k\) ，然后改成对 \(k\) 求和得到
我们改为 \(\epsilon\) $\epsilon$ 是为了方便我们做 \(\texttt{Dirichlet}\) 卷积。
所以我们可以转化为
不过别忘了，这个式子只是 \(\displaystyle \sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\)
我们推导得
用之前那个式子对蓝色的 \(\displaystyle{\color{blue}{\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)}}\) 部分做数论分块求解。
代码：
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cmath>
#include<queue>
#include<map>
#include<stack>
//#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define INL inline
//Tosaka Rin Suki~
using namespace std;
const int N=2000005;
INL int read()
{
int x=0;int w=1;
char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')
{x=(x<<1)+(x<<3)+ch-48,ch=getchar();}
return x*w;
}
INL int mx(int a,int b){return a>b?a:b;}
INL int mn(int a,int b){return a<b?a:b;}
int mu[N];
ll pri[N],cnt;
bool vis[N];
int n;
void Mu()
{
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!vis[i])pri[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pri[j]<=n;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j])mu[i*pri[j]]=-mu[i];
else {mu[i*pri[j]]=0;break;}
}
}
for(int i=1;i<=n;i++)mu[i]+=mu[i-1];
}
int main()
{
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
n=read();
ll ans=0;
Mu();
for(int i=1;i<=n;i++)
{
ll sum=0;
for(ll l=1,r;l<=n/i;l=r+1)
{
r=(n/i)/((n/i)/l);
sum+=(ll)(mu[r]-mu[l-1])*(ll)((n/i)/l)*(ll)((n/i)/l);
}
ans+=sum*i*1ll;
}
ans-=(ll)((ll)n*(ll)(n+1))>>1ll;
ans>>=1;
printf("%lld\n",ans);
return 0;
}
LibreOJ #6539. 奇妙数论题
大家先回顾一下上一题，这题的部分分就是上一题。
题如其名，一道巧妙莫反题。
给定长为 \(n\) 的排列 \(a\) 求
我们有部分分是 \(a_i=i\) 的，所以一直到 \(70\) 分的部分分就是求：
就上一题稍微改一下就有了：
正解的瓶颈其实就是排列变成无序后无从下手。
所以我们先将式子转换成：
这里就是根据上一题的套路，\(\gcd(i,j)\) 看成 \(\gcd(i,j)\times 1\) 于是可以在数论分块的部分直接带入对应下标 \(a_i\) 的 \(\gcd\)。
但是求 \(\gcd(a_i,a_j)\) 也就是 \(\gcd(a_{idx},a_{jdx})\) 的复杂度降不下来，考虑怎么做。
设：
我们观察这个函数，其实就是要求在 \(a\) 这个排列中下标分别为 \(x\) 的倍数的元素彼此之间 \(\gcd\) 的和。
我们设 \(\mathbb S_x\) 是大小为 \(x\) 的倍数的下标的集合。
于是有：
我们进一步抽象，我们将 \(i\in\mathbb S_x\) 映射到每一个 \(a_i\)，\(\forall i\in \mathbb S_x,\) 设 \(a_i\) 的集合为 \(\mathbf S_x\)。
现在有：
\(i,j\) 的集合与范围发生了变化，显然我们不能再将 \([\gcd(i,j)=d]\) 化为 \([\gcd(i,j)=1]\) 去做。
我们有另一种方法求 \([x=d]\) 这一类的式子。
我们观察 \([x=d]\) 长得很像 \(\epsilon(x)=[x=1]\)。
我们有 \(\mu*1=\epsilon\) 所以有 \(\displaystyle \epsilon(x)=\sum_{d|x}\mu(d)\)。
那么我们的 \([x=d]\) 可以化为 \(\displaystyle[\frac{T}{d}=1\space\mathrm{and}\space T|d]\)。
类似的我们展开得 \(\displaystyle [x=d]=\sum_{d|T,T|x}\mu(\frac{T}{d})\)
得到（这个地方看官方题解看了很久才懂的，感觉很 tricky）：
我们知道有 \(\mu * \mathrm{id}=\varphi\)。
观察上面蓝色部分的 \(\displaystyle{\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}}\) 就是 \(\varphi(T)\)
最后得到：
我们回到 \(\texttt{ans}\) 的部分。
有：
具体实现方法：
- 首先我们筛出得到 \(\forall x\in[1,n]\) 每个 \(\mathbf S_x\)。
- 依次统计对于每个 \(\mathbf S_x\) 每个元素的的因子的贡献。
- 得到 \(f(x)\) 并筛出 \(\mu(x)\) 后直接求和即可。
复杂度分析的话，看到 loj 官方题解的评论区中有这样的结论：
我们设对于 \(x\)，\(x\) 的因子个数为 \(d(x)\)。
As Vaclav Kotesovec said Aug 30 2018, \(\displaystyle \sum_{k\leq n}d^2(k)\)is asymptotic to \(\Theta(\frac{1}{\pi^2}n\log^3 n+n\log^2 n)\).
代码（其实就是出题人的写法，好像多数采用了埃筛）：
#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define INL inline
#define Re register
//Tosaka Rin Suki~
INL int read()
{
int x=0,w=1;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x*w;
}
const int N=1e5+5,MOD=1e9+7;
std::vector <int> factors[N];
std::vector <int> pack;
int mu[N],g[N],f[N],d[N],n,a[N],ans;
int main()
{
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
n=read();for(int i=1;i<=n;i++)a[i]=read();
mu[1]=1;
for(int i=1;i<=n;i++)
{
factors[i].push_back(i);
for(int j=i*2;j<=n;j+=i)
mu[j]-=mu[i],factors[j].push_back(i);
}
for(int i=1;i<=n;i++)
for(int j=i;j<=n;j+=i)
g[j]=(g[j]+mu[j/i]*i)%MOD;
for(int i=1;i<=n;i++)
{
for(int j=i;j<=n;j+=i)
for(int k=0;k<factors[a[j]].size();k++)
++d[factors[a[j]][k]],pack.push_back(factors[a[j]][k]);
for(int j=0;j<pack.size();j++)
if(d[pack[j]])
f[i]=(f[i]+1ll*d[pack[j]]*d[pack[j]]*g[pack[j]])%MOD,d[pack[j]]=0;
pack.clear();
}
for(int d=1;d<=n;d++)
{
int s=0;
for(int x=1;x<=n/d;x++)s=(s+1ll*mu[x]*f[d*x])%MOD;
ans=(ans+1ll*s*d)%MOD;
}
printf("%d\n",ans);
return 0;
}
Luogu P3704 [SDOI2017]数字表格
在机房复盘 SDOI 2017 的模拟赛第一次场 A 黑题写题解纪念一下。
这个题要你求：
其中 \(\mathrm{fib_{i}}\) 表示斐波拉契数列第 \(i\) 项。
\(\prod\) 连乘运算做不了 \(\sum\) 的一些卷积公式（尤其是 \(\texttt{Dirichlet}\) 卷积）。我们考虑把这个式子换成可以用 \(\sum\) 推的情况。
首先一般性地，我们保证 \({\color{red}{n\leq m}}\)。
我们枚举 \(\displaystyle\sum_i\sum_j f(i,j)\) 这一类式子的时候会习惯性地枚举一个 \(k\) 让其满足和式内部的条件。
这题即使是 \(\prod\) 也可以套路化地使用这个方法。
所以原式可以转换为
我们发现每个 \(\mathrm{fib}_d\) 都会与自己相乘很多次，于是我们可以将每个 \(d\) 的贡献考虑写成幂的形式。
显然每个 \(\mathrm{fib}_d\) 的幂的指数就是其等于 \(\gcd(i,j)\) 对应的数对 \((i,j)\) 的个数。
现在我们有
我们把指数特地标记出来。因为我们发现他是一个老朋友。
我们在做形似 \(\displaystyle\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\) 的一类莫反题是会将式子转换成 \(\displaystyle \sum_{k=1}^n\bigg(k{\color{blue}{\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=k]}}\bigg)\)
我们考虑怎么将这个式子转化成可做的形式（以下皆是莫反题老套路，跟不上请先去做一些基础莫反题，不是幽灵乐团，别信 CYJian）。
这个东西已经可以 \(\Theta(\sqrt n)\) 级别的数论分块做了，但是我们要把其代入原式继续看。
现在原式变成了：
我们设 \(T=k\times d\)。
我们发现这样方便让我们将 \(\displaystyle\sum_{d=1}^{\lfloor \frac{n}{k}\rfloor}\) 与 \(\mathrm{fib}_k\) 以幂的形式结合。所有变量用 \(T\) 表示出来后贡献到一个 \(\prod\) 中。
写出来就是：
蓝色部分我们可以通过类似筛法预处理的方法求出：
- 线筛 \(\mu\) 函数。
- 对于每一个 \(k\) 枚举其倍数 \(T\) ，乘上这个 \(\mathrm{fib}_k\) 算入这个 \(T\) 的贡献。
- 因为 \(\mu\) 只有三种取值，有影响的只有 \(\{1,-1\}\) 两种。
- 所以问题就变成了乘上 \(\mathrm {fib}_{k}\) 或除以 \(\mathrm {fib}_k\)。这里就变成 \(\mathrm{inv fib}\)（逆元）和 \(\mathrm{fib}\) 两种情况。
外层我们做数论分块的时候，每次求内部只需要一次 \(\Theta (\log p)\)（\(p\) 是模数）和 \(\Theta(\log(\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor))\)的快速幂分别求逆元和与外层指数的幂。
预处理复杂度是 \(\Theta(n\log n+n\log p)\)，数论分块复杂度是 \(\Theta(t\sqrt n (\log p+\log \lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor))\) 小写 \(t\) 是数据组数。
实现看代码（考试时忘了写逆元卡了好久导致 T2 LCT 都没写完）：
#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define INL inline
#define Re register
//Tosaka Rin Suki~
INL int read()
{
int x=0,w=1;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x*w;
}
const int N=1e6+5;
const ll MOD=1e9+7;
ll f[N],finv[N],prime[N],mu[N],cnt,prod[N];
bool vis[N];
INL ll fpow(ll x,ll p)
{
ll res=1;
while(p)
{
if(p&1)res=(1ll*res*x)%MOD;
x=(1ll*x*x)%MOD;p>>=1;
}
return res;
}
INL void get_mu(ll n)
{
f[0]=0,f[1]=1;
for(int i=2;i<=n;i++)f[i]=(f[i-1]+f[i-2])%MOD;
for(int i=0;i<=n;i++)finv[i]=fpow(f[i],MOD-2),prod[i]=1;
vis[1]=1,vis[0]=1,mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!vis[i])prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
{
vis[prime[j]*i]=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++)
{
for(int j=i;j<=n;j+=i)
{
int now=j/i;
if(mu[now]==-1)prod[j]=1ll*prod[j]*finv[i]%MOD;
else if(mu[now]==1)prod[j]=1ll*prod[j]*f[i]%MOD;
}
}
for(int i=1;i<=n;i++)prod[i]=1ll*prod[i-1]*prod[i]%MOD;
}
int main()
{
//freopen("product.in","r",stdin);
//freopen("product.out","w",stdout);
int T=read();
get_mu(1000000);
while(T--)
{
ll n=read(),m=read();
if(n<m)std::swap(n,m);
ll l=1,r,ans=1;
while(l<=m)
{
r=std::min(n/(n/l),m/(m/l));
ans=1ll*ans*fpow(1ll*prod[r]*fpow(prod[l-1],MOD-2)%MOD,1ll*(n/l)*(m/l)%(MOD-1))%MOD;
l=r+1;
}
printf("%lld\n",(ans%MOD+MOD)%MOD);
}
return 0;
}
