N번째 페리 수열의 K번째 항을 빠르게 구하기

문제

0 이상 1 이하인 분모가 $N$이하의 기약분수를 크기 순으로 나열한 수열을 $N$번째 페리 수열 (Farey Sequence)라고 한다. $N$과 $K$가 주어졌을 때, $N$번째 페리 수열의 $K$번째 원소를 구하여라.

 

BOJ 1882와 동일한 문제고, 올해 ICPC 예선에도 출제된 문제이다. 이 글에서는 $N$이 고정이고 $K$가 쿼리로 주어질 때 전처리 $O(N^{\frac{2}{3}})$, 쿼리 $O(\sqrt{N}{\log }^{1.5}N)$ 시간에 이 문제를 해결하는 방법에 대해 다룬다.

 

문제 변형하기

$n$번째 페리 수열의 길이는 $\mathrm{\Theta }(n^2)$임이 널리 알려져있다. (서로소인 쌍의 개수의 비율이 $\frac{6}{{\pi }^2}$에 수렴한다.) 따라서, 페리 수열을 직접 구하는 방식으로는 subquadratic 알고리즘을 얻을 수 없다. 그 대신, 다음과 같은 문제를 풀자.

 

0 이상 1 이하인 실수 $x$가 주어졌을 때, $N$번째 페리 수열의 원소 중 $x$이하인 원소의 개수를 구하여라.

 

위 문제를 $T(N)$의 시간에 풀 수 있다면, 이분탐색을 통해 기존 문제를 $T(N)\log N$ 시간에 해결할 수 있다.

 

포함과 배제의 원리

포함과 배제의 원리를 이용하면, 뫼비우스 함수를 이용해서 $x$이하 원소의 개수를 표현할 수 있다. 다음과 같은 함수 $f$를 생각해 보자.

