我们熟知的FFT算法实际上是将一个多项式在2n个单位根处展开,将其点值对应相乘,并进行逆变换。然而,由于单位根具有“旋转”的特征(即$w_{m}^{j}=w_{m}^{j+m}$),若多项式次数大于二分之长度,FFT将进行一次长度为2n的循环卷积。bluestein的算法是为了解决在任意长度上的循环卷积问题。

我们知道,任何一个n次多项式都可以被n+1个点值进行表示,因此如果我们选取所有形如$w_{n+1}^{i}$的单位根并带入多项式,进行类似于FFT的变化(这里没有证明),理应得到正确的结果。

设多项式A为$\sum_{i=0}^{n}{a_i*x^i}$,$F_k$为$A(w_{n+1}^{k})$,则有:

$F_k=\sum_{i=0}^{n}{a_i*w_{n+1}^{ik}}$

考虑ik的另外一种组合含义,即有两个盒子,每个盒子分别有i个球和k个球,求有多少种随机拿出两个球且分别属于两个盒子的方法,因此$ik=\tbinom{i+k}{2}-\tbinom{i}{2}-\tbinom{k}{2}$。它的意义在下面推导中可见。

因此$F_k=\sum_{i=0}^{n}{a_i*w_{n+1}^{\tbinom{i+k}{2}-\tbinom{i}{2}-\tbinom{k}{2}}}$

$=w_{n+1}^{-\tbinom{k}{2}}\sum_{i=0}^{n}{a_i*w_{n+1}^{-\tbinom{i}{2}}*w_{n+1}^{\tbinom{i+k}{2}}}$

注意到(i+k)-(i)=k,令$A_{-i}=a_i*w_{n+1}^{-\tbinom{i}{2}}$,$B_i=w_{n+1}^{\tbinom{i}{2}}$。因此,A和B的卷积的第k项即为$F_k$。由于A的下标为负数,我们将A的下标集体加上n。于是,一次bluestein操作花了三次长度为4n的FFT操作。

将多项式转化为点值表达后,我们依葫芦画瓢地将对应位置相乘、进行相应的逆变换(即取单位根的共轭)。而此部分正确性的证明过程是与FFT类似的。

例题:poj2821

 1 // 2821
2 #include<cstdio>
3 #include<math.h>
4 #include<cstring>
5 #include<iomanip>
6 #define mod 998244353
7 using namespace std;
8 typedef double ld;
9 const int maxn=(1<<19)+5;
10 const int LIMIT=1<<19;
11 const ld pi=acos(-1);
12 struct com
13 {
14 ld x,y;
15 com(ld a=0,ld b=0):x(a),y(b){}
16 com operator+(const com&A){return com(x+A.x,y+A.y);}
17 com operator-(const com&A){return com(x-A.x,y-A.y);}
18 com operator*(const com&A){return com(x*A.x-y*A.y,x*A.y+y*A.x);}
19 com operator/(const ld&d){return com(x/d,y/d);}
20 com operator/(const com&A){return com(x,y)*com(A.x,-A.y)/(A.x*A.x+A.y*A.y);}
21 void operator/=(const ld&d){x/=d,y/=d;}
22 };
23 int r[maxn];
24 inline void DFT(com*A,int limit,int type)
25 {
26 for(int i=1;i<limit;++i)
27 {
28 r[i]=(r[i>>1]>>1)|((i&1)?(limit>>1):0);
29 if(i<r[i])
30 swap(A[i],A[r[i]]);
31 }
32 for(int len=2;len<=limit;len<<=1)
33 {
34 com w;
35 if(type==1)
36 w=com(cos(pi*2/len),sin(pi*2/len));
37 else
38 w=com(cos(pi*2/len),-sin(pi*2/len));
39 for(int i=0;i<limit;i+=len)
40 {
41 com d(1,0);
42 for(int j=0,p1=i,p2=i+len/2;j<len/2;++j,++p1,++p2)
43 {
44 com a=A[p1],b=A[p2]*d;
45 A[p1]=a+b;
46 A[p2]=a-b;
47 d=d*w;
48 }
49 }
50 }
51 }
52 com tmp1[maxn],tmp2[maxn];
53
54 inline void bluestein(com*A,int n,int type) // n already stands for the number of terms
55 {
56 int limit=1;
57 while(limit<4*n) // 4 times !!!!!!!
58 limit<<=1;
59 for(int i=0;i<limit;++i)
60 tmp1[i]=tmp2[i]=0;
61 for(int i=0;i<n;++i)
62 tmp1[i]=A[i]*com(cos(pi*i*i/n),type*sin(pi*i*i/n));
63 for(int i=0;i<n*2;++i)
64 tmp2[i]=com(cos(pi*(i-n)*(i-n)/n),-type*sin(pi*(i-n)*(i-n)/n));
65 DFT(tmp1,limit,1);
66 DFT(tmp2,limit,1);
67 for(int i=0;i<limit;++i)
68 tmp1[i]=tmp1[i]*tmp2[i];
69 DFT(tmp1,limit,-1);
70 for(int i=0;i<n;++i)
71 A[i]=tmp1[i+n]*com(cos(pi*i*i/n),type*sin(pi*i*i/n))/limit; // dont forget this !!!
72 }
73 com A[maxn],B[maxn],C[maxn];
74 int n;
75 int main()
76 {
77 scanf("%d",&n);
78 --n;
79 for(int i=0;i<=n;++i)
80 scanf("%lf",&A[i].x);
81 for(int i=0;i<=n;++i)
82 scanf("%lf",&B[i].x);
83 bluestein(A,n+1,1);
84 bluestein(B,n+1,1);
85 for(int i=0;i<n+1;++i)
86 A[i]=B[i]/A[i];
87 bluestein(A,n+1,-1);
88 for(int i=0;i<=n;++i)
89 A[i].x/=(n+1);
90 for(int i=0;i<=n;++i)
91 printf("%.4f\n",A[i].x);
92 return 0;
93 }

