#267931 (C++14) No.660 家を通り過ぎないランダムウォーク問題

提出ソース

結果

問題	No.660 家を通り過ぎないランダムウォーク問題
ユーザー	anta
提出日時	2018-06-22 19:00:12
言語	C++14 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	60 ms / 2,000 ms
コード長	16,782 bytes
コンパイル時間	1,774 ms
コンパイル使用メモリ	183,752 KB
実行使用メモリ	4,384 KB
最終ジャッジ日時	2023-09-13 08:03:37
合計ジャッジ時間	3,644 ms
ジャッジサーバーID （参考情報）	judge11 / judge15

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	1 ms 4,376 KB
testcase_01	AC	2 ms 4,380 KB
testcase_02	AC	1 ms 4,384 KB
testcase_03	AC	2 ms 4,376 KB
testcase_04	AC	1 ms 4,376 KB
testcase_05	AC	1 ms 4,376 KB
testcase_06	AC	2 ms 4,376 KB
testcase_07	AC	1 ms 4,376 KB
testcase_08	AC	2 ms 4,376 KB
testcase_09	AC	1 ms 4,380 KB
testcase_10	AC	1 ms 4,376 KB
testcase_11	AC	2 ms 4,380 KB
testcase_12	AC	1 ms 4,380 KB
testcase_13	AC	1 ms 4,376 KB
testcase_14	AC	1 ms 4,380 KB
testcase_15	AC	1 ms 4,376 KB
testcase_16	AC	1 ms 4,376 KB
testcase_17	AC	1 ms 4,376 KB
testcase_18	AC	2 ms 4,376 KB
testcase_19	AC	1 ms 4,380 KB
testcase_20	AC	1 ms 4,376 KB
testcase_21	AC	1 ms 4,376 KB
testcase_22	AC	1 ms 4,380 KB
testcase_23	AC	2 ms 4,380 KB
testcase_24	AC	1 ms 4,376 KB
testcase_25	AC	1 ms 4,376 KB
testcase_26	AC	2 ms 4,376 KB
testcase_27	AC	1 ms 4,376 KB
testcase_28	AC	2 ms 4,376 KB
testcase_29	AC	2 ms 4,380 KB
testcase_30	AC	3 ms 4,376 KB
testcase_31	AC	2 ms 4,384 KB
testcase_32	AC	4 ms 4,376 KB
testcase_33	AC	4 ms 4,376 KB
testcase_34	AC	4 ms 4,380 KB
testcase_35	AC	10 ms 4,380 KB
testcase_36	AC	11 ms 4,376 KB
testcase_37	AC	11 ms 4,380 KB
testcase_38	AC	26 ms 4,380 KB
testcase_39	AC	27 ms 4,376 KB
testcase_40	AC	28 ms 4,380 KB
testcase_41	AC	46 ms 4,376 KB
testcase_42	AC	46 ms 4,380 KB
testcase_43	AC	54 ms 4,376 KB
testcase_44	AC	60 ms 4,376 KB

権限があれば一括ダウンロードができます

ソースコード

raw source code

#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;
	else
		return 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;
	else
		ss << 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";
				else
					cout << "{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());
		}
	}
}

yukicoder

結果

テストケース

ソースコード