| 267 | } |
| 268 | |
| 269 | std::vector<long long> convolution_ll(const std::vector<long long>& a, |
| 270 | const std::vector<long long>& b) { |
| 271 | int n = int(a.size()), m = int(b.size()); |
| 272 | if (!n || !m) return {}; |
| 273 | |
| 274 | static constexpr unsigned long long MOD1 = 754974721; // 2^24 |
| 275 | static constexpr unsigned long long MOD2 = 167772161; // 2^25 |
| 276 | static constexpr unsigned long long MOD3 = 469762049; // 2^26 |
| 277 | static constexpr unsigned long long M2M3 = MOD2 * MOD3; |
| 278 | static constexpr unsigned long long M1M3 = MOD1 * MOD3; |
| 279 | static constexpr unsigned long long M1M2 = MOD1 * MOD2; |
| 280 | static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3; |
| 281 | |
| 282 | static constexpr unsigned long long i1 = |
| 283 | internal::inv_gcd(MOD2 * MOD3, MOD1).second; |
| 284 | static constexpr unsigned long long i2 = |
| 285 | internal::inv_gcd(MOD1 * MOD3, MOD2).second; |
| 286 | static constexpr unsigned long long i3 = |
| 287 | internal::inv_gcd(MOD1 * MOD2, MOD3).second; |
| 288 | |
| 289 | static constexpr int MAX_AB_BIT = 24; |
| 290 | static_assert(MOD1 % (1ull << MAX_AB_BIT) == 1, "MOD1 isn't enough to support an array length of 2^24."); |
| 291 | static_assert(MOD2 % (1ull << MAX_AB_BIT) == 1, "MOD2 isn't enough to support an array length of 2^24."); |
| 292 | static_assert(MOD3 % (1ull << MAX_AB_BIT) == 1, "MOD3 isn't enough to support an array length of 2^24."); |
| 293 | assert(n + m - 1 <= (1 << MAX_AB_BIT)); |
| 294 | |
| 295 | auto c1 = convolution<MOD1>(a, b); |
| 296 | auto c2 = convolution<MOD2>(a, b); |
| 297 | auto c3 = convolution<MOD3>(a, b); |
| 298 | |
| 299 | std::vector<long long> c(n + m - 1); |
| 300 | for (int i = 0; i < n + m - 1; i++) { |
| 301 | unsigned long long x = 0; |
| 302 | x += (c1[i] * i1) % MOD1 * M2M3; |
| 303 | x += (c2[i] * i2) % MOD2 * M1M3; |
| 304 | x += (c3[i] * i3) % MOD3 * M1M2; |
| 305 | // B = 2^63, -B <= x, r(real value) < B |
| 306 | // (x, x - M, x - 2M, or x - 3M) = r (mod 2B) |
| 307 | // r = c1[i] (mod MOD1) |
| 308 | // focus on MOD1 |
| 309 | // r = x, x - M', x - 2M', x - 3M' (M' = M % 2^64) (mod 2B) |
| 310 | // r = x, |
| 311 | // x - M' + (0 or 2B), |
| 312 | // x - 2M' + (0, 2B or 4B), |
| 313 | // x - 3M' + (0, 2B, 4B or 6B) (without mod!) |
| 314 | // (r - x) = 0, (0) |
| 315 | // - M' + (0 or 2B), (1) |
| 316 | // -2M' + (0 or 2B or 4B), (2) |
| 317 | // -3M' + (0 or 2B or 4B or 6B) (3) (mod MOD1) |
| 318 | // we checked that |
| 319 | // ((1) mod MOD1) mod 5 = 2 |
| 320 | // ((2) mod MOD1) mod 5 = 3 |
| 321 | // ((3) mod MOD1) mod 5 = 4 |
| 322 | long long diff = |
| 323 | c1[i] - internal::safe_mod((long long)(x), (long long)(MOD1)); |
| 324 | if (diff < 0) diff += MOD1; |
| 325 | static constexpr unsigned long long offset[5] = { |
| 326 | 0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3}; |