[转][数论]Farey Sequence深入学习

Farey Sequence深入学习

本文所要解决的问题：

1、如何生成N阶Farey序列

2、如何用Stern Brocot tree查找有理数或逼近无理数

3、如何求N阶Farey序列中第K个分数

一、基础知识：

1、Stern Brocot tree

如图是一棵Stern Brocot tree

Stern Brocot tree是一棵二叉树，初始为(0/1, 1/0)，每次在m/n和m'/n'之间插入

(m + m')/(n+n')，如此不断生成下一层。

Stern Brocot tree有很多神奇而优美的性质：

         1）Stern Brocot tree树可以生成所有的有理数

         2）Stern Brocot tree生成的分数都是不可约的，即gcd(m,n) = 1

         3）任意两个在构造时是相邻的分数m/n和m'/n'，有mn' – m'n = 1

         4）中序遍历Stern Brocot tree树可以得到有序的序列

以上性质的证明可见《具体数学》第四章，即参考[1]

2、Farey Sequence

         Farey序列是Stern Brocot tree的左子树，即所有的真分数

二、N阶Farey序列的构造

         1、递归构造

由于中序遍历Stern Brocot tree可以得到有序的序列，故最直接的方法就是一边生成Stern Brocot tree的左子树一边进行中序遍历

Make_farey(int m1, int n1, int m2, int n2) {

         If(n1 + n2 > N)

                   return;

         Make_farey(m1, n1, m1 + m2, n1 + n2)

         printf("%d/%d", m1 + m2, n1 + n2);

         Make_farey(m1 + m2, n1 + n2, m2, n2);

}

