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

2020年08月28日 阅读数:35
这篇文章主要向大家介绍杜教筛 以及积性函数的前世此生 --算法竞赛专题解析(4),主要内容包括基础应用、实用技巧、原理机制等方面,希望对大家有所帮助。

本系列文章将于2021年整理出版,书名《算法竞赛专题解析》。
前驱教材:《算法竞赛入门到进阶》 清华大学出版社 2019.8
网购:京东 当当      做者签名书php

若有建议,请加QQ 群:567554289,或联系做者QQ:15512356html


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

0 杜教筛简介

0.1 杜教筛的核心内容

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

0.2 杜教筛之起源

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

给定一个正整数 n n 2 31 1 n,n≤2^{31} −1 ,求: a n s 1 = i = 1 n ϕ ( i ) a n s 2 = i = 1 n μ ( i ) ans1=\sum_{i=1}^n\phi(i),ans2=\sum_{i=1}^n\mu(i)

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

0.3 本文的结构

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

1 积性函数

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

1.1 积性函数定义

  定义1:ide

  (1)定义在全部正整数上的函数称为算术函数(或数论函数)。
  (2)若是算术函数f对任意两个互质(互素)的正整数 p p q q ,均有 f ( p q ) = f ( p ) f ( q ) f(pq)=f(p)f(q) ,称为积性函数(multiplicative functions,或译为乘性函数)2svg

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

1.2积性函数的基本问题

  (1)计算积性函数 f f 的第 n n f ( n ) f(n)
  (2)计算积性函数 f f 1 n 1~n 范围内全部项: f ( 1 ) f ( 2 ) f ( n ) f(1)、f(2)、\cdots 、f(n)
  (3)计算积性函数 f f n n 项的和: i = 1 n f ( i ) \sum_{i=1}^nf(i) ,即前缀和。这就是杜教筛的目标。
  后面以欧拉函数和莫比乌斯函数为例,解答了这几个问题。

1.3常见积性函数

  下面的积性函数,在狄利克雷卷积中经常用到。
   I ( n ) I(n) :恒等函数, I ( n ) = 1 I(n)=1
   i d ( n ) id(n) :单位函数, i d ( n ) = n id(n)=n
   I k ( n ) I_k(n) :幂函数, I k ( n ) = n k I_k(n)=n^k
   ϵ ( n ) \epsilon(n) :元函数, ϵ ( n ) = { 1 , n=1 0 , n>1 \epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}
   σ ( n ) \sigma(n) :因子和函数, σ ( n ) = d n d \sigma(n)=\sum_{d|n}d
   d ( n ) d(n) :约数个数, d ( n ) = d n 1 d(n)=\sum_{d|n}1

2 欧拉函数

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

2.1 欧拉函数的定义和性质

  定义:设n是一个正整数,欧拉函数 ϕ ( n ) \phi(n) 定义为不超过 n n 且与 n n 互质的正整数的个数。
  定义用公式表示: ϕ ( n ) = i = 1 n [ g c d ( i , n ) = 1 ] \phi(n)=\sum_{i=1}^n[gcd(i,n)=1]
  例如:
  n = 4时,一、3这2个数与4互质, ϕ ( 4 ) = 2 \phi(4)=2
  n = 9时,一、二、四、五、七、8这6个数与9互质, ϕ ( 9 ) = 6 \phi(9)=6
【定理2.1】3 设p和q是互质的正整数,那么 ϕ ( p q ) = ϕ ( p ) ϕ ( q ) \phi(pq)=\phi(p)\phi(q)

  这个定理说明欧拉函数是一个积性函数。有如下推理:
  若 n = p 1 a 1 × p 2 a 2 × × p s a s n=p_1^{a_1}\times p_2^{a_2}\times \cdots\times p_s^{a_s} ,其中 p 1 p 2 p s p_一、p_二、\cdots、p_s 互质, a 1 a 2 a s a_一、a_二、\cdots、a_s 是它们的幂,则 ϕ ( n ) = ϕ ( p 1 a 1 ) × ϕ ( p 2 a 2 ) × × ϕ ( p s a s ) ) \phi(n)=\phi(p_1^{a_1}) \times \phi(p_2^{a_2}) \times \cdots\times\phi(p_s^{a_s}) )
  若是 p 1 p 2 . . . p s p_一、p_二、...、p_s 互质,那么 p 1 a 1 p 2 a 2 p s a s p_1^{a_1}、p_2^{a_2}、\cdots、p_s^{a_s} 也是互质的。
  以 n = 84 n = 84 为例:
   84 = 2 2 × 3 × 7 84=2^2 \times 3 \times 7
   ϕ ( 84 ) = ϕ ( 2 2 × 3 × 7 ) = ϕ ( 2 2 ) × ϕ ( 3 ) × ϕ ( 7 ) ) \phi(84)=\phi(2^2 \times 3 \times 7)=\phi(2^2) \times \phi(3) \times \phi(7))
