fork download
copy 
  1. /**
  2. Input:
  3. Dòng 1: Số testcase
  4. Dòng 2: Gồm hai số n, m
  5. Dòng 3: Gồm n+1 hệ số của đa thức A bậc n (bao gồm cả hệ số của hạng tử x^0)
  6. Dòng 4: Gồm m+1 hệ số của đa thức B bậc m (bao gồm cả hệ số của hạng tử x^0)
  7.  
  8. Output:
  9. Với mỗi testcase in ra 1 dòng duy nhất gồm n+m+1 hệ số của đa thức tích C = A x B
  10.  
  11. Giới hạn:
  12. 1 <= n, m <= (1e6-1)
  13. -1e6 <= Hệ số hạng tử <= 1e6
  14. Tổng n+m trên tất cả các testcase không vượt quá 2e6
  15.  
  16. Ví dụ:
  17.  
  18. Input:
  19. 2
  20. 2 4
  21. 1 2 3
  22. -4 1 2 5 -2
  23. 3 2
  24. 1 4 -4 -1
  25. 1 3 5
  26.  
  27. Output:
  28. -4 -7 -8 12 14 11 -6
  29. 1 7 13 7 -23 -5
  30. **/
  31.  
  32. #include <bits/stdc++.h>
  33. #define up(i,a,b) for (int i = (int)a; i <= (int)b; i++)
  34. using namespace std;
  35. using i32 = int32_t;
  36. using i64 = int64_t;
  37. using u32 = uint32_t;
  38. using u64 = uint64_t;
  39.  
  40. // p = c*2^k + 1, deu < 2^31 (Montgomery32 can mod < 2^31 de bat dau am bang bit 31)
  41. const u32 P1 = 2013265921, R1 = 31;   // 15*2^27+1
  42. const u32 P2 = 2130706433, R2 = 3;    // 127*2^24+1
  43. const i64 M  = (i64)P1 * P2;          // ~4.29e18 > 2*1e18
  44. const int MAXN = 1 << 21;             // n+m+1 <= 2e6+1 <= 2^21
  45.  
  46. u64 power_plain(u64 a, u64 e, u64 mod){
  47.     a %= mod;
  48.     u64 res = 1;
  49.     while (e){
  50.         if (e & 1) res = res * a % mod;
  51.         a = a * a % mod;
  52.         e >>= 1;
  53.     }
  54.     return res;
  55. }
  56.  
  57. // ================== MONTGOMERY (R = 2^32) ==================
  58. struct Montgomery32 {
  59.     u32 mod, inv_mod, r2;
  60.     void set_mod(u32 m){
  61.         mod = m;
  62.         inv_mod = 1;
  63.         for (int i = 0; i < 5; i++) inv_mod *= 2 - mod * inv_mod;
  64.         u64 r = ((u64)1 << 32) % mod;
  65.         r2 = (u32)(r * r % mod);
  66.     }
  67.     // REDC: 0 <= x < mod * 2^32
  68.     u32 reduce(u64 x) const {
  69.         u32 q = (u32)x * inv_mod;
  70.         u64 m = (u64)q * mod;
  71.         u32 y = (u32)((x - m) >> 32);
  72.         return (y >> 31) ? y + mod : y;
  73.     }
  74.     u32 to_mont(u32 a)    const { return reduce((u64)a * r2); }
  75.     u32 from_mont(u32 a)  const { return reduce((u64)a); }
  76.     u32 mul(u32 a, u32 b) const { return reduce((u64)a * b); }
  77.     u32 add(u32 a, u32 b) const { u32 s = a + b; return s >= mod ? s - mod : s; }
  78.     u32 sub(u32 a, u32 b) const { return a >= b ? a - b : a + mod - b; }
  79. };
  80. // =============================================================
  81.  
  82. // Mot bo NTT cho 1 modulo: mt + bang root/root_inv (dang Montgomery)
  83. struct NTT {
  84.     Montgomery32 mt;
  85.     u32 mod;
  86.     vector<u32> root, root_inv;
  87.  
  88.     void init(u32 p, u32 g, int n){
  89.         mod = p;
  90.         mt.set_mod(p);
  91.         root.resize(n); root_inv.resize(n);
  92.         root[1] = root_inv[1] = mt.to_mont(1);
  93.         u32 g_inv = (u32)power_plain(g, p - 2, p);
  94.         for (int k = 2; k * 2 <= n; k <<= 1){
  95.             u32 w_mont     = mt.to_mont((u32)power_plain(g,     (p - 1) / (2 * k), p));
  96.             u32 w_inv_mont = mt.to_mont((u32)power_plain(g_inv, (p - 1) / (2 * k), p));
  97.             for (int j = k / 2; j < k; j++){
  98.                 root[j * 2]         = root[j];
  99.                 root[j * 2 + 1]     = mt.mul(root[j], w_mont);
  100.                 root_inv[j * 2]     = root_inv[j];
  101.                 root_inv[j * 2 + 1] = mt.mul(root_inv[j], w_inv_mont);
  102.             }
  103.         }
  104.     }
  105.  
  106.     // DIF (Gentleman-Sande): vao tu nhien -> ra bit-reversed
  107.     void forward(vector<u32>& a) const {
  108.         int n = a.size();
  109.         for (int len = n; len >= 2; len >>= 1){
  110.             int half = len / 2;
  111.             for (int i = 0; i < n; i += len){
  112.                 for (int j = 0; j < half; j++){
  113.                     u32 w = root[half + j];
  114.                     u32 u = a[i + j], v = a[i + j + half];
  115.                     a[i + j]        = mt.add(u, v);
  116.                     a[i + j + half] = mt.mul(mt.sub(u, v), w);
  117.                 }
  118.             }
  119.         }
  120.     }
  121.  
  122.     // DIT (Cooley-Tukey): vao bit-reversed -> ra tu nhien
  123.     void inverse(vector<u32>& a) const {
  124.         int n = a.size();
  125.         for (int len = 2; len <= n; len <<= 1){
  126.             int half = len / 2;
  127.             for (int i = 0; i < n; i += len){
  128.                 for (int j = 0; j < half; j++){
  129.                     u32 w = root_inv[half + j];
  130.                     u32 u = a[i + j];
  131.                     u32 v = mt.mul(a[i + j + half], w);
  132.                     a[i + j]        = mt.add(u, v);
  133.                     a[i + j + half] = mt.sub(u, v);
  134.                 }
  135.             }
  136.         }
  137.         u32 n_inv = mt.to_mont((u32)power_plain(n, mod - 2, mod));
  138.         for (auto& x : a) x = mt.mul(x, n_inv);
  139.     }
  140.  
  141.     // tra ve he so tich modulo `mod`, o mien thuong (da from_mont)
  142.     vector<u32> multiply(const vector<i32>& a, const vector<i32>& b, int bound) const {
  143.         vector<u32> A(bound), B(bound);
  144.         up(i, 0, (int)a.size() - 1) A[i] = mt.to_mont(a[i] < 0 ? (u32)(a[i] + (i64)mod) : (u32)a[i]);
  145.         up(i, 0, (int)b.size() - 1) B[i] = mt.to_mont(b[i] < 0 ? (u32)(b[i] + (i64)mod) : (u32)b[i]);
  146.         forward(A);
  147.         forward(B);
  148.         up(i, 0, bound - 1) A[i] = mt.mul(A[i], B[i]);
  149.         inverse(A);
  150.         for (auto& x : A) x = mt.from_mont(x);
  151.         return A;
  152.     }
  153. };
  154.  
  155. NTT ntt1, ntt2;
  156. u32 inv_P1_mod_P2;
  157.  
  158. // ================= Garner CRT: 2 modulus =================
  159. // x = r1 + P1 * t2, tra ve gia tri co dau (|x| < M/2)
  160. i64 garner2(u32 r1, u32 r2){
  161.     u64 t2 = (u64)(r2 + (u64)P2 - r1) % P2 * inv_P1_mod_P2 % P2;
  162.     i64 x = (i64)r1 + (i64)P1 * (i64)t2;
  163.     if (x > M / 2) x -= M;
  164.     return x;
  165. }
  166.  
  167. void solve(){
  168.     int n, m; cin >> n >> m;
  169.     vector<i32> a(n + 1), b(m + 1);
  170.     for (auto& x : a) cin >> x;
  171.     for (auto& x : b) cin >> x;
  172.  
  173.     int need = (int)a.size() + (int)b.size() - 1;
  174.     int bound = 1;
  175.     while (bound < need) bound <<= 1;
  176.  
  177.     vector<u32> r1 = ntt1.multiply(a, b, bound);
  178.     vector<u32> r2 = ntt2.multiply(a, b, bound);
  179.     up(i, 0, need - 1) cout << garner2(r1[i], r2[i]) << " \n"[i == need - 1];
  180. }
  181.  
  182. signed main(){
  183.     ios_base::sync_with_stdio(false);
  184.     cin.tie(0);
  185.     #define Task "A"
  186.     if (fopen(Task".inp", "r")){
  187.         freopen(Task".inp", "r", stdin);
  188.         freopen(Task".out", "w", stdout);
  189.     }
  190.  
  191.     ntt1.init(P1, R1, MAXN);
  192.     ntt2.init(P2, R2, MAXN);
  193.     inv_P1_mod_P2 = (u32)power_plain(P1 % P2, P2 - 2, P2);
  194.  
  195.     int tt; cin >> tt;
  196.     while (tt--) solve();
  197. }
Success #stdin #stdout 0.02s 36000KB
comments (?)
stdin
copy 
2
2 4
1 2 3
-4 1 2 5 -2
3 2
1 4 -4 -1
1 3 5
stdout
copy 
-4 -7 -8 12 14 11 -6
1 7 13 7 -23 -5
https://ideone.com/TR6Fey
ntt with 2 modulos
language:
C++14 (gcc 8.3)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!