Processing math: 93%

算法竞赛专题解析(4):杜教筛--以及积性函数的前世今生

本系列是这本算法教材的扩展:《算法竞赛入门到进阶》(


PDF下载地址:https://github.com/luoyongjun999/code 其中的“补充资料”
如有建议,请联系:(1)QQ 群,567554289;(2)作者QQ,15512356

  本文详细讲解了与杜教筛有关的各种知识,包括积性函数、欧拉函数、整除分块、狄利克雷卷积、莫比乌斯反演、杜教筛等。以及它们的:理论知识、典型例题、关键代码、练习题。
  本文的目标是:让一名数论小白,也能最终看懂和应用杜教筛。

0 杜教筛简介

0.1 杜教筛的核心内容

  杜教筛用途:在低于线性时间里,高效率求一些积性函数的前缀和。
  杜教筛算法 = 整除分块 + 狄利克雷卷积 + 线性筛。
  杜教筛公式g(1)S(n)=ni=1h(i)ni=2g(i)S(ni)

0.2 杜教筛之起源

  “杜教筛”,一个亲切的算法名字,从2016年后流行于中国的OI圈。
   杜教,即信息学竞赛的大佬杜瑜皓。本文作者好奇“杜教筛”这个名称的来源,在QQ上联系到他,他说,虽然算法的原理不是他的发明,但确实是他最早应用在OI竞赛中。至于为什么被称为“杜教筛”并广为传播,可能是一个迷了。
  这个算法的原始出处,据skywalkert(唐靖哲,唐老师)的博客说,“这种黑科技大概起源于Project Euler这个网站,由xudyh(杜瑜皓)引入中国的OI、ACM界,目前出现了一些OI模拟题、OJ月赛题、ACM赛题是需要这种技巧在低于线性时间的复杂度下解决一类积性函数的前缀和问题。”
  Project Euler网站也是一个OJ,题库里面有很多高难度的数学题。不过这个OJ并不需要提交代码,只需要提交答案即可。这个网站被推崇的一大原因是:如果提交的答案正确,就有权限看别人对这一题的解题方法,以及上传的代码。
  经典的杜教筛题目,例如洛谷P4213,求数论函数的前缀和。

给定一个正整数nn2311,求: $$ans1=\sum_{i=1}n\phi(i),ans2=\sum_{i=1}n\mu(i)$$

  题目中的ϕ(i)是欧拉函数,μ(i)是莫比乌斯函数。欧拉函数、莫比乌斯函数在算法竞赛中很常见,它们都是积性函数。题目求这2个积性函数的前缀和,由于给的n很大,显然不能直接用暴力的、复杂度O(n)的线性方法算,需要用杜教筛,复杂度是O(n23)
  读者看到O(n23),可能感觉和O(n)相差不大,但实际上是一个很大的优化,两者相差n13倍。在上面的例题中,n=21亿,n23=160万,两者相差1300倍。
  关于“杜教筛”的理论知识介绍,影响较大的有:
  (1)2016年信息学奥林匹克中国国家队候选队员论文集,“积性函数求和的几种方法 任之洲”,其中提到杜教筛:“这个算法在国内信息学竞赛中一般称为杜教筛,可以高效地完成一些常见数论函数的计算。”
  (2)2016年,skywalkert(唐老师)的博客

0.3 本文的结构

  由于杜教筛涉及的知识比较多,不容易讲清楚。所以本文将循循善诱地完成这一困难的任务,本文的顺序是:
  (1)积性函数的定义和一些性质。积性函数的常见计算。
  (2)以欧拉函数为例,讲解积性函数的计算。欧拉函数有什么用、如何求、它的求解和积性函数有什么关系、为什么用到杜教筛。
  (3)整数分块。
  (4)狄利克雷卷积。
  (5)莫比乌斯函数和莫比乌斯反演。
  (6)杜教筛。杜教筛公式、复杂度证明、欧拉函数和莫比乌斯函数的杜教筛实现。

1 积性函数

  如果读者没有深入学过数论,第一次看到这些词语:积性函数、欧拉函数、狄利克雷卷积、莫比乌斯函数、莫比乌斯反演、杜教筛,是不是感到高深莫测、头皮发怵?
  即使读者学过一些数论,但是还没有系统读过初等数论的书,那么在阅读本文之前,最好也读一本,比如《初等数论及其应用》Kenneth H.Rosen著,夏鸿刚译,机械工业出版社。这本书很易读,非常适合学计算机的人阅读:(1)概览并证明了初等数论的理论知识;(2)理论知识的各种应用,可以直接用在算法题目里;(3)大量例题和习题;(4)与计算机算法编程有很多结合。花两天时间通读很有益处。
  本文所有的数学定义都引用自这本书。

1.1 积性函数定义

  定义[《初等数论及其应用》第7章乘性函数]

  (1)定义在所有正整数上的函数称为算术函数(或数论函数)。
  (2)如果算术函数f对任意两个互质(互素)的正整数pq,均有f(pq)=f(p)f(q),称为积性函数(multiplicative functions,或译为乘性函数)[有乘性函数,也有加性函数,参考《初等数论及其应用》182页]

  (3)如果对任意两个正整数pq,均有f(pq)=f(p)f(q),称为完全积性函数。
  性质:
【定理1.1】 积性函数的和函数也是积性函数。如果f是积性函数,那么f的和函数F(n)=d|nf(d)也是积性函数。
  其中d|n表示d整除n
  和函数在杜教筛中起到了关键的作用

1.2积性函数的基本问题

  (1)计算积性函数f的第nf(n)
  (2)计算积性函数f1n范围内所有项:f(1)f(2)f(n)
  (3)计算积性函数fn项的和:ni=1f(i),即前缀和。这就是杜教筛的目标。
  后面以欧拉函数和莫比乌斯函数为例,解答了这几个问题。

1.3常见积性函数

  下面的积性函数,在狄利克雷卷积中常常用到。
  I(n):恒等函数,I(n)=1
  id(n):单位函数,id(n)=n
  Ik(n):幂函数,Ik(n)=nk
  ϵ(n):元函数,ϵ(n)={1,n=10,n>1
  σ(n):因子和函数,σ(n)=d|nd
  d(n):约数个数,d(n)=d|n1

2 欧拉函数

  为了彻底了解积性函数的运算,本节以欧拉函数为例进行讲解,这是一个经典的积性函数。并用这个例子引导出为什么要使用杜教筛。

2.1 欧拉函数的定义和性质

  定义:设n是一个正整数,欧拉函数ϕ(n)定义为不超过n且与n互质的正整数的个数。
  定义用公式表示:ϕ(n)=ni=1[gcd(i,n)=1]
  例如:
  n = 4时,1、3这2个数与4互质,ϕ(4)=2
  n = 9时,1、2、4、5、7、8这6个数与9互质,ϕ(9)=6
【定理2.1】 [有关欧拉函数的所有证明,见《初等数论及其应用》176-180页]: 设p和q是互质的正整数,那么ϕ(pq)=ϕ(p)ϕ(q)

  这个定理说明欧拉函数是一个积性函数。有以下推理:
  若n=pa11×pa22××pass,其中p1p2ps互质,a1a2as是它们的幂,则ϕ(n)=ϕ(pa11)×ϕ(pa22)××ϕ(pass))
  如果p1p2...ps互质,那么pa11pa22pass也是互质的。
  以n=84为例:
  84=22×3×7
  ϕ(84)=ϕ(22×3×7)=ϕ(22)×ϕ(3)×ϕ(7))
【定理2.2】n为正整数,那么n=d|nϕ(d)
  其中d|n表示d整除n,上式对n的正因子求和。 例如n=12,那么d1234612,有:12=ϕ(1)+ϕ(2)+ϕ(3)+ϕ(4)+ϕ(6)+ϕ(12)
  这个定理说明了n与的ϕ(n)关系:n的正因子(包括1和自身)的欧拉函数之和等于n
  有很多方法可以证明,这里给出一种直接易懂的证明。
  证明:分数1n,2n,,nn,共有n个,且互不相等。化简这些分数,得到新的n个分数,它们的分母和分子互质,形如add|nad互质。在所有n个分数中,分母为d的分数,数量是ϕ(d)。所有不同分母的分数,其总数是d|nϕ(d),所以n=d|nϕ(d),得证。
  定理2.2是欧拉函数特别重要的一个性质,这一性质被用到了杜教筛中。

2.2求欧拉函数的通解公式

  欧拉函数有什么用呢?利用它可以推出欧拉定理[欧拉定理:设m是一个正整数,a是一个整数且(a,m)=1,那么aϕ(m)1(mod m)],欧拉定理可用于:求大数的模、求解线性同余方程等,在密码学中有重要应用。所以求ϕ(n)是一个常见的计算。
  从定理2.1可以推出定理2.3,定理2.3是计算ϕ(n)的通解公式。
【定理2.3】 设:n=pa11×pa22××pass为正整数n的质幂因子分解,那么:

ϕ(n)=n(11p1)(11p2)(11pk)=nki=1(11pi)

  上述公式,有两个特殊情况:
  (1)若n是质数,ϕ(n)=n1。这一点很容易理解,请读者思考。
  (2)若n=pkp是质数,有:ϕ(n)=ϕ(pk)=pkpk1=pk1(p1)=pk1ϕ(p)
  这两种情况将用于2.3.2节编程求欧拉函数。
  所以,欧拉函数ϕ(n)的求解,归结到了分解质因子这个问题。分解质因子有一个古老的方法,试除法[试除法很低效,有很多快速分解质因子的方法,参考《初等数论及其应用》93页]:求n的质因子时,逐个检查从2n的所有质数,如果它能整除n,就是一个因子。试除法的复杂度是O(n),效率很低。在密码学中,有高达百位以上的十进制数分解质因子,因此发明了很多高效率的方法。不过,在算法竞赛中,数据规模不大,所以一般就用试除法。
  下面是求欧拉函数的代码。先用试除法分解质因子,再用公式ϕ(n)=nki=1(11pi)求得ϕ(n)。总复杂度仍然是O(n)

求单个欧拉函数的代码
int euler(int n){
    int ans = n;
    for(int p = 2; p*p <= n; ++ p){ //试除法:检查从2到sqrt(n)的每个数
        if(n%p == 0){               //能整除,p是一个因子,而且是质因子,请思考
            ans = ans/p*(p-1);      //求欧拉函数的通式
            while(n%p == 0)         //去掉这个因子的幂,并使得下一个p是质因子
                n /= p;             //减小了n
        }
    }
    if(n != 1)                      //情况(1):n是一个质数,没有执行上面的分解
        ans = ans/n*(n-1);
    return ans;
}

2.3 线性筛(欧拉筛)求1~n内所有的欧拉函数

  现在要求1n范围内所有的欧拉函数,而不是只求一个欧拉函数。前面求一个欧拉函数的复杂度是O(n),如果一个一个地求1n范围内所有的欧拉函数,那么n个的总复杂度是O(nn)。效率太低。
  这个过程可以优化:
  (1)利用积性函数的性质ϕ(p1p2)=ϕ(p1)ϕ(p2),求得前面的ϕ(p1)ϕ(p2)后,可以递推出后面的ϕ(p1p2)
  (2)求1n范围内的质数,用于上一步骤的递推。
  编程时,可以这样操作:
  (1)求得1n内每个数的最小质因子。所谓最小质因子,例如84=22×3×784的最小质因子是2。这步操作可以用欧拉筛,复杂度是O(n)的,不可能更快了。
  (2)在第(1)步基础上,递推求1n内每个数的欧拉函数。这一步的复杂度也是O(n)的。
  两者相加,总复杂度是O(n)的。

2.3.1 欧拉筛

  欧拉筛是一种线性筛,也就是说,它能在O(n)的线性时间里求得1n内所有的质数。
  欧拉筛是对埃氏筛法的改进。埃氏筛法的原理是:每个合数都是多个质数的积;那么从最小的质数2开始,用每一个质数去筛比它大的数,就能筛掉合数。埃氏筛法低效的原因是,一个合数会被它的多个质因子重复筛。请读者自己详细了解埃氏筛法。
  欧拉筛法的原理是:一个合数肯定有一个最小质因子;让每个合数只被它的最小质因子筛选一次,以达到不重复筛的目的。
  具体操作步骤:
  (1)逐一检查2n的所有数。第一个检查的是2,它是第一个质数。
  (2)当检查到第i个数时,利用已经求得的质数去筛掉对应的合数x,而且是用x的最小质因子去筛。
  下面是代码。代码很短很精妙。

欧拉筛求素数
int prime[MAXN];              //保存质数
int vis[MAXN];                //记录是否被筛
int euler_sieve(int n){       //欧拉筛。返回质数的个数。
	int cnt = 0;              //记录质数个数
	memset(vis,0,sizeof(vis));
	memset(prime,0,sizeof(prime));
	for(int i=2;i<=n;i++){             //检查每个数,筛去其中的合数
		if(!vis[i]) prime[cnt++]=i;    //如果没有筛过,是质数,记录。第一个质数是2
		for(int j=0; j<cnt; j++){      //用已经得到的质数去筛后面的数
			if(i*prime[j] >n)  break;  //只筛小于等于n的数
			vis[i*prime[j]]=1;         //关键1。用x的最小质因数筛去x
			if(i%prime[j]==0)  break;  //关键2。如果不是这个数的最小质因子,结束
		}
	}
	return cnt;                        //返回小于等于n的素数的个数
}

  代码中难懂的是后面2个关键行。在“关键1”, 其中的prime[j]是最小质因子,vis[iprime[j]]=1;表示用最小质因子筛去了它的倍数。在“关键2”,及时跳出,避免重复筛。
  下面用图解来说明,请读者对照代码和图解自己理解。特别注意图4,它是“关键2”。

图1. 初始化
图2. i=2,用质数2筛去4(2是4的最小质因子)
图3. i=3,用质数2筛去6,用质数3筛去9
图4. i=4,用质数2筛去8,注意质数3没有用到,循环j被跳出了
图5. i=5,用质数2筛去10,用3筛去15,用5筛去25

  可以发现,每个数x只被筛了一次,而且是被x的最小质因子筛去的。这说明筛法是线性的,复杂度 O(n)
  现在用欧拉筛求1n内每个数的最小质因子。只要简单修改上述程序,直接用vis[MAXN]记录最小质因子就行。代码修改如下,修改了带注释的后2行。

欧拉筛求质数+最小质因子
int prime[MAXN];      //记录质数
int vis[MAXN];      //记录最小质因子
int euler_sieve(int n){     
	int cnt=0;
	memset(vis,0,sizeof(vis));
	memset(prime,0,sizeof(prime));
	for(int i=2;i<=n;i++){               
		if(!vis[i]){ vis[i]=i; prime[cnt++]=i;}    //vis[]记录x的最小质因子
		for(int j=0; j<cnt; j++){       
			if(i*prime[j] >n)  break;       
			vis[i*prime[j]]=prime[j];              //vis[]记录x的最小质因子
			if(i%prime[j]==0)
                break;                  
		}
	}
	return cnt;
}

2.3.2 求1~n内所有的欧拉函数

  下面用一个模板题给出模板代码,它就是“欧拉筛+欧拉函数公式”。
  例题:洛谷P2158 仪仗队 https://www.luogu.com.cn/problem/P2158
  这是欧拉函数的模板题。get_phi()是计算欧拉函数的模板,其中判断了3种情况,细节见注释。
  情况(1):若n是质数,ϕ(n)=n1
  情况(2):若n=pkp是质数,有ϕ(n)=pk1ϕ(p)
  情况(3):若n=p1p2,用ϕ(n)=ϕ(p1)ϕ(p1)递推。

洛谷P2158代码
#include<bits/stdc++.h>
using namespace std;

const int maxn = 50000;
int vis[maxn];   //记录是否被筛;或者用于记录最小质因子
int prime[maxn]; //记录质数
int phi[maxn];   //记录欧拉函数
int sum[maxn];   //计算欧拉函数的和

void get_phi(){  //模板:求1~maxn范围内的欧拉函数
    phi[1]=1;
    int cnt=0;
    for(int i=2;i<maxn;i++) {
        if(!vis[i]) {
            vis[i]=i; //vis[i]=1; //二选一:前者记录最小质因子,后者记录是否被筛
            prime[cnt++]=i;       //记录质数
            phi[i]=i-1;           //情况(1):i是质数,它欧拉函数值=i-1
        }
        for(int j=0;j<cnt;j++) {
            if(i*prime[j] > maxn)  break;
            vis[i*prime[j]] = prime[j]; //vis[i*prime[j]]=1; 
                                 //二选一:前者记录最小质因子,后者记录是否被筛
            if(i%prime[j]==0){          //prime[j]是最小质因子
                phi[i*prime[j]] = phi[i]*prime[j];
                                 //情况(2):i是prime[j]的k次方
                break;
            }
            phi[i*prime[j]] = phi[i]*phi[prime[j]];  
                                //情况(3):i和prime[j]互质,递推出i*prime[j]
        }
    }
}

int main(){
    get_phi();        //计算所有的欧拉函数
    sum[1]=1;
    for(int i=2;i<=maxn;i++)   //打表计算欧拉函数的和
        sum[i]=sum[i-1]+phi[i];
    int n;   scanf("%d",&n);
    if(n==1) printf("0\n");
    else     printf("%d\n",2*sum[n-1]+1);
    return 0;
}

  本节的例子是用欧拉筛求欧拉函数,读者可以认识到:积性函数都可以用欧拉筛来求。后面还有一个例子:用欧拉筛求解积性函数莫比乌斯函数。

2.4 欧拉函数前缀和

  在欧拉函数问题的最后,回到杜教筛的模板题,求:ni=1ϕ(i)
  简单的方法,是用2.3.2节计算出每个欧拉函数,再求和,复杂度是O(n)
  更快的方法,就是本篇的主题“杜教筛”,复杂度是O(n23)
  “杜教筛”需要整除分块、狄利克雷卷积等前置知识。

3 整除分块(数论分块)

  整除分块是杜教筛算法的基本思路。
  整除分块是为了解决一个整除的求和问题:ni=1ni
  其中ni表示n除以i,并向下取整[除法运算,向上取整,英文Ceiling,用符号 表示;向下取整,英文Floor,用符号 表示]n是给定的一个整数。

  如果直接暴力计算,复杂度是O(n)的,若n很大会超时。但是,用整除分块算法来求解,复杂度是O(n)的。
  下面以n=8为例,列出每个情况:

i 1 2 3 4 5 6 7 8
8i 8 4 2 2 1 1 1 1

  对第2行求和,结果是:8i=18i=20
  观察上面第28i的值,发现是逐步变小的,并且很多相等,这是整除操作的规律。例如:
  8/1=8
  8/2=4
  8/3=8/4=2
  8/5=8/6=8/7=8/8=1
  为了对这些整除的结果进行快速求和,自然能想到,把它们分成一块块的,其中每一块的8i相同,把这一块一起算,就快多了。
  设每一块的左右区间是[l,r],上面的例子可以分成4块:[1,1][2,2][3,4][5,8]
  每一块内部的求和计算,只要用数量乘以值就行了,例如[5,8]区间的整除和:(85+1)1=4。这个计算的复杂度是O(1)的。
  最后还有两个问题:
  (1)把ni按相同值分块,一共需要分成几块?或者说,ni有几种取值?这决定了算法的复杂度。
  【答】分块少于2n种,证明如下:
  (a)in时,ni的值有:{n1,n2,n3nn}nin,共n个,此时nin种取值;
  (b)i>n时,有ni<n,此时ni也有n种取值。
两者相加,共2n种。所以,整除分块的数量是O(n)的。
  (2)如何计算?或者说,给定每个块的第一个数l,能推出这个块的最后一个数r吗?
  【答】可以证明:r=n/(n/l)。证明略[证明细节]

  下面是标程。

