Advertisement
Guest User

Untitled

a guest
Nov 21st, 2016
886
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. typedef long double ld;
  3.  
  4. #define pb push_back
  5. #define mp make_pair
  6. #define eprintf(...) fprintf(stderr, __VA_ARGS__)
  7. #define sz(x) ((int)(x).size())
  8. #define TASKNAME "text"
  9.  
  10. const ld pi = acos((ld) -1);
  11.  
  12. //BEGIN ALGO
  13. namespace FFT {
  14.   struct com {
  15.     ld x, y;
  16.  
  17.     com(ld _x = 0, ld _y = 0) : x(_x), y(_y) {}
  18.  
  19.     inline com operator + (const com &c) const {
  20.       return com(x + c.x, y + c.y);
  21.     }
  22.     inline com operator - (const com &c) const {
  23.       return com(x - c.x, y - c.y);
  24.     }
  25.     inline com operator * (const com &c) const {
  26.       return com(x * c.x - y * c.y, x * c.y + y * c.x);
  27.     }
  28.     inline com conj() const {
  29.       return com(x, -y);
  30.     }
  31.   };
  32.  
  33.   const static int maxk = 21, maxn = (1 << maxk) + 1;
  34.   com ws[maxn];
  35.   int dp[maxn];
  36.   com rs[maxn];
  37.   int n, k;
  38.   int lastk = -1;
  39.  
  40.   void fft(com *a, bool torev = 0) {
  41.     if (lastk != k) {
  42.       lastk = k;
  43.       dp[0] = 0;
  44.  
  45.       for (int i = 1, g = -1; i < n; ++i) {
  46.         if (!(i & (i - 1))) {
  47.           ++g;
  48.         }
  49.         dp[i] = dp[i ^ (1 << g)] ^ (1 << (k - 1 - g));
  50.       }
  51.  
  52.       ws[1] = com(1, 0);
  53.       for (int two = 0; two < k - 1; ++two) {
  54.         ld alf = pi / n * (1 << (k - 1 - two));
  55.         com cur = com(cos(alf), sin(alf));
  56.  
  57.         int p2 = (1 << two), p3 = p2 * 2;
  58.         for (int j = p2; j < p3; ++j) {
  59.           ws[j * 2 + 1] = (ws[j * 2] = ws[j]) * cur;
  60.         }
  61.       }
  62.     }
  63.     for (int i = 0; i < n; ++i) {
  64.       if (i < dp[i]) {
  65.         swap(a[i], a[dp[i]]);
  66.       }
  67.     }
  68.     if (torev) {
  69.       for (int i = 0; i < n; ++i) {
  70.         a[i].y = -a[i].y;
  71.       }
  72.     }
  73.     for (int len = 1; len < n; len <<= 1) {
  74.       for (int i = 0; i < n; i += len) {
  75.         int wit = len;
  76.         for (int it = 0, j = i + len; it < len; ++it, ++i, ++j) {
  77.           com tmp = a[j] * ws[wit++];
  78.           a[j] = a[i] - tmp;
  79.           a[i] = a[i] + tmp;
  80.         }
  81.       }
  82.     }
  83.   }
  84.  
  85.   com a[maxn];
  86.   int mult(int na, int *_a, int nb, int *_b, long long *ans) {
  87.     if (!na || !nb) {
  88.       return 0;
  89.     }
  90.     for (k = 0, n = 1; n < na + nb - 1; n <<= 1, ++k) ;
  91.     assert(n < maxn);
  92.     for (int i = 0; i < n; ++i) {
  93.       a[i] = com(i < na ? _a[i] : 0, i < nb ? _b[i] : 0);
  94.     }
  95.     fft(a);
  96.     a[n] = a[0];
  97.     for (int i = 0; i <= n - i; ++i) {
  98.       a[i] = (a[i] * a[i] - (a[n - i] * a[n - i]).conj()) * com(0, (ld) -1 / n / 4);
  99.       a[n - i] = a[i].conj();
  100.     }
  101.     fft(a, 1);
  102.     int res = 0;
  103.     for (int i = 0; i < n; ++i) {
  104.       long long val = (long long) round(a[i].x);
  105.       assert(abs(val - a[i].x) < 1e-1);
  106.       if (val) {
  107.         assert(i < na + nb - 1);
  108.         while (res < i) {
  109.           ans[res++] = 0;
  110.         }
  111.         ans[res++] = val;
  112.       }
  113.     }
  114.     return res;
  115.   }
  116. };
  117. //END ALGO
Advertisement
RAW Paste Data Copied
Advertisement