程序博客网 > 淘宝网店工作室怎么开

FFT-NTT

来源：互联网 发布：淘宝网店工作室怎么开 编辑：程序博客网 时间：2024/04/26 08:47

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

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

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

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

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

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


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

rader算法代码如下：

[cpp] view plain copy
 print?
  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
 print?
  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 点击打开链接

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


大数运算NTT：

[cpp] view plain copy
 print?
  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

点击打开链接

[cpp] view plain copy
 print?
  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. }  
0 0
淘宝网店工作室怎么开 淘宝网店工作室怎么开
原创粉丝点击
热门IT博客
鸟哥的linux私房菜同类鸟哥的linux私房菜同类
腾讯视频会员破解软件腾讯视频会员破解软件
视频融合软件
linux安装java
cs算法
百度云签到源码
矢量地图数据下载
淘宝防蓝光眼镜有用吗
中国网络视听协会
手机打印软件免费版
淘宝管控记录得扣分吗
nginx 多站点
淘宝首页怎么建热点
安徽联通云计算
韦德历届总决赛数据
淘宝交易快照怎么查看
windows系统ios调试
mac 迅雷
知天下图
eve美女捏脸数据
热门问题 老师的惩罚 人脸识别 我在镇武司摸鱼那些年 重生之率土为王 我在大康的咸鱼生活 盘龙之生命进化 天生仙种 凡人之先天五行 春回大明朝 姑娘不必设防，我是瞎子 胸罩的价格 蓝色胸罩 脱胸罩 胸罩图片 如何选胸罩 明星胸罩 胸罩大 胸罩怎么买 不穿胸罩 红色胸罩 怎么穿胸罩 欣姿芳内衣 柔美内衣 性感蕾丝胸罩 哪个品牌的胸罩好 如何穿胸罩 胸罩有多少种 胸罩怎么用 黛安芬女聚拢文胸 无痕硅胶隐形文胸 夏娃之秀硅胶文胸 废柴兄弟1兰菲没戴文胸 废柴兄弟兰菲没戴文胸 木村纱织没有戴文胸 小胸适合穿哪种文胸 平胸适合穿什么文胸 有钢圈和无钢圈文胸的区别 文能 中国文脉 文脉 文脉教育 中国文脉摘抄及赏析 中国文脉 余秋雨 文臣 文艺时代 文艺小清新 文艺句子 文艺大明星 文艺 文艺生活 文艺青年

程序博客网，程序员的互联网技术博客家园。csdn论坛精品 msdn技术资料都在这里