#include #include #include #include using namespace std; const int mod = 1e9 + 7; struct Modint { int n; Modint(int n = 0) : n(n) {} }; Modint operator+(Modint a, Modint b) { return Modint((a.n += b.n) >= mod ? a.n - mod : a.n); } Modint operator-(Modint a, Modint b) { return Modint((a.n -= b.n) < 0 ? a.n + mod : a.n); } Modint operator*(Modint a, Modint b) { return Modint(1LL * a.n * b.n % mod); } Modint &operator+=(Modint &a, Modint b) { return a = a + b; } Modint &operator-=(Modint &a, Modint b) { return a = a - b; } Modint &operator*=(Modint &a, Modint b) { return a = a * b; } Modint modinv(Modint a) { if (a.n == 1) return 1; return modinv(mod % a.n) * (mod - mod / a.n); } Modint operator/(Modint a, Modint b) { return a * modinv(b); } std::vector berlekamp_massey(std::vector s) { using K = Modint; const int N = s.size(); std::vector C(N); std::vector B(N); C[0] = 1; B[0] = 1; int L = 0; int m = 1; K b = 1; for (int n = 0; n < N; n++) { K d = s[n]; for (int i = 1; i <= L; i++) d += C[i] * s[n - i]; if (d.n == 0) { m++; } else if (2 * L <= n) { auto T = C; for (int i = 0; i + m < N; i++) C[i + m] -= B[i] * (d / b); L = n + 1 - L; B = T; b = d; m = 1; } else { for (int i = 0; i + m < N; i++) C[i + m] -= B[i] * (d / b); m++; } } C.resize(L + 1); reverse(C.begin(), C.end()); return C; } vector poly_mod(vector a, const vector &m) { const int n = m.size(); for (int i = a.size() - 1; i >= m.size(); i--) { for (int j = 0; j < m.size(); j++) { a[i - n + j] += a[i] * m[j]; } } a.resize(m.size()); return a; } // a*b mod m(x) vector poly_mul(const vector &a, const vector &b, const vector &m) { vector ret(a.size() + b.size() - 1); for (int i = 0; i < a.size(); i++) { for (int j = 0; j < b.size(); j++) { ret[i + j] += a[i] * b[j]; } } return poly_mod(ret, m); } // a^b mod m(x) vector poly_pow(vector a, long long b, const vector &m) { vector ret(1); ret[0] = 1; while (b > 0) { if (b & 1) ret = poly_mul(ret, a, m); a = poly_mul(a, a, m); b >>= 1; } return poly_mod(ret, m); } Modint linear_recurrence_relation(vector a, long long k) { auto m = berlekamp_massey(a); m.pop_back(); for (Modint &a : m) { a *= mod - 1; } auto p = poly_pow({0, 1}, k, m); Modint res = 0; for (int i = 0; i < p.size(); i++) { res += a[i] * p[i]; } return res; } bool g[3][50]; Modint memo[50][8 * 8]; bool vis[50][8 * 8]; Modint f(int n, int k) { if (k == n * 3) { bool any = false; for (int i = 0; i < 3; i++) { for (int j = 0; j < n; j++) { any |= i + 1 < 3 && !g[i][j] && !g[i + 1][j]; any |= j + 1 < n && !g[i][j] && !g[i][j + 1]; } } return !any; } int y = k % 3; int x = k / 3; if (y == 0 && x >= 2) { bool any = false; for (int i = 0; i < 3; i++) { for (int j = 0; j < x; j++) { any |= i + 1 < 3 && !g[i][j] && !g[i + 1][j]; any |= j + 1 < x && !g[i][j] && !g[i][j + 1]; } } if (any) return 0; int bit = 0; bit = bit << 1 | g[0][x - 1]; bit = bit << 1 | g[1][x - 1]; bit = bit << 1 | g[2][x - 1]; bit = bit << 1 | g[0][x]; bit = bit << 1 | g[1][x]; bit = bit << 1 | g[2][x]; if (vis[x][bit]) return memo[x][bit]; } Modint res = f(n, k + 1); if (y + 1 < 3 && !g[y][x] && !g[y + 1][x]) { g[y][x] = g[y + 1][x] = true; res += f(n, k + 1); g[y][x] = g[y + 1][x] = false; } if (x + 1 < n && !g[y][x] && !g[y][x + 1]) { g[y][x] = g[y][x + 1] = true; res += f(n, k + 1); g[y][x] = g[y][x + 1] = false; } if (y == 0 && x >= 2) { int bit = 0; bit = bit << 1 | g[0][x - 1]; bit = bit << 1 | g[1][x - 1]; bit = bit << 1 | g[2][x - 1]; bit = bit << 1 | g[0][x]; bit = bit << 1 | g[1][x]; bit = bit << 1 | g[2][x]; vis[x][bit] = true; memo[x][bit] = res; } return res; } int main() { vector a; for (int i = 0; i < 40; i++) { memset(vis, false, sizeof(vis)); a.emplace_back(f(i, 0)); } long long n; cin >> n; cout << linear_recurrence_relation(a, n).n << endl; }