<div id="readability-page-1" class="page"><div>Bostan Mori 알고리즘은 2020년 8월 Alin Bostan과 Ryuhei Mori가 작성한 <a href="https://arxiv.org/abs/2008.08822" target="_blank" rel="noopener">이 논문</a>에 소개되어 있는 선형점화식을 가지는 수열의 $N$ 번째 항을 빠르게 구하는 알고리즘이다. Bostan Mori 알고리즘 이외에 선형점화식을 가지는 수열의 $N$ 번째 항을 빠르게 구하는 방법은 <a href="https://blog.myungwoo.kr/148" target="_blank" rel="noopener">이 글</a>을 참고하자.

$D_0, D_1, \cdots, D_{k-1}$과 $c_1, c_2, \cdots, c_k$가 주어졌을 때, $i \ge k$인 $D_i$를 다음과 같은 선형점화식으로 구할 수 있다고 하자.
$$D_i = \sum_{j=1}^{k}{c_jD_{i-j}} = c_1D_{i-1} + c_2D_{i-2} + \cdots + c_kD_{i-k}$$
이러한 선형점화식이 주어졌을 때, $D_N$을 Bostan Mori 알고리즘을 통해 구하는 방법을 알아보자.

$D$의 각 원소를 계수로 가지는 생성함수 $F$와 점화식을 바탕으로 만들어진 다항함수 $Q$가 있다.
$$F(x) = \sum_{n=0}^{\infty}{D_nx^n}\\Q(x) = 1-c_1x-c_2x^2-\cdots-c_kx^k$$
여기서, $F(x) = \frac{P(x)}{Q(x)}$를 만족하는 다항함수 $P$가 존재한다는 것과 $P$의 차수가 $k$ 보다 작다는 것을 보일 수 있다. $P(x) = F(x)Q(x)$이고, $i \ge k$인 $i$에 대해 $P(x)$의 $x^i$ 항의 계수를 구해보면 $D_i - c_1D_{i-1} - c_2D_{i-2} - \cdots - c_kD_{i-k} = 0$이다. 즉, $P(x)$의 차수는 $k$ 보다 작다.

구하고자 하는 것은 $N$ 번째 항인 $D_N$이다. 이는 $F(x) = \frac{P(x)}{Q(x)}$의 $x^N$ 항의 계수이다. 이를 다음과 같은 수식으로 표현하자.
$$\left[x^N\right]\frac{P(x)}{Q(x)}$$
위 분수의 분자와 분모 모두에 $Q(-x)$를 곱해보자. 즉, $\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)}$이다.
이때, $Q(x)$의 차수는 $k$이다. $Q(x)Q(-x)$는 짝함수이기 때문에 $V(x^2)$으로 나타낼 수 있으며 $V$의 차수는 $k$ 이하다. $P(x)Q(-x) = U(x)$로 쓰자. $P(x)$의 차수는 $k$ 보다 작기 때문에 $U(x)$의 차수는 $2k$ 보다 작다. $U(x)$를 홀수차항과 짝수차항으로 $U(x) = U_e(x^2) + xU_o(x^2)$처럼 나누자. 여기서 $U_e(x)$와 $U_o(x)$의 차수는 $k$ 보다 작다. 이를 통해 식을 다음처럼 정리할 수 있다.
$$\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)} = \frac{U_e(x^2)}{V(x^2)} + x\frac{U_o(x^2)}{V(x^2)}$$
유리함수 $P/Q$의 $x^N$ 항의 계수는 다음과 같이 정리된다.
$$\left[ x^N \right]\frac{P(x)}{Q(x)} = \begin{cases} \left[ x^{\frac{N}{2}} \right] \frac{U_e(x)}{V(x)}, &amp; \text{if $N$ is even,}\\ \left[ x^{\frac{N-1}{2}} \right] \frac{U_o(x)}{V(x)}, &amp; \text{otherwise} &nbsp;\end{cases}$$
즉, 분자의 차수가 $k$보다 작고, 분모의 차수가 $k$ 이하인 유리함수 $P/Q$의 $x^N$ 항의 계수를 구하기 위해 같은 조건을 가지고 있는 또다른 유리함수의 $x^{\lfloor N/2 \rfloor}$의 계수를 구하도록 차원을 줄였다. 이처럼 두 번의 다항식 곱셈으로 구하고자 하는 항의 차수를 반으로 줄일 수 있다. 이를 상수항의 계수를 구할 때까지 반복한다. 유리함수 $P/Q$의 상수 계수는 $P(0)/Q(0)$으로 바로 구할 수 있다.