时间复杂度是O(N ^ 2），空间复杂度是O(N)，顺便提一下，N阶Farey序列所包含的元素有(3/π^ 2)O(N ^ 2)+ O(NlgN)

         2、非递归构造

                   对于N阶Farey序列中连续的三个分数m/n，m'/n'，m''/n''

                   有如下递推公式([1]习题4.61)，[ ]表示取下整：

                            m'' = [(n + N) / n']m' - m

                            n'' = [(n + N) / n']n'– n

                   证明：由于这三个数是连续的，故n' + n'' > N，否则m'/n'与m''/n''之间还有分数

                            则N = ((n + N) / n')n'– n >= n'' > N – n' = ((n + N) / n' – 1)n'– n

                            n'' = [(n + N) / n']n' – n

                            又n'm''– m'n'' = 1，得m'' = (m'n'' + 1) / n' = (([(n + N) / n']n' – n)m' + 1) / n'

                            m'' = [(n + N) / n']m' – (nm'– 1) / n' = [(n + N) / n']m' – m

                            得证

                   时间复杂度是O(N ^ 2)，空间复杂度是O(1)

三、在Stern Brocot tree树上查找有理数或逼近无理数

因为Stern Brocot tree是一棵二叉树，所以在它上面查找相当于二分查找（但略有不同），但也可能遇到一些情况使得查找效率退化到线性，这一部分将介绍查找的一般方法以及提高查找效率的方法。

事实上，标题的说法只是为了严谨，计算机里没有无理数，小数在计算机里也不能等同于一个分数，我们真正要解决的问题是如何查找一个分数和如何逼近一个小数。当然，在分子或分母有限制的情况下，查找一个分数也会变成逼近一个分数，不管输入的是分数还是小数，处理的方法都是类似的，所以下面的问题以分数为例。

在具体数学中给出了一种矩阵的表示方法：

我们把分子写在下面，分母写在上面，初始的矩阵为，M(I) = ，是一个单位矩阵。

向矩阵L = 表示向左走，矩阵R = 表示向右走

设M(S) = ，则M(SL) = ，M(SR) = 

定义f(S) = (m + m') / (n + n')

比如M(LRLR) = 对应的分数就是f(LRLR) = 5/8

这与我们在Stern Brocot tree上的搜索是相符合的，事实上，矩阵给出的是一个区间，我们每次修改的是区间的下界或上界。

这样我们对于一个分数就得到了它在Stern Brocot tree树上的字符串表示，形成了一一映射，注意矩阵不满足交换律。

很自然地，我们就得到了如下搜索方式：

S = I

while(f(S) != a/b) {// a/b为要查找的分数

         if(a/b < f(S))

S = SL;

         else

              S = SR       

}

这个方法有很好的可扩展性，稍加修改就可以处理对分母或分子有限制的问题或者是对于小数的逼近，当然实际的程序中我们可能不必真的写成矩阵。

不过这种方法的效率取决于查找的深度，且Stern Brocot tree树不同于二分的地方就是m''/n'' 不等于 (m/n + m'/n') / 2，它的疏密是不同的，若给定一个1/N，则查找的深度就是N了。为了避免退化，我们需要想办法计算连续朝一个方面走的步数，这样我们就不需要一步一步地走了。同时我们可以看到如果遇到经常变换方向的情况，那么这个区间一定会缩小得很快，也就不用担心搜索的深度会太深。

为了计算连续朝一个方向走的步数，我们只需要考虑如下不等式，以向左走为例：

(m' + step * m) / (n' + step * n) <= a / b

n' + step * n <= N

分别求解两个不等式的step，取较小的那个就是连续向左走的步数了。

现在我们再来分析一下矩阵的表示方法，由此引出一个新的结论。
设搜索结束后得到了一个串S
S可能以R或L开头，当S = RS’时
得f(RS’) = f(S’) + 1
即f(RS’) = m/n, f(S’) = (m - n) / n, m > n
同理可得，当S = LS’时
f(LS’) = m/n, f(S’) = m/(n - m), m < n
也就是说我们可以改变m和n的值而不用维护矩阵的状态：
while(m != n) {
if(m < n) {
S = LS;
n = m – n;
}
else {
S = RS;
m = m – n;
}
}
对于这个程序是否很眼熟？
没错！这个就是Euclid算法减法的实现，当然我们一般是用取模来写的，这样一来Farey Sequence就与经典的数论联系起来了。其实Farey Sequence最基本的性质mn’ – nm’ = 1就启发我们已知Farey Sequence中的一项，可以用Extend Euclid求它的前后项。关于Farey Sequence与Euclid算法的关系的讨论，可以见[3]。
这个算法只对分数查找有效，对于小数的逼近还是要用第一种方法，但Euclid算法给出了第一种方法的时间复杂度：O(lgN)。

四、求N阶Farey序列中第K个分数

       此部分全部来自参考[2]

1、 基本的算法

假设我们有了一个子程序，这个程序对于给定的一个实数x，可以返回N阶Farey序列中小于等于x的分数的个数，我们称之为rank(x)。注意：这个子程序非常关键，后面的改进算法都是在这个子程序上做文章。

现在我们可以二分求a，0 < a < N，使得第K个分数在区间[a/N, (a + 1)/N)内，a不要求与N互素。由于对任意q < N，1/q > 1/N，故对任意分母为q，q < N的分数，最多只有1个落在这个区间内。

对于一个给定的q，q < N，由a/N <= p/q < (a + 1)/N得分子p = [((a + 1)q – 1)/N]

这样我们就可以枚举q，并同时检验p/q是否落在这个区间内。

在找第K个分数时我们可以不必存下在这个区间内所有的最简分数然后排序，事实上只需要O(N)的时间复杂度和O(1)的空间复杂度就可以做到：枚举q，存下区间内最小的那个分数，由于这个分数和a/N是紧挨着的，我们可以用第一部分中的递推式生成接下来的分数，由于[a/N, (a + 1)/N)之间最多有N个分数，时间复杂度为O(N)。（也可以用Extend Euclid算法来求a/N的后项）

现在问题的关键就是那个子程序了：

设A[q]:分母为q的最简分数p/q的个数，且p/q <= x，于是有如下的递归式：

在计算A[q]时，如果是枚举q的因子，那还是比较麻烦，可以用类似筛法的方式来求：

for(q = 1; q <= N; ++q)

A[q] = q * a / N;

for(q = 1, ans = 0; q <= N; ++q) {

         ans += A[q];

           for(t = q * 2; t <= N; t += q)

                    A[t] -= A[q];

}

复杂度是（[1]6.59）

因为总共需要lgN次调用该子程序

故时间复杂度是O(N(lgN)^2)，空间复杂度是O(N)

2、 时间上的改进

正如前面所说，我们要改进的是子程序的效率，改进的办法是做一遍预处理，以后只需要用线性的时间计算rank(x)，那么如何预处理呢？

我们观察一下rank(x)，它是[xq]的线性组合，1<= q <=N，而[xq]的系数是独立于[xq]的，也就是说我们可以预先计算出所有[xq]的系数

[xq]最先出现在A[q]里，它的系数是1，然后对于q的倍数t，A[t]中减去了A[q]，即在A[t]中[xt]的系数是1，[xq]的系数是-1，因此这一部分的[xq]的系数是与[xt]是相反的，且对于q的所有的倍数t都是一样

设C[q]:[xq]的系数，则有：

对应于求A[q]的方法，计算C[q]时，我们只要相应地让q从N到1就可以了

预处理需要O(NlgN)，以后每次调用rank(x)需要O(N)

总时间复杂度降为O(NlgN)

3、 空间上的改进

考虑如下的事实：

若[N/q] = [N/q']，则C[q] = C[q']

证明：[xq]最先出现在A[q]，系数为1，然后对于q的倍数m1q，A[m1q]都减去了A[q]，共有[N/q]个m1，然后对于某个m1，m1的倍数m2m1，A[m2m1q]减去了A[m1q]，[x/q]在A[m1q]里的系数是-1，共有[N/(m1q)]个m2，

依次类推，最终是截止在，等价于

，故C[q]只与[N/q]的值有关

因此，对于q > sqrt(N)，实际上只有sqrt(N)个不同的C[q]值

对于q <= sqrt(N)，我们求系数C[q]

对于q >sqrt(N)，我们设D[r]:C[q] = D[r] = D[[N / q]]，对C[q]的表达式稍做变形

将r= [N/q]，D[r] = C[q]，C[qt] = D[[N/(qt)]]代入得：

计算时，先算出D[r]，r<=sqrt(N)再计算C[q]，q <= sqrt(N)

计算D[r]的时间是O(sqrt(N) ^ 2) = O(N)，计算C[q]的时间是O(sqrt(N)lg(sqrt(N)))

但是预处理完后调用rank(x)的时间复杂度还是O(N)

故时间复杂度不变，空间复杂度降到O(sqrt(N))

4、Mobius反演求系数

将上面求D[r]的公式稍做变形，可得

利用Mobius反转公式（见[1]4.61）得：

也就是说D[r]事实上就是Mobius函数的和函数，最终能得到这么优美的结论真是太神奇了！关于Mobius和函数的计算有低于线性的算法，不过对最终的时间复杂度没有提升，这里就不再讨论了。

五、题目汇总：

1、基本概念：

PKU2478

2、生成N阶Farey序列

TOJ2798

PKU3374

SWJTU1481

3、分数或小数的逼近

UVA10077

       ECNU2649

       PKU1650

UVA4275

4、查找N阶Farey序列的第K个分数

       UVA10408

       HDU2432

六、参考：

[1]R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics: A Foundation for Computer Science, 2nd edition, Addison-Wesley, 1994.

[2]Corina E. P tras,cu and Mihai P tras,cu. Computing order statistics in the Farey sequence. In Algorithmic Number Theory, volume 3076 ofLNCS, pages 358_366. Springer, 2004.

[3]杨哲，一类分数问题的研究

---

转自：http://hi.baidu.com/%BC%FB%CF%B0yy/blog/item/d3f10727aa928e1c8a82a1b0.html