SB树和Farey序列

csuzhucong

已于 2024-05-24 23:09:12 修改

阅读量1.6k

点赞数

文章标签：算法

于 2021-04-05 11:49:33 首次发布

本文链接：https://blog.csdn.net/nameofcsdn/article/details/115440090

版权

一，Farey序列和Stern-Brocot树

POJ 2478 Farey Sequence

HDU 4556 Stern-Brocot Tree

Fn其他性质分析

UVA 12995 Farey Sequence

UVALive 7098 Farey Sums

Fn任意项计算

UVA 10408 Farey sequences

HDU 2432 Farey Sequence Again

一，Farey序列和Stern-Brocot树

1，pair序列

考虑一系列由pair组成的序列：

$A_0=\{\frac{0}{1},\frac{1}{0}\}$

$A_{n+1}$ 是 $A_{n}$ 在每2个数m/n和m'/n'之间插入一个数(m+m')/(n+n')得到，

于是， $A_1=\{\frac{0}{1},\frac{1}{1},\frac{1}{0}\},A_2=\{\frac{0}{1},\frac{1}{2},\frac{1}{1},\frac{2}{1},\frac{1}{0}\}$

显然，An中有2^n+1个元素。

性质：

（1）任意一个数列中，任取相邻的2个pair，m/n和m'/n'，有m'n-mn'=1

证明：用数学归纳法，显然成立。

（2）对于任意一个数列的任意一个pair，m/n，有gcd(m,n)=1，即这是最简形式。

根据（1）即可证得

（3）任意一个数列中，任取相邻的2个pair，m/n和m'/n'，有m/n<m'/n'

证明：用数学归纳法，显然成立。

（4）任意一个数列中除了 $\frac{0}{1},\frac{1}{0}$ （0和正无穷）之外，剩下的都是最简形式的正有理数，且递增出现。

（5）对于任意正有理数，一定出现在某个 $A_{n}$ 中。

（6）对于任意连续3项，m/n,a/b,m'/n'，都有，a/b是(m+m')/(n+n')的最简形式。如A2中的1/2,1/1,2/1

2，Stern-Brocot树

可以把上述pair序列之间的关系，表示成一棵二叉树：

其中每一行对应的就是 $A_{n+1}$ 比 $A_{n}$ 多出来的那些元素。

取二叉树前n行，加上2个触角，按照从左往右的顺序，即可得到An

性质：

3，Farey序列

$F_n$ 是把集合 $\{\frac{a}{b}|(a,b)=1,0\leq a\leq b\leq n\}$ 按照递增顺序来排列得到的序列。

F1 = {0/1,1/1}
F2 = {0/1,1/2,1/1}
F3 = {0/1,1/3, 1/2, 2/3,1/1}
F4 = {0/1,1/4, 1/3, 1/2, 2/3, 3/4,1/1}
F5 = {0/1,1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5,1/1}
F6= {0/1,1/6,1/5,1/4,1/3,2/5,1/2,3/5,2/3,3/4,4/5,5/6,1/1}

性质：

（1）从0开始，到1结束，递增序列

（2）0到1之间的最简分数一定出现在某个Fn中

（3）Fn是An的左边一半的子数列

（4）Fn中任取相邻的2个数，m/n和m'/n'，有m'n-mn'=1

（5）由 $F_{n-1}$ 中插入若干个x/n可以得到 $F_{n}$ ，其中x代表1到n-1中和n互质的所有数。

（6）对于Fn中任意连续3项，a/b,c/d,e/f，都有，c/d是(a+e)/(b+f)的最简形式。更准确的，(b+f)/d=(b+n)/d，左边是刚好整除，右边是带余除法，舍弃余数。这一条可作为任意Fn序列根据前2项推出后面所有项的方法。

4，SB树中的路径计算

5，Farey序列计算

（1）根据相邻两项算出下一项

//Fn中根据相邻两项a0/b0和a1/b1算出后一项a2/b2
void fareyNext(int n, int a0,int b0,int a1,int b1,int &a2,int &b2)
{
	int k = (n + b0) / b1;
	a2 = k * a1 - a0;
	b2 = k * b1 - b0;
}

这个递推式利用了Farey序列的性质，相邻三项之间的关系。

（2）根据任意项算出后一项

//求任意x,y，ax-by=gcd(a,b)
long long gcd(long long a, long long b, long long &x, long long &y)
{
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	auto g = gcd(b, a%b, x, y);
	auto t = -a / b * y - x;
	x = -y, y = t;
	return g;
}

//Fn中根据任意项a/b算出后一项c/d
void fareyNext(int n, int a, int b, int &c, int &d)
{
	long long x, y;
	auto g = gcd(b, a, x, y);
	while (y <= n)x += a / g, y += b / g;
	c = x - a / g, d = y - b / g;
}

这个递推式只利用了相邻两项之间的关系，但是需要调用拓展欧几里得算法。

二，OJ实战

Fn中的元素个数

POJ 2478 Farey Sequence

题目：

Description