#include<bits/stdc++.h>
using namespace std;
int main(){
    long long n,ans=0;
    cin >> n;
    for(long long l=1,r;l<=n;l=r+1){
        r = n/(n/l);            //计算r,让分块右移
        ans += (r-l+1)*(n/l);   //求和
        cout << l <<""<< r <<": "<< n/r << endl;  //打印分块
    }
    cout << ans;               //打印和
}

  整除分块习题:洛谷P1829、洛谷P2261、洛谷P3935。

4 狄利克雷卷积

  狄利克雷卷积 Dirichlet Convolution,是计算求和问题的有用工具。狄利克雷卷积是杜教筛用到的核心技术之一。
  定义:
  设fg是算术函数,记fg的狄利克雷卷积是fg,定义为[参考]
(fg)(n)=d|nf(d)g(nd)

  它是一个求和公式,其中d|n表示d整除n,卷积是对正因子求和。
  这个卷积是什么含义呢?现在给出一个特殊的例子。定义恒等函数I(n)=n、常数函数1(n)=1,它们的卷积是:

(I1)(n)=d|nI(d)1(nd)=d|nd1=d|nd=σ(n)

  把它记为:I1=σ
  σ(n)是“因子和函数”的符号。例如n=12,那么σ(n)=1+2+3+4+6+12=28σ(n)是积性函数。
  狄利克雷卷积的性质:
  (1)交换律:fg=gi
  (2)结合律:(fg)h=f(g)h
  (3)分配律:f(g+h)=(fg)+(fh)
  (4)元函数:ϵ(n)={1,n=10,n>1。对任何f,有ef=fe=f
  (5)两个积性函数的狄利克雷卷积仍然是积性函数。
  有一些积性函数算术函数在狄利克雷卷积中常常用到,1.3节已经给出,这里再次列出:
  I(n):恒等函数,I(n)=1
  id(n):单位函数,id(n)=n
  Ik(n):幂函数,Ik(n)=nk
  ϵ(n):元函数,ϵ(n)={1,n=10,n>1
  σ(n):因子和函数,σ(n)=d|nd
  d(n):约数个数,d(n)=d|n1

5 莫比乌斯函数和莫比乌斯反演

  如果读者困惑莫比乌斯函数有什么用,可以先看下一节“5.2 莫比乌斯函数的由来”。

5.1 莫比乌斯函数的定义和性质

  莫比乌斯函数以数学家Möbius的名字命名[https://brilliant.org/wiki/mobius-function/]

  莫比乌斯函数是狄利克雷卷积的重要工具,它也是数论和组合数学的重要的积性函数。
  莫比乌斯函数μ(n)定义为[《初等数论及其应用》200页]

μ(n)={1,如果n=1(1)r,如果n=p1p2...prpi0,其他情况

  一些例子:μ(1)=1,μ(2)=1,μ(3)=1,μ(4)=μ(22)=0,μ(5)=1,μ(6)=μ(2×3)=1,μ(7)=1,μ(8)=μ(23)=0;μ(330)=μ(2×3×5×11)=(1)4=1,μ(660)=μ(22×3×5×11)=0
  莫比乌斯函数μ(n)是积性函数。它有一个和欧拉函数的n=d|nϕ(d)类似的定理:
【定理5.1】莫比乌斯函数的和函数在整数n处的值,满足:
d|nμ(d)={1,n=1 0,n>1
  证明:
  (1)n=1时,显然有F(1)=d|nμ(1)=1
  (2)n>1时。根据积性函数的定义,有F(n)=F(Pa11)F(Pa22)F(Patt),其中n=Pa11Pa22Patt是质因子分解。如果能证明F(Pk)=0即有F(n)=0。因为当i2时,μ(pi)=0,有:

F(pk)=d|pkμ(d)=μ(1)+μ(p)+μ(p2)++μ(pk)=1+(1)+0++0=0

  下面是莫比乌斯函数线性筛代码,求1n内的莫比乌斯函数,和欧拉函数线性筛的代码几乎一样,此处不做解释。

bool vis[maxn];
int prime[maxn];
int Mob[maxn];
void Mobius_sieve(){
     int cnt = 0;
     vis[1] = 1;
     Mob[1] = 1;
     for(int i = 2; i <= maxn; i++){
        if(!vis[i])
            prime[cnt++] = i, Mob[i] = - 1;
        for(int j = 0; j < cnt && 1LL * prime[j] * i <= maxn; j++){
            vis[prime[j] * i] = 1;
            Mob[i * prime[j]] = (i % prime[j] ? -Mob[i]: 0);
            if(i % prime[j] == 0)
                break;
        }
    }
}

5.2 莫比乌斯函数的由来

  看到莫比乌斯函数的定义,读者是不是感到奇怪?它是怎么来的?有什么用?
  本节内容,引用《初等数论及其应用》199页,它给出了莫比乌斯函数的来龙去脉。
  设f为算术函数,f的和函数F的值为F(n)=d|nf(d),显然,Ff的值决定。这种关系可以反过来吗?也就是说,是否存在一种用F求出f的值的简便方法?本节给出这样的公式,看什么样的公式可行。
  若f是算术函数,F是它的和函数F(n)=d|nf(d),按定义展开F(n)n=1,2,...,8,有:
  F(1)=f(1)
  F(2)=f(1)+f(2)
  F(3)=f(1)+f(3)
  F(4)=f(1)+f(2)+f(4)
  F(5)=f(1)+f(5)
  F(6)=f(1)+f(2)+f(3)+f(6)
  F(7)=f(1)+f(7)
  F(8)=f(1)+f(2)+f(4)+f(8)
  根据上面方程,可以解得:
  f(1)=F(1)
  f(2)=F(2)F(1)
  f(3)=F(3)F(1)
  f(4)=F(4)F(2)
  f(5)=F(5)F(1)
  f(6)=F(6)F(3)F(2)+F(1)
  f(7)=F(7)F(1)
  f(8)=F(8)F(4)
  注意到f(n)等于形式为±F(nd)的一些项之和,其中d|n。从这一结果中,可能有这样一个等式,形式是:

f(n)=d|nμ(d)F(nd)

  其中μ是算术函数。如果等式成立,计算得到:
  μ(1)=1,μ(2)=1,μ(3)=1,μ(4)=0,μ(5)=1,μ(6)=1,μ(7)=1,μ(8)=0
  (读者已经发现和前一节莫比乌斯函数的值一样。)
  又F(p)=f(1)+f(p),得f(p)=F(p)F(1),其中p是质数。则μ(p)=1
  又因为F(p2)=f(1)+f(p)+f(p2)
  有:f(p2)=F(p2)(F(p)F(1))F(1)=F(p2)F(p)
  这要求对任意质数p,有μ(p2)=0。类似的推理得出对任意质数p及整数k>1,有μ(pk)=0。如果猜想μ是积性函数,则μ的值就由质数幂处的值决定。
  这就是莫比乌斯函数的由来。

5.3 莫比乌斯反演

【定理5.2】莫比乌斯反演(Möbius inversion)公式。若f是算术函数,Ff的和函数,对任意正整数n,满足F(n)=d|nf(d),则对任意正整数n,有f(n)=d|nμ(d)F(nd)
  从定理5.2可以推出下面的定理5.3。
【定理5.3】f是算术函数,它的和函数为F(n)=d|nf(d),那么如果F是积性函数,则f也是积性函数。
  莫比乌斯反演,不要求f是积性函数。

6 杜教筛

6.1 杜教筛公式

6.1.1 杜教筛公式的推导

  杜教筛要解决的是这一类问题:设f(n)是一个数论函数,计算S(n)=ni=1f(i)
  由于n很大,要求以低于O(n)的复杂度求解,所以用前面提到的线性筛是不够的。如果能用到整除分块,每一块的值相等一起算,就加快了速度,也就是构造形如ni=1ni的整除分块。
  杜教筛的思路,简单地说,是把f(n)构造成能利用整除分块的新的函数,这个构造用到了卷积。
  记S(n)=ni=1f(i)。根据f(n)的性质,构造一个S(n)关于S(ni)的递推式,下面是构造方法[http://fzl123.top/archives/898]

  构造2个积性函数hg,满足卷积:h=gf
  根据卷积的定义,有h(i)=d|ng(d)f(id),求和:

ni=1h(i)=ni=1id|ig(d)f(id)

=nd=1g(d)nd|if(id)

=nd=1g(d)ndi=1f(i)

=nd=1g(d)s(nd)

  注意其中第2行变换ni=1id|ig(d)f(id)=nd=1g(d)nd|if(id),网上昵称“杜教筛变换” [<https://blog.csdn.net/myjs999/article/details/78906549]

  最后一步得到:ni=1h(i)=nd=1g(d)s(nd)
  为计算S(n),把式子右边的第一项提出来:

ni=1h(i)=g(1)S(n)+nd=2g(d)s(nd)

g(1)S(n)=ni=1h(i)nd=2g(d)s(nd)

  式中的id无关,为方便看,改写为:

g(1)S(n)=ni=1h(i)ni=2g(i)s(ni)

  这就是杜教筛公式。

6.1.2 杜教筛算法和复杂度

  一个积性函数求和问题,如果分析得到了杜教筛公式,工作已经完成了一大半,剩下的是编程实现。
  公式g(1)S(n)=ni=1h(i)ni=2g(i)s(ni)S(n)的递归形式,编程时可以递归实现;递归的时候加上记忆化,让每个S[i]只需要算一次。
  其中的h(i)g(i),如果构造得好,就容易计算,见后面欧拉函数和莫比乌斯函数的例子。为方便分析计算复杂度,下面略去h(i)g(i),分析简化后的式子S(n)=ni=2s(ni)的复杂度。
  设S(n)的复杂度是T(n)
  递归展开第一层:S(n)=ni=2s(ni)=s(n2)+s(n3)++s(nn)
  根据整除分块的原理,右边共有O(n)s(ni)。另外,把它们加起来的时候,还有O(n)次合并。总时间是T(n)=O(n)+O(ni=2T(ni))
  再展开一层:T(ni)=O(ni)+O(nik=2T(nik))=O(ni),中间第2项是高阶小量,把它略去。
  合起来:

T(n)=O(n)+O(ni=2T(ni))

=O(n)+O(ni=2ni)

  后一项最大,只考虑它就够了:

T(n)=ni=2O(ni)O(n0nxdx)=O(n34)

  上述计算还能优化,即预处理一部分s(n)s(n)是一个积性函数的前缀和,先用线性筛计算一部分,假设已经预处理了前k个正整数的s(n),且kn。因为S(1)S(k)已经求出,那么这个式子中还没有求出的是k<nin,也就是要算1<ink的这一部分。那么复杂度变为:

T(n)=nki=2O(ni)O(nk0nxdx)=O(nk)

  k取多大合适呢?线性筛计算也是要花时间的,而且是线性的O(k),那么线性筛计算加上杜教筛公式计算,总时间是:

T(n)=O(k)+T(n)=O(k)+O(nk)

  差不多k=n23时,它有最小值,即:O(n23)+O(nn23)=O(n23)
  它比O(n34)要好。
  以上就是杜教筛算法的全部。
  综合起来,杜教筛算法可以概况为:用狄利克雷卷积构造递推式,编程时用整除分块和线性筛优化,算法复杂度达到了O(n23)
  杜教筛算法 = 狄利克雷卷积 + 整除分块 + 线性筛

  应用杜教筛时,狄利克雷卷积是基本的工具。观察卷积的定义(fg)(n)=d|nf(d)g(nd),那么,如果求解的函数有形如d|nf(d))这样的性质,把这个性质与卷积的定义对照,就容易找到gh。下面用欧拉函数和莫比乌斯函数为例说明。算法的具体实现,见6.4节的模板代码。

6.2 欧拉函数前缀和

  求欧拉函数前缀和:ni=1ϕ(i)
  求欧拉函数前缀和的杜教筛公式,需要用到欧拉函数性质:n=d|nϕ(d)
  (1)直接套用杜教筛公式
  按卷积定义h=fg=d|nf(d)g(nd),把n=d|nϕ(d)改写为:id=ϕI,那么有h=idg=I。代入杜教筛公式g(1)S(n)=ni=1h(i)ni=2g(i)s(ni),得:

S(n)=ni=1ini=2s(ni)

=n(n+1)2ni=2s(ni)

  (2)自己推导
  读者如果有兴趣,也可以自己按部就班地以欧拉函数为例,做一次杜教筛公式推导的练习,推导过程是[引用]

n=ϕ(n)+d|n,d<nϕ(d)

ϕ(n)=nd|n,d<nϕ(d)

ni=1ϕ(i)=ni=1(id|i,d<iϕ(d))

ni=1ϕ(i)=n(n+1)2ni=1d|i,d<iϕ(d)

  现在改变枚举变量:枚举id的值,即枚举id的倍数,因为id,所以从2开始:

ni=1ϕ(i)=n(n+1)2idnid=2dnidd=1ϕ(d)

  设S(n)=ni=1ϕ(i),并把id写成i

S(n)=\frac{n(n+1)}{2}-\sum_{i’=2}^ns({\lfloor\frac{n}{i‘}\rfloor})

  这样就得到了欧拉函数的“杜教筛公式”,和(1)一样。

6.3 莫比乌斯函数前缀和

  求莫比乌斯函数前缀和:\sum_{i=1}^n\mu(i)
  求莫比乌斯函数前缀和的杜教筛公式,方法和上一节欧拉函数几乎一样。
  需要用到莫比乌斯函数性质:\sum_{d|n}\mu(d)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases},为方便书写,改写成\sum_{d|n}\mu(d)=[n=1]
  (1)套用杜教筛公式
  按卷积定义h=f*g= \sum_{d|n}f(d)g(\frac{n}{d}),把\sum_{d|n}\mu(d)=[n=1]改写为:\epsilon=\mu*I,那么有h=\epsilon,g=I。代入杜教筛公式g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor}),得:

S(n)=\sum_{i=1}^n\epsilon(i)-\sum_{i=2}^nS({\lfloor\frac{n}{i}\rfloor})

\qquad\qquad=1-\sum_{i=2}^nS({\lfloor\frac{n}{i}\rfloor})

  (2)自己推导[引用]

[n=1]=\mu(n)+\sum_{d|n,d<n}\mu(d)

\mu(n)=[n=1]-\sum_{d|n,d<n}\mu(d)

\sum_{i=1}^n\mu(i)=1-\sum_{i=1}^n\sum_{d|i,d<i}\mu(d)

\sum_{i=1}^n\mu(i)=1- \sum_{\frac id =2}^ {\frac id \leq n} \sum_{d=1}^{d\leq {\lfloor\frac{n}{\frac id}\rfloor} } \sum_{i=1}^n\mu(d)

  令S(n)=\sum_{i=1}^n\mu(i),并把\dfrac id写成i',得:

S(n)=1-\sum_{i‘=2}^nS({\lfloor\frac{n}{i’}\rfloor})

6.4 模板代码

  下面给出洛谷P4213的代码[改写自]

  注释中说明了杜教筛的三个技术:整除分块、线性筛、狄利克雷卷积。其中整除分块、线性筛的代码和前面有关的代码一样。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int maxn = 5e6+7; //超过n^2/3,够大了

int prime[maxn];   //记录质数
bool vis[maxn];    //记录是否被筛;
int mu[maxn];                   //莫比乌斯函数值
ll phi[maxn];                   //欧拉函数值
unordered_map<int,int> summu;   //莫比乌斯函数前缀和
unordered_map<int,ll> sumphi;   //欧拉函数前缀和
inline void init(){             //线性筛预计算一部分答案
    int cnt = 0;
	vis[0] = vis[1] = 1;
	mu[1] = phi[1] = 1;
	for(int i=2;i<maxn;i++){
		if(!vis[i]){
            prime[cnt++] = i;
            mu[i] = -1;
            phi[i] = i-1;
        }
		for(int j=0;j<cnt && i*prime[j]<maxn;j++){
			vis[i*prime[j]] = 1;
			if(i%prime[j]){
				mu[i*prime[j]] = -mu[i];
				phi[i*prime[j]] = phi[i]*phi[prime[j]];
			}
			else{
				mu[i*prime[j]] = 0;
				phi[i*prime[j]] = phi[i]*prime[j];
				break;
			}
		}
	}
	for(int i=1;i<maxn;i++){  //最后,mu[]和phi[]改为记录1~maxn的前缀和。
        mu[i] += mu[i-1];
        phi[i] += phi[i-1];
    }
}
inline int gsum(int x){         // g(i)的前缀和
	return x;
}
inline int getsmu(int x){
	if(x < maxn) return mu[x];   //预计算
	if(summu[x]) return summu[x];  //记忆化
	int ans = 1;                //杜教筛公式中的 1
	for(int l=2,r;l<=x;l=r+1){  //用整除分块计算杜教筛公式
		r = x/(x/l);
		ans -= (gsum(r)-gsum(l-1))*getsmu(x/l);
	}
	return summu[x] = ans/gsum(1);
}
inline ll getsphi(int x){
	if(x < maxn) return phi[x];
	if(sumphi[x]) return sumphi[x];  //记忆化,每个sumphi[x]只用算一次
	ll ans = 1LL*x*(x+1)/2;      //杜教筛公式中的 n(n+1)/2
	for(int l=2,r;l<=x;l=r+1){   //用整除分块计算杜教筛公式,这里算 sqrt(x)次
		r = x/(x/l);
		ans -= (gsum(r)-gsum(l-1))*getsphi(x/l);
	}
	return sumphi[x] = ans/gsum(1);
}
int main(){
	init();  //用线性筛预计算一部分
	int t;
	scanf("%d",&t);
	while(t--){
		int n;
		scanf("%d",&n);
		printf("%lld %d\n",getsphi(n),getsmu(n));
	}
	return 0;
}

6.5 题目

  这里列出网上的一些题目。
  (1)求\sum_{i=1}^ni\phi(i)
  题解:http://fzl123.top/archives/898
  (2)求\sum_{i=1}^n\sum_{j=1}^n\sigma(ij),其中n≤10^9
  题解:https://www.cnblogs.com/zhugezy/p/11312301.html
  (3)hdu 6706
  http://acm.hdu.edu.cn/showproblem.php?pid=6706
  (4)这里列出了一些可提交的题目:
  https://www.cnblogs.com/TSHugh/p/8361040.html
  (4)skywalkert(唐老师)的博客给出了很多习题:
  https://blog.csdn.net/skywalkert/article/details/50500009

致谢

  杜瑜皓:对本文草稿提出细致的修改意见。
  华东理工大学队员:傅志凌(<oj.ecustacm.cn>管理员)、赵李洋、李震、彭博,讨论算法复杂度。

参考文献

[1] 积性函数求和的几种方法,任之洲,2016年信息学奥林匹克中国国家队候选队员论文集
[2] 唐靖哲:https://blog.csdn.net/skywalkert/article/details/50500009
[3] 吉如一:http://jiruyi910387714.is-programmer.com/posts/195270.html
[4] 傅志凌:http://fzl123.top/archives/898

posted on   罗勇军999  阅读(2393)  评论(0编辑  收藏  举报

编辑推荐:
· DeepSeek 解答了困扰我五年的技术问题
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
阅读排行:
· DeepSeek 解答了困扰我五年的技术问题。时代确实变了!
· PPT革命!DeepSeek+Kimi=N小时工作5分钟完成?
· What?废柴, 还在本地部署DeepSeek吗?Are you kidding?
· 赶AI大潮:在VSCode中使用DeepSeek及近百种模型的极简方法
· DeepSeek企业级部署实战指南:从服务器选型到Dify私有化落地

导航

统计

点击右上角即可分享
微信分享提示