FFT-NTT

灰太狼 2022-07-16 00:40 143阅读 0赞

看了很久文档，觉得自己只学会了套模板的能力，理解的代码是怎么写，还有一点原理，看完现在来推一下原理估计又不会了！![难过][sad.gif]

学这个的原因是因为codechef的一道题目，可惜现在还是没有解决![难过][sad.gif]，谁会了求教[  点击打开链接][Link 1]

看了[  ACdream][ACdream]   的博客，觉得不够详细，而且看了之后根本看不懂代码里面写的什么鬼，之后找了[     一份题解 ][Link 2]

然后在百度文库找到了一篇讲得很详细的[     文档    ][Link 3]，看玩总算理解了那么一点点！![难过][sad.gif]

ps：这里我只能写一点代码要做什么，原理不敢亵渎！

1、首先是反转置换，将一个2的n次幂长度的分成两个，奇数是一组，偶数是一组，然后一次递归下去，直到长度为1，这样每个元素的位置会去到他二进制反转之后变成是十进制的位置，做置换可以采用rader算法线性实现！一下引用上面链接博客的一张图：

![Center][]

其实为什么这样想想也很明显，将位置pos编号换成二进制，初始化ans\_pos=0，ans依次做位运算右移，低位为0的时候走向二叉树的左儿子，否则是右儿子，而ans\_pos依次左移，向二叉树左儿子走的时候末尾补0，否则补1，这样ans\_pos的二进制就和刚开始ans的二进制是相反的

rader算法代码如下：

**\[cpp\]**  [view plain][]   [copy][view plain]

[print][view plain] [?][view plain]

1.  ///雷德算法，2^M=len，将第i位的数与“i的二进制反转之后的位”的数交换  
2.  **void** rader(complex \*F,**int** len)  
3.  \{  
4.      **int** j=len/2;///模拟二进制反转进位的的位置  
5.      **for**(**int** i=1;i<len-1;i++)  
6.      \{  
7.          **if**(i<j)swap(F\[i\],F\[j\]);///该出手时就出手  
8.          **int** k=len/2;  
9.          **while**(j>=k)  
10.         \{  
11.             j-=k;  
12.             k>>=1;  
13.         \}  
14.         **if**(j<k)j+=k;  
15.     \}  
16. \}

然后抄了份模板，具体解释就看上面链接ppt的解释，非常详细：

**\[cpp\]**  [view plain][]   [copy][view plain]

[print][view plain] [?][view plain]

1.  **const** **int** N = 500005;  
2.  **const** **double** pi = acos(-1.0);  
3.  **struct** complex\{ ///复数类  
4.      **double** real,imag;  
5.      complex (**double** r=0.0,**double** i=0.0)  
6.      \{  
7.          real=r;  
8.          imag=i;  
9.      \}  
10.     complex operator + (**const** complex &x)  
11.     \{  
12.         **return** complex(x.real+real,x.imag+imag);  
13.     \}  
14.     complex operator - (**const** complex &x)  
15.     \{  
16.         **return** complex(real-x.real,imag-x.imag);  
17.     \}  
18.     complex operator \* (**const** complex &x)  
19.     \{  
20.         **return** complex(x.real\*real-x.imag\*imag,imag\*x.real+real\*x.imag);  
21.     \}  
22. \}vara\[N\],varb\[N\];  
23. **void** FFT(complex \*F,**int** len,**int** t)  
24. \{  
25.     rader(F,len);  
26.     **for**(**int** h=2;h<=len;h<<=1)///L from 1 to M,2^M=len  
27.     \{  
28.         ///J\_product\_two\_power\_of\_M\_subtract\_L  
29.         ///M\_subtract\_L\_product\_h\_equals\_to\_len  
30.         complex wn(cos(t\*2\*pi/h),sin(t\*2\*pi/h));///公比  
31.         **for**(**int** j=0;j<len;j+=h)  
32.         \{  
33.             complex E(1,0);///螺旋因子  
34.             **for**(**int** k=j;k<j+h/2;k++)  
35.             \{  
36.                 complex u=F\[k\];///蝶型操作  
37.                 complex v=E\*F\[k+h/2\];  
38.                 F\[k\]=u+v;///前半部分  
39.                 F\[k+h/2\]=u-v;///后半部分  
40.                 E=E\*wn;  
41.             \}  
42.         \}  
43.     \}  
44.     **if**(t==-1)///IDFT  
45.     **for**(**int** i=0;i<len;i++)  
46.         F\[i\].real/=len;  
47. \}