The Farey Sequence Fn for any integer n with n >= 2 is the set of irreducible rational numbers a/b with 0 < a < b <= n and gcd(a,b) = 1 arranged in increasing order. The first few are
F2 = {1/2}
F3 = {1/3, 1/2, 2/3}
F4 = {1/4, 1/3, 1/2, 2/3, 3/4}
F5 = {1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5}

You task is to calculate the number of terms in the Farey sequence Fn.
Input

There are several test cases. Each test case has only one line, which contains a positive integer n (2 <= n <= 10 6). There are no blank lines between cases. A line with a single 0 terminates the input.
Output

For each test case, you should output one line, which contains N(n) ---- the number of terms in the Farey sequence Fn.
Sample Input

2
3
4
5
0
Sample Output

1
3
5
9

题意：计算Fn中的元素个数

这个题目无非就是求欧拉函数的前n项和呗

因为查询不知道有多少个，所以要先打表。

代码：

#include<iostream>
#include<stdio.h>
using namespace std;
 
long long phi[1000001];		//前i个数的phi之和
 
void get_phi()
{
	for (int i = 1; i <= 1000000; i++)phi[i] = i;
	for (int i = 2; i <= 1000000; i++)
	{
	if (phi[i] == i)for (int j = i; j <= 1000000; j += i)phi[j] = phi[j] / i*(i - 1);
	phi[i] += phi[i - 1];	//去掉这一行就是求欧拉函数
	}
}
 
int main()
{
	get_phi();
	int n;
	while (scanf("%d",&n))
	{
		if (n == 0)break;
		printf("%llu\n", phi[n] - 1);
	}
	return 0;
}

HDU 4556 Stern-Brocot Tree

Description

　　
　　
　　上图是一棵Stern-Brocot树，其生成规则如下：
　　从第1行到第n行，每行相邻两数a/b和c/d，产生中间数(a+c)/(b+d)，置于下一行中。将一行的分数（包括0/1,1/0），进行约分简化，则每一行（包括0/1,1/0,1/1），不会出现两个相同的分数。若分子或者分母大于n，则去掉该分数，将剩下的分数，从小到大排序，得到数列F。
　　现在请您编程计算第n行的数列F的个数。
Input

　　输入包含多组测试用例，每组输入数据是一个正整数n（n<=1000000）。
Output

　　对于每组的测试数据n，请输出第n行的数列F的个数。
Sample Input

1
2
4
6
Sample Output

3
5
13
25

这个题目就是想说明，SB树和Farey序列的关系。

代码就几乎不用再写了，直接把POJ - 2478 Farey Sequence的代码略改即可。

代码：

#include<iostream>
#include<stdio.h>
using namespace std;

long long phi[1000001];	

void get_phi()
{
	for (int i = 1; i <= 1000000; i++)phi[i] = i;
	for (int i = 2; i <= 1000000; i++)
	{
	if (phi[i] == i)for (int j = i; j <= 1000000; j += i)phi[j] = phi[j] / i*(i - 1);
	phi[i] += phi[i - 1];	
	}
}

int main()
{
	get_phi();
	int n;
	while (scanf("%d",&n)!=-1)printf("%llu\n", phi[n]*2+ 1);
	return 0;
}

Fn其他性质分析

UVA 12995 Farey Sequence

题目：

题目的意思是说，在Fn里面，有多少对(mi/ni,mj/nj)使得，他们在Fn里面不相邻，但是满足mj*ni-mi*nj=1

首先注意到，满足mj*ni-mi*nj=1等价于(mi/ni,nj/nj)在某个Farey里面是相邻的。

比如，2/5和1/2虽然在F7里面不相邻，但是在F5里面是相邻的。这个性质，很重要。

根据这个性质，假设F(n-1)里面有g(n-1)组满足mj*ni-mi*nj=1的二元组，F(n)的长度是l(n)

那么答案可以表示为，g(n-1)+2*phi(n)-l(n)+1，其中phi是欧拉phi函数

然后分别计算，g(n-1)=1+（ phi(2)+phi(3)+......+phi(n-1) ）*2

l(n)=2+phi(2)+phi(3)+......+phi(n-1) +phi(n)

所以答案是phi(2)+phi(3)+......+phi(n-1) +phi(n)，即l(n)-2

真巧，POJ - 2478 就是直接要你求出l(n)-2

所以，这个题目，把 POJ - 2478的代码直接复制过来就行了，一个字符都不需要修改。

UVALive 7098 Farey Sums

题目：

这个题目考的就是1个对称性。

在a和b之间插入a+b，在b和a之间插入a+b

那么a/b+b/a就变成了a/(a+b)+(a+b)/b+b/(a+b)+(a+b)/a=a/b+b/a+3

增量是3，整个序列的增量是若干个3的和，这样的3的个数是n的欧拉函数的一半。

所以表达式很容易求出来，先求出前n个数的欧拉函数之和phi[n]，然后答案便是（phi[n] * 3 -1）/ 2

代码：

