結果
| 問題 | No.660 家を通り過ぎないランダムウォーク問題 |
| ユーザー |
 anta
|
| 提出日時 | 2018-06-22 19:00:12 |
| 言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 59 ms / 2,000 ms |
| コード長 | 16,782 bytes |
| コンパイル時間 | 1,980 ms |
| コンパイル使用メモリ | 187,644 KB |
| 実行使用メモリ | 6,948 KB |
| 最終ジャッジ日時 | 2024-06-30 18:03:18 |
| 合計ジャッジ時間 | 3,605 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 45 |
ソースコード
#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);}};struct 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, uint64_t maxRecurse = numeric_limits<uint64_t>::max()) {vector<vector<Num>> best;int bestNum = numeric_limits<int>::max();uint64_t numRecurse = 0;for (int KD = 1; ; ++ KD) {cerr << "KD = " << KD << "..." << endl;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 {++ numRecurse;if (numRecurse % 10000 == 0) cerr << "checking " << numRecurse << " times... (bestNum = " << bestNum << ")" << endl;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 (sol[0] != mint()) {auto inv = sol[0].inverse();for (auto &x : sol) x *= inv;}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;cerr << "bestNum updated: " << bestNum << endl;}return true;} else {return false;}};function<void(int, int)> rec = [&](int i, int j) {if (bestNum <= leastNum) return;if (i == K) {check();++ numRecurse;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 (numRecurse > maxRecurse) return;}if (e == 0)++ currentNum;else-- leastNum;}};if (check())rec(0, 0);if (numRecurse > maxRecurse) {cerr << "(suboptimal result)" << endl;break;}}if (!best.empty()) break;if (numRecurse > maxRecurse) {break;}}if (!best.empty()) {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 x = mintToSigned(a), d = 1;for (int dd = 1; dd <= 60; ++ dd) {int xx = mintToSigned(a * dd);if (abs(x) > abs(xx)) x = xx, d = dd;}int g = gcd(abs(x), abs(d));if (d / g < 0) g *= -1;x /= g, d /= g;stringstream ss;if (d == 1)ss << x;elsess << x << "/" << d;return ss.str();}#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))vector<mint> fact, factinv;void nCr_computeFactinv(int N) {N = min(N, mint::Mod - 1);fact.resize(N + 1); factinv.resize(N + 1);fact[0] = 1;rer(i, 1, N) fact[i] = fact[i - 1] * i;factinv[N] = fact[N].inverse();for (int i = N; i >= 1; i --) factinv[i - 1] = factinv[i] * i;}mint nCr(int n, int r) {if (n >= mint::Mod)return nCr(n % mint::Mod, r % mint::Mod) * nCr(n / mint::Mod, r / mint::Mod);return r > n ? 0 : fact[n] * factinv[n - r] * factinv[r];}mint catalan_number(int n) {return n == 0 ? 1 : nCr(2 * n, n) - nCr(2 * n, n - 1);}int digitsMod(const char s[], int n, int m) {int x = 0; long long pow10 = 1;for (int i = n - 1; i >= 0; i --) {if ((x += pow10 * (s[i] - '0') % m) >= m) x -= m;pow10 = pow10 * 10 % m;}return x;}template<int MOD> ModInt<MOD> operator^(ModInt<MOD> a, unsigned long long k) {ModInt<MOD> r = 1;while (k) {if (k & 1) r *= a;a *= a;k >>= 1;}return r;}int main() {if (0) {vector<mint> seq;if (1) {seq = { 1,1,3,4,19,26,144,197,1171,1597,9878,13432,85216,115597,746371,1010504,6609043,8933858,59008563,79662593,530279894,715116833,790262320,454109424,458522675,507912745,683985782,349356902,641159448,914253017,666123528,256877025,283012377,259341341,985344830,99096838,997431143,359782436,337868019,275073093,823041379,299193173,618332669,623673376,788988015,348012107,960452357,132459813,484572463,876182618,173753254,441164007,288391002,290937936,415692760,211813919,329314918,283643646,565011970,845017018,390800948,209282283,329903021,700520426,397612109,649057223,946472136,142317513,639359324,720717675,983275264,507698034,651899763,234922533,86387948,840663423,815571175,144003918,212527976,763824514,542610091,452847843,921495990,979234007,906307312,656853418,297047315,193463125,69677636,636423896,799010307,909263327,147379495,774893131,945207456,461924132,658630324,333324822,851825875,36247927,769534543,614653765,155284112,370877182,608767657,471252465,903152244,327024448,257349836,431774297,714479094,114462919,350021247,139934784,801799813,539551848,41576353,802846678,659340919,672408792,575487541,368702253,575350782,619200170,562982523,429350007,150624942,945196520,6319297,894502025,374999342,42453831,221527997,3616961,580046947,341012885,276733823,579878644,670787005,162408190,615289155,917206778,768258447,411299351,832373932,18107510,835563777,233896612,738886121,393809140,644927895,625120199,336985781,516279761,519173913,217509375,754521244,424721008,850139747,438485991,920679348,219842245,442535106,128849044,269625477,753265916,440813281,928041043,42439991,348131531,381592048,654276173,47700473,43991743,548001157,730540323,887685216,932597280,801983201,142570249,471547353,382906298,745576798,286187856,510059281,818604713,507091460,180496747,177275620,665329553,200027308,550501620,759618434,249162203,361024171,165602025,282541684,651186107,424876182,377693857,567673772,859949936,529948870,981890830,720454551,978293949,769899874,526887448,623249023,930481714,794783876,957457862,564435488,398080331,34219860,458598094,396377234,642505402,706573841,542803667,314649460,302333235,577578191,692607196,984717758,51551222,571746567,181633622,552463582,538957298,753601027,442293302,778981625,146456018,791162799,991018889,504501545,952054904,231445357,40511975,512943594,225912874,965164521,550144404,476038888,439900254,975615018,195831342,689809368,815767322,935245854,688904403,993516754,439936142,668910186,937167584,657694071,25151922,724216878,897387781,13684183,263175987,781446329,542562472,747813284,240316688,653466082,475627298,753619876,473430472,816466673,551717822,299623712,258591955,913624315,779342388,313872566,912716451,968634939,645165489,119540208,327353173,70136686,329097227,867149612,458779459,300118018,5895168,924317,136386920,430435362,268428795,484792722,744629423,976335117,799648412,269161423,809591377,443320513,337653395,591911355,627672258,275133634,701549043,952539877,22665917 };} else if (1) {nCr_computeFactinv(1000000);for (int n = 0; n < 1000; ++ n) {mint sum;for (int i = 0; i <= n; ++ i)sum += nCr(n * 3, i);//nCr(n * 7, i * 3);seq.push_back(sum);}} else {int n; string s;while (cin >> n >> s) {if (n != seq.size()) abort();seq.push_back(digitsMod(s.c_str(), (int)s.size(), mint::Mod));}}int N = (int)seq.size();cerr << "N = " << N << endl;PolynomiallyRecursiveSequence prs;prs.polynomialBasisKind = PolynomiallyRecursiveSequence::MonominalBasis;auto solution = prs.findMinimalSolution(seq, 100000);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 << ")";if (t > 0) cout << " ";}if (hi >= 0) {cout << "a_";if (i == 0)cout << "n";elsecout << "{n-" << i << "}";}if (i == 0) {cout << " = ";if (K == 1) cout << "0";}}cout << endl;mint multiplier = 1;//12;for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j)solution[i][j] *= multiplier;for (int j = D - 1; j >= 0; -- j) {int c = mintToSigned(solution[0][j]);if (c != 0) {if (c < 0) {for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j)solution[i][j] *= -1;}break;}}/*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]);cout << " },\n";}cout << "} };\n";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 << "vector<mint> seq(init.begin(), init.end());\n";cout << "seq.resize(N + 1);\n";cout << "for (int n = K - 1; n <= N; ++ n) {\n";cout << " mint ";for (int j = 1; j < D; ++ j) {if (j != 1) cout << ", ";cout << "n" << j << " = ";if (j == 1) cout << "n";else cout << "n" << (j - 1) << " * n1";}cout << ";\n";auto outputSum = [&](int i, bool negate) {int t = 0;for (int j = 0; j < D; ++ j)t += solution[i][j] != mint();if (t == 0) {cout << "0";return;}if (t > 1) cout << "(";int k = 0;for (int j = D - 1; j >= 0; -- j) {int c = mintToSigned(solution[i][j] * (negate ? -1 : 1));if (c != 0) {if (k ++ > 0) {if (c < 0) {cout << " - ";c = -c;} else {cout << " + ";}}if (c < 0) {cout << "-";c = -c;}if (j == 0) {cout << c;} else {cout << "n" << j;if (c != 1)cout << " * " << c;}}}if (t > 1) cout << ")";};cout << " mint sum;\n";for (int i = 1; i < K; ++ i) {cout << " sum ";bool negative = false;for (int j = D - 1; j >= 0; -- j) {int c = mintToSigned(-solution[i][j]);if (c != 0) {negative = c < 0;break;}}cout << (negative ? "-=" : "+=");cout << " seq[n - " << i << "] * ";outputSum(i, !negative);cout << ";\n";}cout << " seq[n] = sum / ";outputSum(0, false);cout << ";\n";cout << "}\n";} else {int N;while (~scanf("%d", &N)) {const int K = 8;const array<mint, K - 1> init = { 1, 1, 3, 4, 19, 26, 144 };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, n3 = n2 * n1, n4 = n3 * n1;mint sum;sum += seq[n - 1] * (n4 * 52269863 + n3 * 305931736 + n2 * 258365537 + n1 * 276782898 - 253220885);sum -= seq[n - 2] * (n4 * 242461620 + n3 * 27603152 + n2 * 139205790 + n1 * 267856601 + 277902371);sum -= seq[n - 3] * (n4 * 317850736 + n3 * 33389931 + n2 * 29110388 + n1 * 483840402 + 262098585);sum -= seq[n - 4] * (n4 * 446541231 + n3 * 262847746 + n2 * 474399727 - n1 * 329634661 + 301985484);sum += seq[n - 5] * (n4 * 179688158 - n3 * 122495876 + n2 * 117668564 + n1 * 427842771 + 117888425);sum -= seq[n - 6] * (n4 * 319613082 - n3 * 141115369 + n2 * 359231985 + n1 * 179539766 - 342110626);sum -= seq[n - 7] * (n4 * 76661316 - n3 * 336327801 + n2 * 475890910 + n1 * 251448872 - 347026197);seq[n] = sum / (n4 * 420387432 + n3 * 79612576 - n2 * 73247243 + n1);}mint ans = seq[N];printf("%d\n", ans.get());}}}