博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
Luogu4191 [CTSC2010]性能优化【多项式,循环卷积】
阅读量:5349 次
发布时间:2019-06-15

本文共 2975 字,大约阅读时间需要 9 分钟。

题目描述:设$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_{\frac{n}{p}}^i)$$

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

$$F_r(x)=\sum_{i}a_{ip+r}x^i$$

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

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

1 #include
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

 

转载于:https://www.cnblogs.com/AThousandMoons/p/10991952.html

你可能感兴趣的文章
【VS开发】ATL辅助COM组件开发
查看>>
FlatBuffers In Android
查看>>
《演说之禅》I &amp; II 读书笔记
查看>>
thinkphp3.2接入支付宝支付接口(PC端)
查看>>
response和request
查看>>
【转】在Eclipse中安装和使用TFS插件
查看>>
回到顶部浮窗设计
查看>>
C#中Monitor和Lock以及区别
查看>>
【NOIP2017】奶酪
查看>>
$ 一步一步学Matlab(3)——Matlab中的数据类型
查看>>
5.6.3.7 localeCompare() 方法
查看>>
Linux下好用的简单实用命令
查看>>
常用web字体的使用指南
查看>>
描绘应用程序级的信息
查看>>
poj2406-Power Strings
查看>>
2018/12/18 JS会像Linux一样改变编程
查看>>
php环境搭建脚本
查看>>
FTP主动模式与被动模式说明
查看>>
php 编译常见错误
查看>>
MES架构
查看>>