Farey序列

今天的 Farey 竟然挂掉了，实在咽不下这口气，晚上又整理了一下资料。

以供大家一起学习。

 

Farey序列

       Fn = {a/b | gcd(a,b)=1 && 0<=a,b<=n};

       即由小于或等于 n 的整数所组成的不可再约分数的递增序列，并满足分子分母互质。

       如：

       F1 = {0/1, 1/1}

       F2 = {0/1, 1/2, 1/1}

       F3 = {0/1, 1/3, 1/2, 2/3, 1/1}

 

性质

       除了 F1 ，其余 Farey 序列都有奇数个元素，并且中间值是 1/2 。

       Farey 序列是一个对称序列，头尾之和为 1 。

       假如序列中有三个连续元素 x1/y1, x2/y2, x3/y3 ，则有 x2 = x1+x3; y2 = y1+y3;

并且有 x1*y2 – x2*y1 = 1 。这条性质保证构造出来的分式肯定是不可约分式。

      

构造

       从第三个性质我们可以得出它的构造方法。求 N 阶 Farey 序列：

       Procedure make_farey(x1,y1,x2,y2 : integer)

              If x1+x2>N or y1+y2>N then Return

              make_farey(x1,y1,x1+x2,y1+y2)

              inc(total)

              farey [total] = {x1+x2,y1+y2}

make_farey(x1+x2,y1+y2,x2,y2)

       End Procedure.

       就这么简单，当时竟然没想到，唉。

       这个方法适合随机给定 N ，求 Farey 序列。

 

       如果求出连续的 F1,F2,F3,F4…Fn 的话，朴素构造方法更好。

       Procedure make_farey(n : integer)

              farey[1] = {0/1, 1/1}

total = 2

For i=2 to n

       farey [i] = {0/1}

       For j=2 to total

              If farey [i-1][j].denominator + farey [i][j]. denominator = N then

                     farey[i] += { farey [i-1][j] + farey [i][j] }

End if

              farey [i] += {farey [i-1][j]}

End for

End for

End Procedure.

当时思维局限在这个构造算法上，导致超时。

 

所谓的 Stern-Brocot 树，其实已经在第一个构造算法里面隐含了。

Stern-Brocot 扩展了 Farey 序列，它能构造出小于某个分式的所有分式，当然这些分式是无穷的。所以它能从另一方面证明素数是无穷的。

 

查找第K大元素

今天 pku 3374 Cake Share 有较多查询，得预先构造出 Fn 序列。

如果查询对不可预知 Fn 的时候，我们可以简单的改造算法 1 ，当 total 达到 k 时输出后，立即跳出。时间复杂度 O(K) 。 K 最大就是 |Fn| 。

Fn 的序列大小是可以递推出来的，有一个近似公式，可以让我们大致了解下 Fn 的大小程度。

|Fn| = 0.304*N^2 。 如 |F5000| ≈ 7600000 ，实际有 7600459 个。

那可以看出这个算法复杂度是 O(N^2) 的。

黑书上的解决方案是二分加上统计。没有详细描述。

今天比赛时，想了一下，可以如下操作：

规定另一个操作计算出小于给定分式的不可约分式数目，要求 <=O(NlogN) 。

再对 X/N 的 X 进行二分枚举，找出区间 X/N ~ (X+1)/N 。

在枚举出改区间的所有不可约分式，最多也就 N 个。

统计再输出答案。

总过程时间复杂度是 O(N (logN)^2) 的。关键就在于实现 O(NlogN) 的操作。

比赛时一直想不通它的实现，以后有想法再说了。

 

推荐题目

PKU 2478 Farey Sequence

TJU 2798 Farey Sequence

PKU 3374 Cake Share

// PKU 3374 Accepted 30120K 248MS 
#include  < cstdio >
#include  < string >
#include  < algorithm >
using   namespace  std;
const   int  MAX  =   4000000 ;
int  total;
int  n,k;
int  farey[ 2 ][MAX];

void  make_farey_seq( int  x1, int  y1, int  x2,  int  y2)  ... {
    if(x1+x2 > n || y1+y2 > n) return;
    make_farey_seq(x1, y1,x1+x2, y1+y2);
    total ++;
    farey[0][total] = x1+x2;
    farey[1][total] = y1+y2;
    make_farey_seq(x1+x2, y1+y2,x2,y2);
}

int  main()  ... {
    int t;
    scanf("%d %d", &n, &t);
    if(n == 1) ...{
        while(t --) ...{
            scanf("%d", &k);
            if(k == 1) puts("0/1");
            else if(k == 2) puts("1/1");
            else puts("No Solution");
        }
        return 0;
    }
    total = 1;
    farey[0][1] = 0;
    farey[1][1] = 1;
    make_farey_seq(0,1,1,2);
    farey[0][total+1] = 1;
    farey[1][total+1] = 2;
    total ++;
    int all = 2*total;
    while(t --) ...{
        scanf("%d", &k);
        if(k >= all) puts("No Solution");
        else if(k <= total) printf("%d/%d ", farey[0][k], farey[1][k]);
        else if(k > total) printf("%d/%d ", farey[1][all-k] - farey[0][all-k], farey[1][all-k]);
    }
}