Project Euler 73: Counting fractions in a range

Posted 2021-03-24 metaquant

对于分数(n/d)，其中(n,d)均为正整数，如果(n<d)且两者的最大公约数(HCF(n,d)=1)，这个分数就称之为简分数。如果我们把(dle8)的所有简分数以从小到大的顺序排列，则有：
[ 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}, ]
可以看到(1/3)和(1/2)之间有三个分数。当(dle12000)时，在这个从小到大排列的简分数集合中有多少个分数处于(1/3)和(1/2)之间？

分析：此题至少有两种解法：第一种解法是利用之间提到的Farey Sequence的性质，相对容易理解，但效率不够高；第二种解法则是利用和七十二题中类似的推导方法，将这个问题转化一个动态规划问题，这种方法效率很高，但是理解起来难度也更大。下面分别简述两种方法：

方法一：题目中所列的Farey Sequence满足以下两个性质：

• 两个简分数(a/b)和(c/d)在Farey Sequence中相邻的必要充分条件是(bc-ad=1)
• Farey Sequence中连续三个分数(a/b,c/d,e/f)满足条件(c/d=(a+e)/(b+f))

题目要求在(dle12000)时，Farey Sequence中处于(1/3)和(1/2)之间的分数有多少个？则我们可以如下思考：第一步，首先确定在(dle12000)的条件下，最紧邻(1/3)的分数是多少？根据上面列出的性质一，假设这个分数是(p/q)，则(3p-q=1Rightarrow3p=q+1)，即(q+1)应是三的倍数，则(q)最大只能取(11999)，据此得到(p=12000/3=4000)，即在(dle12000)时，最紧邻(1/3)的分数应该是(4000/11999)。第二步，在已知(1/3)和(4000/11999)紧邻的情况下，我们要推算下一个紧邻(4000/11999)的分数，如此不断推算下去并记录已经产生的分数个数，直到产生的下一个分数成为(1/2)时停止，这时记录的分数个数即为题目所求。根据以上列出的性质二，因为(c/d)是最简分数的形式，则有(kd=b+fRightarrow f=kd-b)，据题意有(kd-ble12000)，在满足此条件时(k)的最大值为(lfloor (12000+b)/d floor)，则紧邻(a/b)和(c/d)的下一个分数(e/f)满足条件：
[ e=lfloor (12000+b)/d floor cdot c-a, f= lfloor (12000+b)/d floor cdot b-b ]
如此，我们就可以按照从小到大的顺序依次列出Farey Sequence中(1/3)以后的所有分数，直到这个分数达到(1/2)截止，也即(f=2)时截止，此时返回分数个数即为题目所求。此种方法的时间复杂度约为(O(n^2))，所以在处理小规模问题尚可应付，而在处理大规模问题时则会耗时很久。下面介绍第二种方法：

方法二：首先定义如下两个函数：
[ egin{aligned}F(N)&=cardigg{frac{k}{n}:frac{1}{3}<frac{k}{n}<frac{1}{2},nle N igg}\R(N)&=cardigg{frac{k}{n}:frac{1}{3}<frac{k}{n}<frac{1}{2},nle N,gcd(k,n)=1 igg}end{aligned} ]
其中(card)表示集合的基数，也即集合中元素的个数，则我们可以明显看出(R(N))即为我们要求的值。同时根据和七十二题中相似的思路，我们有：
[ F(N)=sum_{m=1}^N R(lfloorfrac{N}{m} floor)Rightarrow F(N)=R(N)+sum_{m=2}^N R(lfloorfrac{N}{m} floor) ]
则有：
[ R(N)=F(N)-sum_{m=2}^N R(lfloorfrac{N}{m} floor) ]
显然我们得到了一个和七十二题中非常类似的递归式，因此可以用动态规划解决，并且可以用上整除分块的技巧。现在的问题是我们需要确定(F(N))，只有当(F(N))容易计算时，上面的递归推导才是有意义的。如果我们设(N=6q+r,0le r<6)，则：
[ Fleft( N ight) =sum_{i=1}^N{left( igglfloor frac{i-1}{2} igg floor -igglfloor frac{i}{3} igg floor ight) =qleft( 3q-2+r ight) +egin{cases} 1& r=5\ 0& r e 5\end{cases}} ]
根据上式我们可以很容易的计算出(F(N))。方法二的时间复杂度约为(O(n^{3/4}))，属于亚线性算法，在处理大规模数据时相对于方法一有很大优势。

两种方法的代码如下：

# approach 1, time cost = 2.97 s ± 6.42 ms

def farey_seq(n=12000):
    b, d = 3, 11999
    ans = 0
    while d!=2:
        ans += 1
        k = int((n + b) / d)
        b, d = d, k * d - b
    return ans


# approach 2, time cost = 101 ns ± 0.171 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

def fractions_count(n):
    q,r = divmod(n,6)
    i = 1 if r == 5 else 0
    ans = q*(3*q-2+r)+i
    return ans

@lru_cache(maxsize=2048)
def reduced_fractions_count(n=12000):
    if n == 1:
        return 0
    else:
        return fractions_count(n) - number_theory_block(reduced_fractions_count,n,2)