【定理2.2】 n n 为正整数,那么 n = d n ϕ ( d ) n=\sum_{d|n}\phi(d)
  其中 d n d|n 表示 d d 整除 n n ,上式对 n n 的正因子求和。 例如 n = 12 n=12 ,那么 d d 1 2 3 4 6 12 一、二、三、四、六、12 ,有: 12 = ϕ ( 1 ) + ϕ ( 2 ) + ϕ ( 3 ) + ϕ ( 4 ) + ϕ ( 6 ) + ϕ ( 12 ) 12=\phi(1)+\phi(2)+\phi(3)+\phi(4)+\phi(6)+\phi(12)
  这个定理说明了n与的 ϕ ( n ) \phi(n) 关系: n n 的正因子(包括 1 1 和自身)的欧拉函数之和等于 n n
  有不少方法能够证实,这里给出一种直接易懂的证实。
  证实:分数 1 n , 2 n , , n n \frac 1n,\frac 2n,\cdots,\frac nn ,共有 n n 个,且互不相等。化简这些分数,获得新的 n n 个分数,它们的分母和分子互质,形如 a d \frac ad d n d|n a a d d 互质。在全部 n n 个分数中,分母为 d d 的分数,数量是 ϕ ( d ) \phi(d) 。全部不一样分母的分数,其总数是 d n ϕ ( d ) \sum_{d|n}\phi(d) ,因此 n = d n ϕ ( d ) n=\sum_{d|n}\phi(d) ,得证。
  定理2.2是欧拉函数特别重要的一个性质,这一性质被用到了杜教筛中。

2.2求欧拉函数的通解公式

  欧拉函数有什么用呢?利用它能够推出欧拉定理4,欧拉定理可用于:求大数的模、求解线性同余方程等,在密码学中有重要应用。因此求 ϕ ( n ) \phi(n) 是一个常见的计算。

  从定理2.1能够推出定理2.3,定理2.3是计算 ϕ ( n ) \phi(n) 的通解公式。
