莫比乌斯反演与数论函数

command_block · 2018-07-20 09:01:35 · 算法·理论

本博客的第一篇博文哦!纪念。

---- 莫比乌斯反演，本质是利用莫比乌斯函数与其他函数间卷积关系，对函数做一系列简化，从而更高效的解决问题。

---- 数论函数,更广阔的天地,因为\mu∈\text{数论函数}

0.写在前面

这里面式子可能比较多,大家多看看证明,尤其是关键证明以及题目推导过程,硬记下的结论越少,越不容易忘记。

(看不清式子可以放大网页)

由于写的时候匆忙,公式可能有错误,请私信或者在下方评论。

每个新东西提出后都会做题或者证明,权当熟悉。

不是一个下午就能看懂的,千万不要直接弃疗

• 给大家简要介绍一下符号的意思。

有时候以(a,b)简写gcd(a,b)。

有\sum\limits_{d|6}d=1+2+3+6 (枚举6的因数)

1.狄利克雷卷积与数论函数

——介绍了狄利克雷卷积,数论函数(积性函数)的综合

换句话说就是\large{(f*g)(n)=\sum\limits_{x*y=n}f(x)g(y)}

(不过一般不这样写,因为不方便化式子。)

比如说计算(f*g)(10)

可得f(1)g(10)+f(2)g(5)+f(5)g(2)+f(10)g(1)

(对于某个新东西,能做到的话,理解的最好方式就是手动模拟)

狄利克雷卷积满足以下运算规律(显而易见,不证)：

交换律——f*g=g*f;

结合律——(f*g)*h=f*(g*h);

(狄利克雷卷积是一个对称的结构)

下面的推导都基于狄利克雷卷积,有必要好好理解。

介绍几个简单的数论函数,让你有个基本概念。

附上两个后面会讲的稍微复杂的例子:

这里的p_1...p_m是指x的m个素因子,也就是把x分解的结果。

把\sum合并(根据约数集合的合并),由于a与b互质,所以t与d互质。

(根据积性可以把F_1(d)F_1(t)化为F_1(dt),F_2类似)

考虑数学归纳法(跟构造方法有点像)

设e=G*F,且F是积性函数。

我们要证明对于任意的(a,b)=1的a,b,都满足G(ab)=G(a)G(b)

因为e=\sum\limits_{d|n}F(d)G(n/d)

• 1) a或b为1

  当a=b=1时,e(1)=1=F(1)G(1)

  由于积性,根据上文F(1)=1,得G(1)=1

  所以当a或b=1时,结论显然成立。

• 2) ab>1

  我们考虑已证明了所有的a'<a;\ b'<b时,结论成立。

  G(ab)=\sum\limits_{d|ab}F(d)G(ab/d)-\sum\limits_{d|ab,d≠1}F(d)G(ab/d)

  由e=G*F,得\sum\limits_{d|ab}F(d)G(ab/d)=e(ab),由于ab>1所以e(ab)=0。

  =-\sum\limits_{d|ab,d≠1}F(d)G(ab/d) =-\sum\limits_{i|a,j|b,ij≠1}F(ij)G(ab/ij)

  根据所有的a'<a;\ b'<b时,积性成立得

  =-\sum\limits_{i|a,j|b,ij≠1}F(i)G(j)F(a/i)G(b/j) =-\sum\limits_{i|a,j|b,ij≠1}F(i)G(a/i)F(j)G(b/j)

  下面拆分求和号(乘法分配律),这一步有点难。

  =-\sum\limits_{i|a}F(i)G(a/i)\sum\limits_{j|b,ij≠1}F(j)G(b/j)

  要求ij≠1,即排除i=j=1的情况。

  =F(1)F(1)G(a)G(b)-\sum\limits_{i|a}F(i)G(a/i)\sum\limits_{j|b}F(j)G(b/j)

  由e=G*F,得\sum\limits_{i|a}F(i)G(a/i)\sum\limits_{j|b}F(j)G(b/j)=e(a)e(b)=0得

  =F(1)F(1)G(a)G(b)

  由于积性,根据上文F(1)=1,得

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