NTT也是搞一发模板：

[ACdream 点击打开链接][ACdream]

[点击打开链接][Link 4]这个一份好的解释看完可以了解一二了

大数运算NTT：

**\[cpp\]**  [view plain][]   [copy][view plain]

[print][view plain] [?][view plain]

1.  **const** **int** N = 400005;  
2.  **const** **int** g=3;  
3.  **const** **long** **long** MOD=(479<<21)+1;  
4.  **int** len;  
5.  **char** A\[N\],B\[N\];  
6.  **long** **long** a\[N\],b\[N\],qp\[30\];  
7.  **long** **long** q\_pow(**long** **long** x,**long** **long** y,**long** **long** P)  
8.  \{  
9.      **long** **long** ans=1;  
10.     **while**(y>0)  
11.     \{  
12.         **if**(y&1)ans=ans\*x%P;  
13.         x=x\*x%P;  
14.         y>>=1;  
15.     \}  
16.     **return** ans;  
17. \}  
18. **void** init()  
19. \{  
20.     len=1;  
21.     **for**(**int** i=0;i<21;i++)  
22.     \{  
23.         **int** t=1<<i;  
24.         qp\[i\]=q\_pow(g,(MOD-1)/t,MOD);  
25.     \}  
26.     **int** len1=strlen(A);  
27.     **int** len2=strlen(B);  
28.     **while**(len<=2\*len1||len<=2\*len2)len<<=1;  
29.     **for**(**int** i=0;i<len1;i++)a\[len1-i-1\]=A\[i\]-'0';  
30.     **for**(**int** i=len1;i<len;i++)a\[i\]=0;  
31.     **for**(**int** i=0;i<len2;i++)b\[len2-i-1\]=B\[i\]-'0';  
32.     **for**(**int** i=len2;i<len;i++)b\[i\]=0;  
33. \}  
34. ///雷德算法，2^M=len，将第i位的数与“i的二进制反转之后的位”的数交换  
35. **void** rader(**long** **long** F\[\],**int** len)  
36. \{  
37.     **int** j=len/2;///模拟二进制反转进位的的位置  
38.     **for**(**int** i=1;i<len-1;i++)  
39.     \{  
40.         **if**(i<j)swap(F\[i\],F\[j\]);///该出手时就出手  
41.         **int** k=len/2;  
42.         **while**(j>=k)  
43.         \{  
44.             j-=k;  
45.             k>>=1;  
46.         \}  
47.         **if**(j<k)j+=k;  
48.     \}  
49. \}  
50. **void** NTT(**long** **long** F\[\],**int** len,**int** t)  
51. \{  
52.     **int** id=0;  
53.     rader(F,len);  
54.     **for**(**int** h=2;h<=len;h<<=1)  
55.     \{  
56.         id++;  
57.         **for**(**int** j=0;j<len;j+=h)  
58.         \{  
59.             **long** **long** E=1;///原根次幂  
60.             **for**(**int** k=j;k<j+h/2;k++)  
61.             \{  
62.                 **long** **long** u=F\[k\];///蝶型操作  
63.                 **long** **long** v=(E\*F\[k+h/2\])%MOD;  
64.                 F\[k\]=(u+v)%MOD;///前半部分  
65.                 F\[k+h/2\]=((u-v)%MOD+MOD)%MOD;///后半部分  
66.                 E=(E\*qp\[id\])%MOD;  
67.             \}  
68.         \}  
69.     \}  
70.     //p(F);  
71.     **if**(t==-1)///插值  
72.     \{  
73.         **for**(**int** i=1;i<len/2;i++)swap(F\[i\],F\[len-i\]);///i+lne-i=i;  
74.         **long** **long** inv=q\_pow(len,MOD-2,MOD);///逆元  
75.         **for**(**int** i=0;i<len;i++)F\[i\]=(F\[i\]%MOD\*inv)%MOD;  
76.     \}  
77.     //p(F);  
78. \}  
79. **void** work()///卷积，点乘，插值  
80. \{  
81.     NTT(a,len,1);  
82.     NTT(b,len,1);  
83.     **for**(**int** i=0;i<len;i++)  
84.         a\[i\]=(a\[i\]\*b\[i\])%MOD;  
85.     NTT(a,len,-1);  
86. \}  
87. **void** pushup()///进位  
88. \{  
89.     **for**(**int** i=0;i<len;i++)  
90.     \{  
91.         **if**(a\[i\]>=10)  
92.         \{  
93.             a\[i+1\]+=a\[i\]/10;  
94.             a\[i\]=a\[i\]%10;  
95.         \}  
96.     \}  
97. \}  
98. **void** print()///输出  
99. \{  
100.     **int** high=0;  
101.     **for**(**int** i=len-1;i>=0;i--)  
102.     \{  
103.         **if**(a\[i\])  
104.         \{  
105.             high=i;  
106.             **break**;  
107.         \}  
108.     \}  
109.     **for**(**int** i=high;i>=0;i--)printf("%lld",a\[i\]);  
110.     puts("");  
111. \}

