072. 计算分数个数(Counting fractions)

对于分数\(n/d\),其中\(n,d\)均为正整数,如果\(n<d\)且两者的最大公约数\(HCF(n,d)=1\),这个分数就称之为简分数。如果我们把\(d\le8\)的所有简分数以从小到大的顺序排列,则有:

$$ \frac{1}{8},\frac{1}{7},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{2}{7},\frac{1}{3},\frac{3}{8},\frac{2}{5},\frac{3}{7},\frac{1}{2},\frac{4}{7},\frac{3}{5},\frac{5}{8},\frac{2}{3},\frac{5}{7},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{6}{7},\frac{7}{8}, $$
可以看到这个集合中包含的分数有二十一个。求当\(d\le10^6\)时这个最简分数集合中包含有多少个分数?

分析:题目要求简分数的个数,实际上求与\(d\)互质的数的个数,其中\(d\in[2,10^6]\),很明显这就是求这个区间上每个整数的欧拉函数值之和,即求:

$$ \sum_{i=2}^{n}\varphi(i)=\sum_{i=1}^{n}\varphi(i)-\varphi(1)=\sum_{i=1}^{n}\varphi(i)-1 $$

\(S(n)=\sum_{i=1}^{n}\varphi(i)\),则我们只需要求出\(S(n)\),求\(S(n)\)至少有三类方法:第一种方法根据我们之前导出的计算欧拉函数值的公式,计算这个区间中每个整数的欧拉函数值,由于计算欧拉函数值需要先对整数进行质因数分解,因此当数据规模比较大时,这个方法的效率较低。第二种方法是对传统的埃拉托斯特尼筛进行扩展,使用其可以更快的计算多个数的欧拉函数值,对于这个题的数据规模来说,这种方法的表现已经很不错了,对这种方法的具体介绍大家可以在做完题后参见题目的overview,里面讲的很详细了,这里我不再做介绍。

我这里要讲的是第三种思路,它把这里的求和问题转化为一个动态规划的问题,从而可以大大提高算法的效率,简要叙述如下:要求\(S(n)\),即求满足\(1\le a\le b\le n\)\(gcd(a,b)=1\)的数对\((a,b)\)的个数,易知这样的数对\((a,b)\)共有\(n(n+1)/2\)个,但数对的最大公约数\((a,b)\)可以取值的范围为\([1,n]\),则我们只需要从总数\(n(n+1)/2\)中减去所有\(gcd(a,b)>1\)的数对的个数,剩下的就是\(gcd(a,b)=1\)的数对个数,即为题目所求。假设\(gcd(a,b)=m\),则\(gcd(a/m,b/m)=1\),则要求在\(1\le a \le b \le n\)\(gcd(a,b)=m\)的数对个数,可以转化为求\(1\le a' \le b' \le \lfloor n/m\rfloor\)\(gcd(a',b')=1\)的数对的个数,而后者即等于\(S(\lfloor n/m \rfloor)\),则有:

$$ S(n)=\frac{1}{2}n(n+1)-\sum_{m=2}^n S(\big\lfloor \frac{n}{m} \big\rfloor) $$

以上即是我们想要得到的递归式,从而可以用动态规划的方法加以解决。可以看到,递归式中出现了形如\(f(\rfloor n/m \rfloor)\)的形式,从而可以用整除分块或者数论分块的技巧进一步提高效率,关于数论分块的介绍,可以参见这篇文章

这道题目中要求的欧拉函数的前\(N\)项和的问题,实际是更一般的求解积性函数前缀和这类问题的特殊情况。积性函数的前缀和问题在信息学竞赛中经常会遇到,且已经经过较为深入的研究,存在几种通用的算法,比如这道题目实际上可以用杜教筛的方法加以解决。要理解此类方法需要比较多的前置知识,包括数论函数、积性函数、狄利克雷卷积、莫比乌斯反演等等,如果大家感兴趣,可以谷歌"积性函数前缀和",便可以看到很多相关文章,我这里不再做进一步的介绍。

# time cost = 225 ns ± 0.296 ns

from functools import lru_cache

def number_theory_block(f,n,i=1):
    ans = 0
    while i <= n:
        j = n//(n//i)
        ans += (j-i+1)*f(n//i)
        i = j + 1
    return ans

@lru_cache(maxsize=2048)
def sum_of_euler_phi(n):
    if n == 1:
        return 1
    else:
        return n*(n+1)//2 - number_theory_block(sum_of_euler_phi,n,2)

def main(N=10**6):
    return sum_of_euler_phi(N)-1