2.莫反公式的推导:

——推导了\mu函数,莫比乌斯反演定理

有两个单变量函数 F 与 f 。

设函数 F 与 f 有如下关系：

理解了上述狄利克雷卷积的定义后,不难发现。

我们只要想办法求出I^{-1}就好了。

大佬Möbius把I^{-1}命名为\mu

问题在于这个函数里面是什么呢?

\mu是个积性函数,因为积性函数的逆还是积性函数。

这里有个小技巧,研究一个积性函数,先研究其在质数的幂时的表现。

比如说\mu(p^k)(p是质数)。

k=0时,显然有\mu(1)=1

k=1时,由于\sum\limits_{d|p}I(d)\mu(p/d)=e(p)=0

注意到p为素数,得到I(p)\mu(1)+I(1)\mu(p)=e(p)=0

因为\mu(1)=1得1+\mu(p)=0,即\large{\mu(p)=-1}

k>1时,由于\sum\limits_{d|p^k}I(d)\mu(p^k/d)=e(p^k)=0

注意到p为素数,得到I(p^k)\mu(1)+I(p^{k-1})\mu(p)+...+I(1)\mu(p^k)=e(p^k)=0

所有的I都为1,得到\mu(1)+\mu(p)+\mu(p^2)...+\mu(p^k)=0

代入\mu(1)=1,\mu(p)=-1得到\mu(p^2)...+\mu(p^k)=0

考虑上式k=2的情况,得到\mu(p^2)=0

考虑上式k=3的情况,得到mu(p^2)+mu(p^3)=0也就是mu(p^3)=0

……,得到k=任意更大的数时mu(p^k)=0

利用数学归纳法,证明了k>1时,\mu(p^k)=0。

大家缓口气哦,\mu函数的终极定义马上就证明出来了。

由于上文积性函数性质2: \mu(n=p_1^{k_1}p_2^{k_2}...p_m^{k_m})=\mu(p_1^{k_1})\mu(p_2^{k_2})...\mu(p_m^{k_m})

因为k>1时,\mu(p^k)=0

可以想象到当某一个k>1时,\mu(n)=0

不然的话所有的k=1。

可以想象到\mu(n=p_1p_2...p_m)=(-1)^{m}

mu 函数的\color{blue}\text{定义}浮出水面。

现在来介绍莫比乌斯函数是怎么用来反演的。

• 反演公式:

  • 嵌入式莫比乌斯反演：

  由\mu*I=e得\sum\limits_{d|n}\mu(d)=[n=1]

  注意到[n|m][n/m=1]=[n=m]

  则有[n|m]\sum\limits_{d|(n/m)}\mu(d)=[n=m]

  初学者可以只掌握这一种。

  • 约数的莫比乌斯反演：

  若：\large{F(n)=∑\limits_{d|n}f(d)}

  则：\large{f(n)=∑\limits_{d|n}μ(d)F(n/d)}

  又作f(n)=∑\limits_{d|n}μ(n/d)F(d)(和上面那个式子等价)

  根据\mu=I^{-1}有f*I=F\rightarrow f=F*\mu。

  如果你们喜欢和式的话:证明2

  • 倍数的莫比乌斯反演：

  若：\large{F(n)=∑\limits_{n|d}f(d)}

  则：\large{f(n)=∑\limits_{n|d}μ(d/n)F(d)}

  （这个比较特殊,并不是狄利克雷卷积的形式）

  又作F(n)=\sum\limits_{k}^∞f(kn)与f(n)=\sum\limits_{k}^∞\mu(k)F(kn)

  \color{blue}\text{证明:} ∑\limits_{n|d}μ(n/d)F(d)

  将F(n)的定义代入。

  =∑\limits_{n|d}μ(n/d)\sum\limits_{d|t}f(t) =∑\limits_{k}^{∞}μ(k)\sum\limits_{nk|t}f(t)

  第一个求和号对k没有限制,我们交换求和号,考虑枚举t,看有什么对应的k:

  =∑\limits_{t}f(t)\sum\limits_{nk|t}\mu(k)=∑\limits_{t}f(t)\sum\limits_{k|(t/n),n|t}\mu(k)

  根据嵌入式反演能得到\sum\limits_{k|(t/n),n|t}\mu(k)=[t=n]

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

