[题解] [Mivik 的萌新赛 & Chino 的比赛 2020] Ame 和 Gura 的奇妙探险

给定 MT19937 的 10 个参数 $N$、$M$、$A$、$U$、$S$、$B$、$T$、$C$、$L$ 和 $F$，以及 $N$ 个引擎刚被初始化后生成的随机数，要求推断出初始化 MT19937 使用的种子。

$10\le M<N\le 2\times 10^5$，$0\le A,B,C<2^{32}$，$1\le U,S,T,L\le 31$，$1\le F<2^{32}$，保证 $F$ 是奇数

题目链接

---

思路

首先来拆解一下随机数生成器生成 $N$ 个随机数的过程。

1. 调用 init(uint32_t seed)，index 置为 $N$。
2. 第一次调用 gen() 时由于 index 为 $N$ 会先调用 twist() 函数，从而将 index 置为 0。
3. 总共 $N$ 次调用 gen()，第 $i$ 次调用会依次返回 temper(mt[i - 1])。

然后我们考虑设计出 temper(uint32_t x) 的反函数。我们观察这个函数：

static inline uint32_t temper(uint32_t x) {
	x ^= (x >> U);
	x ^= (x << S) & B;
	x ^= (x << T) & C;
	x ^= (x >> L);
	return x;
}

那么我们只需要实现下面两种操作：

• 根据 (x ^ (x >> L)) 还原出 x
• 根据 (x ^ ((x << L) & V)) 还原出 x

就可以反向进行 temper 操作了。

我们以第一种反向操作为例。我们注意到 (x ^ (x >> L)) 二进制下最高的 L 位和 x 是相同的，于是我们就可以立即得到 x 二进制下的最高 L 位。接下来，由于我们知道了 x 最高的 p1 位，我们可以通过异或还原出 x 最高的 2L 位，以此类推。下图简单解释了这一流程（L = 6）：

AAAAAABBBBBB--------------------|       --> x
      AAAAAABBBBBB--------------|------ --> (x >> L)

我们第一步能知道 A 的部分，然后通过异或就可以得出 B 的部分，以此类推。下面是简要的代码：

inline uint32_t inv_shift_right(uint32_t v, uint8_t bits) {
	if (bits == 32) return v;
	// cur 表示已经得出的部分，k 表示当前已经推导了前几位。
	uint32_t mask = (-1U) << (32 - bits), cur = v & mask;
	for (uint8_t k = bits; k < 32; k += bits)
		cur |= (v ^ (cur >> bits)) & (mask >>= bits);
	return cur;
}

同样地，第二种反向操作也可以通过相似的方式实现：

inline uint32_t inv_shift_left(uint32_t v, uint8_t bits, uint32_t m) {
	if (bits == 32) return v;
	uint32_t mask = (1U << bits) - 1, cur = v & mask;
	for (uint8_t k = bits; k < 32; k += bits)
		cur |= (v ^ ((cur << bits) & m)) & (mask <<= bits);
	return cur;
}

于是我们就可以反向 temper 操作来得到刚调用 twist() 函数后的 mt 数组中的 $N$ 个数。我们观察这个 twist 函数：

// 得到下一状态
inline void twist() {
	for (int i = 0; i < N; ++i) {
		uint32_t tmp = (mt[i] & (1U << 31)) | (mt[keep_in_range(i + 1)] & 0x7fffffffU);
		tmp = (tmp & 1)? ((tmp >> 1) ^ A): (tmp >> 1);
		mt[i] = mt[keep_in_range(i + M)] ^ tmp;
	}
	index = 0;
}

我们发现不行了：我们没法倒推出 tmp，因为我们并不知道 tmp 一开始到底是不是奇数，如果暴力枚举是 $O(2^N)$ 的。怎么办呢？

我们考虑从其它地方切入，观察这个 init() 函数：

// 使用种子初始化
inline void init(uint32_t seed) {
	mt[0] = seed;
	for (int i = 1; i < N; ++i)
		// 注意这里的无符号整型乘法溢出为定义行为，即得到的结果与真实结果在模 2 ^ 32 意义下同余
		mt[i] = F * (mt[i - 1] ^ (mt[i - 1] >> 30)) + i;
	index = N;
}

由于 $F$ 是奇数，那么 $F$ 在模 $2^{32}$ 意义下有逆元。使用这个逆元，我们可以通过 mt[i] 解出 (mt[i - 1] ^ (mt[i - 1] >> 30))，从而用上面提到的反向操作得出 mt[i - 1] 的值。也就是说我们只要知道了 mt[N - 1] 就可以推出整个调用 twist() 之前的 mt 数组，我们只需要枚举一下 mt[N - 1] 的值并检查一下就好了。时间复杂度 $O(N)$。

代码

mivik.h

代码