일반적인 방법으로 다항식의 곱셈을 구현하면 $O(k^2)$ 시간복잡도가 되고, FFT나 NTT를 통해 다항식의 곱셈을 구현하면 $O(k \lg k)$ 시간복잡도가 된다. 따라서, 전체 시간복잡도는 $O(k^2 \lg N)$ 혹은 $O(k \lg k \lg N)$이다.

연습 문제
<a href="https://boj.kr/15572" target="_blank" rel="noopener">https://boj.kr/15572</a>
문제의 MOD 값이 1999이므로 NTT를 사용할 수 없다. $k$의 최댓값이 1000으로 크지 않으므로 $O(k^2 \lg N)$ 시간복잡도로 해결할 수 있다.
<figure id="og_1669791672456" contenteditable="false" data-ke-type="opengraph" data-ke-align="alignCenter" data-og-type="website" data-og-title="15572번: 블록 4" data-og-description="여러 가지 블록들을 이용하여 직사각형 모양을 만들려고 한다. 우리에게는&nbsp;1 × N&nbsp;블록,&nbsp;2 × N&nbsp;블록, ...,&nbsp;N × N&nbsp;블록이 무한하게 있다. 이 블록들을 사용하여&nbsp;N × M&nbsp;모양을 만들고 " data-og-host="www.acmicpc.net" data-og-source-url="https://boj.kr/15572" data-og-url="https://www.acmicpc.net/problem/15572" data-og-image="https://scrap.kakaocdn.net/dn/bXjQVj/hyQJAXnYOB/kwhsNZ1htBkIY2lnQva7Ik/img.png?width=2834&amp;height=1480&amp;face=0_0_2834_1480"><a href="https://boj.kr/15572" target="_blank" rel="noopener" data-source-url="https://boj.kr/15572">

<div>
15572번: 블록 4
여러 가지 블록들을 이용하여 직사각형 모양을 만들려고 한다. 우리에게는&nbsp;1 × N&nbsp;블록,&nbsp;2 × N&nbsp;블록, ...,&nbsp;N × N&nbsp;블록이 무한하게 있다. 이 블록들을 사용하여&nbsp;N × M&nbsp;모양을 만들고
www.acmicpc.net
</div>
</a><img src="https://scrap.kakaocdn.net/dn/bXjQVj/hyQJAXnYOB/kwhsNZ1htBkIY2lnQva7Ik/img.png?width=2834&amp;height=1480&amp;face=0_0_2834_1480"></figure>
<a href="https://boj.kr/13725" target="_blank" rel="noopener">https://boj.kr/13725</a>
MOD 값이 $104857601 = 25 \times 2^{22}+1$ 로 NTT를 쓸 수 있다. 원시근은 2부터 차례대로 구해보면 되는데, 3이 원시근이 된다. $k$의 최댓값이 30000으로 크기 때문에 $O(k \lg k \lg N)$ 시간복잡도로 문제를 해결해야 한다.
<figure id="og_1669791817117" contenteditable="false" data-ke-type="opengraph" data-ke-align="alignCenter" data-og-type="website" data-og-title="13725번: RNG" data-og-description="첫째 줄에 k와 N (1 ≤ k ≤ 30,000, 1 ≤ N ≤ 1018)이 주어진다. 둘째 줄에는 A1, A2, ..., Ak가 셋째 줄에는 C1, C2, ..., Ck가 주어진다. (0 ≤ Ai, Ci < 104857601)" data-og-host="www.acmicpc.net" data-og-source-url="https://boj.kr/13725" data-og-url="https://www.acmicpc.net/problem/13725" data-og-image="https://scrap.kakaocdn.net/dn/bBRPcL/hyQJA4chmU/Dnyz5lOG9CmvrlhoJdt5m1/img.png?width=2834&amp;height=1480&amp;face=0_0_2834_1480"><a href="https://boj.kr/13725" target="_blank" rel="noopener" data-source-url="https://boj.kr/13725">