3.知识储备:

• 分解质因数求法,单个mu(n)

直接按照定义,分解因数暴力求。

代码被和谐
//Time:O(sqrt(n))

• 线性筛法,mu(1~n) (正宗)

利用\mu函数的积性性质。

想深入学习线性筛法的看这里。

#include<iostream>
#include<cstring>
#define MaxNum 10000100
using namespace std;
bool e[MaxNum];
int p[MaxNum],tn,mu[MaxNum],n;
void Mobius()
{
  e[1]=1;mu[1]=1;
  for (int i=2;i<=n;i++){
    if (!e[i]){p[++tn]=i;mu[i]=-1;}
    for (int j=1;j<=tn;j++){
      if (p[j]*i>n)break;
      mu[p[j]*i]=i%p[j]==0 ? 0 : -mu[i];
      e[p[j]*i]=1;
      if (i%p[j]==0)break;
    }
  }
}
int main()
{
  cin>>n; 
  Mobius();
  for(int i=1;i<=n;++i)
   printf("%d: %d\n",i,mu[i]);
  return 0;
}//Time:O(n)

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

4.第一道题:整除分块

(除了特殊说明,除法均为整除)

讲了那么多,那么问题来了:我们问什么要学莫比乌斯反演?这有什么用呢?

来看一道题目。

让你求\large{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^m[gcd(i,j)=1]}，也就是互质数对数。(m\leq n)

大力两重循环+欧几里得O(n^2logn)?

悄悄告诉你,正解O(\sqrt{n})!

• 即$[d|(i,j)]=[d|i][d|j]$, 人话就是 : 即是$i$的约数又是$j$的约数的数是$i,j$的公约数。

则有[(i,j)=1]=\sum\limits_{d|i,d|j}\mu(d),这比较重要。

观察后面的\sum\limits_{d|i}^{n}\sum\limits_{d|j}^m1=(\sum\limits_{d|i}^{n}1)(\sum\limits_{d|j}^m1)

而\sum\limits_{d|i}^{n}1相当于n以内d的倍数的个数,显然为\lfloor n/d\rfloor

就得到ANS=\sum\limits_{d=1}^n\mu(d)\lfloor n/d\rfloor\lfloor m/d\rfloor

这样子就可以O(n)计算了。

看到和式后面的\lfloor n/d\rfloor \lfloor m/d\rfloor,是可以整除分块的。

问题是还乘了个\mu(d),根据套路弄个前缀和搞定。

(如果这句听不懂,把入门小记再看一遍)

即\mu([l,r])的和乘上\lfloor n/d\rfloor \lfloor m/d\rfloor。

Code:

#include<iostream>
#define MaxNum 100100
using namespace std;
int p[MaxNum/8],tn,mu[MaxNum],n,m;
bool e[MaxNum];
long long calc(int n,int m)
{
  long long ans=0;
  for (int l=1,r=0;l<=min(n,m);l=r+1){
    r=min(n/(n/l),m/(m/l));// 分段
    ans+=1ll*(mu[r]-mu[l-1])*(n/l)*(m/l);
  }return ans; 
}
int main()
{
    e[1]=1;mu[1]=1;
    for (int i=2;i<=MaxNum;i++){
      if (!e[i]){p[++tn]=i;mu[i]=-1;}
      for (int j=1;j<=tn;j++){
        if (p[j]*i>MaxNum)break;
        mu[p[j]*i]=i%p[j]==0 ? 0 : -mu[i];
        e[p[j]*i]=1;
        if (i%p[j]==0)break;
      }
    }for (int i=2;i<=MaxNum;i++)mu[i]+=mu[i-1];
    //筛法弄出mu前缀和,参看前文
    cin>>n>>m;
    cout<<calc(n,m)<<endl;
    return 0;
}