bluestein算法的更多相关文章

  1. 算法系列:FFT 002

    转载自http://blog.jobbole.com/58246/ 快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一.没有正规计算机科学课程背景 ...

  2. 快速傅里叶变换(FFT)算法【详解】

    快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一.我打开一本老旧的算法书,欣赏了JW Cooley 和 John Tukey 在1965年的文章 ...

  3. B树——算法导论(25)

    B树 1. 简介 在之前我们学习了红黑树,今天再学习一种树--B树.它与红黑树有许多类似的地方,比如都是平衡搜索树,但它们在功能和结构上却有较大的差别. 从功能上看,B树是为磁盘或其他存储设备设计的, ...

  4. 分布式系列文章——Paxos算法原理与推导

    Paxos算法在分布式领域具有非常重要的地位.但是Paxos算法有两个比较明显的缺点:1.难以理解 2.工程实现更难. 网上有很多讲解Paxos算法的文章,但是质量参差不齐.看了很多关于Paxos的资 ...

  5. 【Machine Learning】KNN算法虹膜图片识别

    K-近邻算法虹膜图片识别实战 作者:白宁超 2017年1月3日18:26:33 摘要:随着机器学习和深度学习的热潮,各种图书层出不穷.然而多数是基础理论知识介绍,缺乏实现的深入理解.本系列文章是作者结 ...

  6. 红黑树&mdash;&mdash;算法导论(15)

    1. 什么是红黑树 (1) 简介     上一篇我们介绍了基本动态集合操作时间复杂度均为O(h)的二叉搜索树.但遗憾的是,只有当二叉搜索树高度较低时,这些集合操作才会较快:即当树的高度较高(甚至一种极 ...

  7. 散列表(hash table)——算法导论(13)

    1. 引言 许多应用都需要动态集合结构,它至少需要支持Insert,search和delete字典操作.散列表(hash table)是实现字典操作的一种有效的数据结构. 2. 直接寻址表 在介绍散列 ...

  8. 虚拟dom与diff算法 分析

    好文集合: 深入浅出React(四):虚拟DOM Diff算法解析 全面理解虚拟DOM,实现虚拟DOM

  9. 简单有效的kmp算法

    以前看过kmp算法,当时接触后总感觉好深奥啊,抱着数据结构的数啃了一中午,最终才大致看懂,后来提起kmp也只剩下“奥,它是做模式匹配的”这点干货.最近有空,翻出来算法导论看看,原来就是这么简单(先不说 ...

  10. 神经网络、logistic回归等分类算法简单实现

    最近在github上看到一个很有趣的项目,通过文本训练可以让计算机写出特定风格的文章,有人就专门写了一个小项目生成汪峰风格的歌词.看完后有一些自己的小想法,也想做一个玩儿一玩儿.用到的原理是深度学习里 ...

随机推荐

  1. VisualVM通过jstatd方式远程监控远程主机

    配置好权限文件 [root@test bin]# cd $JAVA_HOME/bin [root@test bin]# vim jstatd.all.policy grant codebase &qu ...

  2. OC中的字典

    // ********************不可变最字典***************** /* NSDictionary * dic = [NSDictionary dictionaryWithO ...

  3. js传递json格式对象到服务器

    var message = new Object();message.event = "test";message.params = new Object();message.pa ...

  4. Ceph剖析:数据分布之CRUSH算法与一致性Hash

    作者:吴香伟 发表于 2014/09/05 版权声明:可以任意转载,转载时务必以超链接形式标明文章原始出处和作者信息以及版权声明 数据分布是分布式存储系统的一个重要部分,数据分布算法至少要考虑以下三个 ...

  5. poj3295Tautology

    http://poj.org/problem?id=3295 这几天补一补poj之前落下的题吧 枚举 #include <iostream> #include<cstdio> ...

  6. oracle 过程函数,包的区别和联系

    一.过程与函数区别 1.过程可以有0~N个返回参数,通过OUT or IN OUT参数返回:函数有且仅有1个返回值,通过return语句返回. 2.调用过程时,可做为单独的语句执行:调用函数时,函数必 ...

  7. 【Beta】阶段 第三次Daily Scrum Meeting

    每日任务 ·1.本次会议为第三次 Meeting 会议 ·2.本次会议在周三上午9:40召开,会议时间为10分钟 一.今日站立式会议照片 二.每个人的工作(有work item的ID) 三.工作中遇到 ...

  8. linux下使用crontab实现定时PHP计划任务失败的原因分析

    这篇文章主要介绍了linux下使用crontab实现定时PHP计划任务失败的原因分析,需要的朋友可以参考下   很多人在linux下使用crontab实现PHP执行定时任务却未能成功,不能生成缓存.本 ...

  9. Java实现AES加密,异常java.security.InvalidKeyException: Illegal key size 的解决

    Java实现AES加密,抛出异常如下:java.security.InvalidKeyException: Illegal key size 代码参考 http://my.oschina.net/Ja ...

  10. mybatis源码- 反射模块一(跟着MyBatis学反射):类级别信息的封装

    目录 1 JavaBean 规范 2 Reflector和ReflectorFactory 2.1 Reflector 属性 2.1.1 属性 2.1.2 Invoker 接口 2.2 Reflect ...