#include using namespace std; // Fast Fourier transform // https://cp-algorithms.com/algebra/fft.html // https://drive.google.com/file/d/1B9BIfATnI_qL6rYiE5hY9bh20SMVmHZ7/view using cpx = complex; const double PI = acos(-1); vector roots = {{0, 0}, {1, 0}}; void ensure_capacity(int min_capacity) { for (int len = roots.size(); len < min_capacity; len *= 2) { for (int i = len >> 1; i < len; i++) { roots.emplace_back(roots[i]); double angle = 2 * PI * (2 * i + 1 - len) / (len * 2); roots.emplace_back(cos(angle), sin(angle)); } } } void fft(vector &z, bool inverse) { int n = z.size(); assert((n & (n - 1)) == 0); ensure_capacity(n); for (unsigned i = 1, j = 0; i < n; i++) { int bit = n >> 1; for (; j >= bit; bit >>= 1) j -= bit; j += bit; if (i < j) swap(z[i], z[j]); } for (int len = 1; len < n; len <<= 1) { for (int i = 0; i < n; i += len * 2) { for (int j = 0; j < len; j++) { cpx root = inverse ? conj(roots[j + len]) : roots[j + len]; cpx u = z[i + j]; cpx v = z[i + j + len] * root; z[i + j] = u + v; z[i + j + len] = u - v; } } } if (inverse) for (int i = 0; i < n; i++) z[i] /= n; } vector multiply_bigint(const vector &a, const vector &b, int base) { int need = a.size() + b.size(); int n = 1; while (n < need) n <<= 1; vector p(n); for (size_t i = 0; i < n; i++) { p[i] = cpx(i < a.size() ? a[i] : 0, i < b.size() ? b[i] : 0); } fft(p, false); // a[w[k]] = (p[w[k]] + conj(p[w[n-k]])) / 2 // b[w[k]] = (p[w[k]] - conj(p[w[n-k]])) / (2*i) vector ab(n); cpx r(0, -0.25); for (int i = 0; i < n; i++) { int j = (n - i) & (n - 1); ab[i] = (p[i] * p[i] - conj(p[j] * p[j])) * r; } fft(ab, true); vector result(need); long long carry = 0; for (int i = 0; i < need; i++) { long long d = (long long)(ab[i].real() + 0.5) + carry; carry = d / base; result[i] = d % base; } return result; } vector multiply_mod(const vector &a, const vector &b, int m) { int need = a.size() + b.size() - 1; int n = 1; while (n < need) n <<= 1; vector A(n); for (size_t i = 0; i < a.size(); i++) { int x = (a[i] % m + m) % m; A[i] = cpx(x & ((1 << 15) - 1), x >> 15); } fft(A, false); vector B(n); for (size_t i = 0; i < b.size(); i++) { int x = (b[i] % m + m) % m; B[i] = cpx(x & ((1 << 15) - 1), x >> 15); } fft(B, false); vector fa(n); vector fb(n); for (int i = 0, j = 0; i < n; i++, j = n - i) { cpx a1 = (A[i] + conj(A[j])) * cpx(0.5, 0); cpx a2 = (A[i] - conj(A[j])) * cpx(0, -0.5); cpx b1 = (B[i] + conj(B[j])) * cpx(0.5, 0); cpx b2 = (B[i] - conj(B[j])) * cpx(0, -0.5); fa[i] = a1 * b1 + a2 * b2 * cpx(0, 1); fb[i] = a1 * b2 + a2 * b1; } fft(fa, true); fft(fb, true); vector res(need); for (int i = 0; i < need; i++) { long long aa = (long long)(fa[i].real() + 0.5); long long bb = (long long)(fb[i].real() + 0.5); long long cc = (long long)(fa[i].imag() + 0.5); res[i] = (aa % m + (bb % m << 15) + (cc % m << 30)) % m; } return res; }