后面你会看到,莫比乌斯反演不分块复杂度就是一堆垃圾。

分块是肯定要分块的,这辈子都要分块的

恭喜你,已经在莫比乌斯反演的路上迈出了第一步!

附送可以利用上述代码AC的练手题Link1+Link2

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

5.技巧进阶:莫反常用结论以及练习

莫比乌斯反演是极具技巧性的。

主要的思维难度就建立在反演和分块两部分。

看题吧。

第一题P2568 GCD

求\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=prime]

在上一题中,我们解决了求\sum\limits_{i=1}^{n}\sum\limits_{j=1}^n[(i,j)=1]的问题,而且我们已经能够做到O(\sqrt{n})求解了。

我们尝试求\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=d]

• 注意到,[(i,j)=d]蕴含d|i且d|j.

• \color{blue}\text{技巧:}$ 既然$i,j$都是$d$的倍数,我们不妨把$i,j$都除以$d

  此时,后面原有的i需要还原成id,原有的j要还原成jd

• [(id,jd)=d]=[d(i,j)=d]=[(i,j)=1]

对于每个素数求解的复杂度是O\Big(\sum\limits_{p\in Prime}^{n}{\sqrt{n/p}}\Big)

将素数视为平均分布,可估计额\frac{1}{\log n}\sum\limits_{p=1}^{n}\sqrt{1/p}

接下来需要一点微积分知识 : \sum\limits_{p=1}^{n}\sqrt{1/p}=\sqrt{n}

总的复杂度就是O(n/\log n)

对于不会计算复杂度的同学,可以暴力统计一下式子的值,看看能不能过就好了。

Code:

#include<iostream>
#include<cstdio>
#define MaxNum 10000010
using namespace std;
long long p[MaxNum/10],tn,mu[MaxNum+50],anss,ttn,n,aaa;
bool e[MaxNum+50];
void getmu()
{
  e[1]=1;mu[1]=1;
  for (long long i=2;i<=MaxNum;i++){
    if (!e[i])mu[p[++tn]=i]=-1;
    for (int j=1;j<=tn&&i*p[j]<MaxNum;j++){
      e[i*p[j]]=1;
      mu[i*p[j]]=i%p[j]? -mu[i] : 0;
      if (!i%p[j])break;
    }
  }
} 
//sum(1,N/n)sum(1,N/n)[gcd(i,j)==1]
long long gg(int n){
  long long ans=0;
  for (int l=1,r=0;l<=n;l=r+1){
    r=n/(n/l);          // 分段
    ans+=1ll*(mu[r]-mu[l-1])*(n/l)*(n/l);
  }return ans;
}
int main()
{
  getmu();
  for (int i=2;i<=MaxNum;i++)mu[i]+=mu[i-1];
  cin>>n;
  long long ans=0;
  for (int i=2;i<=n;i++)
   if(!e[i])ans+=gg(n/i);
  cout<<ans;
  return 0;
}

P2257 YY的GCD

貌似和上一题相同呢？发现T<=10000。

上一题的O(n/\log n)还要乘上数据组数,肯定是跑不过去的。

我们将上面的最终式子整理:

明显的,d和p都是k的约数,我们再枚举p,得到d=k/p。（把枚举k的和式放在最前面）

我们发现只有μ(k/p)一项和p有关,所以我们可以把μ(k/p)连着∑\limits_{p∈prime,p|k}一起放到最后面。

得到\sum\limits_{k=1}^{min(n,m)}\lfloor n/k\rfloor\lfloor m/k\rfloor∑\limits_{p∈prime,p|k}μ(k/p)

前面的\sum\limits_{k=1}^{min(n,m)}\lfloor n/k\rfloor\lfloor m/k\rfloor就可以整除分块了。

至于后面的∑\limits_{p∈prime,p|k}μ(k/p)还是套路,维护一个有关于k的前缀和。

怎么弄出前缀和,就是考验数论口胡基本功的时候啦,这里也介绍一下。

设g(k)=∑\limits_{p∈prime,p|k}μ(k/p)