【定理2.3】 设: n = p 1 a 1 × p 2 a 2 × × p s a s n=p_1^{a_1}\times p_2^{a_2}\times \cdots\times p_s^{a_s} 为正整数 n n 的质幂因子分解,那么:
ϕ ( n ) = n ( 1 1 p 1 ) ( 1 1 p 2 ) ( 1 1 p k ) = n i = 1 k ( 1 1 p i ) \phi(n)=n(1-\frac 1{p_1})(1-\frac 1{p_2})\cdots(1-\frac 1{p_k})=n \prod_{i=1}^k(1-\frac 1{p_i})
  上述公式,有两个特殊状况:
  (1)若 n n 是质数, ϕ ( n ) = n 1 \phi(n)=n-1 。这一点很容易理解,请读者思考。
  (2)若 n = p k n=p^k p p 是质数,有: ϕ ( n ) = ϕ ( p k ) = p k p k 1 = p k 1 ( p 1 ) = p k 1 ϕ ( p ) \phi(n)=\phi(p^k)=p^k-p^{k-1}=p^{k-1}(p-1)=p^{k-1}\phi(p)
  这两种状况将用于2.3.2节编程求欧拉函数。
  因此,欧拉函数 ϕ ( n ) \phi(n) 的求解,归结到了分解质因子这个问题。分解质因子有一个古老的方法,试除法5:求 n n 的质因子时,逐个检查从 2 2 n \sqrt{n} 的全部质数,若是它能整除 n n ,就是一个因子。试除法的复杂度是 O ( n ) O(\sqrt{n}) ,效率很低。在密码学中,有高达百位以上的十进制数分解质因子,所以发明了不少高效率的方法。不过,在算法竞赛中,数据规模不大,因此通常就用试除法。

  下面是求欧拉函数的代码。先用试除法分解质因子,再用公式 ϕ ( n ) = n i = 1 k ( 1 1 p i ) \phi(n)=n \prod_{i=1}^k(1-\frac 1{p_i}) 求得 ϕ ( n ) \phi(n) 。总复杂度仍然是 O ( n ) O(\sqrt{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内全部的欧拉函数

  如今要求 1 n 1~n 范围内全部的欧拉函数,而不是只求一个欧拉函数。前面求一个欧拉函数的复杂度是 O ( n ) O(\sqrt{n}) ,若是一个一个地求 1 n 1~n 范围内全部的欧拉函数,那么 n n 个的总复杂度是 O ( n n ) O(n\sqrt{n}) 。效率过低。
  这个过程能够优化:
  (1)利用积性函数的性质 ϕ ( p 1 p 2 ) = ϕ ( p 1 ) ϕ ( p 2 ) \phi(p_1p_2)=\phi(p_1)\phi(p_2) ,求得前面的 ϕ ( p 1 ) \phi(p_1) ϕ ( p 2 ) \phi(p_2) 后,能够递推出后面的 ϕ ( p 1 p 2 ) \phi(p_1p_2)
  (2)求 1 n 1~n 范围内的质数,用于上一步骤的递推。
  编程时,能够这样操做:
  (1)求得 1 n 1~n 内每一个数的最小质因子。所谓最小质因子,例如 84 = 2 2 × 3 × 7 84=2^2 \times 3 \times 7 84 84 的最小质因子是 2 2 。这步操做能够用欧拉筛,复杂度是 O ( n ) O(n) 的,不可能更快了。
  (2)在第(1)步基础上,递推求 1 n 1~n 内每一个数的欧拉函数。这一步的复杂度也是 O ( n ) O(n) 的。
  二者相加,总复杂度是 O ( n ) O(n) 的。

2.3.1 欧拉筛

  欧拉筛是一种线性筛,也就是说,它能在 O ( n ) O(n) 的线性时间里求得 1 n 1~n 内全部的质数。
  欧拉筛是对埃氏筛法的改进。埃氏筛法的原理是:每一个合数都是多个质数的积;那么从最小的质数2开始,用每个质数去筛比它大的数,就能筛掉合数。埃氏筛法低效的缘由是,一个合数会被它的多个质因子重复筛。请读者本身详细了解埃氏筛法。
  欧拉筛法的原理是:一个合数确定有一个最小质因子;让每一个合数只被它的最小质因子筛选一次,以达到不重复筛的目的。
  具体操做步骤:
  (1)逐一检查 2 n 2~n 的全部数。第一个检查的是 2 2 ,它是第一个质数。
  (2)当检查到第 i i 个数时,利用已经求得的质数去筛掉对应的合数 x x ,并且是用 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”, 其中的 p r i m e [ j ] prime[j] 是最小质因子, v i s [ i p r i m e [ j ] ] = 1 ; vis[i*prime[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 只被筛了一次,并且是被 x x 的最小质因子筛去的。这说明筛法是线性的,复杂度 O ( n ) O(n)
  如今用欧拉筛求 1 n 1~n 内每一个数的最小质因子。只要简单修改上述程序,直接用 v i s [ M A X N ] 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 是质数, ϕ ( n ) = n 1 \phi(n)=n-1
  状况(2):若 n = p k n=p^k p p 是质数,有 ϕ ( n ) = p k 1 ϕ ( p ) \phi(n)=p^{k-1}\phi(p)
  状况(3):若 n = p 1 p 2 n=p_1p_2 ,用 ϕ ( n ) = ϕ ( p 1 ) ϕ ( p 1 ) \phi(n)=\phi(p_1)\phi(p_1) 递推。

洛谷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 欧拉函数前缀和

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

3 整除分块(数论分块)

  整除分块是杜教筛算法的基本思路。
  整除分块是为了解决一个整除的求和问题: i = 1 n n i \sum_{i=1}^n\lfloor\dfrac{n}{i}\rfloor
  其中 n i \lfloor\dfrac{n}{i}\rfloor 表示 n n 除以 i i ,并向下取整6 n n 是给定的一个整数。

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

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

  对第2行求和,结果是: i = 1 8 8 i = 20 \sum_{i=1}^8\lfloor\frac{8}{i}\rfloor=20
  观察上面第 2 2 8 i \dfrac 8i 的值,发现是逐步变小的,而且不少相等,这是整除操做的规律。例如:
   8 / 1 = 8 8/1 = 8
   8 / 2 = 4 8/2 = 4
   8 / 3 = 8 / 4 = 2 8/3 = 8/4 = 2
   8 / 5 = 8 / 6 = 8 / 7 = 8 / 8 = 1 8/5 = 8/6 = 8/7 = 8/8 = 1
  为了对这些整除的结果进行快速求和,天然能想到,把它们分红一块块的,其中每一块的 8 i \dfrac 8i 相同,把这一块一块儿算,就快多了。
  设每一块的左右区间是 [ l , r ] [l, r] ,上面的例子能够分红4块: [ 1 , 1 ] [ 2 , 2 ] [ 3 , 4 ] [ 5 , 8 ] [1,1]、[2,2]、[3,4]、[5,8]
  每一块内部的求和计算,只要用数量乘以值就好了,例如 [ 5 , 8 ] [5,8] 区间的整除和: ( 8 5 + 1 ) 1 = 4 (8-5+1)*1=4 。这个计算的复杂度是 O ( 1 ) O(1) 的。
  最后还有两个问题:
  (1)把 n i \lfloor\dfrac{n}{i}\rfloor 按相同值分块,一共须要分红几块?或者说, n i \lfloor\dfrac{n}{i}\rfloor 有几种取值?这决定了算法的复杂度。
  【答】分块少于 2 n 2\sqrt{n} 种,证实以下:
  (a) i n i\leq\sqrt{n} 时, n i \dfrac{n}{i} 的值有: { n 1 , n 2 , n 3 n n } \lbrace \frac{n}{1},\frac{n}{2},\frac{n}{3}\cdots \frac{n}{\sqrt{n}} \rbrace n i n \frac{n}{i}\geq\sqrt{n} ,共 n \sqrt{n} 个,此时 n i \lfloor\dfrac{n}{i}\rfloor n \sqrt{n} 种取值;
  (b) i > n i>\sqrt{n} 时,有 n i < n \dfrac{n}{i}<\sqrt{n} ,此时 n i \lfloor\dfrac{n}{i}\rfloor 也有 n \sqrt{n} 种取值。
二者相加,共 2 n 2\sqrt{n} 种。因此,整除分块的数量是 O ( n ) O(\sqrt{n}) 的。
  (2)如何计算?或者说,给定每一个块的第一个数 l l ,能推出这个块的最后一个数 r r 吗?
  【答】能够证实: r = n / ( n / l ) r = n/(n/l) 。证实略7

  下面是标程。

#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;               //打印和
}

  整除分块习题:洛谷P182九、洛谷P226一、洛谷P3935。

4 狄利克雷卷积

  狄利克雷卷积 Dirichlet Convolution,是计算求和问题的有用工具。狄利克雷卷积是杜教筛用到的核心技术之一。
  定义:
  设 f g f,g 是算术函数,记 f f g g 的狄利克雷卷积是 f g f*g ,定义为8
( f g ) ( n ) = d n f ( d ) g ( n d ) (f*g)(n)= \sum_{d|n}f(d)g(\frac{n}{d})

  它是一个求和公式,其中 d n d|n 表示 d d 整除 n n ,卷积是对正因子求和。
  这个卷积是什么含义呢?如今给出一个特殊的例子。定义恒等函数 I ( n ) = n I(n)=n 、常数函数 1 ( n ) = 1 1(n)=1 ,它们的卷积是:
( I 1 ) ( n ) = d n I ( d ) 1 ( n d ) = d n d 1 = d n d = σ ( n ) (I*1)(n)=\sum_{d|n}I(d)1(\frac{n}{d})=\sum_{d|n}{d} \cdot1=\sum_{d|n}d=\sigma(n)
  把它记为: I 1 = σ I*1=\sigma
   σ ( n ) \sigma(n) 是“因子和函数”的符号。例如n=12,那么 σ ( n ) = 1 + 2 + 3 + 4 + 6 + 12 = 28 \sigma(n)=1+2+3+4+6+12=28 σ ( n ) \sigma(n) 是积性函数。
  狄利克雷卷积的性质:
  (1)交换律: f g = g f f*g=g*f
  (2)结合律: ( f g ) h = f ( g h ) (f*g)*h=f*(g*h)
  (3)分配律: f ( g + h ) = ( f g ) + ( f h ) f*(g+h)=(f*g)+(f*h)
  (4)元函数: ϵ ( n ) = { 1 , n=1 0 , n>1 \epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases} 。对任何 f f ,有 e f = f e = f e*f=f*e=f
  (5)两个积性函数的狄利克雷卷积仍然是积性函数。
  有一些积性函数算术函数在狄利克雷卷积中经常用到,1.3节已经给出,这里再次列出:
   I ( n ) I(n) :恒等函数, I ( n ) = 1 I(n)=1
   i d ( n ) id(n) :单位函数, i d ( n ) = n id(n)=n
   I k ( n ) I_k(n) :幂函数, I k ( n ) = n k I_k(n)=n^k
   ϵ ( n ) \epsilon(n) :元函数, ϵ ( n ) = { 1 , n=1 0 , n>1 \epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}
   σ ( n ) \sigma(n) :因子和函数, σ ( n ) = d n d \sigma(n)=\sum_{d|n}d
   d ( n ) d(n) :约数个数, d ( n ) = d n 1 d(n)=\sum_{d|n}1

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

  若是读者困惑莫比乌斯函数有什么用,能够先看下一节“5.2 莫比乌斯函数的由来”。

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

  莫比乌斯函数以数学家Möbius的名字命名9

  莫比乌斯函数是狄利克雷卷积的重要工具,它也是数论和组合数学的重要的积性函数。
  莫比乌斯函数 μ ( n ) \mu(n) 定义为10

μ ( n ) = { 1 , 若是n=1 ( 1 ) r , 若是 n = p 1 p 2 . . . p r p i 0 , 其余状况 \mu(n)= \begin{cases} 1, & \text {若是n=1}\\ (-1)^r, & \text{若是$n = p_1 p_2...p_r,其中p_i是不一样的质数$} \\ 0, & \text{其余状况} \end{cases}
  一些例子: μ ( 1 ) = 1 , μ ( 2 ) = 1 , μ ( 3 ) = 1 , μ ( 4 ) = μ ( 2 2 ) = 0 , μ ( 5 ) = 1 , μ ( 6 ) = μ ( 2 × 3 ) = 1 , μ ( 7 ) = 1 , μ ( 8 ) = μ ( 2 3 ) = 0 ; μ ( 330 ) = μ ( 2 × 3 × 5 × 11 ) = ( 1 ) 4 = 1 , μ ( 660 ) = μ ( 2 2 × 3 × 5 × 11 ) = 0 \mu(1)=1,\mu(2)=-1,\mu(3)=-1,\mu(4)=\mu(2^2)=0,\mu(5)=-1,\mu(6)=\mu(2\times3)=1,\mu(7)=-1,\mu(8)=\mu(2^3)=0; \mu(330)=\mu(2\times3\times5\times11)=(-1)^4=1,\mu(660)=\mu(2^2\times3\times5\times11)=0。
  莫比乌斯函数 μ ( n ) \mu(n) 是积性函数。它有一个和欧拉函数的 n = d n ϕ ( d ) n=\sum_{d|n}\phi(d) 相似的定理:
【定理5.1】莫比乌斯函数的和函数在整数 n n 处的值,知足:
d n μ ( d ) = { 1 , n=1 0 , n>1 \sum_{d|n}\mu(d)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}
  证实:
  (1) n = 1 n=1 时,显然有 F ( 1 ) = d n μ ( 1 ) = 1 F(1)=\sum_{d|n}\mu(1)=1
  (2) n > 1 n>1 时。根据积性函数的定义,有 F ( n ) = F ( P 1 a 1 ) F ( P 2 a 2 ) F ( P t a t ) F(n)=F(P_1^{a_1})F(P_2^{a_2})\cdots F(P_t^{a_t}) ,其中 n = P 1 a 1 P 2 a 2 P t a t n=P_1^{a_1}P_2^{a_2} \cdots P_t^{a_t} 是质因子分解。若是能证实 F ( P k ) = 0 F(P^k)=0 即有 F ( n ) = 0 F(n)=0 。由于当 i 2 i≥2 时, μ ( p i ) = 0 \mu(p^i)=0 ,有:
F ( p k ) = d p k μ ( d ) = μ ( 1 ) + μ ( p ) + μ ( p 2 ) + + μ ( p k ) = 1 + ( 1 ) + 0 + + 0 = 0 F(p^k)=\sum_{d|p^k}\mu(d)=\mu(1)+\mu(p)+\mu(p^2)+\cdots+\mu(p^k)=1+(-1)+0+\cdots+0=0
  下面是莫比乌斯函数线性筛代码,求 1 n 1~n 内的莫比乌斯函数,和欧拉函数线性筛的代码几乎同样,此处不作解释。

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 的和函数 F F 的值为 F ( n ) = d n f ( d ) F(n)=\sum_{d|n}f(d) ,显然, F F f f 的值决定。这种关系能够反过来吗?也就是说,是否存在一种用 F F 求出 f f 的值的简便方法?本节给出这样的公式,看什么样的公式可行。
  若 f f 是算术函数, F F 是它的和函数 F ( n ) = d n f ( d ) F(n)=\sum_{d|n}f(d) ,按定义展开 F ( n ) n = 1 , 2 , . . . , 8 F(n),n=1,2,...,8 ,有:
   F ( 1 ) = f ( 1 ) F(1) = f(1)
   F ( 2 ) = f ( 1 ) + f ( 2 ) F(2) = f(1) + f(2)
   F ( 3 ) = f ( 1 ) + f ( 3 ) F(3) = f(1) + f(3)
   F ( 4 ) = f ( 1 ) + f ( 2 ) + f ( 4 ) F(4) = f(1) + f(2) + f(4)
   F ( 5 ) = f ( 1 ) + f ( 5 ) F(5) = f(1) + f(5)
   F ( 6 ) = f ( 1 ) + f ( 2 ) + f ( 3 ) + f ( 6 ) F(6) = f(1) + f(2) + f(3) + f(6)
   F ( 7 ) = f ( 1 ) + f ( 7 ) F(7) = f(1) + f(7)
   F ( 8 ) = f ( 1 ) + f ( 2 ) + f ( 4 ) + f ( 8 ) F(8) = f(1) + f(2) + f(4) + f(8)
  根据上面方程,能够解得:
   f ( 1 ) = F ( 1 ) f(1) = F(1)
   f ( 2 ) = F ( 2 ) F ( 1 ) f(2) = F(2) - F(1)
   f ( 3 ) = F ( 3 ) F ( 1 ) f(3) = F(3) - F(1)
   f ( 4 ) = F ( 4 ) F ( 2 ) f(4) = F(4) - F(2)
   f ( 5 ) = F ( 5 ) F ( 1 ) f(5) = F(5) - F(1)
   f ( 6 ) = F ( 6 ) F ( 3 ) F ( 2 ) + F ( 1 ) f(6) = F(6) - F(3) - F(2) + F(1)
   f ( 7 ) = F ( 7 ) F ( 1 ) f(7) = F(7) - F(1)
   f ( 8 ) = F ( 8 ) F ( 4 ) f(8) = F(8) - F(4)
  注意到 f ( n ) f(n) 等于形式为 ± F ( n d ) \pm F(\frac nd) 的一些项之和,其中 d n d|n 。从这一结果中,可能有这样一个等式,形式是:
f ( n ) = d n μ ( d ) F ( n d ) f(n)=\sum_{d|n}\mu(d)F(\frac nd)
  其中 μ \mu 是算术函数。若是等式成立,计算获得: μ ( 1 ) = 1 , μ ( 2 ) = 1 , μ ( 3 ) = 1 , μ ( 4 ) = 0 , μ ( 5 ) = 1 , μ ( 6 ) = 1 , μ ( 7 ) = 1 , μ ( 8 ) = 0 \mu(1)=1,\mu(2)=-1,\mu(3)=-1,\mu(4)=0,\mu(5)=-1,\mu(6)=1,\mu(7)=-1,\mu(8)=0 。(读者已经发现和前一节莫比乌斯函数的值同样。)
  又 F ( p ) = f ( 1 ) + f ( p ) F(p)=f(1)+f(p) ,得 f ( p ) = F ( p ) F ( 1 ) f(p)=F(p)-F(1) ,其中 p p 是质数。则 μ ( p ) = 1 \mu(p)=-1
  又由于 F ( p 2 ) = f ( 1 ) + f ( p ) + f ( p 2 ) F(p^2)=f(1)+f(p)+f(p^2)
  有: f ( p 2 ) = F ( p 2 ) ( F ( p ) F ( 1 ) ) F ( 1 ) = F ( p 2 ) F ( p ) f(p^2)=F(p^2)-(F(p)-F(1))-F(1)=F(p^2)-F(p)
  这要求对任意质数 p p ,有 μ ( p 2 ) = 0 \mu(p^2)=0 。相似的推理得出对任意质数 p p 及整数 k > 1 k>1 ,有 μ ( p k ) = 0 \mu(p^k)=0 。若是猜测 μ \mu 是积性函数,则 μ \mu 的值就由质数幂处的值决定。
  这就是莫比乌斯函数的由来。

5.3 莫比乌斯反演

【定理5.2】莫比乌斯反演(Möbius inversion)公式。若 f f 是算术函数, F F f f 的和函数,对任意正整数 n n ,知足 F ( n ) = d n f ( d ) F(n)=\sum_{d|n}f(d) ,则对任意正整数 n n ,有 f ( n ) = d n μ ( d ) F ( n d ) f(n)=\sum_{d|n}\mu(d)F(\frac nd)
  从定理5.2能够推出下面的定理5.3。
【定理5.3】 f f 是算术函数,它的和函数为 F ( n ) = d n f ( d ) F(n)=\sum_{d|n}f(d) ,那么若是 F F 是积性函数,则 f f 也是积性函数。
  莫比乌斯反演,不要求 f f 是积性函数。

6 杜教筛

6.1 杜教筛公式

6.1.1 杜教筛公式的推导

  杜教筛要解决的是这一类问题:设 f ( n ) f(n) 是一个数论函数,计算 S ( n ) = i = 1 n f ( i ) S(n)=\sum_{i=1}^nf(i)
  因为 n n 很大,要求以低于 O ( n ) O(n) 的复杂度求解,因此用前面提到的线性筛是不够的。若是能用到整除分块,每一块的值相等一块儿算,就加快了速度,也就是构造形如 i = 1 n n i \sum_{i=1}^n\lfloor\frac{n}{i}\rfloor 的整除分块。
  杜教筛的思路,简单地说,是把f(n)构形成能利用整除分块的新的函数,这个构造用到了卷积。
  记 S ( n ) = i = 1 n f ( i ) S(n)=\sum_{i=1}^nf(i) 。根据 f ( n ) f(n) 的性质,构造一个 S ( n ) S(n) 关于 S ( n i ) S(\lfloor\frac{n}{i}\rfloor) 的递推式,下面是构造方法11

  构造2个积性函数 h h g g ,知足卷积: h = g f h=g*f
  根据卷积的定义,有 h ( i ) = d i g ( d ) f ( i d ) h(i)= \sum_{d|i}g(d)f(\frac{i}{d}) ,求和:
i = 1 n h ( i ) = i = 1 n d i i g ( d ) f ( i d ) \sum_{i=1}^nh(i)=\sum_{i=1}^n\sum_{d|i}^ig(d)f(\frac{i}{d})
= d = 1 n g ( d ) d i n f ( i d ) \qquad\qquad =\sum_{d=1}^ng(d)\sum_{d|i}^nf(\frac{i}{d})
= d = 1 n g ( d ) i = 1 n d f ( i ) \qquad\qquad =\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)
= d = 1 n g ( d ) s ( n d ) \qquad\qquad =\sum_{d=1}^ng(d)s({\lfloor\frac{n}{d}\rfloor})
  注意其中第2行变换 i = 1 n d i i g ( d ) f ( i d ) = d = 1 n g ( d ) d i n f ( i d ) \sum_{i=1}^n\sum_{d|i}^ig(d)f(\frac{i}{d})=\sum_{d=1}^ng(d)\sum_{d|i}^nf(\frac{i}{d}) ,网上昵称“杜教筛变换” 11

  最后一步获得: i = 1 n h ( i ) = d = 1 n g ( d ) s ( n d ) \sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)s({\lfloor\frac{n}{d}\rfloor})
  为计算 S ( n ) S(n) ,把式子右边的第一项提出来:
i = 1 n h ( i ) = g ( 1 ) S ( n ) + d = 2 n g ( d ) s ( n d ) \sum_{i=1}^nh(i)=g(1)S(n)+\sum_{d=2}^ng(d)s({\lfloor\frac{n}{d}\rfloor})
g ( 1 ) S ( n ) = i = 1 n h ( i ) d = 2 n g ( d ) s ( n d ) g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)s({\lfloor\frac{n}{d}\rfloor})
  式中的 i i d d 无关,为方便看,改写为:
g ( 1 ) S ( n ) = i = 1 n h ( i ) i = 2 n g ( i ) s ( n i ) g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor})
  这就是杜教筛公式。

