#267807 (C++14) No.93 ペガサス

提出ソース

結果

問題	No.93 ペガサス
ユーザー	anta
提出日時	2018-06-22 13:48:19
言語	C++14 (gcc 13.2.0 + boost 1.83.0)
結果	AC
実行時間	2 ms / 5,000 ms
コード長	10,664 bytes
コンパイル時間	1,855 ms
コンパイル使用メモリ	180,492 KB
実行使用メモリ	4,384 KB
最終ジャッジ日時	2023-09-13 07:55:40
合計ジャッジ時間	3,261 ms
ジャッジサーバーID （参考情報）	judge12 / judge13

このコードへのチャレンジ（β）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 4,380 KB
testcase_01	AC	1 ms 4,380 KB
testcase_02	AC	1 ms 4,380 KB
testcase_03	AC	2 ms 4,376 KB
testcase_04	AC	1 ms 4,380 KB
testcase_05	AC	1 ms 4,376 KB
testcase_06	AC	2 ms 4,380 KB
testcase_07	AC	2 ms 4,376 KB
testcase_08	AC	1 ms 4,380 KB
testcase_09	AC	1 ms 4,380 KB
testcase_10	AC	1 ms 4,376 KB
testcase_11	AC	1 ms 4,384 KB
testcase_12	AC	2 ms 4,376 KB
testcase_13	AC	1 ms 4,380 KB
testcase_14	AC	2 ms 4,376 KB
testcase_15	AC	2 ms 4,380 KB
testcase_16	AC	2 ms 4,376 KB
testcase_17	AC	1 ms 4,380 KB
testcase_18	AC	1 ms 4,376 KB
testcase_19	AC	2 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);
	}
};

// idea: use falling factorial basis for polynomials

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

yukicoder

結果

テストケース

ソースコード