先线筛出\mu,然后枚举p,对于 p|k 的k将g[k]+=\mu(k/p)。(贡献模式,好好理解)

这样的话对于一个素数p,复杂度为O(n/p)

复杂度和埃氏筛相同,为O(n\log\log n)

每个询问O(\sqrt{n}),总的复杂度O(t\sqrt{n}+n\log\log n)

亮出代码:

#include<algorithm>
#include<cstdio>
using namespace std;
int t,n,m,tn,p[1000500],mu[10000500],g[10000500];
bool e[10000500];
void getmu()
{
  e[1]=1;mu[1]=1;
  for (int i=2;i<=10000100;i++){
    if (!e[i]){
      p[++tn]=i;
      mu[i]=-1;
    }for (int j=1;j<=tn;j++){
      if (1ll*i*p[j]>10000100)break;
      e[i*p[j]]=1;
      mu[i*p[j]]=i%p[j] ? -mu[i] : 0;
      if (!i%p[j])break;
    }
  }//线筛出mu 
}
long long calc(int n,int m)
{
  long long ans=0;
  int l=1,r;
  for (;l<=min(n,m);l=r+1){
    r=min(n/(n/l),m/(m/l));//整除分块 
    ans+=1ll*(n/l)*(m/l)*(g[r]-g[l-1]);
  }return ans;
}
int main()
{
scanf("%d",&t);
getmu();
for (int i=1;i<=tn;i++)
 for (int j=p[i];j<=10000100;j+=p[i])
  g[j]+=mu[j/p[i]];
for (int i=2;i<=10000100;i++)g[i]+=g[i-1];
//预处理前缀和 
while(t--){

  scanf("%d%d",&n,&m);
  printf("%lld\n",calc(n,m)); 

}return 0;}

P5221 Product & My Sol

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

6.探索新知:其他数论函数

想必经过上面的学习，大家对\mu已经有了一定的了解，运用反演定理也更熟练了些，那么我们一起向数论函数这个广阔的天地进发吧!

设(n,m)=1,

得\varphi(n)\varphi(m)=\sum\limits_{i=1}^n[(n,i)=1]\sum\limits_{j=1}^m[(m,j)=1]

即\sum\limits_{i=1}^n\sum\limits_{j=1}^m[(m,j)=1\ {\bf and}\ (n,i)=1]

由于互质,[(m,j)=1\ {\bf and}\ (n,i)=1]等价于[(nm,im+jn)=1]。

因为\gcd\Big((im+jn)\bmod n,n\Big)=\gcd(im,n)=1,且\gcd\Big((im+jn)\bmod m,m\Big)=\gcd(jn,m)=1

我们还要证明所有的(in+jm)\bmod {nm}恰好组成[1...nm]的集合。

首先,这些二元组的大小为nm,和目标集合大小相同,我们只需要证明这些元素两两不同即可。

前置芝士 : ca=cb\pmod p且(c,p)=1,则有a=b. (可见同余系基本定理1)

• 采用反证法,假如有am+bn=cm+dn\pmod {nm}

  我们将两边对m取模得到bn=dn\pmod {m},由引理得b=d,同理有a=c,矛盾。

那么(in+jm)\bmod {nm}的集合与[1...nm]的集合一一对应,且贡献模式相同。

于是上式等价于\sum\limits_{i=1}^{nm}[(nm,i)=1]=\varphi(nm)

如何求积性函数的前n个值呢?

线性筛!

对于线性筛,我们要解决\varphi(prime^k)类的函数值。(下面的p均为素数)

则np=p^{(k+1)}*m,那么\varphi(np)=\varphi(m)*\varphi(p^{(k+1)})=\varphi(m)*(p-1)p^k

原命题明显成立。

Code:

bitset<MaxN> e;
int p[MaxN/10],tn;
long long ans,phi[MaxN];
void getphi()
{
  phi[1]=1;
  for (int i=2;i<=n;i++){
    if (!e[i]){
      p[++tn]=i;
      phi[i]=i-1;
    }for (int j=1,t;j<=tn&&(t=i*p[j])<=n;j++){
      e[t]=1;
      phi[t]=phi[i]*(i%p[j]?p[j]-1:p[j]);
      //此处注意
      if (i%p[j]==0)break;
    }
  }
}

• 代码里"此处注意":

  如果(i,p)=1,则\varphi(ip)=\varphi(i)*\varphi(p)=(p-1)\varphi(i)

  如果p|i,根据上文的定理,\varphi(ip)=\varphi(i)*p

• 附: 有些时候我们需要求出单个的\varphi(n),如果按照定义大力求,至少得O(n)次\gcd。

  我们把n分解得p_1^{k_1}*p_2^{k_2}*...*p_m^{k_m}

  那么\varphi(n)=\varphi(p_1^{k_1})*\varphi(p_2^{k_2})*...*\varphi(p_m^{k_m})=\dfrac{(p_1-1)p_1^k}{p_1}*\dfrac{(p_2-1)p_2^k}{p_2}*...\dfrac{(p_m-1)p_m^k}{p_m}

  =(p_1^kp_2^{k_2}...p_m^{k_m})*\left(\dfrac{p_1-1}{p_1}*\dfrac{p_2-1}{p_2}*...*\dfrac{p_m-1}{p_m}\right)=n\left(\dfrac{p_1-1}{p_1}*\dfrac{p_2-1}{p_2}*...*\dfrac{p_m-1}{p_m}\right) \color{blue}\text{性质:}$ $\varphi(n=p_1^kp_2^{k_2}...p_m^{k_m})=n\left(\dfrac{p_1-1}{p_1}*\dfrac{p_2-1}{p_2}*...*\dfrac{p_m-1}{p_m}\right)

  根据上式,把n质因数分解即可求出\varphi(n),复杂度O(\sqrt{n})

  Code:

int phi(int n)
{
  int ans=n;
  for (int i=2;i*i<=n;i++)
   if (n%i==0){
     ans=ans/i*(i-1);
     //根据上文式子计算
     while(n%i==0)n/=i;
   }
  if (n>1)ans=ans/n*(n-1);
  //如果分解不尽,那么肯定还剩一个素数
  return ans;
}

我们来看看\varphi的经典应用:

题目

相信你能直接想到一个莫比乌斯反演的做法:

• 设f(n)是\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=1]

  考虑枚举gcd那么答案是\sum\limits_{d=1}^nd\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=d]

  后面的\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=d]=\sum\limits_{i=1}^{(n/d)}\sum\limits_{j=1}^{(n/d)}[(i,j)=1]=f(n/d)

  答案化为\sum\limits_{d=1}^nd*f(n/d)。

  对f(n)反演,……即为∑\limits_{d=1}^{n}μ(d)(n/d)^2,用整除分块求一次f时间复杂度O(\sqrt{n})

  按\sum\limits_{d=1}^nd*f(n/d)计算即可O(n\sqrt{n})。

  考虑对f(n/d)的n/d整除分块,即可O(n)-O(n^{3/4})(方法参照上文,证明需用积分)。

其实利用\varphi可以很方便地解决上述问题,复杂度是O(n)-O(\sqrt{n})。

设f(n)是\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=1]

我们知道\varphi(n)的定义是:小于等于n的数中,与n互质的数的个数。

我们把\varphi(1)+\varphi(2)+...+\varphi(n)加起来,则可以得到每个数与自己更小的数互质的次数和。

即\sum\limits_{i=1}^n\sum\limits_{j=1}^i[(i,j)=1]=\sum\limits_{i=1}^n\varphi(i)

现在，每个数只能统计比自己小的互质产生贡献，也就是，只统计了数对(x,y)中x>=y的部分。

我们考虑(x,y)中x<=y的部分贡献,显然和上式相同。

（想像一个邻接矩阵）

乘以2就好了，但是由于(1,1)=1会多算一次,所以答案减1。

我们设S(n)=\sum\limits_{i=1}^n\varphi(i)(前缀和),上式变为f(n)=2*s(n)-1