6.1.2 杜教筛算法和复杂度

  一个积性函数求和问题,若是分析获得了杜教筛公式,工做已经完成了一大半,剩下的是编程实现。
  公式 g ( 1 ) S ( n ) = i = 1 n h ( i ) i = 2 n g ( i ) s ( n i ) g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor}) S ( n ) S(n) 的递归形式,编程时能够递归实现;递归的时候加上记忆化,让每一个 S [ i ] S[i] 只须要算一次。
  其中的 h ( i ) h(i) g ( i ) g(i) ,若是构造得好,就容易计算,见后面欧拉函数和莫比乌斯函数的例子。为方便分析计算复杂度,下面略去 h ( i ) h(i) g ( i ) g(i) ,分析简化后的式子 S ( n ) = i = 2 n s ( n i ) S(n)=\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor}) 的复杂度。
  设 S ( n ) S(n) 的复杂度是 T ( n ) T(n)
  递归展开第一层: S ( n ) = i = 2 n s ( n i ) = s ( n 2 ) + s ( n 3 ) + + s ( n n ) S(n)=\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})=s({\lfloor\frac{n}{2}\rfloor})+s({\lfloor\frac{n}{3}\rfloor})+\cdots+s({\lfloor\frac{n}{n}\rfloor})
  根据整除分块的原理,右边共有 O ( n ) O(\sqrt{n}) s ( n i ) s({\lfloor\frac{n}{i}\rfloor}) 。另外,把它们加起来的时候,还有 O ( n ) O(\sqrt{n}) 次合并。总时间是 T ( n ) = O ( n ) + O ( i = 2 n T ( n i ) ) T(n)=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}T(\frac ni))
  再展开一层: T ( n i ) = O ( n i ) + O ( k = 2 n i T ( n i k ) ) = O ( n i ) T(\frac ni)=O(\sqrt{\frac ni})+O(\sum_{k=2}^{\sqrt {\frac ni}}T(\frac {\frac ni}k))=O(\sqrt{\frac ni}) ,中间第2项是高阶小量,把它略去。
  合起来:
