#ifndef ___CLASS_MODINT
#define ___CLASS_MODINT

#include <cstdint>

template<std::uint32_t mod>
class modint {
private:
	std::uint32_t n;
public:
	modint() : n(0) {};
	modint(std::uint64_t n_) : n(n_ % mod) {};
	bool operator==(const modint& m) const { return n == m.n; }
	bool operator!=(const modint& m) const { return n != m.n; }
	std::uint32_t get() const { return n; }
	modint& operator+=(const modint& m) { n += m.n; n = (n < mod ? n : n - mod); return *this; }
	modint& operator-=(const modint& m) { n += mod - m.n; n = (n < mod ? n : n - mod); return *this; }
	modint& operator*=(const modint& m) { n = std::uint64_t(n) * m.n % mod; return *this; }
	modint operator+(const modint& m) const { return modint(*this) += m; }
	modint operator-(const modint& m) const { return modint(*this) -= m; }
	modint operator*(const modint& m) const { return modint(*this) *= m; }
	modint binpow(std::uint64_t b) const {
		modint ans = 1, m = modint(*this);
		while (b) {
			if (b & 1) ans *= m;
			m *= m;
			b >>= 1;
		}
		return ans;
	}
	modint inv() { return (*this).binpow(mod - 2); }
};

#endif // ___CLASS_MODINT

