Luogu4191 [CTSC2010]性能优化多项式，循环卷积

Posted 2021-12-30 athousandmoons

题目描述：设$A,B$为$n-1$次多项式，求$A*B^C$在系数模$n+1$，长度为$n$的循环卷积。

数据范围：$n\leq 5*10^5,C\leq 10^9$，且$n$的质因子不超过7，$n+1$为质数。

---

这就是一个循环卷积，在$n=2^k$的情况下可以直接使用FFT/NTT，但是这里不行。

由于$n$的质因子很小，设$n=2^k_1*3^k_2*5^k_3*7^k_4$，我们考虑将$n$分治为$p$份（NTT就是$p=2$的情况，这里$p=2,3,5,7$）

$$F(\omega_n^i)=\sum_i=0^p-1\omega_n^iF_i(\omega_\fracnp^i)$$

设$a_i=F(x)[x^i]$

$$F_r(x)=\sum_ia_ip+rx^i$$

而$rev$数组其实就是把模$p$同余的放在一起。

NTT的逆变换就是把次数上的$i$改成$-i$，不过代码实现里面我直接写了对$A[1:n-1]$进行reverse（这里说的就是反序）

 1 #include<bits/stdc++.h>
 2 #define Rint register int
 3 using namespace std;
 4 typedef long long LL;
 5 const int N = 500003;
 6 int n, C, mod, pri[N], tot, Wn[N];
 7 inline void add(int &a, int b)a += b; if(a >= mod) a -= mod;
 8 inline int kasumi(int a, int b)
 9     int res = 1;
10     while(b)
11         if(b & 1) res = (LL) res * a % mod;
12         a = (LL) a * a % mod;
13         b >>= 1;
14     
15     return res;
16 
17 inline void factor(int n)
18     for(Rint i = 2;i * i <= n;i ++)
19         if(!(n % i)) pri[++ tot] = i, n /= i, -- i;
20     if(n > 1) pri[++ tot] = n;
21 
22 inline int primitive()
23     for(Rint i = 2;;i ++)
24         bool flag = true;
25         for(Rint j = 1;j <= tot && flag;j ++)
26             if(kasumi(i, n / pri[j]) == 1) flag = false;
27         if(flag) return i;
28     
29 
30 int a[N], b[N], tmp[N];
31 inline void Rev(int *A)
32     for(Rint i = tot, block = n;i;block /= pri[i], i --)
33         for(Rint num = 0, j = 0;j < n;j += block)
34             for(Rint k = 0;k < pri[i];k ++)
35                 for(Rint l = 0;l < block;l += pri[i])
36                     tmp[num ++] = A[j + k + l];
37         for(Rint i = 0;i < n;i ++) A[i] = tmp[i];
38     
39 
40 inline void NTT(int *A, int type)
41     Rev(A);
42     for(Rint i = 1, block = 1;i <= tot;i ++)
43         int mid = block, wi = Wn[n / (block *= pri[i])];
44         for(Rint j = 0;j < n;j ++) tmp[j] = 0;
45         for(Rint j = 0;j < n;j += block)
46             int wk = 1;
47             for(Rint k = 0;k < block;k ++)
48                 for(Rint l = k % mid, w = 1;l < block;l += mid, w = (LL) w * wk % mod)
49                     add(tmp[j + k], (LL) w * A[j + l] % mod);
50                 wk = (LL) wk * wi % mod;
51             
52         
53         for(Rint j = 0;j < n;j ++) A[j] = tmp[j];
54     
55     if(type == -1)
56         std :: reverse(A + 1, A + n);
57         for(Rint i = 0;i < n;i ++)
58             A[i] = (LL) A[i] * n % mod;
59     
60 
61 int main()
62     scanf("%d%d", &n, &C); mod = n + 1;
63     for(Rint i = 0;i < n;i ++) scanf("%d", a + i);
64     for(Rint i = 0;i < n;i ++) scanf("%d", b + i);
65     factor(n);
66     Wn[0] = 1; Wn[1] = primitive();
67     for(Rint i = 2;i <= n;i ++) Wn[i] = (LL) Wn[i - 1] * Wn[1] % mod;
68     NTT(a, 1); NTT(b, 1);
69     for(Rint i = 0;i < n;i ++)
70         a[i] = (LL) a[i] * kasumi(b[i], C) % mod;
71     NTT(a, -1);
72     for(Rint i = 0;i < n;i ++)
73         printf("%d\n", a[i]);
74 

Luogu4191