T ( n ) = O ( n ) + O ( i = 2 n T ( n i ) ) T(n)=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}T(\frac ni))
= O ( n ) + O ( i = 2 n n i ) \qquad\qquad=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}\sqrt{\frac ni})
  后一项最大,只考虑它就够了:
T ( n ) = i = 2 n O ( n i ) O ( 0 n n x d x ) = O ( n 3 4 ) T(n)=\sum_{i=2}^{\sqrt {n}}O(\sqrt{\frac ni})\approx O(\int_0^{\sqrt{n}} {\sqrt{\frac nx}} \,{\rm d}x)=O(n^\frac 34)
  上述计算还能优化,即预处理一部分 s ( n ) s(n) s ( n ) s(n) 是一个积性函数的前缀和,先用线性筛计算一部分,假设已经预处理了前 k k 个正整数的 s ( n ) s(n) ,且 k n k\geq \sqrt{n} 。由于 S ( 1 ) S(1) S ( k ) S(k) 已经求出,那么这个式子中尚未求出的是 k < n i n k<{\lfloor\frac{n}{i}\rfloor}\leq{n} ,也就是要算 1 < i n k 1<i\leq{\frac nk} 的这一部分。那么复杂度变为:
T ( n ) = i = 2 n k O ( n i ) O ( 0 n k n x d x ) = O ( n k ) T(n)=\sum_{i=2}^{\frac nk}O(\sqrt{\frac ni})\approx O(\int_0^{\frac nk} {\sqrt{\frac nx}} \,{\rm d}x)=O(\frac n{\sqrt k})
   k k 取多大合适呢?线性筛计算也是要花时间的,并且是线性的 O ( k ) O(k) ,那么线性筛计算加上杜教筛公式计算,总时间是:
T ( n ) = O ( k ) + T ( n ) = O ( k ) + O ( n k ) T'(n)=O(k)+T(n)=O(k)+O(\frac n{\sqrt k})
  差很少 k = n 2 3 k=n^{\frac 23} 时,它有最小值,即: O ( n 2 3 ) + O ( n n 2 3 ) = O ( n 2 3 ) O(n^{\frac 23})+O(\frac n{\sqrt{n^{\frac 23}}})=O(n^{\frac 23})
  它比 O ( n 3 4 ) O(n^{\frac 34}) 要好。
  以上就是杜教筛算法的所有。
  综合起来,杜教筛算法能够概况为:用狄利克雷卷积构造递推式,编程时用整除分块和线性筛优化,算法复杂度达到了 O ( n 2 3 ) O(n^{\frac 23})
  杜教筛算法 = 狄利克雷卷积 + 整除分块 + 线性筛

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

6.2 欧拉函数前缀和

  求欧拉函数前缀和: i = 1 n ϕ ( i ) \sum_{i=1}^n\phi(i)
  求欧拉函数前缀和的杜教筛公式,须要用到欧拉函数性质: n = d n ϕ ( d ) n=\sum_{d|n}\phi(d)
  (1)直接套用杜教筛公式
  按卷积定义 h = f g = d n f ( d ) g ( n d ) h=f*g= \sum_{d|n}f(d)g(\frac{n}{d}) ,把 n = d n ϕ ( d ) n=\sum_{d|n}\phi(d) 改写为: i d = ϕ I id=\phi*I ,那么有 h = i d g = I h=id,g=I 。代入杜教筛公式 g ( 1 ) S ( n ) = i = 1 n h ( i ) i = 2 n g ( i ) s ( n i ) g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor}) ,得:
