BZOJ-3557. [Ctsc2014]随机数

露露、花花和萱萱最近对计算机中的随机数产生了兴趣。为大家所熟知的是，由计算机生成的随机数序列并非真正的随机数，而是由一定法则生成的伪随机数。

某一天，露露了解了一种生成随机数的方法，称为 Mersenne twister。给定初始参数 $m \in \mathbb Z^+$，$x \in \mathbb Z \cap [0, 2^m)$ 和初值 $M_0 \in \mathbb Z \cap [0, 2^m)$，它通过下列递推式构造伪随机数列 ${M_n}$：

[ \begin{equation} M_n = \left { \begin{aligned} &2M_{n-1}& ~2M_{n-1} < 2^m \ (2M_{n-1} &-2^m) \oplus x & ~2M_{n-1} \geq 2^m \ \end{aligned} \right. \end{equation} ]

其中 $\oplus$ 是二进制异或运算（C/C++ 中的 ^ 运算）。而参数 $x$ 的选取若使得该数列在长度趋于无穷时，近似等概率地在 $\mathbb Z \cap (0, 2^m)$ 中取值，就称 $x$ 是好的。例如，在 $m > 1$ 时 $x=0$ 就显然不是好的。

在露露向伙伴介绍了 Mersenne twister 之后，花花想用一些经典的随机性测试来检验它的随机性强度。为此，花花使用计算机计算了一些 $M_k$。

但细心的萱萱注意到，花花在某次使用二进制输入 $k$ 时，在末尾多输入了 $l$ 个 $0$。她正想告诉花花这个疏忽，然而花花已经计算并记录了错误的 $M_k$ 值而没有记录 $k$ 的值。虽然这其实不是什么致命的问题，但是在萱萱告诉花花她的这个疏漏时，作为完美主义者的花花还是恳求萱萱帮她修正 $M_k$ 的值。萱萱便把这个任务交给了他的 AI ——你！

【输入格式】

输入文件 random.in 的第一行包含一个正整数 $m$，第二行为二进制表示的 $x$（共 $m$ 个数，从低位到高位排列），第三行为二进制表示的 $M_0$（排列方式同 $x$），第四行包含一个整数 $type$。

接下来分为两种可能的情况：

$type = 0$（萱萱记下了花花的输入）：则第五行包含一个整数，表示萱萱记下来的正确的 $k$ 值。
$type = 1$（萱萱未能记下花花的输入）：则第五行为 $l$，第六行输入花花计算出错误的二进制表示的 $M_k$。

【输出格式】

输出文件 random.out 仅一行，为一个 $m$ 位的 01 串，表示你求得的正确的 $M_k$（同样要求从低位到高位排列）

【数据范围】

对于 $type = 0$ 的部分，$m, k \leq 10^6$

对于 $type = 1$ 的部分，$m \leq 10^3, k \leq 10^{18}, l \leq 10$，$x$ 是“好的”

每个数据点时限 20s，并且开启编译器优化

首先对于 Mersenne twister 这个生成方法，由于是二进制的移位，异或，我们可以认为一个数是 $GF(2)$ 上的多项式（或者说系数是在 $\bmod 2$ 意义下的多项式）

这样对于一个数 $a$ 就对应一个多项式 $A(x)$，并且移位和异或这两个操作都可以对应到多项式的运算

$2a \Leftrightarrow xA(x)$
$a \oplus b \Leftrightarrow A(x) \pm B(x)$

再看题目中的部分，对于 $2M_{n-1}<2^m$ 的情况，对应的就是 $M_n(x) = xM_{n-1}(x)$

对于另外的一种情况，第一步也是乘 $x$，这时得到的多项式最高项次数是 $m$，我们可以考虑将其放在 $\bmod x^m$ 下，这样 $2M_{n-1} - 2^m \Leftrightarrow xM_{n-1}(x) \bmod x^m$，那剩下的异或就相当于是要减去一个 $x$ 对应的多项式 $X(x)$，这只要改成把多项式放在 $\bmod (x^m + X(x))$ 意义下就可以了，然后对于第一种情况，由于最高项次数小于 $m$，模不对其产生影响

因此整个操作就可以统一规约成 $M_n(x) = xM_{n-1}(x) \bmod (x^m + X(x))$

具体的解释？我们来看看样例

当 $m = 10, M_0 = (1100000111)_2, x = (0011100111)_2$ 时，可以得到对应的多项式

\[\begin{eqnarray*} M_0(x) &=& x^9 + x^8 + x^2 + x + 1 \\ X(x) &=& x^7 + x^6 + x^5 + x^2 + x + 1 \end{eqnarray*}\]

然后来模拟一下

\[\begin{eqnarray*} M_1(x) &=& xM_0(x) &\pmod {x^{10} + X(x)}\\ &=& x^{10} + x^9 + x^3 + x^2 + x &\pmod {x^{10} + X(x)}\\ &=& x^9 + x^7 + x^6 + x^5 + x^3 + 1 \end{eqnarray*}\]

好，现在来看 $type=0$ 的部分，就是给定 $M_0$ 求 $M_k$

这相当于求 $x^kM_0(x) \bmod (x^m + X(x))$，因此利用快速傅里叶变换可以在 $\mathcal O(m\log m \log k)$ 的复杂度内计算完（关于多项式除法和求模…… 我过一段时间再写好了）
然后是 $type = 1$ 的部分，已知 $l, X(x), M_0(x), M_{k2^l}(x)$，要求求出 $M_k(x)$

因为已经知道了 $x^{k2^l}M_0(x) \equiv M_{k2^l}(x) \pmod {x^m + X(x)}$

那么两边同时乘上 $M_0(x)$ 的逆元就可以得到
\[x^{k2^l} \equiv M_{k2^l}(x)M_0^{-1}(x) \pmod {x^m + X(x)}\]
然后由于 $x$ 是“好的”，也就是说前 $2^m - 1$ 个数字中间 $\mathbb Z \cap (0, 2^m)$ 都会出现过一次（否则最后不会近似等概率），这样就可以转换为 $ x \equiv x^{2^m} \pmod {x^m + X(x)} $，然后可以得到
\[\begin{eqnarray*} x^{k2^l2^{m-l}} &\equiv& \left(M_{k2^l}(x)M_0^{-1}(x)\right )^{2^{m-l}} &\pmod {x^m + X(x)} \\ x^k &\equiv& \left(M_{k2^l}(x)M_0^{-1}(x)\right )^{2^{m-l}} &\pmod {x^m + X(x)} \end{eqnarray*}\]
剩下最后只要在左右两边乘上 $M_0(x)$ 就可以得到答案，复杂度是 $\mathcal O(m^2\log m)$

至于如何求在 $\bmod (2^m + X(x))$ 下的逆元？这可以用扩展欧几里德算法（我还是过一段时间写）然后你需要不断地优化常数，另外为了不出现精度误差，你可以选一个质数做 FFT，但要求在 $GF(2)$ 下，因此每次计算好了乘法之后要转换到这里（也就是模 2，注意正负）否则会爆出去，还有要注意的就是在求逆元的时候不能两次连续的乘法而不进行转换，原因也是因为这个

然后下面是有点长的代码

/* BZOJ-3557: [Ctsc2014]随机数
 *   多项式乘、除、求模及GCD  */
#include <cstdio>
#include <algorithm>

using std::swap;
using std::fill;
using std::copy;
using std::reverse;
using std::reverse_copy;

const int MaxL = 21;
typedef int value_t;
typedef long long calc_t;
typedef long long power_t;
const value_t mod_base = 479, mod_exp = 21;
const value_t mod_v = (mod_base << mod_exp) + 1;
const value_t primitive_root = 3;

int fft_tot, type;
value_t eps[1 << MaxL], inv_eps[1 << MaxL];

void mod2(value_t& x)
{
	if(x >= (mod_v >> 1))
		x = (mod_v - x) & 1;
	else x &= 1;
}

value_t pow(value_t x, power_t p)
{
	value_t v = 1;
	for(; p; p >>= 1, x = (calc_t)x * x % mod_v)
		if(p & 1) v = (calc_t)x * v % mod_v;
	return v;
}

value_t inc(value_t x, value_t delta) 
{
	x += delta;
	return x >= mod_v ? x - mod_v : x;
}

value_t dec(value_t x, value_t delta)
{
	x -= delta;
	return x < 0 ? x + mod_v : x;
}

void init_eps(int tot)
{
	fft_tot = tot;
	value_t base = pow(primitive_root, (mod_v - 1) / tot);
	value_t inv_base = pow(base, mod_v - 2);
	eps[0] = inv_eps[0] = 1;
	for(int i = 1; i != tot; ++i)
	{
		eps[i] = (calc_t)eps[i - 1] * base % mod_v;
		inv_eps[i] = (calc_t)inv_eps[i - 1] * inv_base % mod_v;
	}
}

void transform(int n, value_t *x, value_t *w = eps)
{
	for(int i = 0, j = 0; i != n; ++i)
	{
		if(i > j) swap(x[i], x[j]);
		for(int l = n >> 1; (j ^= l) < l; l >>= 1);
	}

	for(int i = 2, t = 1; i <= n; t = i, i <<= 1)
	{
		int r = fft_tot / i;
		for(int p = 0; p < n; p += i)
		{
			for(int l = 0, s = 0; l != t; ++l, s += r)
			{
				value_t z = (calc_t)x[p + l + t] * w[s] % mod_v;
				x[p + l + t] = dec(x[p + l], z);
				x[p + l] = inc(x[p + l], z);
			}
		}
	}
}

void inverse_transform(int n, value_t *x)
{
	transform(n, x, inv_eps);
	value_t inv = pow(n, mod_v - 2);
	for(int i = 0; i != n; ++i)
		x[i] = (calc_t)x[i] * inv % mod_v;
//	for(int i = 0; i != n; ++i) mod2(x[i]);
}

// A(x)B(x) = 1 (mod x^n)
void polynomial_inverse(int n, value_t *A, value_t *B)
{
	if(n == 1)
	{
		B[0] = pow(A[0], mod_v - 2);
		return;
	}

	static value_t tmp[1 << MaxL];
	int p = 1;
	while(p < n << 1) p <<= 1;

	int half = (n + 1) >> 1;
	polynomial_inverse(half, A, B);
	fill(B + half, B + p, 0);
	transform(p, B);

	copy(A, A + n, tmp);
	fill(tmp + n, tmp + p, 0);
	transform(p, tmp);

	for(int i = 0; i != p; ++i)
		tmp[i] = (calc_t)B[i] * tmp[i] % mod_v;
	if(n > 100000)
	{
		inverse_transform(p, tmp);
		for(int i = 0; i != p; ++i) mod2(tmp[i]);
		transform(p, tmp);
	}
	for(int i = 0; i != p; ++i)
		B[i] = dec(2, (calc_t)B[i] * tmp[i] % mod_v);
	inverse_transform(p, B);
	for(int i = 0; i != n; ++i) mod2(B[i]);
	fill(B + n, B + p, 0);
}

// A(x) = B(x)D(x) + R(x)
void polynomial_division(int n, int m, value_t *A, value_t *B, value_t *D, value_t *R)
{
	static bool mod_solved = false;
	static value_t Z[1 << MaxL];
	static value_t A0[1 << MaxL], B0[1 << MaxL];

	int p = 1, t = n - m + 1;
	while(p < t << 1) p <<= 1;

	if(!mod_solved || type)
	{
		if(!type) mod_solved = true;
		fill(A0, A0 + p, 0);
		reverse_copy(B, B + m, A0);
		polynomial_inverse(t, A0, Z);
		fill(Z + t, Z + p, 0);
		for(int i = 0; i != t; ++i) mod2(Z[i]);
		transform(p, Z);
	}

	reverse_copy(A, A + n, A0);
	fill(A0 + t, A0 + p, 0);
	transform(p, A0);

	for(int i = 0; i != p; ++i)
		A0[i] = (calc_t)A0[i] * Z[i] % mod_v;
	inverse_transform(p, A0);
	reverse(A0, A0 + t);
	for(int i = 0; i != t; ++i) mod2(A0[i]);
	if(D) copy(A0, A0 + t, D);
	
	if(!R) return;
	for(p = 1; p < n; p <<= 1);
	if(t < p) fill(A0 + t, A0 + p, 0);
	transform(p, A0);

	copy(B, B + m, B0);
	fill(B0 + m, B0 + p, 0);
	transform(p, B0);

	for(int i = 0; i != p; ++i)
		A0[i] = (calc_t)A0[i] * B0[i] % mod_v;

	inverse_transform(p, A0);
	for(int i = 0; i != m; ++i)
	{
		R[i] = dec(A[i], A0[i]);
		mod2(R[i]);
	}

	fill(R + m, R + p, 0);
}

int mod_len;
value_t MOD[1 << MaxL], M0[1 << MaxL];
void multiply(int n, value_t *A, value_t *B, value_t *C)
{
	static value_t B0[1 << MaxL];
	int p = 1;
	while(p < n) p <<= 1;
	if(A == B)
	{
		int pos = p - 1;
		for(; pos && !A[pos]; --pos);
		for(p = 1; p <= pos << 1; p <<= 1);

		copy(A, A + p, C);
		transform(p, C);
		for(int i = 0; i != p; ++i)
			C[i] = (calc_t)C[i] * C[i] % mod_v;
		inverse_transform(p, C);
	} else {
		int pos1 = p - 1, pos2 = p - 1;
		for(; pos1 && !A[pos1]; --pos1);
		for(; pos2 && !B[pos2]; --pos2);
		for(p = 1; p <= pos1 + pos2; p <<= 1);

		copy(B, B + p, B0);
		copy(A, A + p, C);
		transform(p, C);
		transform(p, B0);
		for(int i = 0; i != p; ++i)
			C[i] = (calc_t)C[i] * B0[i] % mod_v;
		inverse_transform(p, C);
	}

	for(int i = 0; i != p; ++i) mod2(C[i]);
	if(p < n) fill(C + p, C + n, 0);
	bool flag = false;
	for(int i = mod_len - 1; i < p; ++i)
	{
		if(C[i]) 
		{
			flag = true;
			break;
		}
	}

	if(flag) polynomial_division(n, mod_len, C, MOD, 0, C);
}

// AX + BY = 1
typedef std::pair<int, int> len_t;
len_t polynomial_exgcd(int n, int m, value_t *A, value_t *B, value_t *X, value_t *Y, value_t *R)
{
	while(n && !A[n - 1]) --n;
	while(m && !B[m - 1]) --m;
	if(!m)
	{
		X[0] = 1, Y[0] = 0;
		return len_t(1, 1);
	}

	if(n < m) 
	{
		len_t len = polynomial_exgcd(m, n, B, A, Y, X, R);
		std::swap(len.first, len.second);
		return len;
	}

	int p = 1;
	while(p < (n - m + 1) + m) p <<= 1;
	value_t *D = new value_t[p << 1];

	polynomial_division(n, m, A, B, D, R);
	len_t l = polynomial_exgcd(m, m, B, R, X, Y, A);

	while(p < (n - m + 1) + l.second) p <<= 1;

	len_t ret_len;
	copy(X, X + l.first, A);
	fill(A + l.first, A + p, 0);
	copy(Y, Y + l.second, X);
	ret_len.first = l.second;

	fill(D + (n - m + 1), D + p, 0);
	transform(p, D);
	fill(Y + l.second, Y + p, 0);
	transform(p, Y);

	for(int i = 0; i != p; ++i)
		Y[i] = (calc_t)D[i] * Y[i] % mod_v;
	delete[] D;

	inverse_transform(p, Y);
	for(int i = 0; i != p; ++i)
	{
		Y[i] = dec(A[i], Y[i]);
		mod2(Y[i]);
	}

	int len = p;
	while(len && !Y[len - 1]) --len;
	ret_len.second = len;
	return ret_len;
}

value_t temp[5][1 << MaxL];
void solve0(int m)
{
	power_t k;
	std::scanf("%lld", &k);
	mod_len = m + 1;

	int p = 1;
	while(p < mod_len << 1) p <<= 1;
	init_eps(p);

	value_t *V = temp[0], *B = temp[1];
	int flag = 1, lenv = 0, lenb = 1;

	while(k)
	{
		if(!flag)
		{
			if(k & 1) multiply(mod_len << 1, B, V, V);
			if(k >> 1) multiply(mod_len << 1, B, B, B);
		} else {
			if(k & 1) lenv += lenb;
			lenb <<= 1;

			if(lenb > m >> 1 || lenv > m >> 1)
			{
				flag = 0;
				B[lenb] = 1;
				V[lenv] = 1;
			}
		} 

		k >>= 1;
	}

	if(flag) V[lenv] = 1;
	multiply(mod_len << 1, V, M0, M0);
	for(int i = 0; i != m; ++i)
		std::printf("%d", M0[i]);
}

void solve1(int m)
{
	int l;
	std::scanf("%d", &l);
	value_t *Mk = temp[0], *K = temp[1];
	value_t *X = temp[2], *Y = temp[3];
	mod_len = m + 1;

	int p = 1;
	while(p < std::max(2 << l, m << 1)) p <<= 1;
	init_eps(p << 1);

	copy(MOD, MOD + p, K);
	copy(M0, M0 + p, Mk);
	len_t len = polynomial_exgcd(m, m + 1, Mk, K, X, Y, temp[4]);

	for(int i = 0; i != m; ++i) std::scanf("%d", Mk + i);
	fill(Mk + m, Mk + p, 0);

	multiply(len.first + m, Mk, X, Mk);

	type = 0;
	for(int i = 0; i != m - l; ++i)
		multiply(mod_len << 1, Mk, Mk, Mk);

	multiply(mod_len << 1, M0, Mk, M0);
	for(int i = 0; i != m; ++i)
		std::printf("%d", M0[i]);
}

int main()
{
	int m;
	std::scanf("%d", &m);
	for(int i = 0; i != m; ++i)
		std::scanf("%d", MOD + i);
	MOD[m] = 1;
	for(int i = 0; i != m; ++i)
		std::scanf("%d", M0 + i);
	std::scanf("%d", &type);
	if(type == 0) solve0(m);
	else solve1(m);
	return 0;
}

本文遵守 CC BY-NC 4.0 许可协议。