#include <queue>
#include <vector>
#include <iostream>
#include <algorithm>
#include <functional>
using namespace std;
const int mod = 1000000009;
const int inf = 1012345678;
using modulo = modint<mod>;
const vector<int> dx = { 0, 1, 0, -1 };
const vector<int> dy = { 1, 0, -1, 0 };
vector<vector<modulo> > calc(vector<string> S, bool will_cut) {
	// S must be a square grid
	int N = S.size();
	vector<vector<int> > dist(N, vector<int>(N, inf)); dist[0][0] = 0;
	vector<vector<modulo> > ans(N, vector<modulo>(N)); ans[0][0] = 1;
	queue<pair<int, int> > que; que.push(make_pair(0, 0));
	while (!que.empty()) {
		pair<int, int> u = que.front(); que.pop();
		int d = dist[u.first][u.second];
		for (int i = 0; i < 4; ++i) {
			int tx = u.first + dx[i], ty = u.second + dy[i];
			if (0 <= tx && tx < N && 0 <= ty && ty < N && S[tx][ty] == '.') {
				if (dist[tx][ty] == inf) {
					dist[tx][ty] = d + 1;
					que.push(make_pair(tx, ty));
				}
				if (dist[tx][ty] == d + 1) {
					ans[tx][ty] += ans[u.first][u.second];
				}
			}
		}
	}
	int mindiff = inf / 2;
	for (int i = 0; i < N; ++i) {
		for (int j = 0; j < N; ++j) {
			if (i != N - 1 && j != N - 1) {
				ans[i][j] = 0;
			}
			else {
				mindiff = min(mindiff, dist[i][j] - i - j);
			}
		}
	}
	if (will_cut) {
		for (int i = 0; i < N; ++i) {
			for (int j = 0; j < N; ++j) {
				if ((i == N - 1 || j == N - 1) && dist[i][j] - i - j != mindiff) {
					ans[i][j] = 0;
				}
			}
		}
	}
	return ans;
}
modulo solve_easy(int N, int M, vector<int> X, vector<int> Y) {
	vector<string> S(N, string(N, '.'));
	for (int i = 0; i < M; ++i) {
		S[X[i]][Y[i]] = '#';
	}
	vector<vector<modulo> > vans = calc(S, false);
	modulo ans = vans[N - 1][N - 1];
	return ans;
}
modulo solve(int N, int M, vector<int> X, vector<int> Y) {
	vector<modulo> fact(2 * N), inv(2 * N), factinv(2 * N);
	fact[0] = 1; inv[1] = 1; factinv[0] = 1;
	for (int i = 1; i < 2 * N; ++i) fact[i] = fact[i - 1] * i;
	for (int i = 2; i < 2 * N; ++i) inv[i] = inv[mod % i] * (mod - mod / i);
	for (int i = 1; i < 2 * N; ++i) factinv[i] = factinv[i - 1] * inv[i];
	function<modulo(int, int)> comb = [&](int x, int y) {
		if (y < 0 || x < y) return modulo(0);
		return fact[x] * factinv[y] * factinv[x - y];
	};
	modulo ans;
	if (N <= 3 || N <= 2 * M + 1) {
		return solve_easy(N, M, X, Y);
	}
	else {
		int lim = max(M, 1);
		vector<string> bs(lim, string(lim, '.')), ts(lim, string(lim, '.'));
		vector<int> RX, RY;
		vector<bool> bd(lim), br(lim), tu(lim), tl(lim);
		for (int i = 0; i < M; ++i) {
			if (X[i] < lim && Y[i] < lim) bs[X[i]][Y[i]] = '#';
			else if (X[i] >= N - lim && Y[i] >= N - lim) ts[N - X[i] - 1][N - Y[i] - 1] = '#';
			else {
				if (X[i] == lim && Y[i] < lim) bd[Y[i]] = true;
				if (Y[i] == lim && X[i] < lim) br[X[i]] = true;
				if (X[i] == N - lim - 1 && Y[i] >= N - lim) tu[N - Y[i] - 1] = true;
				if (Y[i] == N - lim - 1 && X[i] >= N - lim) tl[N - X[i] - 1] = true;
				RX.push_back(X[i]);
				RY.push_back(Y[i]);
			}
		}
		int R = RX.size();
		for (int i = 0; i < R; ++i) {
			for (int j = 0; j < i; ++j) {
				if (RX[j] > RX[i] || (RX[j] == RX[i] && RY[j] > RY[i])) {
					swap(RX[j], RX[i]);
					swap(RY[j], RY[i]);
				}
			}
		}
		vector<vector<modulo> > bans = calc(bs, true);
		vector<vector<modulo> > tans = calc(ts, true);
		vector<int> BX, BY, TX, TY;
		for (int i = 0; i < lim; ++i) {
			if (!bd[i] && bans[lim - 1][i].get() > 0) {
				BX.push_back(lim);
				BY.push_back(i);
			}
			if (!br[i] && bans[i][lim - 1].get() > 0) {
				BX.push_back(i);
				BY.push_back(lim);
			}
			if (!tu[i] && tans[lim - 1][i].get() > 0) {
				TX.push_back(N - lim - 1);
				TY.push_back(N - i - 1);
			}
			if (!tl[i] && tans[i][lim - 1].get() > 0) {
				TX.push_back(N - i - 1);
				TY.push_back(N - lim - 1);
			}
		}
		RX.insert(RX.begin(), BX.begin(), BX.end()); RX.insert(RX.end(), TX.begin(), TX.end());
		RY.insert(RY.begin(), BY.begin(), BY.end()); RY.insert(RY.end(), TY.begin(), TY.end());
		int BS = BX.size(), TS = TX.size();
		vector<modulo> dp(BS + R + TS);
		for (int i = 0; i < BS; ++i) {
			if (RX[i] == lim) dp[i] = modulo(0) - bans[lim - 1][RY[i]];
			if (RY[i] == lim) dp[i] = modulo(0) - bans[RX[i]][lim - 1];
		}
		for (int i = BS; i < BS + R + TS; ++i) {
			for (int j = 0; j < i && j < BS + R; ++j) {
				dp[i] -= dp[j] * comb((RX[i] - RX[j]) + (RY[i] - RY[j]), RX[i] - RX[j]);
			}
		}
		for (int i = 0; i < TS; ++i) {
			modulo sub;
			if (TX[i] == N - lim - 1) sub = tans[lim - 1][N - TY[i] - 1];
			if (TY[i] == N - lim - 1) sub = tans[N - TX[i] - 1][lim - 1];
			ans += dp[BS + R + i] * sub;
		}
	}
	return ans;
}
int main() {
	int N, M;
	cin >> N >> M;
	vector<int> X(M), Y(M);
	for (int i = 0; i < M; ++i) {
		cin >> X[i] >> Y[i]; --X[i], --Y[i];
	}
	modulo ans = solve(N, M, X, Y);
	cout << ans.get() << endl;
	return 0;
}

### Prompt

Please provide a CPP coded solution to the problem described below:

There is an NxN grid.

<image>
(Figure when N = 3)

It costs 1 to move to an adjacent square. However, it cannot be moved diagonally. Also, there are M squares with obstacles, and you are not allowed to enter those squares. How many ways are there to move from mass (1, 1) to mass (N, N) at the lowest cost? The answer can be big, so take the mod at 1000000009 (= 109 + 9) and output it.

Constraints

* 2 ≤ N ≤ 106

* 0 ≤ M ≤ 50

* 1 ≤ Xi, Yi ≤ N

* If i ≠ j, then (Xi, Yi) ≠ (Xj, Yj)

* (Xi, Yi) ≠ (1, 1)