S ( n ) = i = 1 n i i = 2 n s ( n i ) S(n)=\sum_{i=1}^ni-\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})
= n ( n + 1 ) 2 i = 2 n s ( n i ) \qquad\qquad=\frac{n(n+1)}{2}-\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})
  (2)本身推导
  读者若是有兴趣,也能够本身循序渐进地以欧拉函数为例,作一次杜教筛公式推导的练习,推导过程是12

n = ϕ ( n ) + d n , d < n ϕ ( d ) n=\phi(n)+\sum_{d|n,d<n}\phi(d)
ϕ ( n ) = n d n , d < n ϕ ( d ) \phi(n)=n-\sum_{d|n,d<n}\phi(d)
i = 1 n ϕ ( i ) = i = 1 n ( i d i , d < i ϕ ( d ) ) \sum_{i=1}^n\phi(i)=\sum_{i=1}^n(i-\sum_{d|i,d<i}\phi(d))
i = 1 n ϕ ( i ) = n ( n + 1 ) 2 i = 1 n d i , d < i ϕ ( d ) \sum_{i=1}^n\phi(i)=\frac{n(n+1)}{2}-\sum_{i=1}^n\sum_{d|i,d<i}\phi(d)
  如今改变枚举变量:枚举 i d \dfrac id 的值,即枚举 i i d d 的倍数,由于 i d i≠d ,因此从2开始:
