fork download
copy 
  1. /**
  2.  * Đầu vào: Hai đa thức có độ dài n, m. 1 <= n, m <= 1e6-1 thể hiện bậc cao nhất của đa thức là n+1, m+1
  3.  Hệ số hạng tử nằm trong khoảng [-1e9; 1e9]
  4.  
  5.  * Đầu ra: Đa thức kết quả ở dạng đầy đủ có giá trị trong khoảng [-1e24; 1e24], không mod
  6.  **/
  7.  
  8. #include <bits/stdc++.h>
  9. #define up(i,a,b) for (int i = (int)a; i <= (int)b; i++)
  10. using namespace std;
  11. using i32 = int32_t;
  12. using i64 = int64_t;
  13. using u32 = uint32_t;
  14. using u64 = uint64_t;
  15. using i128 = __int128;
  16. using u128 = unsigned __int128;
  17.  
  18. // ================= Fast IO (USACO style) =================
  19. const int BUF_SZ = 1 << 20;
  20.  
  21. inline namespace Input {
  22.     char buf[BUF_SZ];
  23.     int pos, len;
  24.     inline char next_char() {
  25.         if (pos == len) {
  26.             pos = 0;
  27.             len = (int)fread(buf, 1, BUF_SZ, stdin);
  28.             if (!len) return EOF;
  29.         }
  30.         return buf[pos++];
  31.     }
  32.     inline i64 read_int() {
  33.         char ch;
  34.         while ((ch = next_char()) != '-' && !isdigit(ch));
  35.         bool neg = false;
  36.         if (ch == '-') { neg = true; ch = next_char(); }
  37.         i64 x = ch - '0';
  38.         while (isdigit(ch = next_char())) x = x * 10 + (ch - '0');
  39.         return neg ? -x : x;
  40.     }
  41. }
  42.  
  43. inline namespace Output {
  44.     char buf[BUF_SZ];
  45.     int pos;
  46.     inline void flush_out() {
  47.         fwrite(buf, 1, pos, stdout);
  48.         pos = 0;
  49.     }
  50.     inline void write_char(char c) {
  51.         if (pos == BUF_SZ) flush_out();
  52.         buf[pos++] = c;
  53.     }
  54.     inline void write_int(i64 x) {
  55.         static char num_buf[24];
  56.         if (x < 0) { write_char('-'); x = -x; }
  57.         int len = 0;
  58.         do { num_buf[len++] = (char)('0' + (x % 10)); x /= 10; } while (x);
  59.         while (len) write_char(num_buf[--len]);
  60.     }
  61.     // in i128 co dau: tach 1 lan theo 10^18 (1 phep chia 128-bit), phan con lai dung u64
  62.     inline void write_i128(i128 x) {
  63.         if (pos + 48 > BUF_SZ) flush_out();          // dam bao du cho (dau + <=39 chu so)
  64.         u128 ux;
  65.         if (x < 0) { buf[pos++] = '-'; ux = (u128)0 - (u128)x; } else ux = (u128)x;
  66.         const u64 E18 = 1000000000000000000ULL;
  67.         u64 hi = (u64)(ux / E18);
  68.         u64 lo = (u64)(ux - (u128)hi * E18);
  69.         char t[24];
  70.         if (hi) {
  71.             // hi co the dai toi 21 chu so neu |x| ~ 1e39; voi bai nay hi < 4e7
  72.             int l = 0;
  73.             do { t[l++] = (char)('0' + hi % 10); hi /= 10; } while (hi);
  74.             while (l) buf[pos++] = t[--l];
  75.             for (int d = 17; d >= 0; d--) { buf[pos + d] = (char)('0' + lo % 10); lo /= 10; }
  76.             pos += 18;                                // lo dem du 18 chu so
  77.         } else {
  78.             int l = 0;
  79.             do { t[l++] = (char)('0' + lo % 10); lo /= 10; } while (lo);
  80.             while (l) buf[pos++] = t[--l];
  81.         }
  82.     }
  83.     void init_output() { int r = atexit(flush_out); (void)r; }   // khong dung assert: -DNDEBUG se xoa mat atexit
  84. }
  85. // =========================================================
  86.  
  87. // 3 so nguyen to NTT-friendly chuan, cung primitive root g = 3 (deu < 2^31)
  88. const u32 P1 = 998244353;   // 119*2^23+1
  89. const u32 P2 = 167772161;   // 5*2^25+1
  90. const u32 P3 = 469762049;   // 7*2^26+1
  91. const u32 R  = 3;
  92. const i128 M = (i128)P1 * P2 * P3;
  93. const int MAXN = 1 << 21;   // n+m+1 <= 2e6+1 <= 2^21
  94.  
  95. u64 power_plain(u64 a, u64 e, u64 mod){
  96.     a %= mod;
  97.     u64 res = 1;
  98.     while (e){
  99.         if (e & 1) res = res * a % mod;
  100.         a = a * a % mod;
  101.         e >>= 1;
  102.     }
  103.     return res;
  104. }
  105.  
  106. // ================== MONTGOMERY (R = 2^32) ==================
  107. struct Montgomery32 {
  108.     u32 mod, inv_mod, r2;
  109.     void set_mod(u32 m){
  110.         mod = m;
  111.         inv_mod = 1;
  112.         for (int i = 0; i < 5; i++) inv_mod *= 2 - mod * inv_mod;
  113.         u64 r = ((u64)1 << 32) % mod;
  114.         r2 = (u32)(r * r % mod);
  115.     }
  116.     // REDC: 0 <= x < mod * 2^32
  117.     u32 reduce(u64 x) const {
  118.         u32 q = (u32)x * inv_mod;
  119.         u64 m = (u64)q * mod;
  120.         u32 y = (u32)((x - m) >> 32);
  121.         return (y >> 31) ? y + mod : y;
  122.     }
  123.     u32 to_mont(u32 a)    const { return reduce((u64)a * r2); }
  124.     u32 from_mont(u32 a)  const { return reduce((u64)a); }
  125.     u32 mul(u32 a, u32 b) const { return reduce((u64)a * b); }
  126.     u32 add(u32 a, u32 b) const { u32 s = a + b; return s >= mod ? s - mod : s; }
  127.     u32 sub(u32 a, u32 b) const { return a >= b ? a - b : a + mod - b; }
  128. };
  129. // =============================================================
  130.  
  131. // Mot bo NTT cho 1 modulo: mt + bang root/root_inv (dang Montgomery)
  132. struct NTT {
  133.     Montgomery32 mt;
  134.     u32 mod;
  135.     vector<u32> root, root_inv;
  136.  
  137.     void init(u32 p, u32 g, int n){
  138.         mod = p;
  139.         mt.set_mod(p);
  140.         root.resize(n); root_inv.resize(n);
  141.         root[1] = root_inv[1] = mt.to_mont(1);
  142.         u32 g_inv = (u32)power_plain(g, p - 2, p);
  143.         for (int k = 2; k * 2 <= n; k <<= 1){
  144.             u32 w_mont     = mt.to_mont((u32)power_plain(g,     (p - 1) / (2 * k), p));
  145.             u32 w_inv_mont = mt.to_mont((u32)power_plain(g_inv, (p - 1) / (2 * k), p));
  146.             for (int j = k / 2; j < k; j++){
  147.                 root[j * 2]         = root[j];
  148.                 root[j * 2 + 1]     = mt.mul(root[j], w_mont);
  149.                 root_inv[j * 2]     = root_inv[j];
  150.                 root_inv[j * 2 + 1] = mt.mul(root_inv[j], w_inv_mont);
  151.             }
  152.         }
  153.     }
  154.  
  155.     // DIF (Gentleman-Sande): vao tu nhien -> ra bit-reversed
  156.     void forward(vector<u32>& a) const {
  157.         int n = a.size();
  158.         for (int len = n; len >= 2; len >>= 1){
  159.             int half = len / 2;
  160.             for (int i = 0; i < n; i += len){
  161.                 for (int j = 0; j < half; j++){
  162.                     u32 w = root[half + j];
  163.                     u32 u = a[i + j], v = a[i + j + half];
  164.                     a[i + j]        = mt.add(u, v);
  165.                     a[i + j + half] = mt.mul(mt.sub(u, v), w);
  166.                 }
  167.             }
  168.         }
  169.     }
  170.  
  171.     // DIT (Cooley-Tukey): vao bit-reversed -> ra tu nhien
  172.     void inverse(vector<u32>& a) const {
  173.         int n = a.size();
  174.         for (int len = 2; len <= n; len <<= 1){
  175.             int half = len / 2;
  176.             for (int i = 0; i < n; i += len){
  177.                 for (int j = 0; j < half; j++){
  178.                     u32 w = root_inv[half + j];
  179.                     u32 u = a[i + j];
  180.                     u32 v = mt.mul(a[i + j + half], w);
  181.                     a[i + j]        = mt.add(u, v);
  182.                     a[i + j + half] = mt.sub(u, v);
  183.                 }
  184.             }
  185.         }
  186.         u32 n_inv = mt.to_mont((u32)power_plain(n, mod - 2, mod));
  187.         for (auto& x : a) x = mt.mul(x, n_inv);
  188.     }
  189.  
  190.     // tra ve he so tich modulo `mod`, o mien thuong (da from_mont)
  191.     vector<u32> multiply(const vector<i32>& a, const vector<i32>& b, int bound) const {
  192.         vector<u32> A(bound), B(bound);
  193.         // sửa multiply: đưa về [0, mod) đúng cách
  194.         up(i, 0, (int)a.size()-1){ i64 v = a[i] % (i64)mod; if (v < 0) v += mod; A[i] = mt.to_mont((u32)v); }
  195.         up(i, 0, (int)b.size()-1){ i64 v = b[i] % (i64)mod; if (v < 0) v += mod; B[i] = mt.to_mont((u32)v); }
  196.         forward(A);
  197.         forward(B);
  198.         up(i, 0, bound - 1) A[i] = mt.mul(A[i], B[i]);
  199.         inverse(A);
  200.         for (auto& x : A) x = mt.from_mont(x);
  201.         return A;
  202.     }
  203. };
  204.  
  205. NTT ntt1, ntt2, ntt3;
  206. u32 inv_P1_mod_P2, inv_P1P2_mod_P3;
  207.  
  208. // ================= Garner CRT: 3 modulus =================
  209. // x = r1 + P1*t2 + P1*P2*t3, tra ve gia tri co dau (|x| < M/2)
  210. i128 garner3(u32 r1, u32 r2, u32 r3){
  211.     // ghep (r1 mod P1) va (r2 mod P2): t2 = (r2 - r1) * P1^-1 mod P2
  212.     u64 t2 = (u64)(((i64)r2 - (i64)r1) % (i64)P2 + P2) % P2 * inv_P1_mod_P2 % P2;
  213.     u64 x12 = (u64)r1 + (u64)P1 * t2;                    // < P1*P2 ~ 1.7e17
  214.     // ghep x12 (mod P1*P2) voi (r3 mod P3): t3 = (r3 - x12) * (P1*P2)^-1 mod P3
  215.     u64 t3 = (u64)(((i64)r3 - (i64)(x12 % P3)) % (i64)P3 + P3) % P3 * inv_P1P2_mod_P3 % P3;
  216.     i128 x = (i128)x12 + (i128)P1 * P2 * t3;
  217.     if (x > M / 2) x -= M;
  218.     return x;
  219. }
  220.  
  221. void solve(){
  222.     int n = (int)read_int(), m = (int)read_int();
  223.     vector<i32> a(n + 1), b(m + 1);
  224.     for (auto& x : a) x = (i32)read_int();
  225.     for (auto& x : b) x = (i32)read_int();
  226.  
  227.     int need = (int)a.size() + (int)b.size() - 1;
  228.     int bound = 1;
  229.     while (bound < need) bound <<= 1;
  230.  
  231.     vector<u32> r1 = ntt1.multiply(a, b, bound);
  232.     vector<u32> r2 = ntt2.multiply(a, b, bound);
  233.     vector<u32> r3 = ntt3.multiply(a, b, bound);
  234.  
  235.     up(i, 0, need - 1){
  236.         write_i128(garner3(r1[i], r2[i], r3[i]));
  237.         write_char(i == need - 1 ? '\n' : ' ');
  238.     }
  239. }
  240.  
  241. signed main(){
  242.     init_output();
  243.     #define Task "A"
  244.     if (fopen(Task".inp", "r")){
  245.         freopen(Task".inp", "r", stdin);
  246.         freopen(Task".out", "w", stdout);
  247.     }
  248.  
  249.     ntt1.init(P1, R, MAXN);
  250.     ntt2.init(P2, R, MAXN);
  251.     ntt3.init(P3, R, MAXN);
  252.     inv_P1_mod_P2   = (u32)power_plain(P1 % P2, P2 - 2, P2);
  253.     inv_P1P2_mod_P3 = (u32)power_plain((u64)P1 * P2 % P3, P3 - 2, P3);
  254.  
  255.     int tt = (int)read_int();
  256.     while (tt--) solve();
  257. }
  258.  
Success #stdin #stdout 0.03s 52300KB
comments (?)
stdin
copy 
1
7 8
-4 -10 7 10 3 0 -3 -5
2 -9 -8 -3 -7 0 -4 10 9
stdout
copy 
-8 16 136 49 -82 -58 -93 -62 -116 -11 187 155 39 -10 -77 -45
https://ideone.com/lHUXXX
ntt with 3 modulos (big i128 output and USACO fastio)
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!