* (Xi, Yi) ≠ (N, N)

Input

Input is given in the following format

> N M
> X1 Y1
> X2 Y2
> ……
> XM YM
>

Xi, Yi means that there is an obstacle in (Xi, Yi)

Output

Output the total number of routes on one line. However, if there is no route to move from the mass (1, 1) to the mass (N, N), output 0.

Examples

Input

3 0


Output

6


Input

3 1
2 2


Output

2


Input

5 8
4 3
2 1
4 2
4 4
3 4
2 2
2 4
1 4


Output

1


Input

1000 10
104 87
637 366
393 681
215 604
707 876
943 414
95 327
93 415
663 596
661 842


Output

340340391


Input

2 2
1 2
2 1


Output

0

### Response

```cpp
#ifndef ___CLASS_MODINT
#define ___CLASS_MODINT

#include <cstdint>

template<std::uint32_t mod>
class modint {
private:
	std::uint32_t n;
public:
	modint() : n(0) {};
	modint(std::uint64_t n_) : n(n_ % mod) {};
	bool operator==(const modint& m) const { return n == m.n; }
	bool operator!=(const modint& m) const { return n != m.n; }
	std::uint32_t get() const { return n; }
	modint& operator+=(const modint& m) { n += m.n; n = (n < mod ? n : n - mod); return *this; }
	modint& operator-=(const modint& m) { n += mod - m.n; n = (n < mod ? n : n - mod); return *this; }
	modint& operator*=(const modint& m) { n = std::uint64_t(n) * m.n % mod; return *this; }
	modint operator+(const modint& m) const { return modint(*this) += m; }
	modint operator-(const modint& m) const { return modint(*this) -= m; }
	modint operator*(const modint& m) const { return modint(*this) *= m; }
	modint binpow(std::uint64_t b) const {
		modint ans = 1, m = modint(*this);
		while (b) {
			if (b & 1) ans *= m;
			m *= m;
			b >>= 1;
		}
		return ans;
	}
	modint inv() { return (*this).binpow(mod - 2); }
};

#endif // ___CLASS_MODINT

#include <queue>
#include <vector>
#include <iostream>
#include <algorithm>
#include <functional>
using namespace std;
const int mod = 1000000009;
const int inf = 1012345678;
using modulo = modint<mod>;
const vector<int> dx = { 0, 1, 0, -1 };
const vector<int> dy = { 1, 0, -1, 0 };
vector<vector<modulo> > calc(vector<string> S, bool will_cut) {
	// S must be a square grid
	int N = S.size();
	vector<vector<int> > dist(N, vector<int>(N, inf)); dist[0][0] = 0;
	vector<vector<modulo> > ans(N, vector<modulo>(N)); ans[0][0] = 1;
	queue<pair<int, int> > que; que.push(make_pair(0, 0));
	while (!que.empty()) {
		pair<int, int> u = que.front(); que.pop();
		int d = dist[u.first][u.second];
		for (int i = 0; i < 4; ++i) {
			int tx = u.first + dx[i], ty = u.second + dy[i];
			if (0 <= tx && tx < N && 0 <= ty && ty < N && S[tx][ty] == '.') {
				if (dist[tx][ty] == inf) {
					dist[tx][ty] = d + 1;
					que.push(make_pair(tx, ty));
				}
				if (dist[tx][ty] == d + 1) {
					ans[tx][ty] += ans[u.first][u.second];
				}
			}
		}
	}
	int mindiff = inf / 2;
	for (int i = 0; i < N; ++i) {
		for (int j = 0; j < N; ++j) {
			if (i != N - 1 && j != N - 1) {
				ans[i][j] = 0;
			}
			else {
				mindiff = min(mindiff, dist[i][j] - i - j);
			}
		}
	}
	if (will_cut) {
		for (int i = 0; i < N; ++i) {
			for (int j = 0; j < N; ++j) {
				if ((i == N - 1 || j == N - 1) && dist[i][j] - i - j != mindiff) {
					ans[i][j] = 0;
				}
			}
		}
	}
	return ans;
}
modulo solve_easy(int N, int M, vector<int> X, vector<int> Y) {
	vector<string> S(N, string(N, '.'));
	for (int i = 0; i < M; ++i) {
		S[X[i]][Y[i]] = '#';
	}
	vector<vector<modulo> > vans = calc(S, false);
	modulo ans = vans[N - 1][N - 1];
	return ans;
}
modulo solve(int N, int M, vector<int> X, vector<int> Y) {
	vector<modulo> fact(2 * N), inv(2 * N), factinv(2 * N);
	fact[0] = 1; inv[1] = 1; factinv[0] = 1;
	for (int i = 1; i < 2 * N; ++i) fact[i] = fact[i - 1] * i;
	for (int i = 2; i < 2 * N; ++i) inv[i] = inv[mod % i] * (mod - mod / i);
	for (int i = 1; i < 2 * N; ++i) factinv[i] = factinv[i - 1] * inv[i];
	function<modulo(int, int)> comb = [&](int x, int y) {
		if (y < 0 || x < y) return modulo(0);
		return fact[x] * factinv[y] * factinv[x - y];
	};
	modulo ans;
	if (N <= 3 || N <= 2 * M + 1) {
		return solve_easy(N, M, X, Y);
	}
	else {
		int lim = max(M, 1);
		vector<string> bs(lim, string(lim, '.')), ts(lim, string(lim, '.'));
		vector<int> RX, RY;
		vector<bool> bd(lim), br(lim), tu(lim), tl(lim);
		for (int i = 0; i < M; ++i) {
			if (X[i] < lim && Y[i] < lim) bs[X[i]][Y[i]] = '#';
			else if (X[i] >= N - lim && Y[i] >= N - lim) ts[N - X[i] - 1][N - Y[i] - 1] = '#';
			else {
				if (X[i] == lim && Y[i] < lim) bd[Y[i]] = true;
				if (Y[i] == lim && X[i] < lim) br[X[i]] = true;
				if (X[i] == N - lim - 1 && Y[i] >= N - lim) tu[N - Y[i] - 1] = true;
				if (Y[i] == N - lim - 1 && X[i] >= N - lim) tl[N - X[i] - 1] = true;
				RX.push_back(X[i]);
				RY.push_back(Y[i]);
			}
		}
		int R = RX.size();
		for (int i = 0; i < R; ++i) {
			for (int j = 0; j < i; ++j) {
				if (RX[j] > RX[i] || (RX[j] == RX[i] && RY[j] > RY[i])) {
					swap(RX[j], RX[i]);
					swap(RY[j], RY[i]);
				}
			}
		}
		vector<vector<modulo> > bans = calc(bs, true);
		vector<vector<modulo> > tans = calc(ts, true);
		vector<int> BX, BY, TX, TY;
		for (int i = 0; i < lim; ++i) {
			if (!bd[i] && bans[lim - 1][i].get() > 0) {
				BX.push_back(lim);
				BY.push_back(i);
			}
			if (!br[i] && bans[i][lim - 1].get() > 0) {
				BX.push_back(i);
				BY.push_back(lim);
			}
			if (!tu[i] && tans[lim - 1][i].get() > 0) {
				TX.push_back(N - lim - 1);
				TY.push_back(N - i - 1);
			}
			if (!tl[i] && tans[i][lim - 1].get() > 0) {
				TX.push_back(N - i - 1);
				TY.push_back(N - lim - 1);
			}
		}
		RX.insert(RX.begin(), BX.begin(), BX.end()); RX.insert(RX.end(), TX.begin(), TX.end());
		RY.insert(RY.begin(), BY.begin(), BY.end()); RY.insert(RY.end(), TY.begin(), TY.end());
		int BS = BX.size(), TS = TX.size();
		vector<modulo> dp(BS + R + TS);
		for (int i = 0; i < BS; ++i) {
			if (RX[i] == lim) dp[i] = modulo(0) - bans[lim - 1][RY[i]];
			if (RY[i] == lim) dp[i] = modulo(0) - bans[RX[i]][lim - 1];
		}
		for (int i = BS; i < BS + R + TS; ++i) {
			for (int j = 0; j < i && j < BS + R; ++j) {
				dp[i] -= dp[j] * comb((RX[i] - RX[j]) + (RY[i] - RY[j]), RX[i] - RX[j]);
			}
		}
		for (int i = 0; i < TS; ++i) {
			modulo sub;
			if (TX[i] == N - lim - 1) sub = tans[lim - 1][N - TY[i] - 1];
			if (TY[i] == N - lim - 1) sub = tans[N - TX[i] - 1][lim - 1];
			ans += dp[BS + R + i] * sub;
		}
	}
	return ans;
}
int main() {
	int N, M;
	cin >> N >> M;
	vector<int> X(M), Y(M);
	for (int i = 0; i < M; ++i) {
		cin >> X[i] >> Y[i]; --X[i], --Y[i];
	}
	modulo ans = solve(N, M, X, Y);
	cout << ans.get() << endl;
	return 0;
}

```