i = 1 n ϕ ( i ) = n ( n + 1 ) 2 i d = 2 i d n d = 1 d n i d ϕ ( d ) \sum_{i=1}^n\phi(i)=\frac{n(n+1)}{2}-\sum_{\frac id=2}^{\frac id \leq n}\sum_{d=1}^{ d \leq {\lfloor\frac{n}{\frac id}\rfloor} }\phi(d)
  设 S ( n ) = i = 1 n ϕ ( i ) S(n)=\sum_{i=1}^n\phi(i) ,并把 i d \dfrac id 写成 i i'
S ( n ) = n ( n + 1 ) 2 i = 2 n s ( n i ) S(n)=\frac{n(n+1)}{2}-\sum_{i’=2}^ns({\lfloor\frac{n}{i‘}\rfloor})
  这样就获得了欧拉函数的“杜教筛公式”,和(1)同样。

6.3 莫比乌斯函数前缀和

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

  (2)本身推导13

[ n = 1 ] = μ ( n ) + d n , d < n μ ( d ) [n=1]=\mu(n)+\sum_{d|n,d<n}\mu(d)
μ ( n ) = [ n = 1 ] d n , d < n μ ( d ) \mu(n)=[n=1]-\sum_{d|n,d<n}\mu(d)
i = 1 n μ ( i ) = 1 i = 1 n d i , d < i μ ( d ) \sum_{i=1}^n\mu(i)=1-\sum_{i=1}^n\sum_{d|i,d<i}\mu(d)
i = 1 n μ ( i ) = 1 i d = 2 i d n d = 1 d n i d i = 1 n μ ( 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 ) = i = 1 n μ ( i ) S(n)=\sum_{i=1}^n\mu(i) ,并把 i d \dfrac id 写成 i i' ,得:
S ( n ) = 1 i = 2 n S ( n i ) S(n)=1-\sum_{i‘=2}^nS({\lfloor\frac{n}{i’}\rfloor})

6.4 模板代码

  下面给出洛谷P4213的代码14

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

#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)求 i = 1 n i ϕ ( i ) \sum_{i=1}^ni\phi(i)
  题解:http://fzl123.top/archives/898
  (2)求 i = 1 n j = 1 n σ ( i j ) \sum_{i=1}^n\sum_{j=1}^n\sigma(ij) ,其中 n 1 0 9 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


本文的注脚:


  1. 《初等数论及其应用》第7章乘性函数 ↩︎

  2. 有乘性函数,也有加性函数,参考《初等数论及其应用》182页。 ↩︎

  3. 有关欧拉函数的全部证实,见《初等数论及其应用》176-180页。 ↩︎

  4. [欧拉定理:设 m m 是一个正整数, a a 是一个整数且 ( a , m ) = 1 (a, m)=1 ,那么 a ϕ ( m ) 1 ( m o d   m ) a^{\phi(m)}\equiv1(mod\ m) ]。 ↩︎

  5. 试除法很低效,有不少快速分解质因子的方法,参考《初等数论及其应用》93页。 ↩︎

  6. 除法运算,向上取整,英文Ceiling,用符号 \lceil \rceil 表示;向下取整,英文Floor,用符号 \lfloor \rfloor 表示。 ↩︎

  7. 证实细节见:https://blog.csdn.net/qq_41021816/article/details/84842956↩︎

  8. 参考:https://brilliant.org/wiki/dirichlet-convolution/↩︎

  9. https://brilliant.org/wiki/mobius-function/↩︎

  10. 《初等数论及其应用》200页。 ↩︎

  11. https://blog.csdn.net/myjs999/article/details/78906549 ↩︎ ↩︎

  12. 引用自:https://blog.csdn.net/qq_34454069/article/details/79492437 ↩︎

  13. 引用自:https://blog.csdn.net/qq_34454069/article/details/79492437 ↩︎

  14. 改写自:https://blog.csdn.net/KIKO_caoyue/article/details/100061406 ↩︎