$$f(i,j)=\{\begin{array}{ll}gcd(i,j)& \text{if}\frac{i}{j}\le x\\ 0& \text{else}\end{array}$$

구하고자 하는 값은 다음과 같이 쓸 수 있다. (대괄호는 아이버슨 괄호로, 조건이 참이면 1, 거짓이면 0이다.)

$$\sum _{1\le i,j\le N}[f(i,j)=1]$$

$\sum _{d|n}\mu (d)=[n=1]$을 적용해서 식을 변형하자. 포함과 배제의 원리를 생각하면 참이라는 것을 쉽게 알 수 있다.

$$\sum _{1\le i,j\le N}\sum _{d|f(i,j)}\mu (d)$$

$d$를 기준으로 식을 다시 쓰면 다음과 같다.

$$\sum _{d=1}^N\mu (d)\times g(d)$$

$\mu (d)$가 더해진 횟수인 $g(d)$를 알아내야 하는데, 이는 $f(i,j)$가 $d$의 배수인 $(i,j)$의 개수와 동일하다. $i=d{i}^{\prime }$, $j=d{j}^{\prime }$으로 놓고 계산해 보면 이는 $\sum _{k=1}^{\lfloor \frac{N}{d}\rfloor }\lfloor kx\rfloor$임을 알 수 있다. 따라서, 구하려는 값은 다음과 같다.

$$\sum _{d=1}^N(\mu (d)\sum _{k=1}^{\lfloor \frac{N}{d}\rfloor }\lfloor kx\rfloor )$$

$\mu (1),\dots ,\mu (N)$은 체를 통해 $O(N)$에 구할 수 있고, $\lfloor kx\rfloor$의 prefix sum도 $O(N)$에 모두 구할 수 있기 때문에 $x$이하 원소의 개수를 $O(N)$에 계산할 수 있다. 따라서, 전체 문제를 $O(N\log N)$에 해결하였다.

 

식 전개가 복잡하다고 느껴지면 포함과 배제의 원리만 사용해서 결론을 바로 도출해 보자. 표 형태로 $(i,j)$를 나열하고 포함과 배제의 원리를 적용하면 직관적으로 식을 얻을 수 있다.

 

Stern-Brocot tree를 이용한 이분탐색

double 등의 실수 자료형을 이용하여 이분탐색을 해도 되지만, 분모가 $N$이하인 유리수 $x$에 대해 이분탐색을 해도 충분하다. Stern-Brocot tree를 이용하면 정수 자료형만으로 이를 구현할 수 있다.

 

Stern-Brocot tree는 $\frac{a}{b}<\frac{a+c}{b+d}<\frac{c}{d}$임을 이용하여 모든 유리수를 생성하는 트리이다. 현재 $l=\frac{a}{b}$가 true, $r=\frac{c}{d}$가 false인 것을 알고 있을 때, $mid=\frac{a+c}{b+d}$의 참/거짓을 판별한 후 $l$과 $r$ 중 하나를 갱신해주면 된다. 이를 Stern-Brocot tree 위에 그려보면 트리를 한 칸씩 내려가는 형태임을 알 수 있다.

 

그러나, 이를 단순히 구현하면 깊이 $N$까지 들어갈 수 있게 되어 시간복잡도가 최악의 경우 $\mathrm{\Theta }(N)$이 된다는 문제점이 있다. 이를 해결하기 위해 $mid=l$을 반복하다가 $mid=r$을 처음으로 하는 지점(또는 그 반대), 즉 "꺾이는 점"만 계산하게 되면 분모가 최대 $N$일 때 "꺾이는 점"은 $O(\log N)$개임을 증명할 수 있다. 이를 이용하여 꺾이는 점을 이분탐색으로 계산하면 $O({\log }^{2}N)$에 구현할 수 있고, exponential search를 이용하면 $O(\log N)$에 구현할 수 있다. 내 코드는 다음과 같다. (f(x)가 참인 최대 유리수를 구하는 코드이다.)

struct Frac{
	int a, b;
	Frac(){}
	Frac(int _a, int _b): a(_a), b(_b) {}
};

int n;
ll k;

int exp_search(Frac L, Frac R, const auto &ok){
	if (!ok(L, R, 1)) return 0;

	int l = 1, r = 2;
	while(ok(L, R, r)) l <<= 1, r <<= 1;

	while(r-l > 1){
		int mid = (l+r)>>1;
		if (ok(L, R, mid)) l = mid;
		else r = mid;
	}

	return l;
}

Frac search(){
	Frac L(0, 1), R(1, 0);

	auto chk_L = [](Frac L, Frac R, int x){
		return L.b*x + R.b <= n && !f(Frac(L.a*x + R.a, L.b*x + R.b));
	};

	auto chk_R = [](Frac L, Frac R, int x){
		return L.b + R.b*x <= n && f(Frac(L.a + R.a*x, L.b + R.b*x));
	};
	
	while(true){
		int x = exp_search(L, R, chk_L);
		if (x) R = Frac(L.a*x + R.a, L.b*x + R.b);

		int y = exp_search(L, R, chk_R);
		if (y) L = Frac(L.a + R.a*y, L.b + R.b*y);

		if (!x && !y) break;
	}

	return L;
}

 

 

Sublinear 최적화

앞서 구했던 식의 항을 모두 풀어쓰면 다음과 같이 쓸 수 있다.

$$\sum _{ij\le N}\mu (i)\lfloor jx\rfloor$$

$\mu (i)$ 대신 $\lfloor jx\rfloor$를 고정하고 항을 묶으면 다음과 같다. ($M(x)$는 메르텐스 함수이며, $M(x)=\sum _{i=1}^x\mu (x)$로 정의된다.)

$$\sum _{d=1}^N\lfloor dx\rfloor M\left(\frac{N}{d}\right)$$

$\frac{N}{d}$로 가능한 값이 루트개임에 주목하자. $\frac{N}{d}>B$면 그냥 계산하고, $\frac{N}{d}\le B$면 $\frac{N}{d}$값이 같은 $d$를 묶어서 한 번에 계산하자. 이때, $x=\frac{p}{q}$인 $N$이하의 양의 정수 $p$, $q$가 존재하는 경우만 생각하면 충분하기 때문에, $\sum _{i=l}^r\lfloor \frac{pi}{q}\rfloor$꼴의 식을 빠르게 계산할 수 있으면 되고, 이는 $O(\log N)$에 계산할 수 있음이 잘 알려져있다. (ACL 또는 yosupo judge에서 찾아볼 수 있다.)

따라서, $M\left(\frac{N}{d}\right)$의 값을 모두 알고 있으면 $O(N/B+B\log N)$ 시간에 해결할 수 있고, $B=\sqrt{\frac{N}{\log N}}$일 때 $O(\sqrt{N\log N})$으로 최적이다.

 

$M\left(\frac{N}{d}\right)$의 빠른 전처리 방법은 Xudyh's sieve 또는 Mertens trick으로 알려져 있다. $O(N^{\frac{2}{3}})$에 계산할 수 있다. (https://xy-plane.tistory.com/19, https://codeforces.com/blog/entry/54150, https://rkm0959.tistory.com/190)

 

따라서, 전처리 $O(N^{\frac{2}{3}})$, 쿼리 $O(\sqrt{N}{\log }^{1.5}N)$ 시간에 문제를 해결하였다.