初识FFT和NTT
来源:互联网 发布:网络诈骗举报电话多少 编辑:程序博客网 时间:2024/03/29 05:05
看了很久文档,觉得自己只学会了套模板的能力,理解的代码是怎么写,还有一点原理,看完现在来推一下原理估计又不会了!
学这个的原因是因为codechef的一道题目,可惜现在还是没有解决,谁会了求教 点击打开链接
看了 ACdream 的博客,觉得不够详细,而且看了之后根本看不懂代码里面写的什么鬼,之后找了 一份题解
然后在百度文库找到了一篇讲得很详细的 文档 ,看玩总算理解了那么一点点!
ps:这里我只能写一点代码要做什么,原理不敢亵渎!
1、首先是反转置换,将一个2的n次幂长度的分成两个,奇数是一组,偶数是一组,然后一次递归下去,直到长度为1,这样每个元素的位置会去到他二进制反转之后变成是十进制的位置,做置换可以采用rader算法线性实现!一下引用上面链接博客的一张图:
其实为什么这样想想也很明显,将位置pos编号换成二进制,初始化ans_pos=0,ans依次做位运算右移,低位为0的时候走向二叉树的左儿子,否则是右儿子,而ans_pos依次左移,向二叉树左儿子走的时候末尾补0,否则补1,这样ans_pos的二进制就和刚开始ans的二进制是相反的
rader算法代码如下:
///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换void rader(complex *F,int len){ int j=len/2;///模拟二进制反转进位的的位置 for(int i=1;i<len-1;i++) { if(i<j)swap(F[i],F[j]);///该出手时就出手 int k=len/2; while(j>=k) { j-=k; k>>=1; } if(j<k)j+=k; }}
然后抄了份模板,具体解释就看上面链接ppt的解释,非常详细:
const int N = 500005;const double pi = acos(-1.0);struct complex{///复数类 double real,imag; complex (double r=0.0,double i=0.0) { real=r; imag=i; } complex operator + (const complex &x) { return complex(x.real+real,x.imag+imag); } complex operator - (const complex &x) { return complex(real-x.real,imag-x.imag); } complex operator * (const complex &x) { return complex(x.real*real-x.imag*imag,imag*x.real+real*x.imag); }}vara[N],varb[N];void FFT(complex *F,int len,int t){ rader(F,len); for(int h=2;h<=len;h<<=1)///L from 1 to M,2^M=len { ///J_product_two_power_of_M_subtract_L ///M_subtract_L_product_h_equals_to_len complex wn(cos(t*2*pi/h),sin(t*2*pi/h));///公比 for(int j=0;j<len;j+=h) { complex E(1,0);///螺旋因子 for(int k=j;k<j+h/2;k++) { complex u=F[k];///蝶型操作 complex v=E*F[k+h/2]; F[k]=u+v;///前半部分 F[k+h/2]=u-v;///后半部分 E=E*wn; } } } if(t==-1)///IDFT for(int i=0;i<len;i++) F[i].real/=len;}
NTT也是搞一发模板:
ACdream 点击打开链接
点击打开链接这个一份好的解释看完可以了解一二了
大数运算NTT:
const int N = 400005;const int g=3;const long long MOD=(479<<21)+1;int len;char A[N],B[N];long long a[N],b[N],qp[30];long long q_pow(long long x,long long y,long long P){ long long ans=1; while(y>0) { if(y&1)ans=ans*x%P; x=x*x%P; y>>=1; } return ans;}void init(){ len=1; for(int i=0;i<21;i++) { int t=1<<i; qp[i]=q_pow(g,(MOD-1)/t,MOD); } int len1=strlen(A); int len2=strlen(B); while(len<=2*len1||len<=2*len2)len<<=1; for(int i=0;i<len1;i++)a[len1-i-1]=A[i]-'0'; for(int i=len1;i<len;i++)a[i]=0; for(int i=0;i<len2;i++)b[len2-i-1]=B[i]-'0'; for(int i=len2;i<len;i++)b[i]=0;}///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换void rader(long long F[],int len){ int j=len/2;///模拟二进制反转进位的的位置 for(int i=1;i<len-1;i++) { if(i<j)swap(F[i],F[j]);///该出手时就出手 int k=len/2; while(j>=k) { j-=k; k>>=1; } if(j<k)j+=k; }}void NTT(long long F[],int len,int t){ int id=0; rader(F,len); for(int h=2;h<=len;h<<=1) { id++; for(int j=0;j<len;j+=h) { long long E=1;///原根次幂 for(int k=j;k<j+h/2;k++) { long long u=F[k];///蝶型操作 long long v=(E*F[k+h/2])%MOD; F[k]=(u+v)%MOD;///前半部分 F[k+h/2]=((u-v)%MOD+MOD)%MOD;///后半部分 E=(E*qp[id])%MOD; } } } //p(F); if(t==-1)///插值 { for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);///i+lne-i=i; long long inv=q_pow(len,MOD-2,MOD);///逆元 for(int i=0;i<len;i++)F[i]=(F[i]%MOD*inv)%MOD; } //p(F);}void work()///卷积,点乘,插值{ NTT(a,len,1); NTT(b,len,1); for(int i=0;i<len;i++) a[i]=(a[i]*b[i])%MOD; NTT(a,len,-1);}void pushup()///进位{ for(int i=0;i<len;i++) { if(a[i]>=10) { a[i+1]+=a[i]/10; a[i]=a[i]%10; } }}void print()///输出{ int high=0; for(int i=len-1;i>=0;i--) { if(a[i]) { high=i; break; } } for(int i=high;i>=0;i--)printf("%lld",a[i]); puts("");}
hdu1402
点击打开链接
#include <iostream>#include <string.h>#include <algorithm>#include <stdio.h>#include <math.h>using namespace std;const int N = 400005;const double pi = acos(-1.0);const int g=3;//const int len=1<<18;const long long MOD=(479<<21)+1;int len;char A[N],B[N];long long a[N],b[N],qp[30];void p(long long a[]){ for(int i=0;i<len;i++)printf("%lld ",a[i]);puts("");}long long q_pow(long long x,long long y,long long P){ long long ans=1; while(y>0) { if(y&1)ans=ans*x%P; x=x*x%P; y>>=1; } return ans;}void init(){ len=1; for(int i=0;i<21;i++) { int t=1<<i; qp[i]=q_pow(g,(MOD-1)/t,MOD); } int len1=strlen(A); int len2=strlen(B); while(len<2*len1||len<2*len2)len<<=1; for(int i=0;i<len1;i++)a[len1-i-1]=A[i]-'0'; for(int i=len1;i<len;i++)a[i]=0; for(int i=0;i<len2;i++)b[len2-i-1]=B[i]-'0'; for(int i=len2;i<len;i++)b[i]=0;}void rader(long long F[],int len){ int j=len/2; for(int i=1;i<len-1;i++) { if(i<j)swap(F[i],F[j]); int k=len/2; while(j>=k) { j-=k; k>>=1; } if(j<k)j+=k; }}void NTT(long long F[],int len,int t){ int id=0; rader(F,len); for(int h=2;h<=len;h<<=1) { id++; for(int j=0;j<len;j+=h) { long long E=1; for(int k=j;k<j+h/2;k++) { long long u=F[k]; long long v=(E*F[k+h/2])%MOD; F[k]=(u+v)%MOD; F[k+h/2]=((u-v)%MOD+MOD)%MOD; E=(E*qp[id])%MOD; } } } //p(F); if(t==-1) { for(int i=1;i<len/2;i++)swap(F[i],F[len-i]); long long inv=q_pow(len,MOD-2,MOD); for(int i=0;i<len;i++)F[i]=(F[i]%MOD*inv)%MOD; } //p(F);}void work(){ NTT(a,len,1); NTT(b,len,1); for(int i=0;i<len;i++) a[i]=(a[i]*b[i])%MOD; NTT(a,len,-1);}void pushup(){ for(int i=0;i<len;i++) { if(a[i]>=10) { a[i+1]+=a[i]/10; a[i]=a[i]%10; } }}void print(){ int high=0; for(int i=len-1;i>=0;i--) { if(a[i]) { high=i; break; } } for(int i=high;i>=0;i--)printf("%lld",a[i]); puts("");}int main(){ while(~scanf("%s%s",A,B)) { init(); work(); pushup(); print(); } return 0;}
0 0
- 初识FFT和NTT
- fft & ntt
- FFT-NTT
- fft/ntt
- BZOJ2179【FFT】【NTT】
- [UOJ34]FFT && NTT 模板
- FFT、NTT小结
- FFT,NTT学习笔记
- FFT及NTT模板
- hiho1388 FFT/NTT
- FFT & NTT学习心得
- FFT&NTT(草稿)
- FFT&&FWT&&NTT
- HDU4609 NTT||FFT
- FFT & NTT 学习 模板
- 快速傅里叶变换(FFT)和数论变换(NTT)模板
- FFT+NTT 学习资料收集
- 如何调试fft与ntt
- vim 的基本操作
- 《TCP-IP详解卷1:协议》读书笔记四 ARP地址解析协议 .
- 关于未在本地计算机上注册"Microsoft.Ace.Oledb.12.0"提供程序的一个细节
- 算法导论插入排序算法python实现
- JS笔记整理(五)
- 初识FFT和NTT
- Centos 配置svn 版本管理工具
- hduoj2504(又见GCD)
- 1641 - ASCII Area
- 《TyptScript语言手册》第1章-介绍
- leetcode: Palindrome Partitioning
- java实现消费者与生产者队列
- ubuntu14.10中tomcat8设置管理员帐号
- 类成员函数指针