<div>
13725번: RNG
첫째 줄에 k와 N (1 ≤ k ≤ 30,000, 1 ≤ N ≤ 1018)이 주어진다. 둘째 줄에는 A1, A2, ..., Ak가 셋째 줄에는 C1, C2, ..., Ck가 주어진다. (0 ≤ Ai, Ci &lt; 104857601)
www.acmicpc.net
</div>
</a><img src="https://scrap.kakaocdn.net/dn/bBRPcL/hyQJA4chmU/Dnyz5lOG9CmvrlhoJdt5m1/img.png?width=2834&amp;height=1480&amp;face=0_0_2834_1480"></figure>
아래 코드는 BOJ 13725 RNG 문제를 해결하는 코드이며, mod int와 NTT를 이용한 다항식 템플릿을 포함한 코드다.
</div></div>

PS 이야기

Bostan Mori 알고리즘은 2020년 8월 Alin Bostan과 Ryuhei Mori가 작성한 이 논문에 소개되어 있는 선형점화식을 가지는 수열의 $N$ 번째 항을 빠르게 구하는 알고리즘이다. Bostan Mori 알고리즘 이외에 선형점화식을 가지는 수열의 $N$ 번째 항을 빠르게 구하는 방법은 이 글을 참고하자.
 
$D_0, D_1, \cdots, D_{k-1}$과 $c_1, c_2, \cdots, c_k$가 주어졌을 때, $i \ge k$인 $D_i$를 다음과 같은 선형점화식으로 구할 수 있다고 하자.
$$D_i = \sum_{j=1}^{k}{c_jD_{i-j}} = c_1D_{i-1} + c_2D_{i-2} + \cdots + c_kD_{i-k}$$
이러한 선형점화식이 주어졌을 때, $D_N$을 Bostan Mori 알고리즘을 통해 구하는 방법을 알아보자.
 
$D$의 각 원소를 계수로 가지는 생성함수 $F$와 점화식을 바탕으로 만들어진 다항함수 $Q$가 있다.
$$F(x) = \sum_{n=0}^{\infty}{D_nx^n}\\Q(x) = 1-c_1x-c_2x^2-\cdots-c_kx^k$$
여기서, $F(x) = \frac{P(x)}{Q(x)}$를 만족하는 다항함수 $P$가 존재한다는 것과 $P$의 차수가 $k$ 보다 작다는 것을 보일 수 있다. $P(x) = F(x)Q(x)$이고, $i \ge k$인 $i$에 대해 $P(x)$의 $x^i$ 항의 계수를 구해보면 $D_i - c_1D_{i-1} - c_2D_{i-2} - \cdots - c_kD_{i-k} = 0$이다. 즉, $P(x)$의 차수는 $k$ 보다 작다.
 
