模板大全

__int128 mu = -1ull / mod;
inline ll reduce(ll x) {
	ll r = x - (mu * x >> 64) * mod;
	return r >= mod ? r - mod : r;
}

vector<vector<int>> find_rel(const vector<int> &a, int deg) {
	int n = a.size(), B = (n + 2) / (deg + 2), C = B * (deg + 1), R = n - (B - 1), c = 0;
	assert(B >= 2 && R >= C - 1);
	vector<vector<int>> b(R, vector<int>(C));
	for (int i = 0; i < R; i++) {
		for (int j = 0; j < B; j++) {
			for (int k = 0, x = a[i + j]; k <= deg; k++, x = (ll)x * (i + j) % mod) {
				b[i][j * (deg + 1) + k] = x;
			}
		}
	}
	for (c = 0; c < C; c++) {
		int p = -1;
		for (int i = c; i < R; i++) {
			if (b[i][c]) {
				p = i;
				break;
			}
		}
		if (p == -1) break;
		swap(b[p], b[c]);
		int w = qpow(b[c][c]);
		for (int i = c; i < C; i++) {
			b[c][i] = (ll)b[c][i] * w % mod;
		}
		for (int i = c + 1; i < R; i++) if (b[i][c]) {
			int t = b[i][c];
			for (int j = c; j < C; j++) {
				b[i][j] = (b[i][j] + mod - (ll)b[c][j] * t % mod) % mod;
			}
		}
	}
	assert(c != C);
	for (int i = c - 1; ~i; i--) if (b[i][c]) {
		for (int j = i - 1; ~j; j--) {
			b[j][c] = (b[j][c] + mod - (ll)b[i][c] * b[j][i] % mod) % mod;
		}
	}
	int od = c / (deg + 1);
	vector<vector<int>> ret(od + 1, vector<int>(deg + 1));
	ret[0][c % (deg + 1)] = 1;
	for (int i = c - 1; ~i; i--) {
		ret[od - i / (deg + 1)][i % (deg + 1)] = mod - b[i][c];
	}
	for (int i = 0; i <= od; i++) {
		vector<int> tmp(deg + 1);
		for (int k = 0; k <= deg; k++) {
			for (int j = k, s = 1; j <= deg; j++) {
				tmp[k] = (tmp[k] + (ll)s * ret[i][j]) % mod;
				s = (ll)s * (mod - i) % mod * (j + 1) % mod * qpow(j + 1 - k) % mod;
			}
		}
		ret[i] = tmp;
	}
	return ret;
}