答案变为\sum\limits_{d=1}^nd*f(n/d)=\sum\limits_{d=1}^nd*(2*S(n/d)-1)已经可以做到O(n)了，再整除分快一次就可以O(\sqrt{n})了。

然后使用线性筛就得到O(n)-O(\sqrt{n})的算法了。

Code:

#include<cstdio>
#include<bitset>
#define MaxN 10000500
using namespace std;
int n;
bitset<MaxN> e;
int p[MaxN/10],tn;
long long ans,phi[MaxN];
void getphi()
{
  phi[1]=1;
  for (int i=2;i<=n;i++){
    if (!e[i]){
      p[++tn]=i;
      phi[i]=i-1;
    }for (int j=1,t;j<=tn&&(t=i*p[j])<=n;j++){
      e[t]=1;
      phi[t]=phi[i]*(i%p[j]?p[j]-1:p[j]);
      if (i%p[j]==0)break;
    }
  }
}
int main()
{
  scanf("%d",&n);
  getphi();
  for (int i=1;i<=n;i++)phi[i]+=phi[i-1];
  for (int i=1;i<=n;i++)ans+=(phi[n/i]*2-1)*i;
  printf("%lld",ans);
}

这个解法相对莫比乌斯反演有什么好处呢?

代码短啊!

如果多组询问的话,就可以使用数论分块O(\sqrt{n})回答啦。

重申经典结论:\sum\limits_{i=1}^n\sum\limits_{j=1}^n[(i,j)=1]=\left(\sum\limits_{i=1}^n\varphi(i)*2\right)-1

既然欧拉函数这么好用,为啥还要反演呢?

求\sum\limits_{i=1}^n\sum\limits_{j=1}^m(i,j)? 欧拉函数就没法直接按照定义做了。

这么说太过于人类智慧了,下面,我们考虑使用狄利克雷卷积来系统分析。

2)一些狄利克雷卷积结论:

• id=I*\varphi

即\sum\limits_{d|n}\varphi(d)=n

这个东西是比较重要的,相应地,证明也很长。

我们知道F是个积性函数,那么我们只要证明所有的F(prime^k)=prime^k成立。

然后根据积性和唯一分解定理,把结论扩展到全体正整数。

根据\varphi(p^k)=(p-1)p^{(k-1)},\varphi(1)=1

证毕。

根据I^{-1}=\mu,可得id=I*\varphi等价于\mu*id=\varphi

即\sum\limits_{d|n}\mu(d)(n/d)=\varphi(n)

这个结论极其有用,比如说我们求\sum\limits_{i=1}^n\sum\limits_{j=1}^m(i,j)

利用莫比乌斯反演,变形为\sum\limits_{d=1}^nd*∑\limits_{t=1}^{n/d}μ(t)\lfloor n/dt\rfloor\lfloor m/dt\rfloor

设p=dt,上式即为

可以交换求和号,得到

发现中间的\sum\limits_{d|p}d*μ(p/d)就是\mu*id所以等于\varphi(p)

(这个式子是可以整除分块的)

这样就可以解决上文留下的问题了。

3)除数函数相关

• σ=id*I

把n分解得到p_1^{k_1}*p_2^{k_2}*...*p_m^{k_m}

有\mu(ij)=\mu(i)\mu(j)[i\perp j]

这个容易,假如(i,j)大于1,那么一定会导致出现平方因子,式子的值为0.

否则按照积性分解即可。

• 除数函数d(n)

有d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[x\perp y]

比较难,只会从右边推到左边。

把i分解得到p_1^{k_1}*p_2^{k_2}*...*p_m^{k_m},j分解得到p_1^{k_1'}*p_2^{k_2'}*...*p_m^{k_m'}

显然,每个质因子是独立的,我们考虑i=p^a;j=p^b的情况。

大眼观察可得只有 x'取0...a,且y'取0这a+1个,

有\varphi(ij)=\dfrac{\varphi(i)\varphi(j)(i,j)}{\varphi((i,j))}

默认p指质数。

考虑到[p|ij]=[p|i\ or\ p|j]=[p|i]+[p|j]-[p|(i,j)]