BZOJ-3557. [Ctsc2014]随机数

Posted by miskcoo on May 21, 2015

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

某一天,露露了解了一种生成随机数的方法,称为 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$。

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

  1. $type = 0$(萱萱记下了花花的输入):则第五行包含一个整数,表示萱萱记下来的正确的 $k$ 值。
  2. $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)$,并且移位和异或这两个操作都可以对应到多项式的运算

  1. $2a \Leftrightarrow xA(x)$
  2. $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*}\]
  1. 好,现在来看 $type=0$ 的部分,就是给定 $M_0$ 求 $M_k$

    这相当于求 $x^kM_0(x) \bmod (x^m + X(x))$,因此利用快速傅里叶变换可以在 $\mathcal O(m\log m \log k)$ 的复杂度内计算完(关于多项式除法和求模…… 我过一段时间再写好了)

  2. 然后是 $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,注意正负)否则会爆出去,还有要注意的就是在求逆元的时候不能两次连续的乘法而不进行转换,原因也是因为这个

然后下面是有点长的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
/* 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;
}