这是 CTSC2014 一道陈老师的可怕题目!
露露、花花和萱萱最近对计算机中的随机数产生了兴趣。为大家所熟知的是,由计算机生成的随机数序列并非真正的随机数,而是由一定法则生成的伪随机数。
某一天,露露了解了一种生成随机数的方法,称为 Mersenne twister。给定初始参数 ,
和初值
,它通过下列递推式构造伪随机数列
:
其中 是二进制异或运算(C/C++ 中的 ^ 运算)。而参数
的选取若使得该数列在长度趋于无穷时,近似等概率地在
中取值,就称
是好的。例如,在
时
就显然不是好的。
在露露向伙伴介绍了 Mersenne twister 之后,花花想用一些经典的随机性测试来检验它的随机性强度。为此,花花使用计算机计算了一些 。
但细心的萱萱注意到,花花在某次使用二进制输入 时,在末尾多输入了
个
。她正想告诉花花这个疏忽,然而花花已经计算并记录了错误的
值而没有记录
的值。虽然这其实不是什么致命的问题,但是在萱萱告诉花花她的这个疏漏时,作为完美主义者的花花还是恳求萱萱帮她修正
的值。萱萱便把这个任务交给了他的 AI ——你!
【输入格式】
输入文件 random.in 的第一行包含一个正整数 ,第二行为二进制表示的
(共
个数,从低位到高位排列),第三行为二进制表示的
(排列方式同
),第四行包含一个整数
。
接下来分为两种可能的情况:
(萱萱记下了花花的输入):则第五行包含一个整数,表示萱萱记下来的正确的
值。
(萱萱未能记下花花的输入):则第五行为
,第六行输入花花计算出错误的二进制表示的
。
【输出格式】
输出文件 random.out 仅一行,为一个 位的 01 串,表示你求得的正确的
(同样要求从低位到高位排列)
【数据范围】
对于 的部分,
对于 的部分,
,
是“好的”
每个数据点时限 20s,并且开启编译器优化
这是 CTSC2014 一道陈老师的可怕题目!(不要问我为什么要重复一遍,写了一天脑子已经不太清楚了QAQ)
首先对于 Mersenne twister 这个生成方法,由于是二进制的移位,异或,我们可以认为一个数是 上的多项式(或者说系数是在
意义下的多项式)
这样对于一个数 就对应一个多项式
,并且移位和异或这两个操作都可以对应到多项式的运算
再看题目中的部分,对于 的情况,对应的就是
对于另外的一种情况,第一步也是乘 ,这时得到的多项式最高项次数是
,我们可以考虑将其放在
下,这样
,那剩下的异或就相当于是要减去一个
对应的多项式
,这只要改成把多项式放在
意义下就可以了,然后对于第一种情况,由于最高项次数小于
,模不对其产生影响
因此整个操作就可以统一规约成
具体的解释?我们来看看样例
当 时,可以得到对应的多项式
然后来模拟一下
好,现在来看
的部分,就是给定
求
这相当于求 ,因此利用快速傅里叶变换可以在
的复杂度内计算完(关于多项式除法和求模…… 我过一段时间再写好了QAQ)
然后是
的部分,已知
,要求求出
因为已经知道了
那么两边同时乘上 的逆元就可以得到
然后由于 是“好的”,也就是说前
个数字中间
都会出现过一次(否则最后不会近似等概率),这样就可以转换为
,然后可以得到
剩下最后只要在左右两边乘上 就可以得到答案,复杂度是
至于如何求在 下的逆元?这可以用扩展欧几里德算法(我还是过一段时间写QAQ)然后你需要不断地优化常数,另外为了不出现精度误差,你可以选一个质数做 FFT,但要求在
下,因此每次计算好了乘法之后要转换到这里(也就是模 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; } |
好厉害,膜拜学长