前言
题目链接：洛谷。
题意简述
\(n\) 个点，\(m\) 条边。每条边上有商店，经过一条边花费 \(1\) 单位时间，如果在边上的商店购物，额外花费 \(1\) 单位时间。需要购买 \(4\) 种物品，每个商店售出 \(1 \sim 4\) 种物品不等。请问，在 \(T\) 个单位时间内，从 \(1\) 出发购物，得到这 \(4\) 种物品，并回到 \(1\) 的方案数，对 \(5557\) 取模。你在购完并回到 \(1\) 后可以立即停止，不能在其他点逗留。
\(n \leq 25\)，\(m \leq 500\)，\(T \leq 10^9\)。
题目分析
四种物品很少，想到使用状压 DP。记 \(f_{i, u, \mathcal{S}}\) 表示经过了 \(i\) 个单位时间，现在在 \(u\)，目前购买了的物品集合为 \(\mathcal{S}\)，的方案数。初始状态 \(f_{0, 1, \varnothing} = 1\)，答案即为 \(\sum \limits _ {i = 0} ^ T f_{i, 1, \{1, 2, 3, 4\}}\)。
转移其实不难，对于 \((u, v, S')\) 这条边，枚举是否进行购买即可。
- 不购买\[f_{i + 1, v, S} \gets f_{i, u, S} \]
- 购买\[f_{i + 2, v, S \cup S'} \gets f_{i, u, S} \]
其中 \(a \gets b\) 表示将 \(b\) 累加到 \(a\) 上。
直接做是 \(\Theta(2^{|S|}T(n + m))\)，\(T\) 很大，肯定超时。转移这么简单，考虑使用矩阵快速幂优化。
由于 \(f_{i}\) 需要用到 \(f_{i - 1}\) 和 \(f_{i - 2}\)，并且需要额外记录 \(ans\)，所以我们转移是从 \(\displaystyle \begin{bmatrix} f_{i - 2} & f_{i - 1} & ans \end{bmatrix}\) 转移到 \(\displaystyle \begin{bmatrix} f_{i - 1} & f_{i} & ans' \end{bmatrix}\)。
考虑构造转移矩阵，注意到，把状态 \(\lambda\) 的值乘上 \(k\) 累加到 \(\varphi\)，相当于在转移矩阵的第 \(\lambda\) 行，第 \(\varphi\) 列累加一个 \(k\)。（采用转移矩阵放置在左边的形式。）当然，我们这里所有 \(k = 1\)。这一部分预处理是简单的。考虑什么时候累加答案。发现，如果 \(v = 1\)，并且目标状态 \(S\) 或者 \(S \cup S'\) 是全集，不光要累加到 \(f\) 上，也要累加到 \(ans\) 上。\(f_{i - 1}\) 放到答案矩阵左侧这一步，可以在转移矩阵中套一个单位矩阵。
\[\begin{gathered} {\large \left [ \begin{array}{ccc:ccc:c} f_{i-2, 1, \varnothing} & \cdots & f_{i-2, n, \{1, 2, 3, 4\}} & f_{i-1, 1, \varnothing} & \cdots & f_{i-1, n, \{1, 2, 3, 4\}} & ans \\ \end{array} \right ]} \\ {\Huge \times} \\ {\large \left [ \begin{array}{ccc:ccc:c} 0 & \cdots & 0 & & & & \\ \vdots & \ddots & \vdots & & & & \\ 0 & \cdots & 0 & & & & \\ \hdashline 1 & \cdots & 0 & & & & \\ \vdots & \ddots & \vdots & & & & \\ 0 & \cdots & 1 & & & & \\ \hdashline 0 & \cdots & 0 & 0 & \cdots & 0 & 1 \\ \end{array} \right ]} \\ {\Huge \Downarrow} \\ {\large \left [ \begin{array}{ccc:ccc:c} f_{i-1, 1, \varnothing} & \cdots & f_{i-1, n, \{1, 2, 3, 4\}} & f_{i, 1, \varnothing} & \cdots & f_{i, n, \{1, 2, 3, 4\}} & ans' \\ \end{array} \right ]} \end{gathered} \]构造转移矩阵的代码
// a => false: f[i-2], true: f[i-1]
inline int cal(bool a, int u, int st) {
return (u - 1) * (1 << 4) + st + (a ? n * (1 << 4) : 0);
}
B = 2 * n * (1 << 4); // ans 的位置
for (int u = 1; u <= n; ++u)
for (int st = 0; st < 1 << 4; ++st) {
toadd(trans(cal(true, u, st), cal(false, u, st)), 1); // 单位矩阵
for (int o = xym.head[u]; o; o = xym[o].nxt) {
int v = xym[o].to;
toadd(trans(cal(true, u, st), cal(true, v, st)), 1);
if (v == 1 && st == 0b1111)
toadd(trans(cal(true, u, st), B), 1);
toadd(trans(cal(false, u, st), cal(true, v, st | xym[o].len)), 1);
if (v == 1 && (st | xym[o].len) == 0b1111)
toadd(trans(cal(false, u, st), B), 1);
}
}
toadd(trans(B, B), 1); // 右下角 ans 累加
时间复杂度降到了：\(\Theta\Big((n 2^{|S|})^3 \log T\Big) = \Theta(n^316^{|S|}\log T)\)。还是会超时，但是能拿到 \(70\)pts。
考场上，我想到这里就不会做了。事实上，我们发现 \(16^{|S|}\) 显然是算法瓶颈所在，在矩阵乘法立方的加持下，不能记录 \(2^{|S|}\)。那怎么做呢？
我们可以在一开始钦定某些商品不能购买，得到在「钦定意义」下的「至多」购买 \(x\) 个商品的方案 \(f(x)\)。我们求的是「恰好」购买 \(x = 4\) 个商品的方案 \(g(x)\)。使用二项式反演即可，参见我的学习笔记，适用 \(\mathtt{Thm}\ 2.1\)，这里给出公式：
\[f(x) = \sum \limits _ {i = 0} ^ x \dbinom{n - i}{n - x} g(i) \Longleftrightarrow g(x) = \sum _ {i = 0} ^ x (-1)^{x - i} \binom{n - i}{n - x} f(i) \]至于「钦定」后怎么求出 \(f(x)\)，同样使用矩阵快速幂，不过这次矩阵大小是 \(2n + 1\)，总的时间复杂度为 \(\Theta(2^{|S|}n^3\log T)\)。
代码
$70$ 分代码
#include <cstdio>
#include <cstring>
#include <unordered_map>
#include <limits>
#include <iostream>
using namespace std;
const int N = 30, M = 520;
namespace Mod_Int_Class {
template <typename T, typename _Tp>
constexpr bool in_range(_Tp val) {
return std::numeric_limits<T>::min() <= val && val <= std::numeric_limits<T>::max();
}
template <auto mod = 1000000007, typename T = int, typename S = long long>
class Mod_Int {
static_assert(in_range<T>(mod), "mod must in the range of type T.");
static_assert(std::is_integral<T>::value, "type T must be an integer.");
static_assert(std::is_integral<S>::value, "type S must be an integer.");
public:
constexpr Mod_Int() noexcept = default;
template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value>>
constexpr Mod_Int(_Tp v) noexcept: val(0) {
if (0 <= v && v < mod) val = v;
else val = (v % mod + mod) % mod;
}
constexpr T const& raw() const {
return this -> val;
}
static constexpr inline T add(T a, T b) {
return a >= mod - b ? a + b - mod : a + b;
}
static constexpr inline T sub(T a, T b) {
return a < b ? a - b + mod : a - b;
}
static constexpr inline T mul(T a, T b) {
return static_cast<S>(a) * b % mod;
}
static constexpr inline T& toadd(T& a, T b) {
return a = add(a, b);
}
static constexpr inline T& tosub(T& a, T b) {
return a = sub(a, b);
}
static constexpr inline T& tomul(T& a, T b) {
return a = mul(a, b);
}
template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value>>
static constexpr inline T pow(T a, _Tp p) {
T res = 1;
for (; p; p >>= 1, tomul(a, a))
if (p & 1) tomul(res, a);
return res;
}
template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value>>
static constexpr inline bool is_prime(_Tp val) {
if (val < 2) return false;
for (_Tp i = 2; i * i <= val; ++i)
if (val % i == 0)
return false;
return true;
}
static constexpr inline T inv(T a) requires (is_prime(mod)) {
assert(a != 0);
return pow(a, mod - 2);
}
constexpr Mod_Int& operator + () const {
return *this;
}
constexpr Mod_Int operator - () const {
return sub(0, val);
}
constexpr Mod_Int inv() const {
return inv(val);
}
constexpr friend inline Mod_Int operator + (Mod_Int a, Mod_Int b) {
return add(a.val, b.val);
}
constexpr friend inline Mod_Int operator - (Mod_Int a, Mod_Int b) {
return sub(a.val, b.val);
}
constexpr friend inline Mod_Int operator * (Mod_Int a, Mod_Int b) {
return mul(a.val, b.val);
}
constexpr friend inline Mod_Int operator / (Mod_Int a, Mod_Int b) {
return mul(a.val, inv(b.val));
}
template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value>>
constexpr friend inline Mod_Int operator ^ (Mod_Int a, _Tp p) {
return pow(a.val, p);
}
constexpr friend inline Mod_Int& operator += (Mod_Int& a, Mod_Int b) {
return a = add(a.val, b.val);
}
constexpr friend inline Mod_Int& operator -= (Mod_Int& a, Mod_Int b) {
return a = sub(a.val, b.val);
}
constexpr friend inline Mod_Int& operator *= (Mod_Int& a, Mod_Int b) {
return a = mul(a.val, b.val);
}
constexpr friend inline Mod_Int& operator /= (Mod_Int& a, Mod_Int b) {
return a = mul(a.val, inv(b.val));
}
template <typename _Tp, typename = std::enable_if_t<std::is_integral<_Tp>::value>>
constexpr friend inline Mod_Int& operator ^= (Mod_Int& a, _Tp p) {
return a = pow(a.val, p);
}
constexpr friend inline bool operator == (Mod_Int a, Mod_Int b) {
return a.val == b.val;
}
constexpr friend inline bool operator != (Mod_Int a, Mod_Int b) {
return a.val != b.val;
}
constexpr Mod_Int& operator ++ () {
this -> val + 1 == mod ? this -> val = 0 : ++this -> val;
return *this;
}
constexpr Mod_Int& operator -- () {
this -> val == 0 ? this -> val = mod - 1 : --this -> val;
return *this;
}
constexpr Mod_Int operator ++ (int) {
Mod_Int res = *this;
this -> val + 1 == mod ? this -> val = 0 : ++this -> val;
return res;
}
constexpr Mod_Int operator -- (int) {
Mod_Int res = *this;
this -> val == 0 ? this -> val = mod - 1 : --this -> val;
return res;
}
friend std::istream& operator >> (std::istream& is, Mod_Int<mod, T, S>& x) {
T ipt;
return is >> ipt, x = ipt, is;
}
friend std::ostream& operator << (std::ostream& os, Mod_Int<mod, T, S> x) {
return os << x.val;
}
protected:
T val;
};
using mint = Mod_Int<5557, short, int>;
constexpr mint operator ""_m (unsigned long long x) {
return mint(x);
}
constexpr mint operator ""_mod (unsigned long long x) {
return mint(x);
}
}
using namespace Mod_Int_Class;
int n, m, t, B;
struct Graph {
struct node {
int to, len, nxt;
} edge[M];
int tot = 1, head[N];
void add(int u, int v, int w) {
edge[++tot] = {v, w, head[u]};
head[u] = tot;
}
inline node & operator [] (const int x) {
return edge[x];
}
} xym;
inline int cal(bool a, int u, int st) {
return (u - 1) * (1 << 4) + st + (a ? n * (1 << 4) : 0);
}
using umap = unordered_map<short, mint>;
struct Matrix {
umap val[801];
umap & operator [] (int x) {
return val[x];
}
mint & operator () (int a, int b) {
// a 转移到 b
return val[a][b];
}
inline Matrix friend operator * (Matrix &a, Matrix &b) {
static Matrix res;
for (register int i = 0; i <= B; ++i) res[i].clear();
for (register int i = 0; i <= B; ++i)
for (auto& [k, av]: a[i])
for (auto& [j, bv]: b[k])
res[i][j] += av * bv;
return res;
}
} trans, base;
inline int calc(char c) {
switch (c) {
case 'B': return 0;
case 'J': return 1;
case 'M': return 2;
case 'P': return 3;
default: return -1;
}
}
signed main() {
scanf("%d%d", &n, &m);
for (int i = 1, u, v; i <= m; ++i) {
static char str[8];
scanf("%d%d%s", &u, &v, str);
int w = 0;
for (int j = 0; str[j]; ++j) w |= 1 << calc(str[j]);
xym.add(u, v, w);
}
B = 2 * n * (1 << 4);
scanf("%d", &t);
for (int u = 1; u <= n; ++u)
for (int st = 0; st < 1 << 4; ++st) {
++trans(cal(true, u, st), cal(false, u, st));
for (int o = xym.head[u]; o; o = xym[o].nxt) {
int v = xym[o].to;
++trans(cal(true, u, st), cal(true, v, st));
if (v == 1 && st == 0b1111)
++trans(cal(true, u, st), B);
++trans(cal(false, u, st), cal(true, v, st | xym[o].len));
if (v == 1 && (st | xym[o].len) == 0b1111)
++trans(cal(false, u, st), B);
}
}
++trans(B, B);
base[0][cal(true, 1, 0)] = 1;
for (; t; trans = trans * trans, t >>= 1) {
if (t & 1) base = base * trans;
}
printf("%hd", base[0][B].raw());
return 0;
}
#include <cstdio>
#include <cstring>
#include <unordered_map>
#include <limits>
#include <iostream>
using namespace std;
const int N = 30, M = 520;
namespace Mod_Int_Class {
// see 70 pts code
}
using namespace Mod_Int_Class;
int n, m, t, B;
struct Graph {
struct node {
int to, len, nxt;
} edge[M];
int tot = 1, head[N];
void add(int u, int v, int w) {
edge[++tot] = {v, w, head[u]};
head[u] = tot;
}
inline node & operator [] (const int x) {
return edge[x];
}
} xym;
inline int cal(bool a, int u) {
return (u - 1) + (a ? n : 0);
}
struct Matrix {
mint val[51][51];
mint * operator [] (const int x) {
return val[x];
}
mint const * operator [] (const int x) const {
return val[x];
}
mint & operator () (const int a, const int b) {
// a 转移到 b
return val[a][b];
}
Matrix() {
memset(val, 0x00, sizeof (val));
}
inline Matrix friend operator * (const Matrix &a, const Matrix &b) {
static Matrix res;
memset(res.val, 0x00, sizeof (res.val));
for (register int i = 0; i <= B; ++i)
for (register int k = 0; k <= B; ++k)
for (register int j = 0; j <= B; ++j)
res[i][j] += a[i][k] * b[k][j];
return res;
}
} trans, base;
inline int calc(char c) {
switch (c) {
case 'B': return 0;
case 'J': return 1;
case 'M': return 2;
case 'P': return 3;
default: return -1;
}
}
inline bool subset(int a, int b) {
return (a & b) == a;
}
inline mint solve(int S) {
memset(trans.val, 0x00, sizeof (trans.val));
memset(base.val, 0x00, sizeof (base.val));
for (int u = 1; u <= n; ++u){
++trans(cal(true, u), cal(false, u));
for (int o = xym.head[u]; o; o = xym[o].nxt) {
int v = xym[o].to;
++trans(cal(true, u), cal(true, v));
if (v == 1)
++trans(cal(true, u), B);
if (subset(xym[o].len, S)) {
++trans(cal(false, u), cal(true, v));
if (v == 1)
++trans(cal(false, u), B);
}
}
}
++trans(B, B);
base[0][cal(true, 1)] = 1;
for (int x = t; x; trans = trans * trans, x >>= 1) {
if (x & 1) base = base * trans;
}
return base[0][B];
}
mint f[4], ans;
signed main() {
scanf("%d%d", &n, &m);
for (int i = 1, u, v; i <= m; ++i) {
static char str[8];
scanf("%d%d%s", &u, &v, str);
int w = 0;
for (int j = 0; str[j]; ++j) w |= 1 << calc(str[j]);
xym.add(u, v, w);
}
B = 2 * n;
scanf("%d", &t);
for (int st = 0; st < 1 << 4; ++st)
f[__builtin_popcount(st)] += solve(st);
for (int i = 0; i <= 4; ++i) {
// C(4 - i, 4 - 4) = 1
if ((4 - i) & 1)
ans -= f[i];
else
ans += f[i];
}
printf("%hu", ans.raw());
return 0;
}