#include<iostream>
#include<stdio.h>
using namespace std;

int phi[10001];		

void get_phi()
{
	for (int i = 1; i <= 10000; i++)phi[i] = i;
	for (int i = 2; i <= 10000; i++)
	{
		if (phi[i] == i)for (int j = i; j <= 10000; j += i)phi[j] = phi[j] / i*(i - 1);
		phi[i] += phi[i - 1];	
	}
}

int main()
{
	get_phi();
	int p, n;
	scanf("%d", &p);
	for (int i = 1; i <= p; i++)
	{
		scanf("%d%d", &n, &n);
		printf("%d %d/2\n", i, phi[n] * 3 -1);
	}
	return 0;
}

Fn任意项计算

UVA 10408 Farey sequences

题意：

求出Fn的第k项

思路一：一阶递推式

#include<cstdio>

//求任意x,y，ax-by=gcd(a,b)
long long gcd(long long a, long long b, long long &x, long long &y)
{
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	auto g = gcd(b, a%b, x, y);
	auto t = -a / b * y - x;
	x = -y, y = t;
	return g;
}

//Fn中根据任意项a/b算出后一项c/d
void fareyNext(int n, int a, int b, int &c, int &d)
{
	long long x, y;
	auto g = gcd(b, a, x, y);
	while (y <= n)x += a / g, y += b / g;
	c = x - a / g, d = y - b / g;
}

int main()
{
	int n, k;
	while (~scanf("%d%d", &n, &k)) {
		int a0 = 0, b0 = 1, a1, b1;
		while (k--) {
			fareyNext(n, a0, b0, a1, b1);
			a0 = a1, b0 = b1;
		}
		printf("%d/%d\n", a1, b1);
	}
	return 0;
}

思路二：二阶递推式

#include<cstdio>

//Fn中根据相邻两项a0/b0和a1/b1算出后一项a2/b2
void fareyNext(int n, int a0,int b0,int a1,int b1,int &a2,int &b2)
{
	int k = (n + b0) / b1;
	a2 = k * a1 - a0;
	b2 = k * b1 - b0;
}

int main()
{
	int n, k;
	while (~scanf("%d%d", &n, &k)) {
		int a0 = 0, b0 = 1, a1 = 1, b1 = n, a2, b2;
		while (--k) {
			fareyNext(n, a0, b0, a1, b1, a2, b2);
			a0 = a1, a1 = a2;
			b0 = b1, b1 = b2;
		}
		printf("%d/%d\n", a1, b1);
	}
	return 0;
}

HDU 2432 Farey Sequence Again

The Farey sequence Fn for any positive integer n is the set of irreducible rational numbers a/b with 0<a<b<=n and (a, b) = 1 arranged in increasing order. Here (a, b) mean the greatest common divisor of a and b. For example:
            F2 = {1/2}
            F3 = {1/3, 1/2, 2/3}
      Given two positive integers N and K, with K less than N, you need to find out the K-th smallest element of the Farey sequence FN.

Input

The first line of input is the number of test case T, 1<=T<=1000. Then each test case contains two positive integers N and K. 1<=K<N<=10^9.

Output

For each test case output the Kth smallest element of the Farey sequence FN in a single line.

Sample

Inputcopy	Outputcopy
3 2 1 100 68 101 69	1/2 2/83 1/42

思路：

n太大了，不能用直接递推的方法。

注意到k远远小于Fn的元素个数，所以不难推出，Fn的前k项中，分子只有3种情况，即1,2,3（除了第一个是0/1）

于是，我们只需要寻找x，使得不超过6/x的数的数目刚好是k

这题的边界条件太难搞了。。。


#include <iostream>
#include <algorithm>
using namespace std;

// class Bsearch 模板

class FS :public Bsearch<long long>{
public:
	FS(int n,int k) {
		this->n = n, this->k = k;
	}
	bool isOk(long long x) //若isOk(x)且!isOk(y)则必有y<x
	{
		long long s = 1, z = 0;
		s += max(z, num((x + 5) / 6, n, 1));
		s += max(z, num((x + 2) / 3, n, 2));
		s += max(z, num((x + 1) / 2, n, 3));
		return s <= k;
	}
	long long num(long long low, long long high, int d)
	{
		if (d == 1)return high - low + 1;
		low = (low - 1) / d * (d - 1) + (low - 1) % d;
		high = high / d * (d - 1) + high % d;
		return high - low;
	}
	int n, k;
};

int main()
{
	int n, k;
	scanf("%d", &n);
	while (~scanf("%d%d", &n, &k)) {
		long long x = FS(n, k + 1).find(6, 9876543210);
		long long a = (x + 5) / 6;
		long long b = (x + 2) / 3;
		long long c = (x + 1) / 2;
		pair<int, long long> ans = make_pair(1, a);
		if (b <= n && b < a * 2)ans = make_pair(2, b);
		if (c <= n && c * ans.first < ans.second * 3)ans = make_pair(3, c);
		cout << ans.first << "/" << ans.second << endl;
	}
	return 0;
}