구하고자 하는 것은 $N$ 번째 항인 $D_N$이다. 이는 $F(x) = \frac{P(x)}{Q(x)}$의 $x^N$ 항의 계수이다. 이를 다음과 같은 수식으로 표현하자.
$$\left[x^N\right]\frac{P(x)}{Q(x)}$$
위 분수의 분자와 분모 모두에 $Q(-x)$를 곱해보자. 즉, $\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)}$이다.
이때, $Q(x)$의 차수는 $k$이다. $Q(x)Q(-x)$는 짝함수이기 때문에 $V(x^2)$으로 나타낼 수 있으며 $V$의 차수는 $k$ 이하다. $P(x)Q(-x) = U(x)$로 쓰자. $P(x)$의 차수는 $k$ 보다 작기 때문에 $U(x)$의 차수는 $2k$ 보다 작다. $U(x)$를 홀수차항과 짝수차항으로 $U(x) = U_e(x^2) + xU_o(x^2)$처럼 나누자. 여기서 $U_e(x)$와 $U_o(x)$의 차수는 $k$ 보다 작다. 이를 통해 식을 다음처럼 정리할 수 있다.
$$\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)} = \frac{U_e(x^2)}{V(x^2)} + x\frac{U_o(x^2)}{V(x^2)}$$
유리함수 $P/Q$의 $x^N$ 항의 계수는 다음과 같이 정리된다.
$$\left[ x^N \right]\frac{P(x)}{Q(x)} = \begin{cases}
\left[ x^{\frac{N}{2}} \right] \frac{U_e(x)}{V(x)}, & \text{if $N$ is even,}\\
\left[ x^{\frac{N-1}{2}} \right] \frac{U_o(x)}{V(x)}, & \text{otherwise}
 \end{cases}$$
즉, 분자의 차수가 $k$보다 작고, 분모의 차수가 $k$ 이하인 유리함수 $P/Q$의 $x^N$ 항의 계수를 구하기 위해 같은 조건을 가지고 있는 또다른 유리함수의 $x^{\lfloor N/2 \rfloor}$의 계수를 구하도록 차원을 줄였다. 이처럼 두 번의 다항식 곱셈으로 구하고자 하는 항의 차수를 반으로 줄일 수 있다. 이를 상수항의 계수를 구할 때까지 반복한다. 유리함수 $P/Q$의 상수 계수는 $P(0)/Q(0)$으로 바로 구할 수 있다.
 
일반적인 방법으로 다항식의 곱셈을 구현하면 $O(k^2)$ 시간복잡도가 되고, FFT나 NTT를 통해 다항식의 곱셈을 구현하면 $O(k \lg k)$ 시간복잡도가 된다. 따라서, 전체 시간복잡도는 $O(k^2 \lg N)$ 혹은 $O(k \lg k \lg N)$이다.
 
