Advertisement
peltorator

Modulo FFT

Oct 14th, 2019
503
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 1.92 KB | None | 0 0
  1. const ll MOD = 786433;
  2. const ll root = 3, st = 65536;
  3.  
  4. ll binpow(ll a, ll b)
  5. {
  6.     ll ans = 1;
  7.     while (b)
  8.     {
  9.         if (b & 1)
  10.             ans = ans * a % MOD;
  11.         a = a * a % MOD;
  12.         b >>= 1;
  13.     }
  14.     return ans;
  15. }
  16.  
  17. int rev(int x, int l)
  18. {
  19.     int ans = 0;
  20.     for (int i = 0; i < l; i++)
  21.         ans |= (((x >> i) & 1) << (l - i - 1));
  22.     return ans;
  23. }
  24.  
  25. const ll root_inv = binpow(root, MOD - 2);
  26.  
  27. void fft(vector<ll> &a, bool inv)
  28. {
  29.     int n = sz(a);
  30.     int len = 0;
  31.     while ((1 << len) < n)
  32.         len++;
  33.     for (int i = 0; i < n; i++)
  34.     {
  35.         int j = rev(i, len);
  36.         if (j > i)
  37.             swap(a[i], a[j]);
  38.     }
  39.     for (int l = 1; l < n; l *= 2)
  40.     {
  41.         ll w = binpow((inv ? root_inv : root), st / l / 2LL);
  42.         for (int i = 0; i < n; i += 2 * l)
  43.         {
  44.             ll wp = 1;
  45.             for (int j = 0; j < l; j++)
  46.             {
  47.                 ll u = a[i + j], v = a[i + j + l] * wp % MOD;
  48.                 a[i + j] = u + v;
  49.                 if (a[i + j] >= MOD)
  50.                     a[i + j] -= MOD;
  51.                 a[i + j + l] = u - v;
  52.                 if (a[i + j + l] < 0)
  53.                     a[i + j + l] += MOD;
  54.                 wp = wp * w % MOD;
  55.             }
  56.         }
  57.     }
  58.     if (inv)
  59.     {
  60.         ll t = binpow(n, MOD - 2);
  61.         for (int i = 0; i < n; i++)
  62.             a[i] = a[i] * t % MOD;
  63.     }
  64. }
  65.  
  66. vector<ll> mult(const vector<ll> &a, const vector<ll> &b)
  67. {
  68.     int nn = sz(a) + sz(b) + 10;
  69.     int n = 1;
  70.     while (n < 2 * nn)
  71.         n *= 2;
  72.     n = min(n, st);
  73.     vector<ll> aa(n, 0), bb(n, 0);
  74.     for (int i = 0; i < sz(a); i++)
  75.         aa[i] = a[i];
  76.     for (int i = 0; i < sz(b); i++)
  77.         bb[i] = b[i];
  78.     fft(aa, 0);
  79.     fft(bb, 0);
  80.     for (int i = 0; i < n; i++)
  81.         aa[i] = aa[i] * bb[i] % MOD;
  82.     fft(aa, 1);
  83.     while (sz(aa) && aa.back() == 0)
  84.         aa.pop_back();
  85.     return aa;
  86. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement