BZOJ 1919 [Ctsc2010]性能优化

题目：
http://www.lydsy.com/JudgeOnline/problem.php?id=1919

题意：
给出两个长度为$n$的整数序列$a[0..n-1],b[0..n-1]$和非负整数$C$。
对于两个长度为$n$的整数序列，定义$*$运算，结果为一个长度为$n$的整数序列，例如$f*g=h$，则有$h[k]=\sum_{i+j\equiv k\pmod n}f[i]\cdot g[j]$。
求$a*b*b*\cdots*b$每一位模$(n+1)$的值，其中有$C$个$*$运算，$(n+1)$是质数，$n$的质因数大小均不超过$10$。
$n\leq 5\cdot10^5,a[i],b[i],C\leq 10^9$

题解：
由原根的性质可知，长度为$n$的FFT即可支持$*$运算，难点在于$b^C$使得值域过大，即使能够快速计算长度为$n$的FFT，使用复数运算的FFT也很难得到精确的答案。

先考虑如何快速计算长度为$n$的FFT。
当$n=2^k$时，FFT每次是将序列一分为二，然后利用分治的技巧来进行合并。
因此当$n=2^{k_1}\cdot 3^{k_2}\cdot 5^{k_3}\cdot 7^{k_4}$时，FFT每次可能将序列一分为$p(p=2,3,5,7)$，合并时的式子需要重新推导。
不妨设是将$p$个长度为$n$的式子合并成一个长度为$p\cdot n$的式子，即利用$p$组$n$个点值得到$pn$个点值。
由于分裂时将模$p$意义相同的部分放在了一起，所以对于合并后的多项式

$$F(x)=\sum_{0\leq i<pn}{a_ix^i}$$
拆分的$p$个多项式分别为
$$F_r(x)=\sum_{0\leq i<n}{a_{ip+r}x^i}$$
故有

$$F(\omega^{an+b}_{pn})=\sum_{0\leq r<p}{(\omega^{an+b}_{pn})^rF_r(\omega^{b}_{n})}$$
于是可以 $O(p)$合并出每个点的值，而这样的分治层数是 $O(\sum_{i=1}^{4}{k_i})=O(\log n)$，每层的复杂度是 $O(pn)=O(7n)$，因此整体的复杂度是 $O(n\log n)$。
上述方法也可非递归实现，在分裂过程中注意每段之间互不影响，在合并过程中注意存储方式即可，笔者的做法就是迭代的做法。

再考虑解决精度问题，由同余关系的性质，可以使得每次计算相乘时的值域降低到$O(n^2)$，但需要将单位复根映射到模意义下的剩余系中。
由于$(n+1)$是质数，$\varphi(n+1)=n$，所以在模$(n+1)$意义下存在原根$g$，使得$g\equiv\omega_n\pmod{n+1}$，于是利用NTT代替FFT计算即可。
由于模$(n+1)$意义下原根的数量为$\varphi(n)= n\prod_{p_i\;is\;prime,p_i|n}{1-\frac{1}{p_i}}$，而$n$的质因子大小不超过10，所以期望检查$\frac{35}{8}\approx 5$次就可以找到原根了。

上述做法基于$n$是10-smooth number，即Cooley–Tukey FFT algorithm，而对于更强性质的$n$，可以使用Bluestein’s algorithm。

代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const int maxn = 500001;
int n, m, mod, tot, p[maxn], pw[maxn], a[maxn], b[maxn];
int mod_pow(int x, int k)
{
    int ret = 1;
    for( ; k > 0; k >>= 1, x = (LL)x * x % mod)
        if(k & 1)
            ret = (LL)ret * x % mod;
    return ret;
}
void NTT(int x[maxn], int flag)
{
    // go deeper
    static int y[maxn] = {};
    int *cur = x, *nxt = y;
    for(int i = tot - 1, delta = n / p[i]; i > 0; --i, delta /= p[i], swap(cur, nxt))
        for(int j = 0, *np = nxt; j < n; j += delta * p[i])
            for(int k = 0; k < p[i]; ++k)
                for(int l = 0, *cp = cur + j + k; l < delta; ++l, ++np, cp += p[i])
                    *np = *cp;
    // recursion
    for(int i = 0, clen = 1, nlen = p[i]; i < tot; ++i, clen = nlen, nlen *= p[i], swap(cur, nxt))
        for(int j = 0, k = 0, ww = 1, delta = 0; j < n; ++j, k = k + 1 < clen ? k + 1 : 0, ww = (LL)ww * pw[i] % mod, delta = delta + nlen > j ? delta : delta + nlen)
        {
            nxt[j] = 0;
            for(int t = 0, www = 1; t < nlen; t += clen, www = (LL)www * ww % mod)
                nxt[j] = (nxt[j] + (LL)www * cur[delta + t + k]) % mod;
        }
    if(flag == -1)
    {
        reverse(cur + 1, cur + n);
        for(int i = 0; i < n; ++i)
            cur[i] = (LL)cur[i] * n % mod; // n * n mod (n + 1) = 1
    }
    if(cur != x)
        memcpy(x, cur, n * sizeof(int));
}
int main()
{
    int tmp;
    scanf("%d%d", &n, &m);
    mod = n + 1;
    tmp = n;
    m = (m - 1) % n + 1;
    for(int i = 2; i * i <= tmp; ++i)
        for( ; tmp % i == 0; tmp /= i, p[tot++] = i);
    if(tmp > 1)
        p[tot++] = tmp;
    for(int ori = 2; ; ++ori)
    {
        bool flag = 1;
        for(int i = 0; i < tot && flag; ++i)
            if(!i || p[i - 1] != p[i])
                flag &= mod_pow(ori, n / p[i]) != 1;
        if(flag)
        {
            pw[tot - 1] = ori;
            for(int i = tot - 2; i >= 0; --i)
                pw[i] = mod_pow(pw[i + 1], p[i + 1]);
            break;
        }
    }
    for(int i = 0; i < n; ++i)
        scanf("%d", a + i);
    NTT(a, 1);
    for(int i = 0; i < n; ++i)
        scanf("%d", b + i);
    NTT(b, 1);
    for(int i = 0; i < n; ++i)
        a[i] = (LL)a[i] * mod_pow(b[i], m) % mod;
    NTT(a, -1);
    for(int i = 0; i < n; ++i)
        printf("%d\n", a[i]);
    return 0;
}