연습 문제
https://boj.kr/15572
문제의 MOD 값이 1999이므로 NTT를 사용할 수 없다. $k$의 최댓값이 1000으로 크지 않으므로 $O(k^2 \lg N)$ 시간복잡도로 해결할 수 있다.
https://boj.kr/13725
MOD 값이 $104857601 = 25 \times 2^{22}+1$ 로 NTT를 쓸 수 있다. 원시근은 2부터 차례대로 구해보면 되는데, 3이 원시근이 된다. $k$의 최댓값이 30000으로 크기 때문에 $O(k \lg k \lg N)$ 시간복잡도로 문제를 해결해야 한다.
아래 코드는 BOJ 13725 RNG 문제를 해결하는 코드이며, mod int와 NTT를 이용한 다항식 템플릿을 포함한 코드다.
더보기
#include <bits/stdc++.h>
using namespace std;
using lld = long long;
template <int MOD>
struct MINT{
	int v{};
	MINT() = default;
	constexpr MINT(lld v_){
		v = -MOD <= v_ && v_ < MOD ? v_ : v_ % MOD;
		if (v < 0) v += MOD;
	}
	friend istream& operator >> (istream& is, MINT& a){ lld t; is >> t; a = MINT(t); return is; }
	friend ostream& operator << (ostream& os, const MINT& a){ return os << a.v; }
	bool operator == (const MINT& ot)const{ return v == ot.v; }
	bool operator != (const MINT& ot)const{ return v != ot.v; }
	static constexpr MINT pow(MINT a, lld b){
		MINT ret = 1;
		for (;b;b>>=1,a*=a) if (b&1) ret *= a;
		return ret;
	}
	static constexpr MINT inv(const MINT& t){ return pow(t, MOD-2); }
	constexpr MINT operator - ()const{ return MINT(-v); }
	constexpr MINT& operator += (const MINT& ot){ if ((v += ot.v) >= MOD) v -= MOD; return *this; }
	constexpr MINT& operator -= (const MINT& ot){ if ((v -= ot.v) < 0) v += MOD; return *this; }
	constexpr MINT& operator *= (const MINT& ot){ v = (lld)v*ot.v%MOD; return *this; }
	constexpr MINT& operator /= (const MINT& ot){ *this *= inv(ot); return *this; }
	constexpr MINT operator + (const MINT& ot)const{ return MINT(*this) += ot; }
	constexpr MINT operator - (const MINT& ot)const{ return MINT(*this) -= ot; }
	constexpr MINT operator * (const MINT& ot)const{ return MINT(*this) *= ot; }
	constexpr MINT operator / (const MINT& ot)const{ return MINT(*this) /= ot; }
	constexpr operator int32_t ()const{ return v; }
	constexpr operator int64_t ()const{ return v; }
};
template <int prime = (119<<23)+1, int primitive_root = 3>
struct FFT{
	using T = MINT<prime>;
	static constexpr int P = prime;
	static constexpr int PR = primitive_root;
	static constexpr int PRI = T::inv(PR);
	static constexpr int MSZ = (P - 1) & (1 - P);
	static constexpr int W = T::pow(PR, (P - 1) / MSZ);
	static constexpr int WI = T::pow(PRI, (P - 1) / MSZ);
	static void fft(int sz, vector<T>& a, bool inv){
		T w = T::pow(inv ? WI : W, MSZ / sz);
		for (int ssz=sz>>1;ssz;ssz>>=1){
			for (int i=0;i<sz;i+=ssz<<1){
				T wt = 1;
				for (int j=0;j<ssz;++j){
					T& lft = a[i + j];
					T& rgh = a[i + j + ssz];
					T lftold = lft;
					lft += rgh;
					rgh = (lftold - rgh) * wt;
					wt *= w;
				}
			}
			w *= w;
		}
		if (inv){
			T invsz = T::inv(sz);
			for (int i=0;i<sz;i++) a[i] *= invsz;
		}
		for (int j=0,i=1;i<sz;i++){
			int dg = sz>>1;
			for (;j&dg;dg>>=1) j ^= dg;
			j ^= dg;
			if (i < j) swap(a[i], a[j]);
		}
	}
	static vector<T> convolution(const vector<T>& a, const vector<T>& b){
		int sz = a.size() + b.size();
		int fftsz = sz;
		while (fftsz & (fftsz - 1)) fftsz += fftsz & -fftsz;
		vector<T> av(fftsz);
		vector<T> bv(fftsz);
		for (int i=0,ilim=size(a);i<ilim;i++) av[i] = a[i];
		for (int i=0,ilim=size(b);i<ilim;i++) bv[i] = b[i];
		fft(fftsz, av, false);
		fft(fftsz, bv, false);
		for (int i=0;i<fftsz;i++) av[i] *= bv[i];
		fft(fftsz, av, true);
		return av;
	}
};
// alternatives
using FFT1 = FFT<1'012'924'417, 5>; // (483<<21)+1
using FFT2 = FFT<1'004'535'809, 3>; // (479<<21)+1
template <int prime = (119<<23)+1, int primitive_root = 3>
struct Poly{
	using T = MINT<prime>;
	using fft_t = FFT<prime, primitive_root>;
	vector<T> a;
	Poly() = default;
	Poly(T a0): a(1, a0){ normalize(); }
	Poly(const vector<T>& a): a(a) { normalize(); }
	int size()const{ return a.size(); }
	int deg()const{ return a.size()-1; }
	T& operator [](int i){ return a[i]; }
	typename vector<T>::const_iterator begin()const{ return a.begin(); }
	typename vector<T>::const_iterator end()const{ return a.end(); }
	void push_back(const T& v){ a.push_back(v); }
	template<typename... Args> void emplace_back(Args&&... args){ a.emplace_back(args...); }
	void pop_back(){ a.pop_back(); }
	void resize(int n){ a.resize(n); }
	void normalize(){ while (!a.empty() && a.back() == T(0)) a.pop_back(); }
	Poly reversed()const{ auto b = a; reverse(b.begin(), b.end()); return b; }
	Poly trim(int n)const{ return vector<T>(a.begin(), a.begin() + min(n, size())); }
	Poly inv(int n){
		Poly q(T(1) / a[0]);
		for (int i=1;i<n;i<<=1){
			Poly p = Poly(2) - q*trim(i<<1);
			q = (p * q).trim(i<<1);
		}
		return q.trim(n);
	}
	Poly& operator *= (const T& x){ for (T& v: a) v *= x; normalize(); return *this; }
	Poly& operator /= (const T& x){ return *this *= T(1) / x; }
	Poly& operator += (const Poly& ot){
		if (size() < ot.size()) a.resize(ot.size());
		for (int i=0;i<ot.size();i++) a[i] += ot.a[i];
		normalize();
		return *this;
	}
	Poly& operator -= (const Poly& ot){
		if (size() < ot.size()) a.resize(ot.size());
		for (int i=0;i<ot.size();i++) a[i] -= ot.a[i];
		normalize();
		return *this;
	}
	Poly& operator *= (const Poly& ot){
		*this = fft_t::convolution(a, ot.a);
		normalize();
		return *this;
	}
	Poly& operator /= (const Poly& ot){
		if (deg() < ot.deg()) return *this = Poly();
		int sz = deg()-ot.deg()+1;
		Poly ra = reversed().trim(sz), rb = ot.reversed().trim(sz).inv(sz);
		*this = (ra*rb).trim(sz);
		while (size() < sz) emplace_back();
		reverse(a.begin(), a.end());
		normalize();
		return *this;
	}
	Poly& operator %= (const Poly& ot){
		if (deg() < ot.deg()) return *this;
		Poly tmp = *this; tmp /= ot; tmp *= ot;
		*this -= tmp;
		normalize();
		return *this;
	}
	Poly operator * (const T& x)const{ return Poly(*this) *= x; }
	Poly operator / (const T& x)const{ return Poly(*this) /= x; }
	Poly operator + (const Poly& ot)const{ return Poly(*this) += ot; }
	Poly operator - (const Poly& ot)const{ return Poly(*this) -= ot; }
	Poly operator * (const Poly& ot)const{ return Poly(*this) *= ot; }
	Poly operator / (const Poly& ot)const{ return Poly(*this) /= ot; }
	Poly operator % (const Poly& ot)const{ return Poly(*this) %= ot; }
};
constexpr int MOD = 104857601, PR = 3; // (25 << 22) + 1
using mint_t = MINT<MOD>;
using poly_t = Poly<MOD, PR>;
int main(){
	cin.tie(0)->sync_with_stdio(0);
	int K; lld N; cin >> K >> N; N--;
	vector<mint_t> A(K), B(K+1);
	for (auto& v: A) cin >> v;
	for (int i=1;i<=K;i++) cin >> B[i], B[i] *= -1;
	B[0] = 1;
	poly_t F(A), Q(B);
	poly_t P = (F*Q).trim(K);
	for (;N;N>>=1){
		poly_t Q2 = Q;
		for (int i=1;i<Q2.size();i+=2) Q2[i] *= -1;
		poly_t U = P*Q2, V = Q*Q2;
		P.resize(K); U.resize(K<<1); V.resize(K<<1|1); // prevent out of index
		for (int i=N&1;i<(K<<1);i+=2) P[i>>1] = U[i];
		Q.resize(K+1);
		for (int i=0;i<=K;i++) Q[i] = V[i<<1];
	}
	cout << P[0]/Q[0] << '\n';
}

Bostan Mori 알고리즘

댓글