結果
| 問題 | No.93 ペガサス |
| ユーザー |
 anta
|
| 提出日時 | 2018-06-22 13:48:19 |
| 言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 2 ms / 5,000 ms |
| コード長 | 10,664 bytes |
| コンパイル時間 | 1,835 ms |
| コンパイル使用メモリ | 182,936 KB |
| 実行使用メモリ | 5,376 KB |
| 最終ジャッジ日時 | 2024-06-30 17:56:34 |
| 合計ジャッジ時間 | 2,478 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge3 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 16 |
ソースコード
#include "bits/stdc++.h"using namespace std;template<int MOD>struct ModInt {static const int Mod = MOD;unsigned x;ModInt() : x(0) { }ModInt(signed sig) { int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; }ModInt(signed long long sig) { int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; }int get() const { return (int)x; }ModInt &operator+=(ModInt that) { if ((x += that.x) >= MOD) x -= MOD; return *this; }ModInt &operator-=(ModInt that) { if ((x += MOD - that.x) >= MOD) x -= MOD; return *this; }ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; }ModInt &operator/=(ModInt that) { return *this *= that.inverse(); }ModInt operator+(ModInt that) const { return ModInt(*this) += that; }ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }ModInt operator/(ModInt that) const { return ModInt(*this) /= that; }ModInt inverse() const {signed a = x, b = MOD, u = 1, v = 0;while (b) {signed t = a / b;a -= t * b; std::swap(a, b);u -= t * v; std::swap(u, v);}if (u < 0) u += Mod;ModInt res; res.x = (unsigned)u;return res;}bool operator==(ModInt that) const { return x == that.x; }bool operator!=(ModInt that) const { return x != that.x; }ModInt operator-() const { ModInt t; t.x = x == 0 ? 0 : Mod - x; return t; }};typedef ModInt<1000000007> mint;struct GaussianEliminationCore {using Num = mint;using RowVec = vector<Num>;static void multiplySubtract(RowVec &x, const RowVec &y, int m, Num c) {if (c == Num()) return;for (int j = 0; j < m; ++ j)x[j] -= y[j] * c;}static Num dotProduct(const RowVec &x, const RowVec &y, int m) {Num sum = Num();for (int j = 0; j < m; ++ j)sum += x[j] * y[j];return sum;}vector<RowVec> basis;vector<Num> invDiagonal;vector<int> order;void init(int m) {basis.assign(m, RowVec());invDiagonal.assign(m, Num());order.clear();}int size() const { return (int)basis.size(); }void eliminate(RowVec &row, vector<Num> &coeffs) const {for (int i : order) {Num c = row[i] * invDiagonal[i];coeffs[i] = c;multiplySubtract(row, basis[i], size(), c);}}void add(const RowVec &row, int k) {assert(row[k] != Num() && invDiagonal[k] == Num());basis[k] = row;invDiagonal[k] = row[k].inverse();order.push_back(k);}};// idea: use falling factorial basis for polynomialsstruct PolynomiallyRecursiveSequence : private GaussianEliminationCore {enum PolynomalBasisKind {MonominalBasis,FallingFactorialBasis,} polynomialBasisKind = MonominalBasis;void getPolynomialBasis(vector<Num> &xs, int D, int n) {for (int j = 0; j < D; ++ j) {Num x = 1;if (polynomialBasisKind == MonominalBasis) {for (int k = 0; k < j; ++ k)x *= n;} else {for (int k = 0; k < j; ++ k)x *= n - k;}xs[j] = x;}}vector<Num> solve(const vector<Num> &seq, const vector<pair<int, int>> &indices) {int K = 0, D = 0;for (auto ix : indices) {if (K <= ix.first) K = ix.first + 1;if (D <= ix.second) D = ix.second + 1;}if (K == 0 || D == 0) return {};int N = (int)seq.size(), m = (int)indices.size();init(m);vector<Num> newRow(m), coeffs(m);for (int n = K - 1; n < N; ++ n) {vector<Num> ys(K), xs(D);for (int i = 0; i < K; ++ i) ys[i] = seq[n - i];getPolynomialBasis(xs, D, n);newRow.clear();for (auto ix : indices)newRow.push_back(ys[ix.first] * xs[ix.second]);eliminate(newRow, coeffs);for (int k = 0; k < m; ++ k) if (newRow[k] != mint()) {add(newRow, k);break;}if (order.size() == m) return {};}int k = 0;for (; invDiagonal[k] != Num(); ++ k);vector<Num> res(m);res[k] = 1;reverse(order.begin(), order.end());for (int i : order) {Num dp = dotProduct(res, basis[i], m);res[i] -= dp * invDiagonal[i];}return res;}vector<vector<Num>> findMinimalSolution(const vector<Num> &seq) {vector<vector<Num>> best;int bestNum = numeric_limits<int>::max();for (int KD = 1; ; ++ KD) {for (int K = 1; K <= KD; ++ K) if (KD % K == 0) {int D = KD / K;vector<bool> enabled(KD, true);int currentNum = KD, leastNum = 0;auto check = [&]() -> bool {vector<pair<int, int>> indices;for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) if (enabled[i * D + j])indices.emplace_back(i, j);auto sol = solve(seq, indices);if (!sol.empty()) {if (bestNum > currentNum) {best.assign(K, vector<Num>(D));int t = 0;for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) if (enabled[i * D + j])best[i][j] = sol[t ++];bestNum = currentNum;}return true;} else {return false;}};function<void(int, int)> rec = [&](int i, int j) {if (bestNum <= leastNum) return;if (i == K) {check();return;}if (j == D) return rec(i + 1, 0);for (int e = 0; e < 2; ++ e) {enabled[i * D + j] = e != 0;if (e == 0)-- currentNum;else++ leastNum;if (e == 1 || check()) {rec(i, j + 1);}if (e == 0)++ currentNum;else-- leastNum;}};rec(0, 0);}if (!best.empty()) break;}cerr << "bestNum = " << bestNum << endl;return best;}};template<typename T>T gcd(T x, T y) { if (y == 0)return x; else return gcd(y, x%y); }template<int MOD> int mintToSigned(ModInt<MOD> a) {int x = a.get();if (x <= MOD / 2)return x;elsereturn x - MOD;}string mintToSignedRatio(mint a, int den = 36) {int x = mintToSigned(a * den);int g = gcd(abs(x), abs(den));if (den / g < 0) g *= -1;x /= g, den /= g;stringstream ss;if (den == 1)ss << x;elsess << x << "/" << den;return ss.str();}int main() {if(0) {vector<mint> seq = {1, 1, 2, 2, 8, 28, 152, 952, 7208, 62296, 605864, 6522952, 76951496, 986411272, 647501133, 653303042, 170637030, 248109503, 700583494,619914523, 682935856, 443753916, 423068688, 507501942, 315541972, 110825117, 848156395, 798418282, 920964362, 23823302, 114894774,279365223, 992413784, 833179437, 785518302, 524368220, 42214454, 140345871, 188150268, 808714798, 718376249, 732000901, 955005007,139255097, 484615744, 615066955, 726914809, 856989248, 460819998, 321277105, 536397091, 555447300, 597473569, 217709372, 24981477,143561526, 171000806, 137649694, 749333590, 700935246, 916763337, 762367836, 296796066, 236278263, 398507715, 148909632, 568524543,926513708, 163591024, 339393165, 549241395, 548924577, 915489821, 706913104, 380913764, 993919668, 895691202, 628078606, 542382606,735060428, 385303214, 453133962, 470556393, 439972973, 4764973, 459438929, 49172129, 93448766, 14767450, 302365655, 44994640,637650527, 462797839, 174866371, 963824426, 761996745, 999013044, 209330964, 997280223, 561428453, 300321098, 733666220, 501835787,614206349, 790510106, 3598758, 51698818, 133690163, 663340013, 193264830, 514261029, 953300090, 783013838, 841987389, 44621114,487627404};int N = (int)seq.size();cerr << "N = " << N << endl;PolynomiallyRecursiveSequence prs;prs.polynomialBasisKind = PolynomiallyRecursiveSequence::MonominalBasis;auto solution = prs.findMinimalSolution(seq);if (solution.empty()) {cerr << "No solution found" << endl;return 1;}int K = (int)solution.size(), D = (int)solution[0].size();for (int n = K - 1; n < N; ++ n) {vector<mint> ys(K), xs(D);for (int i = 0; i < K; ++ i) ys[i] = seq[n - i];prs.getPolynomialBasis(xs, D, n);mint sum;for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j)sum += solution[i][j] * ys[i] * xs[j];if (sum != mint())cerr << "err" << endl;}for (int i = 0; i < K; ++ i) {auto p = solution[i];if (i == 0) {for (auto &x : p) x = -x;}int t = 0;for(auto x : p) t += x != mint();int hi = D - 1;for (; hi >= 0 && p[hi] == mint(); -- hi);if (hi >= 0 && mintToSignedRatio(p[hi])[0] == '-') {cout << (i > 1 ? " - " : "-");for (auto &x : p) x = -x;} else {if(i > 1) cout << " + ";}if (hi == 0 && p[hi] == 1) {} else {if (t > 1) cout << "(";auto &o = cout;bool first = true;for (int j = D - 1; j >= 0; -- j) {string c = mintToSignedRatio(p[j]);if (c != "0") {if (first && c[0] == '-') o << "-";else if (first) o << "";else if (c[0] == '-') o << " - ";else o << " + ";if (j != 0 && (c == "1" || c == "-1")) o << "";else if (c[0] == '-') o << c.substr(1);else o << c;if (j == 0) o << "";else if (j == 1) o << "n";else {if (prs.polynomialBasisKind == PolynomiallyRecursiveSequence::MonominalBasis) {o << "n^" << j;} else {o << "n";for (int k = 1; k < j; ++ k)o << "(n-" << k << ")";}}first = false;}}if (first) o << "0";if (t > 1) cout << ")";cout << " ";}cout << "a_";if (i == 0)cout << "n";elsecout << "{n-" << i << "}";if (i == 0) {cout << " = ";if (K == 1) cout << "0";}}cout << endl;cout << "const int K = " << K << ";\n";cout << "const array<mint, K - 1> init = { ";for (int j = 0; j < K - 1; ++ j) cout << (j == 0 ? "" : ", ") << seq[j].get();cout << " };\n";cout << "const array<array<mint, " << D << ">, K> coeffs = { {\n";for (int i = 0; i < K; ++ i) {cout << "\t{ ";for (int j = 0; j < D; ++ j) cout << (j == 0 ? "" : ", ") << mintToSigned(solution[i][j] * 2);cout << " },\n";}}int N;while (cin >> N) {const int K = 14;const array<mint, K - 1> init = { 1, 1, 2, 2, 8, 28, 152, 952, 7208, 62296, 605864, 6522952, 76951496 };const array<array<mint, 3>, K> coeffs = { {{ 4, 0, 0 },{ -18, -4, 0 },{ 6, 16, 0 },{ 18, -11, 2 },{ 116, -7, -9 },{ -138, -6, 8 },{ -64, -78, 14 },{ -90, 136, -18 },{ -632, 256, -20 },{ 892, -298, 22 },{ 1674, -354, 18 },{ -1480, 307, -16 },{ -240, 55, -3 },{ 240, -44, 2 },} };const mint invDen = coeffs[0][0].inverse();vector<mint> seq(init.begin(), init.end());seq.resize(N + 1);for (int n = K - 1; n <= N; ++ n) {mint n1 = n, n2 = n1 * n1;mint sum;for (int i = 1; i < K; ++ i)sum += (coeffs[i][0] + coeffs[i][1] * n1 + coeffs[i][2] * n2) * seq[n - i];seq[n] = -sum * invDen;}mint ans = seq[N];printf("%d\n", ans.get());}}