hdu1402

[点击打开链接][Link 5]

**\[cpp\]**  [view plain][]   [copy][view plain]

[print][view plain] [?][view plain]

1.  \#include <iostream>  
2.  \#include <string.h>  
3.  \#include <algorithm>  
4.  \#include <stdio.h>  
5.  \#include <math.h>  
6.    
7.  **using** **namespace** std;  
8.    
9.  **const** **int** N = 400005;  
10. **const** **double** pi = acos(-1.0);  
11. **const** **int** g=3;  
12. //const int len=1<<18;  
13. **const** **long** **long** MOD=(479<<21)+1;  
14. **int** len;  
15. **char** A\[N\],B\[N\];  
16. **long** **long** a\[N\],b\[N\],qp\[30\];  
17. **void** p(**long** **long** a\[\])  
18. \{  
19.     **for**(**int** i=0;i<len;i++)printf("%lld ",a\[i\]);puts("");  
20. \}  
21. **long** **long** q\_pow(**long** **long** x,**long** **long** y,**long** **long** P)  
22. \{  
23.     **long** **long** ans=1;  
24.     **while**(y>0)  
25.     \{  
26.         **if**(y&1)ans=ans\*x%P;  
27.         x=x\*x%P;  
28.         y>>=1;  
29.     \}  
30.     **return** ans;  
31. \}  
32. **void** init()  
33. \{  
34.     len=1;  
35.     **for**(**int** i=0;i<21;i++)  
36.     \{  
37.         **int** t=1<<i;  
38.         qp\[i\]=q\_pow(g,(MOD-1)/t,MOD);  
39.     \}  
40.     **int** len1=strlen(A);  
41.     **int** len2=strlen(B);  
42.     **while**(len<2\*len1||len<2\*len2)len<<=1;  
43.     **for**(**int** i=0;i<len1;i++)a\[len1-i-1\]=A\[i\]-'0';  
44.     **for**(**int** i=len1;i<len;i++)a\[i\]=0;  
45.     **for**(**int** i=0;i<len2;i++)b\[len2-i-1\]=B\[i\]-'0';  
46.     **for**(**int** i=len2;i<len;i++)b\[i\]=0;  
47. \}  
48. **void** rader(**long** **long** F\[\],**int** len)  
49. \{  
50.     **int** j=len/2;  
51.     **for**(**int** i=1;i<len-1;i++)  
52.     \{  
53.         **if**(i<j)swap(F\[i\],F\[j\]);  
54.         **int** k=len/2;  
55.         **while**(j>=k)  
56.         \{  
57.             j-=k;  
58.             k>>=1;  
59.         \}  
60.         **if**(j<k)j+=k;  
61.     \}  
62. \}  
63. **void** NTT(**long** **long** F\[\],**int** len,**int** t)  
64. \{  
65.     **int** id=0;  
66.     rader(F,len);  
67.     **for**(**int** h=2;h<=len;h<<=1)  
68.     \{  
69.         id++;  
70.         **for**(**int** j=0;j<len;j+=h)  
71.         \{  
72.             **long** **long** E=1;  
73.             **for**(**int** k=j;k<j+h/2;k++)  
74.             \{  
75.                 **long** **long** u=F\[k\];  
76.                 **long** **long** v=(E\*F\[k+h/2\])%MOD;  
77.                 F\[k\]=(u+v)%MOD;  
78.                 F\[k+h/2\]=((u-v)%MOD+MOD)%MOD;  
79.                 E=(E\*qp\[id\])%MOD;  
80.             \}  
81.         \}  
82.     \}  
83.     //p(F);  
84.     **if**(t==-1)  
85.     \{  
86.         **for**(**int** i=1;i<len/2;i++)swap(F\[i\],F\[len-i\]);  
87.         **long** **long** inv=q\_pow(len,MOD-2,MOD);  
88.         **for**(**int** i=0;i<len;i++)F\[i\]=(F\[i\]%MOD\*inv)%MOD;  
89.     \}  
90.     //p(F);  
91. \}  
92. **void** work()  
93. \{  
94.     NTT(a,len,1);  
95.     NTT(b,len,1);  
96.     **for**(**int** i=0;i<len;i++)  
97.         a\[i\]=(a\[i\]\*b\[i\])%MOD;  
98.     NTT(a,len,-1);  
99. \}  
100. **void** pushup()  
101. \{  
102.     **for**(**int** i=0;i<len;i++)  
103.     \{  
104.         **if**(a\[i\]>=10)  
105.         \{  
106.             a\[i+1\]+=a\[i\]/10;  
107.             a\[i\]=a\[i\]%10;  
108.         \}  
109.     \}  
110. \}  
111. **void** print()  
112. \{  
113.     **int** high=0;  
114.     **for**(**int** i=len-1;i>=0;i--)  
115.     \{  
116.         **if**(a\[i\])  
117.         \{  
118.             high=i;  
119.             **break**;  
120.         \}  
121.     \}  
122.     **for**(**int** i=high;i>=0;i--)printf("%lld",a\[i\]);  
123.     puts("");  
124. \}  
125. **int** main()  
126. \{  
127.     **while**(~scanf("%s%s",A,B))  
128.     \{  
129.         init();  
130.         work();  
131.         pushup();  
132.         print();  
133.     \}  
134.     **return** 0;  
135. \}

[sad.gif]: /images/20220716/3b18d1ceb450434fb22ea2ad50c32091.png
[Link 1]: /images/20220716/bfe0bb156bb14d54a63d3ff65dcb026a.png
[ACdream]: http://blog.csdn.net/acdreamers/article/details/39026505
[Link 2]: http://www.csdn123.com/html/itweb/20130725/20149_20120_20137.htm
[Link 3]: http://wenku.baidu.com/view/8bfb0bd476a20029bd642d85.html
[Center]: /images/20220716/820db3eea9af49f09bc8cfd6c590da02.png
[view plain]: http://blog.csdn.net/caoguo_app_android/article/details/44067483#
[Link 4]: http://wenku.baidu.com/link?url=aYf2IZcizK5GWzH8xl0Wcrd9OLNOCJ4Ne7vp4uwcEhOI8ooSXJAmAzb4GrENfUnP0ZwZF3HO85FkNFfIPKWTRQNWAfF92RRZU-nSgfw5DCq
[Link 5]: http://acm.hdu.edu.cn/showproblem.php?pid=1402

FFT-NTT

发表评论取消回复

还没有评论，来说两句吧...

相关阅读