[CTSC2010] 性能优化 题解

不难发现就是求 \(A*B^C\) 的长度为 \(n\) 的循环卷积（这都看不出来还配做这个题？）

一种及其 naive 的想法是倍增多项式快速幂模拟循环卷积。

因为值域很小，直接用 FFT 是足够的，复杂度 \(O(n\log n\log C)\)，怎么看都过不了。

我们知道 FFT 本质上就是做的循环卷积。

所以说如果我们能求出长度为 \(n\) 的 DFT 和 IDFT，那么我们直接 DFT 了 \(A,B\)，然后把 \(B\) 的点值求个 \(C\) 次方，然后再 IDFT 回去。套一个 Bluestein 就是 \(O(n\log n)\) 的了。

这样复杂度是 \(O(9n\log n)\) 的，有点点卡常。

但是显然有一个性质我们没有用到，就是 \(n\) 能表示成若干个不超过 \(10\) 的正整数的乘积。

不妨考虑普通 FFT 的本质，是把它分成奇偶两部分 \(A_0(x)=a_0x^0+a_2x^1...,A_1(x)=a_1x^0+a_2x^1...\)，然后递归下去算，最后这个点的 \(A(x)=A_0(x^2)+xA_1(x^2)\)。

类似的，我们可以把 \(A\) 分成 \(p\) 份，\(A_0,A_1...A_p\)，递归下去，回溯的时候 \(O(p)\) 合并。

总复杂度是 \(O(\sum np)\) 的，因为 \(n\) 的质因子只有 \(2,3,5,7\)，所以这个复杂度是可以接受的。

#include <bits/stdc++.h>
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define RI register int
typedef long long LL;

#define FILEIO(name) freopen(name".in", "r", stdin), freopen(name".out", "w", stdout);

using namespace std;

namespace IO {
  char buf[1000000], *p1 = buf, *p2 = buf;
  inline char gc() {
    if (p1 == p2) p2 = (p1 = buf) + fread(buf, 1, 1000000, stdin);
    return p1 == p2 ? EOF : *(p1++);
  }
  template <class T> inline void read(T &n) {
    n = 0; RI ch = gc(), f;
    while ((ch < '0' || ch > '9') && ch != '-') ch = gc();
    f = (ch == '-' ? ch = gc(), -1 : 1);
    while (ch >= '0' && ch <= '9') n = n * 10 + (ch ^ 48), ch = gc();
    n *= f;
  }
  char Of[105], *O1 = Of, *O2 = Of;
  template <class T> inline void print(T n, char ch = '\n') {
    if (n < 0) putchar('-'), n = -n;
    if (n == 0) putchar('0');
    while (n) *(O1++) = (n % 10) ^ 48, n /= 10;
    while (O1 != O2) putchar(*(--O1));
    putchar(ch);
  }
}

using IO :: read;
using IO :: print;

int const MAXN = 5e5 + 5;
int mod;
LL omega[MAXN];
LL a[MAXN], b[MAXN], t[MAXN];

LL qpow(LL a, LL k) {
  LL re = 1;
  for (; k; k >>= 1, a = a * a % mod)
    if (k & 1) re = re * a % mod;
  return re;
}
int Prime[105];
int Getwn() {
  int tmp = mod - 1, cnt = 0, sq = sqrt(mod - 1);
  for (RI i = 2; i <= sq; ++i)
    if (tmp % i == 0) {
      Prime[++cnt] = i;
      while (tmp % i == 0)
        tmp /= i;
    }
  if (tmp != 1) Prime[++cnt] = tmp;
  for (RI i = 2; 666; ++i) {
    int flag = 1;
    for (RI j = 1; flag && j <= cnt; ++j)
      if (qpow(i, (mod - 1) / Prime[j]) == 1)
        flag = 0;
    if (flag) return i;
  }
}
void FFT(LL *a, int len, int l, int r, int op) {
  if (len == 1) return void();
  for (RI i = l; i <= r; ++i) t[i] = a[i];
  int p = -1;
  if (len % 7 == 0) p = 7;
  if (len % 5 == 0) p = 5;
  if (len % 3 == 0) p = 3;
  if (len % 2 == 0) p = 2;
  int L[10] = {0}, R[10] = {0};
  int cur = l, sz = len / p;
  for (RI i = 0; i < p; ++i)
    L[i] = cur, R[i] = cur - 1, cur += sz;
  for (RI i = l; i <= r; ++i)
    a[++R[(i - l) % p]] = t[i];
  for (RI i = 0; i < p; ++i) 
    FFT(a, sz, L[i], R[i], op);
  for (RI i = l; i <= r; ++i) t[i] = a[i];
  LL wn = (mod - 1) / len, w = 0;
  for (RI i = l; i <= r; ++i) {
    LL tmp = 0; a[i] = 0;
    for (RI j = (i - l) % (len / p) + l; j <= r; j += sz) {
      a[i] = (a[i] + t[j] * omega[tmp] % mod) % mod;
      tmp = (tmp + w) % (mod - 1);
    }
    w = (w + op * wn + mod - 1) % (mod - 1);
  }
}

int main() {
  
#ifdef LOCAL
  FILEIO("a");
#endif

  int n, C; read(n), read(C); mod = n + 1;
  LL wn = Getwn(); omega[0] = 1;
  for (RI i = 1; i < n; ++i)
    omega[i] = omega[i - 1] * wn % mod;
  for (RI i = 0; i < n; ++i) read(a[i]);
  for (RI i = 0; i < n; ++i) read(b[i]);
  FFT(a, n, 0, n - 1, 1);
  FFT(b, n, 0, n - 1, 1);
  for (RI i = 0; i < n; ++i)
    a[i] = a[i] * qpow(b[i], C) % mod;
  FFT(a, n, 0, n - 1, -1);
  for (RI i = 0; i < n; ++i)
    print(a[i] * qpow(n, mod - 2) % mod);

  return 0;
}

// created by Daniel yuan
/*
     ________
    /        \
   / /      \ \
  / /        \ \
  \            /
   \  ______  /
    \________/
*/

Daniel_yuan's Blog

[CTSC2010] 性能优化题解