数学思想和算法

欧拉算法

http://blog.sina.com.cn/s/blog_497ab7660101gd4t.html

在计算固体力学中多用Lagrange 列式，计算流体力学用Euler列式，但在解决流体－固体耦合问题时需要一种将两种方法的优点结合起来的算法，即Arbitrary Lagrange-Euler算法，简称为ALE算法。ALE最早是为了解决流体动力学问题而引入的，并使用有限差分方法。

Donea,Belytschko等人分别将ALE法引入有限元当中，用于求解流体于结构相互作用问题。Hughes等人建立了ALE描述的运动学理论，并使用有限元法解决了粘性不可压缩流体和自由表面流动问题。随着ALE技术的不断完善，一些专业计算软件开始加入ALE功能，LS—Dyna

是目前具有较成熟的ALE算法的大型通用有限元程序，程序中最先采用简化ALE，后来发展到多物质ALE，其应用领域主要是流固耦合方面的计算。

原文地址：欧拉算法作者：我会好好的

欧拉算法

　　微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较

低，一般不在工程中单独进行运算。所谓数值求解，就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。对于常微分方程：

　　dy/dx=f(x,y)，x∈[a,b]

　　y(a)=y0

　　可以将区间[a,b]分成n段，那么方程在第xi点有y'(xi)=f(xi,y(xi))，再用向前差商近似代替导数则为：(y(xi+1)-y(xi))/h= f(xi,y(xi))，在这里，h是步长，即相邻两个结点间的距离。因此可以根据xi点和yi点的数值计算出yi+1来：

　　yi+1= yi+h*f(xi ,yi)，i=0,1,2,L

　　这就是欧拉格式，若初值yi+1是已知的，则可依据上式逐步算出数值解y1,y2,L。

　　为简化分析，人们常在yi为准确即yi=y(xi)的前提下估计误差y(xi+1)-yi+1，这种误差称为局部截断误差。

　　如果一种数值方法的局部截断误差为O(h^p+1)，则称它的精度是p阶的，或称之为p阶方法。欧拉格式的局部截断误差为O(h^2)，由此可知欧拉格式仅为一阶方法。

　　欧拉公式：

　　y(xi+1)=yi+h*f(xi,yi)

　　且xi=x0+i*h (i=0,1,2,…,n-1)

　　局部截断误差是O(h^2)

　　

改进的欧拉算法

　　先用欧拉法求得一个初步的近似值，称为预报值，然后用它替代梯形法右端的yi+1再直接计算fi+1，得到校正值yi+1，这样建立的预报-校正系统称为改进的欧拉格式：

　　预报值 y~i+1=yi+1 + h*f(xi,yi)

　　校正值 yi+1 =yi+(h/2)*[f(xi,yi)+f(xi+1,y~i+1)]

　　它有下列平均化形式：

　　yp=yi+h*f(xi,yi)

　　且 yc=yi+h*f(xi+1,yp)

　　且 yi+1=(xp+yc)/2

　　它的局部截断误差为O(h^3)，可见，改进欧拉格式较欧拉格式提高了精度，其截断误差比欧拉格式提高了一阶。

　　注：欧拉法用差商 [y(xi+1)-y(xi)]/h 近似代替y(xi)的导数，局部截断误差较大；改进欧拉法先用欧拉法求出预报值，再利用梯形公式求出校正值，局部截断误差比欧拉法低了一阶，较大程度地提高了计算精度。

  
改进欧拉算法

#include<iostream.h>
#define N 20
void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n)
{
int i;
float yp,yc,x=x0,y=y0,h=(xn-x0)/n;
cout<<"x[0]="<<x<<'t'<<"y[0]"<<y<<endl;
for(i=1;i<=n;i++)
{
yp=y+h*f1(x,y);
x=x0+i*h;
yc=y+h*f1(x,yp);
y=(yp+yc)/2.0;
cout<<"x["<<i<<"]="<<x<<" y["<<i<<"]="<<y<<endl;
}
}
void main()
{
float xn=5.0,x0=0.0,y0=2.0;
float f1(float ,float);
ModEuler(f1,x0,y0,xn,N);
}
float f1(float x,float y)
{
return -x*y*y;
}

========

高斯算法

首项加末项乘以末项数除以2这样的算法称为高斯算法。

一次数学课上，老师让学生练习算数。于是让他们一个小时内算出1+2+3+4+5+6+……+100的得数。全班只有高斯用了不到20分钟给出了答案，因为他想到了用（1+100）+（2+99）+（3+98）……+（50+51）…………一共有50个101，所以50×101就是1加到一百的得数。后

来人们把这种简便算法称作高斯算法。

具体的方法是：

首项加末项乘以末项数除以2

项数的计算方法是末项减去首项除以项差（每两项之间的差）加1.

1+2+3+4+5+······+n

字母表示：n（1+n）/2

等差数列求和公式　Sn=(a1+an)n/2　Sn=n(2a1+(n-1)d)/2; d=公差　Sn=An2+Bn; A=d/2,B=a1-(d/2)

========

高斯模糊的算法

http://www.ruanyifeng.com/blog/2012/11/gaussian_blur.html

通常，图像处理软件会提供"模糊"（blur）滤镜，使图片产生模糊的效果。

"模糊"的算法有很多种，其中有一种叫做"高斯模糊"（Gaussian Blur）。它将正态分布（又名"高斯分布"）用于图像处理。

本文介绍"高斯模糊"的算法，你会看到这是一个非常简单易懂的算法。本质上，它是一种数据平滑技术（data smoothing），适用于多个场合，图像处理恰好提供了一个直观的应用实例。

一、高斯模糊的原理

所谓"模糊"，可以理解成每一个像素都取周边像素的平均值。

上图中，2是中间点，周边点都是1。

"中间点"取"周围点"的平均值，就会变成1。在数值上，这是一种"平滑化"。在图形上，就相当于产生"模糊"效果，"中间点"失去细节。

显然，计算平均值时，取值范围越大，"模糊效果"越强烈。

上面分别是原图、模糊半径3像素、模糊半径10像素的效果。模糊半径越大，图像就越模糊。从数值角度看，就是数值越平滑。

接下来的问题就是，既然每个点都要取周边像素的平均值，那么应该如何分配权重呢？

如果使用简单平均，显然不是很合理，因为图像都是连续的，越靠近的点关系越密切，越远离的点关系越疏远。因此，加权平均更合理，距离越近的点权重越大，距离越远的点权重越小。

二、正态分布的权重

正态分布显然是一种可取的权重分配模式。

在图形上，正态分布是一种钟形曲线，越接近中心，取值越大，越远离中心，取值越小。

计算平均值的时候，我们只需要将"中心点"作为原点，其他点按照其在正态曲线上的位置，分配权重，就可以得到一个加权平均值。

三、高斯函数

上面的正态分布是一维的，图像都是二维的，所以我们需要二维的正态分布。

正态分布的密度函数叫做"高斯函数"（Gaussian function）。它的一维形式是：

其中，μ是x的均值，σ是x的方差。因为计算平均值的时候，中心点就是原点，所以μ等于0。

根据一维高斯函数，可以推导得到二维高斯函数：

有了这个函数 ，就可以计算每个点的权重了。

四、权重矩阵

假定中心点的坐标是（0,0），那么距离它最近的8个点的坐标如下：

更远的点以此类推。

为了计算权重矩阵，需要设定σ的值。假定σ=1.5，则模糊半径为1的权重矩阵如下：

这9个点的权重总和等于0.4787147，如果只计算这9个点的加权平均，还必须让它们的权重之和等于1，因此上面9个值还要分别除以0.4787147，得到最终的权重矩阵。

五、计算高斯模糊

有了权重矩阵，就可以计算高斯模糊的值了。

假设现有9个像素点，灰度值（0-255）如下：

每个点乘以自己的权重值：

得到

将这9个值加起来，就是中心点的高斯模糊的值。

对所有点重复这个过程，就得到了高斯模糊后的图像。如果原图是彩色图片，可以对RGB三个通道分别做高斯模糊。

六、边界点的处理

如果一个点处于边界，周边没有足够的点，怎么办？

一个变通方法，就是把已有的点拷贝到另一面的对应位置，模拟出完整的矩阵。

七、参考文献

* How to program a Gaussian Blur without using 3rd party libraries

========

十、从头到尾彻底理解傅里叶变换算法

http://www.cnblogs.com/v-July-v/archive/2011/02/20/1983676.html

推荐阅读：The Scientist and Engineer's Guide to Digital Signal Processing，By Steven W. Smith, Ph.D。此书地址：http://www.dspguide.com/pdfbook.htm。

博主说明：I、本文中阐述离散傅里叶变换方法，是根据此书：The Scientist and Engineer's Guide to Digital Signal Processing，By Steven W. Smith, Ph.D.而翻译而成的，此书地址：http://www.dspguide.com/pdfbook.htm。II、同时，有相当一部分内容编辑

整理自dznlong的博客，也贴出其博客地址，向原创的作者表示致敬：http://blog.csdn.net/dznlong 。这年头，真正静下心写来原创文章的人，很少了。

------------------------------------

从头到尾彻底理解傅里叶变换算法、上

前言

第一部分、  DFT

第一章、傅立叶变换的由来

第二章、实数形式离散傅立叶变换（Real DFT）

从头到尾彻底理解傅里叶变换算法、下

第三章、复数

第四章、复数形式离散傅立叶变换

 

前言：

“关于傅立叶变换，无论是书本还是在网上可以很容易找到关于傅立叶变换的描述，但是大都是些故弄玄虚的文章，太过抽象，尽是一些让人看了就望而生畏的公式的罗列，让人很难能够从感性上得到理解”---dznlong，

那么，到底什么是傅里叶变换算法列?傅里叶变换所涉及到的公式具体有多复杂列?

傅里叶变换（Fourier transform）是一种线性的积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。

   哦，傅里叶变换原来就是一种变换而已，只是这种变换是从时间转换为频率的变化。这下，你就知道了，傅里叶就是一种变换，一种什么变换列?就是一种从时间到频率的变化或其相互转化。

ok，咱们再来总体了解下傅里叶变换，让各位对其有个总体大概的印象，也顺便看看傅里叶变换所涉及到的公式，究竟有多复杂：

以下就是傅里叶变换的4种变体（摘自，维基百科）

连续傅里叶变换

   一般情况下，若“傅里叶变换”一词不加任何限定语，则指的是“连续傅里叶变换”。连续傅里叶变换将平方可积的函数f（t）表示成复指数函数的积分或级数形式。

这是将频率域的函数F(ω)表示为时间域的函数f（t）的积分形式。

连续傅里叶变换的逆变换 (inverse Fourier transform)为：

即将时间域的函数f（t）表示为频率域的函数F(ω)的积分。

一般可称函数f（t）为原函数，而称函数F(ω)为傅里叶变换的像函数，原函数和像函数构成一个傅里叶变换对（transform pair）。

除此之外，还有其它型式的变换对，以下两种型式亦常被使用。在通信或是信号处理方面，常以来代换，而形成新的变换对:

 或者是因系数重分配而得到新的变换对：

 一种对连续傅里叶变换的推广称为分数傅里叶变换（Fractional Fourier Transform）。分数傅里叶变换(fractional Fourier transform,FRFT)指的就是傅里叶变换(Fourier transform,FT)的广义化。

分数傅里叶变换的物理意义即做傅里叶变换 a 次，其中 a 不一定要为整数；而做了分数傅里叶变换之后，信号或输入函数便会出现在介于时域(time domain)与频域(frequency domain)之间的分数域(fractional domain)。

当f（t）为偶函数（或奇函数）时，其正弦（或余弦）分量将消亡，而可以称这时的变换为余弦变换（cosine transform）或正弦变换（sine transform）.

另一个值得注意的性质是，当f（t）为纯实函数时，F(−ω) = F*(ω)成立.

傅里叶级数

   连续形式的傅里叶变换其实是傅里叶级数 (Fourier series)的推广，因为积分其实是一种极限形式的求和算子而已。对于周期函数，其傅里叶级数是存在的：

其中Fn为复幅度。对于实值函数，函数的傅里叶级数可以写成：

其中an和bn是实频率分量的幅度。

离散时域傅里叶变换

   离散傅里叶变换是离散时间傅里叶变换（DTFT）的特例（有时作为后者的近似）。DTFT在时域上离散，在频域上则是周期的。DTFT可以被看作是傅里叶级数的逆变换。

离散傅里叶变换

   离散傅里叶变换（DFT），是连续傅里叶变换在时域和频域上都离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序

列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。在实际应用中通常采用快速傅里叶变换以高效计算DFT。

   为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换，必须将函数xn定义在离散点而非连续域内，且须满足有限性或周期性条件。这种情况下，使用离散傅里叶变换（DFT），将函数xn表示为下面的求和形式：

其中Xk是傅里叶幅度。直接使用这个公式计算的计算复杂度为O（n*n），而快速傅里叶变换（FFT）可以将复杂度改进为O（n*lgn）。（后面会具体阐述FFT是如何将复杂度降为O（n*lgn）的。）计算复杂度的降低以及数字电路计算能力的发展使得DFT成为在信号处理领域

十分实用且重要的方法。

   下面，比较下上述傅立叶变换的4种变体，

   如上，容易发现：函数在时（频）域的离散对应于其像函数在频（时）域的周期性。反之连续则意味着在对应域的信号的非周期性。也就是说，时间上的离散性对应着频率上的周期性。同时，注意，离散时间傅里叶变换，时间离散，频率不离散，它在频域依然是连续

的。

   如果，读到此，你不甚明白，大没关系，不必纠结于以上4种变体，继续往下看，你自会豁然开朗。（有什么问题，也恳请提出，或者批评指正）

   ok， 本文，接下来，由傅里叶变换入手，后重点阐述离散傅里叶变换、快速傅里叶算法，到最后彻底实现FFT算法，全篇力求通俗易懂、阅读顺畅，教你从头到尾彻底理解傅里叶变换算法。由于傅里叶变换，也称傅立叶变换，下文所称为傅立叶变换，同一个变换，不

同叫法，读者不必感到奇怪。

第一部分、DFT

第一章、傅立叶变换的由来

    要理解傅立叶变换，先得知道傅立叶变换是怎么变换的，当然，也需要一定的高等数学基础，最基本的是级数变换，其中傅立叶级数变换是傅立叶变换的基础公式。

 

一、傅立叶变换的提出

    傅立叶是一位法国数学家和物理学家，原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier于1807年在法国科学学会上发表了一篇论文，论文里描述运用正弦曲线来描述温度分布，论文里有个在当时具有争议性的决断：任何连续周期信号都可以由一组适当

的正弦曲线组合而成。

    当时审查这个论文拉格朗日坚决反对此论文的发表，而后在近50年的时间里，拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率。直到拉格朗日死后15年这个论文才被发表出来。

    谁是对的呢？拉格朗日是对的：正弦曲线无法组合成一个带有棱角的信号。但是，我们可以用正弦曲线来非常逼近地表示它，逼近到两种表示方法不存在能量差别，基于此，傅立叶是对的。

    为什么我们要用正弦曲线来代替原来的曲线呢？如我们也还可以用方波或三角波来代替呀，分解信号的方法是无穷多的，但分解信号的目的是为了更加简单地处理原来的信号。

    用正余弦来表示原信号会更加简单，因为正余弦拥有原信号所不具有的性质：正弦曲线保真度。一个正余弦曲线信号输入后，输出的仍是正余弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的。且只有正余弦曲线才拥有这样的性质，正因如此

我们才不用方波或三角波来表示。

二、傅立叶变换分类

    根据原信号的不同类型，我们可以把傅立叶变换分为四种类别：

1、非周期性连续信号        傅立叶变换（Fourier Transform） 

2、周期性连续信号           傅立叶级数(Fourier Series) 

3、非周期性离散信号        离散时域傅立叶变换（Discrete Time Fourier Transform） 

4、周期性离散信号           离散傅立叶变换(Discrete Fourier Transform) 

       下图是四种原信号图例（从上到下，依次是FT，FS，DTFT，DFT）： 

    这四种傅立叶变换都是针对正无穷大和负无穷大的信号，即信号的的长度是无穷大的，我们知道这对于计算机处理来说是不可能的，那么有没有针对长度有限的傅立叶变换呢？没有。因为正余弦波被定义成从负无穷小到正无穷大，我们无法把一个长度无限的信号组

合成长度有限的信号。

    面对这种困难，方法是：把长度有限的信号表示成长度无限的信号。如，可以把信号无限地从左右进行延伸，延伸的部分用零来表示，这样，这个信号就可以被看成是非周期性离散信号，我们可以用到离散时域傅立叶变换（DTFT）的方法。也可以把信号用复制的方

法进行延伸，这样信号就变成了周期性离散信号，这时我们就可以用离散傅立叶变换方法（DFT）进行变换。本章我们要讲的是离散信号，对于连续信号我们不作讨论，因为计算机只能处理离散的数值信号，我们的最终目的是运用计算机来处理信号的。

 

    但是对于非周期性的信号，我们需要用无穷多不同频率的正弦曲线来表示，这对于计算机来说是不可能实现的。所以对于离散信号的变换只有离散傅立叶变换（DFT）才能被适用，对于计算机来说只有离散的和有限长度的数据才能被处理，对于其它的变换类型只有在

数学演算中才能用到，在计算机面前我们只能用DFT方法，后面我们要理解的也正是DFT方法。

    这里要理解的是我们使用周期性的信号目的是为了能够用数学方法来解决问题，至于考虑周期性信号是从哪里得到或怎样得到是无意义的。

 

    每种傅立叶变换都分成实数和复数两种方法，对于实数方法是最好理解的，但是复数方法就相对复杂许多了，需要懂得有关复数的理论知识，不过，如果理解了实数离散傅立叶变换(real DFT)，再去理解复数傅立叶变换就更容易了，所以我们先把复数的傅立叶变换

放到一边去，先来理解实数傅立叶变换，在后面我们会先讲讲关于复数的基本理论，然后在理解了实数傅立叶变换的基础上再来理解复数傅立叶变换。

 

    还有，这里我们所要说的变换(transform)虽然是数学意义上的变换，但跟函数变换是不同的，函数变换是符合一一映射准则的，对于离散数字信号处理（DSP），有许多的变换：傅立叶变换、拉普拉斯变换、Z变换、希尔伯特变换、离散余弦变换等，这些都扩展了函

数变换的定义，允许输入和输出有多种的值，简单地说变换就是把一堆的数据变成另一堆的数据的方法。

 

三、一个关于实数离散傅立叶变换(Real DFT)的例子

       先来看一个变换实例，下图是一个原始信号图像：

       这个信号的长度是16，于是可以把这个信号分解9个余弦波和9个正弦波（一个长度为N的信号可以分解成N/2+1个正余弦信号，这是为什么呢？结合下面的18个正余弦图,我想从计算机处理精度上就不难理解，一个长度为N的信号，最多只能有N/2+1个不同频率，再

多的频率就超过了计算机所能所处理的精度范围），如下图：

        9个余弦信号：

        9个正弦信号：

       把以上所有信号相加即可得到原始信号，至于是怎么分别变换出9种不同频率信号的，我们先不急，先看看对于以上的变换结果，在程序中又是该怎么表示的，我们可以看看下面这个示例图：

 

    上图中左边表示时域中的信号，右边是频域信号表示方法，

从左向右，-->，表示正向转换(Forward DFT)，从右向左，<--，表示逆向转换(Inverse DFT)，

用小写x[]表示信号在每个时间点上的幅度值数组, 用大写X[]表示每种频率的副度值数组（即时间x-->频率X）, 

因为有N/2+1种频率，所以该数组长度为N/2+1，

    X[]数组又分两种，一种是表示余弦波的不同频率幅度值：Re X[]，

另一种是表示正弦波的不同频率幅度值：Im X[]，

    Re是实数(Real)的意思，Im是虚数(Imagine)的意思，采用复数的表示方法把正余弦波组合起来进行表示，但这里我们不考虑复数的其它作用，只记住是一种组合方法而已，目的是为了便于表达（在后面我们会知道，复数形式的傅立叶变换长度是N，而不是N/2+1）。

如此，再回过头去，看上面的正余弦各9种频率的变化，相信，问题不大了。

第二章、实数形式离散傅立叶变换（Real DFT）

       上一章，我们看到了一个实数形式离散傅立叶变换的例子，通过这个例子能够让我们先对傅立叶变换有一个较为形象的感性认识，现在就让我们来看看实数形式离散傅立叶变换的正向和逆向是怎么进行变换的。在此，我们先来看一下频率的多种表示方法。

 

一、   频域中关于频率的四种表示方法

 

1、序号表示方法，根据时域中信号的样本数取0 ~ N/2，用这种方法在程序中使用起来可以更直接地取得每种频率的幅度值，因为频率值跟数组的序号是一一对应的: X[k]，取值范围是0 ~ N/2；

2、分数表示方法，根据时域中信号的样本数的比例值取0 ~ 0.5: X[ƒ]，ƒ = k/N，取值范围是0 ~ 1/2；

3、用弧度值来表示，把ƒ乘以一个2π得到一个弧度值，这种表示方法叫做自然频率(natural frequency)：X[ω]，ω = 2πƒ = 2πk/N，取值范围是0 ~ π；

4、以赫兹(Hz)为单位来表示，这个一般是应用于一些特殊应用，如取样率为10 kHz表示每秒有10,000个样本数：取值范围是0到取样率的一半。

 

二、   DFT基本函数

 

ck[i] = cos(2πki/N)

sk[i] = sin(2πki/N)

    其中k表示每个正余弦波的频率，如为2表示在0到N长度中存在两个完整的周期，10即有10个周期，如下图：

       上图中至于每个波的振幅(amplitude)值(Re X[k],Im X[k])是怎么算出来的,这个是DFT的核心，也是最难理解的部分，我们先来看看如何把分解出来的正余弦波合成原始信号(Inverse DFT)。

 

三、   合成运算方法(Real Inverse DFT)

 

DFT合成等式（合成原始时间信号，频率-->时间，逆向变换）：

如果有学过傅立叶级数，对这个等式就会有似曾相识的感觉，不错！这个等式跟傅立叶级数是非常相似的：

           当然，差别是肯定是存在的，因为这两个等式是在两个不同条件下运用的，至于怎么证明DFT合成公式，这个我想需要非常强的高等数学理论知识了，这是研究数学的人的工作，对于普通应用者就不需要如此的追根究底了，但是傅立叶级数是好理解的，我们

起码可以从傅立叶级数公式中看出DFT合成公式的合理性。

                                  _            _

       DFT合成等式中的Im X[k]和Re X[k]跟之前提到的Im X[k]和Re X[k]是不一样的，下面是转换方法（关于此公式的解释，见下文）:

       

       但k等于0和N/2时,实数部分的计算要用下面的等式:

              

       上面四个式中的N是时域中点的总数，k是从0到N/2的序号。

       为什么要这样进行转换呢？这个可以从频谱密度(spectral density)得到理解，如下图就是个频谱图：

       

       这是一个频谱图，横坐标表示频率大小，纵坐标表示振幅大小，原始信号长度为N（这里是32），经DFT转换后得到的17个频率的频谱，频谱密度表示每单位带宽中为多大的振幅，那么带宽是怎么计算出来的呢？看上图，除了头尾两个，其余点的所占的宽度是2/N

，这个宽度便是每个点的带宽，头尾两个点的带宽是1/N,而Im X[k]和Re X[k]表示的是频谱密度，即每一个单位带宽的振幅大小，但表示2/N（或1/N）带宽的振幅大小，所以分别应当是Im X[k]和Re X[k]的2/N（或1/N）。

 

频谱密度就象物理中物质密度，原始信号中的每一个点就象是一个混合物，这个混合物是由不同密度的物质组成的，混合物中含有的每种物质的质量是一样的，除了最大和最小两个密度的物质外，这样我们只要把每种物质的密度加起来就可以得到该混合物的密度了，又

该混合物的质量是单位质量，所以得到的密度值跟该混合物的质量值是一样的。

 

       至于为什么虚数部分是负数，这是为了跟复数DFT保持一致，这个我们将在后面会知道这是数学计算上的需要（Im X[k]在计算时就已经加上了一个负号（稍后，由下文，便可知），再加上负号，结果便是正的，等于没有变化）。

 

       如果已经得到了DFT结果，这时要进行逆转换，即合成原始信号，则可按如下步骤进行转换：

1、先根据上面四个式子计算得出的值；

2、再根据DFT合成等式得到原始信号数据。

下面是用BASIC语言来实现的转换源代码：
‘DFT逆转换方法
‘/XX[]数组存储计算结果（时域中的原始信号）
‘/REX[]数组存储频域中的实数分量，IMX[]为虚分量
‘
DIM XX[511]
DIM REX[256]
DIM IMX[256]
‘
PI = 3.14159265
N% = 512
‘
GOSUB XXXX ‘转到子函数去获取REX[]和IMX[]数据
‘
‘
‘
FOR K% = 0 TO 256
REX[K%] = REX[K%] / (N%/2)
IMX[K%] = -IMX[K%] / (N%/2)
NEXT k%
‘
REX[0] = REX[0] / N
REX[256] = REX[256] / N
‘
‘ 初始化XX[]数组
FOR I% = 0 TO 511
XX[I%] = 0
NEXT I%
‘
‘
‘
‘
‘
FOR K% =0 TO 256
FOR I%=0 TO 511
‘
XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%)
XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%)
‘
NEXT I%
NEXT K%
‘
END

上面代码中420至490换成如下形式也许更好理解，但结果都是一样的：
 FOR I% =0 TO 511

   FOR K%=0 TO 256

 ‘

      XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%) 

      XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%)

 ‘

   NEXT I%

 NEXT K%

 

四、   分解运算方法（DFT）

 

      有三种完全不同的方法进行DFT：一种方法是通过联立方程进行求解, 从代数的角度看，要从N个已知值求N个未知值，需要N个联立方程，且N个联立方程必须是线性独立的，但这是这种方法计算量非常的大且极其复杂，所以很少被采用；第二种方法是利用信号的相

关性（correlation）进行计算，这个是我们后面将要介绍的方法；第三种方法是快速傅立叶变换（FFT），这是一个非常具有创造性和革命性的的方法，因为它大大提高了运算速度，使得傅立叶变换能够在计算机中被广泛应用，但这种算法是根据复数形式的傅立叶变换

来实现的，它把N个点的信号分解成长度为N的频域，这个跟我们现在所进行的实域DFT变换不一样，而且这种方法也较难理解，这里我们先不去理解，等先理解了复数DFT后，再来看一下FFT。有一点很重要，那就是这三种方法所得的变换结果是一样的，经过实践证明，当

频域长度为32时，利用相关性方法进行计算效率最好，否则FFT算法效率较高。现在就让我们来看一下相关性算法。

 

利用第一种方法、信号的相关性(correlation)可以从噪声背景中检测出已知的信号，我们也可以利用这个方法检测信号波中是否含有某个频率的信号波：把一个待检测信号波乘以另一个信号波，得到一个新的信号波，再把这个新的信号波所有的点进行相加，从相加的结

果就可以判断出这两个信号的相似程度。如下图：

        上面a和 b两个图是待检测信号波，图a很明显可以看出是个3个周期的正弦信号波，图b的信号波则看不出是否含有正弦或余弦信号，图c和d都是个3个周期的正弦信号波，图e和f分别是a、b两图跟c、d两图相乘后的结果，图e所有点的平均值是0.5，说明信号a含

有振幅为1的正弦信号c，但图f所有点的平均值是0，则说明信号b不含有信号d。这个就是通过信号相关性来检测是否含有某个信号的方法。

 

       第二种方法：相应地，我也可以通过把输入信号和每一种频率的正余弦信号进行相乘（关联操作），从而得到原始信号与每种频率的关联程度（即总和大小），这个结果便是我们所要的傅立叶变换结果，下面两个等式便是我们所要的计算方法：

       第二个式子中加了个负号，是为了保持复数形式的一致，前面我们知道在计算时又加了个负号，所以这只是个形式的问题，并没有实际意义，你也可以把负号去掉，并在计算时也不加负号。

       这里有一点必须明白一个正交的概念：两个函数相乘，如果结果中的每个点的总和为0，则可认为这两个函数为正交函数。要确保关联性算法是正确的，则必须使得跟原始信号相乘的信号的函数形式是正交的，我们知道所有的正弦或余弦函数是正交的，这一点我

们可以通过简单的高数知识就可以证明它，所以我们可以通过关联的方法把原始信号分离出正余弦信号。当然，其它的正交函数也是存在的，如：方波、三角波等形式的脉冲信号，所以原始信号也可被分解成这些信号，但这只是说可以这样做，却是没有用的。

       下面是实域傅立叶变换的BASIC语言代码：

 

       到此为止，我们对傅立叶变换便有了感性的认识了吧。但要记住，这只是在实域上的离散傅立叶变换，其中虽然也用到了复数的形式，但那只是个替代的形式，并无实际意义，现实中一般使用的是复数形式的离散傅立叶变换，且快速傅立叶变换是根据复数离散傅

立叶变换来设计算法的，在后面我们先来复习一下有关复数的内容，然后再在理解实域离散傅立叶变换的基础上来理解复数形式的离散傅立叶变换。更多见下文：十、从头到尾彻底理解傅里叶变换算法、下（July、dznlong）

  

========

几个快速傅立叶变换算法

http://blog.sina.com.cn/s/blog_7090b86701018ryx.html

离散傅里叶变换(DFT)

DFT的的正变换和反变换分别为（1）和（2）式。假设有N个数据，则计算一个频率点需要N次复数乘法和N-1次复数加法，整个DFT需要N*N次复数乘法和N（N-1）次复数加法；由于一次的复数乘法需要进行4次的实数乘法和2次的复数加法，一次的复数加法需要两次的实数

加法，因此整个DFT需要4*N*N次的实数乘法和2*N（N-1）+2*N*N≈4*N*N次的复数加法。当N比较大时，所需的计算工作量相当大，例如N=1024时大约需要400万次乘法运算，对于实时信号处理来说，将对计算设备提出非常苛刻的要求，于是就提出如何能够减少计算DFT的

运算量的问题。

1965年，库力和图基在《计算数学》杂志上发表《机器计算傅立叶级数的一种算法》，此文是最早提出FFT算法的。此后关于DFT的快速算法称为人们研究的热点课题，也正是FFT算法的出现，才使得数字信号处理能够应用于实时场合并进一步推动数字信号处理的发展和应

用。几个快速傅立叶变换算法 <wbr> <wbr>转

大多数的FFT算法都是利用（3）式的周期性、共轭对称性、可压缩和扩展性的特点来压缩计算量。

1）、根据DFT定义进行计算的代码

直接利用DFT的定义进行计算的算法计算量非常大。

//Data为输入数据指针，Log2N=log2(length),flag=-1时为正变换，flag=1时为反变换，变换结果由指针Data指向的原空间返回
void dft(complex*Data,int Log2N,int flag)
{
int i,j,length;
complex wn;
length=1<<Log2N;
complex*temp=new complex[length];
for(i=0;i
{
temp[i]=0;
for(j=0;j
{
wn=complex(cos(2.0*pi/length*i*j),sin(flag*2.0*pi/length*i*j));
temp[i]+=Data[j]*wn;
}
}
if(flag==1)
for(i=0;i
Data[i]=temp[i]/length;
delete[] temp;
}

2）、倒位序重排
基2、基4和分裂基的DIT、DIF都需要进行倒位序重排。DIT输入为倒位序，输出为正常顺序，DIT输入为正常顺序，输出为倒位序。因此使用DIT要对输入先进行倒位序重排，DIT要对输出进行倒位序重排，才能得到正常顺序的结果。倒序序重排的实现主要是利用加1是从高

位开始加，并且进位是向低位进的特点，这刚好和正位序相反，因此可以得到正位序相应的倒位序。还有一种方法，就是得到各个位上的数码，然后倒转再计入权值，就得到倒位序。
（1）、基2的倒位序重排程序

void reverse2(complex *data,int Log2N)
{
int i,j;
int RevNum;
int MaxPos,CurPos,MaxValue;
complex temp;
MaxValue=(1<<Log2N)-1;
MaxPos=1<<(Log2N-1);
RevNum=0;
for(i=1;i
{
CurPos=MaxPos;
while((CurPos&RevNum)!=0)
{
RevNum=RevNum&(~CurPos);
CurPos=CurPos>>1;
}
RevNum=RevNum|CurPos;
if (RevNum
{
temp=data[RevNum];
data[RevNum]=data[i];
data[i]=temp;
}
}
}

（2）、基4的倒位序重排程序

//data为数据指针,Log4N=log4(length)。

void reverse4(complex *Data,int Log4N)

{

 int i,MaxValue,length,MaxPos,CurPos,RevNum;

 complex temp;

 length=1<<(2*Log4N);

 MaxPos=length/4;

 MaxValue=length-1;

 RevNum=0;

 for(i=1;i

 {

  CurPos=MaxPos;

  while(3*CurPos<(RevNum+1))

  {

   RevNum=RevNum-3*CurPos;

   CurPos=CurPos/4;

  }

  RevNum=RevNum+CurPos;

  if(RevNum

  {

   temp=Data[i];

   Data[i]=Data[RevNum];

   Data[RevNum]=temp;

  }

 }

}

3）、基2时间抽选FFT（DIT2）

把时域的数字信号序列按照奇偶进行分组计算，可以进行如下的变换：

几个快速傅立叶变换算法 <wbr> <wbr>转

从变换结果可以知道，一个长度为N的DFT可以变换成长度为N/2的两个子序列的组合。依次类推，可以直到转为N/2个2点的傅立叶变换的组合。由于经过多次的奇偶抽选，输入数据要按照基2倒序的原则进行重排，输出数据为正常顺序。下面首先用递归的形式进行算法的

描述，由于递归方法没有充分利用DIT2算法的优点---原位计算，因此递归形式只是为了描述的清晰，并不是使用的程序。

void dit2rec(complex*InData,complex*OutData,int length,int flag)

{

   complex*EvenData=new complex(length/2);//偶序列组成一个子序列

   complex*OddData  =new complex(length/2);//奇序列组成一个子序列

   complex*EvenResult=new complex(length/2);//偶序列的变换结果

   complex*OddResult=new complex(length/2);//奇序列的变换结果

   int i,j;

   if(length==1)//单个数的变换，不变，也是递归结束的条件

   {

      if (flag==1)

         OutData[0]=InData[0]/length;

      return;

   }

  for(i=0;i

  {

    EvenData[i]=InData[2*i];

    OddData[i]=InData[2*i+1];

  }

  dit2rec(EvenData,EvenResult,length/2,sign);//对奇偶序列分别进行变换

  dit2rec(OddData,EvenResult,length/2,sign);

  for(i=0;i

  {

    OutData[i]=EvenData+OddData*complex(cos(2*pi*i/length),sin(flag*2*pi*i/length));

    OutData[i+length/2]=EvenData- OddData*complex(cos(2*pi*i/length),sin(flag*2*pi*i/length));

  }

  delete[] EvenData,OddData,EvenResult,OddResult;

  return;

}

非递归实现如下（现不考虑输入的倒序数问题）：

//data为数据指针，返回值相同，Log2N=log2(length）,flag=-1时为正变换，flag=1时为逆变换。

complex* dit2(complex*Data,int Log2N,int flag)

{

 int i,j,k;

 int length=1<<Log2N;

 complex wn,temp;

 int step;

 int index0,index1;

 for(i=1;i<=Log2N;i++)

 {

  step=1<<i;

   for(j=0;j

   {

    wn=complex(cos(2*pi/step*j),sin(flag*2*pi/step*j));

    for(k=0;k

    {

     index0=k*step+j;

     index1=k*step+step/2+j;

     temp=Data[index0];

     Data[index0]=temp+wn*Data[index1];

     Data[index1]=temp-wn*Data[index1];

    }

   }

 }

 if(flag==1)//如果为反变换，要除以序列的长度

  for(int i=0;i

   Data[i]/=length;

 return Data;

}

当i=1时，也就是第一次循环并没有必要进行复数运算，因为j只能取1，wn为实数，这个时间可以节省。因此可以改进为:

void dit2(complex*Data,int Log2N,int flag)

{

   int i,j,k,step,length;

   complex wn,temp,deltawn;

  length=1<<Log2N;

   for(i=0;i

   {

      temp=Data[i];

      Data[i]=Data[i]+Data[i+1];

      Data[i+1]=temp-Data[i+1];

   }

   for(i=2;i<=Log2N;i++)

   {

      wn=1;step=1<<i;deltawn=complex(cos(2.0*pi/step),sin(flag*2.0*pi/step));；

      for(j=0;j

      {        

        for(i=0;i

        {

           temp=Data[i*step+step/2+j]*wn;

           Data[i*step+step/2+j]=Data[i*step+j]-temp;

           Data[i*step+j]=Data[i*step+j]+temp;

         }

         wn=wn*deltawn;

      }

   }

   if(flag==1)

   for(i=0;i

     Data[i]/=length;

}

4）、基2频率抽选FFT（DIF2）

几个快速傅立叶变换算法 <wbr> <wbr>转

//DIF2的递归版本实现：

void dif2rec(complex*InData,complex*OutData,int length,int flag)

{

   complex* LData=new complex(length/2);

   complex* LResult=new complex(length/2);

   complex* RData=new complex(length/2);

   complex* RResult=new complex(length/2);

   complex temp;

   int i;

   if(length==1)

   {

       if(flag==1)

          OutData[0]=InData[0]/length;

       return;

   }

   for(i=0;i

   {

     LData[i]=InData[i];

     RData[i]=InData[i+length/2];

   }

     for(i=0;i

   {

     temp=LData[i];

     LData[i]=LData[i]+RData[i];

     RData[i]=(temp-RData[i])* complex(cos(2*pi*i/length),sin(flag*2*pi*i/length))

   }

     dit2rec(LData,LResult,length/2,sign);

     dit2rec(RData,RResult,length/2,sign);

   for(i=0;i

   {

      OutData[2*i]=LResult[i];

      OutData[2*i+1]=RResult[i];

    }

    return;

 

}

 

非递归实现如下（现不考虑输入的倒序数问题）：

//data为数据指针，返回值相同，r=log2(length）,flag=-1时为正变换，flag=1时为逆变换。

complex*dif2(complex*Data,int Log2N,int flag)

{

 int i,j,k,step,length;

 int index0,index1;

 complex wn,temp;

 length=1<<Log2N;

 for(i=1;i<=Log2N;i++)

 {

  step=1<<(Log2N-i+1);

  for(j=0;j

  {

   wn=complex(cos(2*pi/step*j),sin(2*pi/step*j*flag));

   for(k=0;k

   {

    index0=k*step+j;

    index1=k*step+step/2+j;

    temp=Data[index0];

    Data[index0]=temp+Data[index1];

    Data[index1]=(temp-Data[index1])*wn;

   }

  }

 }

 if(flag==1)

  for(i=0;i

   Data[i]/=length;

 return Data;

}

和DIT一样，最外层的最后一个循环可以另外独立出来，因为最后一个循环没有必要进行复数运算，这样可以减少复数运算的次数。

5）、基4时间抽选FFT（DIT4）

基四时间抽选快速傅立叶算法

几个快速傅立叶变换算法 <wbr> <wbr>转 

几个快速傅立叶变换算法 <wbr> <wbr>转

几个快速傅立叶变换算法 <wbr> <wbr>转

几个快速傅立叶变换算法 <wbr> <wbr>转 

DIT4源程序：data为数据指针，返回值相同，Log4N=log4(length）,flag=-1时为正变换，flag=1时为逆变换。

complex*dit4(complex*Data,int Log4N,int flag)

{

 int i,j,k,length,step,index0,index1,index2,index3;

 complex wn1,wn2,wn3,img,temp0,temp1,temp2,temp3;

 length=1<<(2*Log4N);

 img=complex(0,1);

 /////////////////////////////////////////////////////////

 for(i=1;i<=Log4N;i++)

 {

  step=1<<(2*i);

  for(j=0;j

  {

  

   wn1=complex(cos(2*pi/step*j),sin(flag*2*pi/step*j));

   wn2=complex(cos(4*pi/step*j),sin(flag*4*pi/step*j));

   wn3=complex(cos(6*pi/step*j),sin(flag*6*pi/step*j));

   for(k=0;k

   {

    index0=k*step+j;index1=k*step+step/4+j;index2=k*step+step/2+j;index3=k*step+3*step/4+j;

    temp0=Data[index0]+Data[index1]*wn1+Data[index2]*wn2+Data[index3]*wn3;

    temp1=Data[index0]-Data[index1]*wn1*img-Data[index2]*wn2+Data[index3]*wn3*img;

    temp2=Data[index0]-Data[index1]*wn1+Data[index2]*wn2-Data[index3]*wn3;

    temp3=Data[index0]+Data[index1]*wn1*img-Data[index2]*wn2-Data[index3]*wn3*img;

    Data[index0]=temp0;

    Data[index1]=temp1;

    Data[index2]=temp2;

    Data[index3]=temp3;

   }

  }

 }

 if(flag==1)

 {

  for(i=0;i

  {

   Data[i]/=length;

  }

 }

 return Data;

}

6）、基4频率抽选FFT（DIF4）

几个快速傅立叶变换算法 <wbr> <wbr>转

几个快速傅立叶变换算法 <wbr> <wbr>转

几个快速傅立叶变换算法 <wbr> <wbr>转

几个快速傅立叶变换算法 <wbr> <wbr>转

DIF4源程序：data为数据指针，返回值相同，Log4N=log4(length）,flag=-1时为正变换，flag=1时为逆变换。

complex* dif4(complex*Data,int Log4N,int flag)

{

 int i,j,k,step,length,index0,index1,index2,index3;

 complex wn1,wn2,wn3,img,temp0,temp1,temp2,temp3;

 length=1<<(2*Log4N);

 img=complex(0,1);

 

 for(i=1;i<=Log4N;i++)

 {

  step=1<<2*(Log4N-i+1);

  for(j=0;j

  {

   complex wn1=complex(cos(2*pi/step*j),sin(flag*2*pi/step*j));

   complex wn2=complex(cos(2*pi/step*2*j),sin(flag*2*pi/step*2*j));

   complex wn3=complex(cos(2*pi/step*3*j),sin(flag*2*pi/step*3*j));

   for(k=0;k

   {

    index0=k*step+j,index1=k*step+step/4+j;index2=k*step+step/2+j;index3=k*step+step*3/4+j;

    temp0=Data[index0]+Data[index1]+Data[index2]+Data[index3];

    temp1=(Data[index0]-Data[index1]*img-Data[index2]+Data[index3]*img)*wn1;

    temp2=(Data[index0]-Data[index1]+Data[index2]-Data[index3])*wn2;

    temp3=(Data[index0]+Data[index1]*img-Data[index2]-Data[index3]*img)*wn3;

    Data[index0]=temp0;

    Data[index1]=temp1;

    Data[index2]=temp2;

    Data[index3]=temp3;

   }

  }

 }

if(flag==1)

 for(i=0;i

  Data[i]/=length;

 return Data;

}

 

7）、时间抽选分裂基FFT（srfftdit）

基4FFT运算量比基2FFT运算量少，但如果不满足N=4^m就无法用基4FFT进行运算，又由于基2算法中偶序号部分FFT运算所需乘法次数比技术序号部分乘法次数少，因而对偶序数部分用基2FFT算法，而奇序数部分用基4算法，则运算量可以得到降低。1984年，法国的

P.Dohamel和H.Hollmann把基2分解和基4分解融合在一起，提出“分裂基”FFT算法（split-radix FFT或SRFFT），这种算法在目前已知的N=2^m点FFT算法中乘、加法次数最少，已经非常接近理论上所需乘法次数的最小值，而且它的运算结构与基2FFT相似，没有更多的取

指或存储数据的花销，特别适宜用ASIC实现。分裂基算法与混合基算法的不同之处在于，分裂基算法同一级内既有基2部分又有基4部分，而不像组合数FFT中同一级内只用同一种基运算，它们的相似之处在于可以有时间抽取算法，也可以有频率抽取算法。

将输入数据分成三个子序列：

 几个快速傅立叶变换算法 <wbr> <wbr>转

做DFT运算有

 几个快速傅立叶变换算法 <wbr> <wbr>转

又有上式k范围的限制，还应该补充下列式子才能完整计算N个频率点

 几个快速傅立叶变换算法 <wbr> <wbr>转

接着对N/2点DFT和N/4点DFT按以上方法分解成三部分，如此进行到最后不能分解位置。

  

8）、频率抽选分裂基FFT（srfftdif）

  

void srfftdif(complex*Data,int Log2N,int sign)

{

 int i,j,k,len,step2,step4,ix,id,i0,i1,i2,i3;

 complex temp0,temp1,temp2,temp3;

 complex wn1,wn3,img=complex(0,1);

 len=1<<Log2N;

 for(i=1;i

 {

  step2=1<<(Log2N-i+1);

  step4=step2/4;

  for(j=0;j

  {

   wn1=complex(cos(2.0*pi/step2*j),sin(sign*2.0*pi/step2*j));

   wn3=wn1*wn1*wn1;

   ix=j;

   id=2*step2;

   do

   {

    for(k=ix;k

    {

     i0=k;

     i1=k+step4;

     i2=k+2*step4;

     i3=k+3*step4;

     temp0=Data[i0]+Data[i2];

     temp1=Data[i1]+Data[i3];

     temp2=((Data[i0]-Data[i2])-img*(Data[i1]-Data[i3]))*wn1;

     temp3=((Data[i0]-Data[i2])+img*(Data[i1]-Data[i3]))*wn3;

     Data[i0]=temp0;Data[i1]=temp1;Data[i2]=temp2;Data[i3]=temp3;

    }

    ix=2*id-step2+j;

    id=4*id;

   } while(ix<(len-1));

  }

 }

 ix=0;

 id=4;

 do

 {

  for(i=ix;i

  {

   i0=i;

   i1=i+1;

   temp0=Data[i0]+Data[i1];

   temp1=Data[i0]-Data[i1];

   Data[i0]=temp0;

   Data[i1]=temp1;

  }

  ix=2*id-2;

  id=4*id;

 } while(ix<(len-1));

}

9）、高基FFT

10）、混合基FFT

11）、素因子FFT

12）、线性调频Z变换（chirp-z变换）算法

data为数据指针，Log2N=log2(length)

void czt(complex*data,int Log2N,int numofsample,double initfreq,double finalfreq)

{

 int i,len,totallen,Log2Totallen,lenoffft;

 double deltaomiga=2.0*pi*(finalfreq-initfreq)/(numofsample-1.0);

 len=1<<Log2N;totallen=len+numofsample-1;

 ////////////////////////////////////////////

 i=1;Log2Totallen=0;

 while(i

 {

  i=i*2;

  Log2Totallen++;

 }

 lenoffft=i;

 complex *gn=new complex[lenoffft];

 complex *hn=new complex[lenoffft];

 ///////////////////////////////////////////

 

 for(i=0;i

          gn[i]=complex(cos(i*i/2.0*deltaomiga),sin(-i*i/2.0*deltaomiga))*complex(cos(2.0*pi*initfreq*i),sin(-2.0*pi*initfreq*i))*data[i];

 for(i=len;i

 reverse2(gn,Log2Totallen);

 dit2(gn,Log2Totallen,-1);

 for(i=0;i

         hn[i]=complex(cos(i*i/2.0*deltaomiga),sin(i*i/2.0*deltaomiga));

 for(i=totallen;i

 reverse2(hn,Log2Totallen);

 dit2(hn,Log2Totallen,-1);

 for(i=0;i

  data[i]=gn[i]*hn[i];

 reverse2(data,Log2Totallen);

 dit2(data,Log2Totallen,1);

 

 for(i=0;i

 {

         data[i]=complex(cos(deltaomiga*i*i/2.0),sin(-deltaomiga*i*i/2.0))*data[i];

 }

        delete[] hn;

        delete[] gn;

}

========

算法总结：判断一个数是否为素数

http://blog.csdn.net/arvonzhang/article/details/8564836

1.约定

x%y为x取模y，即x除以y所得的余数，当x<y时，x%y=x，所有取模的运算对象都为整数。

x^y表示x的y次方。乘方运算的优先级高于乘除和取模，加减的优先级最低。

见到x^y/z这样，就先算乘方，再算除法。

A/B，称为A除以B，也称为B除A。

若A%B=0，即称为A可以被B整除，也称B可以整除A。

A*B表示A乘以B或称A乘B，B乘A，B乘以A……都一样。

复习一下小学数学

公因数：两个不同的自然数A和B，若有自然数C可以整除A也可以整除B，那么C就是A和B的公因数。

公倍数：两个不同的自然数A和B，若有自然数C可以被A整除也可以被B整除，那么C就是A和B的公倍数。

互质数：两个不同的自然数，它们只有一个公因数1，则称它们互质。

费马是法国数学家，又译“费尔马”，此人巨牛，他的简介请看下面。不看不知道，一看吓一跳。

费马人物简介：http://baike.baidu.com/view/6303430.htm?fromId=5194&redirected=seachword　　

2.费马小定理：

有N为任意正整数，P为素数，且N不能被P整除（显然N和P互质），则有：N^P%P=N(即：N的P次方除以P的余数是N)。

但是我查了很多资料见到的公式都是这个样子：

(N^(P-1))%P=1后来分析了一下，两个式子其实是一样的，可以互相变形得到。

原式可化为：(N^P-N)%P=0(即：N的P次方减N可以被P整除，因为由费马小定理知道N的P次方除以P的余数是N)把N提出来一个，N^P就成了你N*(N^(P-1))，那么(N^P-N)%P=0可化为：

(N*(N^(P-1)-1))%P=0

请注意上式，含义是：N*(N^(P-1)-1)可以被P整除

又因为N*(N^(P-1)-1)必能整除N（这不费话么!）

所以，N*(N^(P-1)-1)是N和P的公倍数，小学知识了^_^

又因为前提是N与P互质，而互质数的最小公倍数为它们的乘积，所以一定存在

正整数M使得等式成立：N*(N^(P-1)-1)=M*N*P

两边约去N，化简之：N^(P-1)-1=M*P

因为M是整数，显然：N^(P-1)-1)%P=0即：N^(P-1)%P=1

============================================

3.积模分解公式

先有一个引理，如果有：X%Z=0，即X能被Z整除，则有：(X+Y)%Z=Y%Z

设有X、Y和Z三个正整数，则必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z

想了很长时间才证出来，要分情况讨论才行：

1.当X和Y都比Z大时，必有整数A和B使下面的等式成立：

X=Z*I+A（1）

Y=Z*J+B（2）

不用多说了吧，这是除模运算的性质！

将（1）和（2）代入(X*Y)modZ得：((Z*I+A)(Z*J+B))%Z乘开，再把前三项的Z提一个出来，变形为：(Z*(Z*I*J+I*A+I*B)+A*B)%Z（3）

因为Z*(Z*I*J+I*A+I*B)是Z的整数倍……晕，又来了。

概据引理，（3）式可化简为：(A*B)%Z又因为：A=X%Z，B=Y%Z，代入上面的式子，就成了原式了。

2.当X比Z大而Y比Z小时，一样的转化：

X=Z*I+A

代入(X*Y)%Z得：

(Z*I*Y+A*Y)%Z

根据引理，转化得：(A*Y)%Z

因为A=X%Z，又因为Y=Y%Z，代入上式，即得到原式。

同理，当X比Z小而Y比Z大时，原式也成立。

3.当X比Z小，且Y也比Z小时，X=X%Z，Y=Y%Z，所以原式成立。

=====================================================

4.快速计算乘方的算法

如计算2^13，则传统做法需要进行12次乘法。

[cpp] view plain copy

/*计算n^p*/  

unsigned power(unsigned n,unsigned p)  

{  

    for(int i=0;i<p;i++) n*=n;  

    return n;  

}  

该死的乘法，是时候优化一下了！把2*2的结果保存起来看看，是不是成了：

4*4*4*4*4*4*2

再把4*4的结果保存起来：16*16*16*2 

一共5次运算，分别是2*2、4*4和16*16*16*2

这样分析，我们算法因该是只需要计算一半都不到的乘法了。

为了讲清这个算法，再举一个例子2^7：2*2*2*2*2*2*2 

两两分开：(2*2)*(2*2)*(2*2)*2 

如果用2*2来计算，那么指数就可以除以2了，不过剩了一个，稍后再单独乘上它。

再次两两分开，指数除以2： ((2*2)*(2*2))*(2*2)*2

实际上最后一个括号里的2 * 2是这回又剩下的，那么，稍后再单独乘上它 现在指数已经为1了，可以计算最终结果了：16*4*2=128

优化后的算法如下：

[cpp]

unsigned Power(unsigned n,unsigned p)   

{  

   unsigned main=n; //用main保存结果  

   unsigned odd=1; //odd用来计算“剩下的”乘积  

   while (p>1)   

   {//一直计算，直到指数小于或等于1  

        if((p%2)!=0)  

        {// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它  

乘到结果中  

            odd*=main; //把“剩下的”乘起来  

        }  

        main*=main; //主体乘方  

        p/=2; //指数除以2  

   }  

   return main*odd; //最后把主体和“剩下的”乘起来作为结果  

}  

够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0!

完美版：

[cpp] view plain copy

unsigned Power(unsigned n, unsigned p)   

{ // 计算n的p次方  

    unsigned odd = 1; //oddk用来计算“剩下的”乘积  

    while (p > 1)  

    { // 一直计算到指数小于或等于1  

       if (( p & 1 )!=0)  

      { // 判断p是否奇数，偶数的最低位必为0  

             odd *= n; // 若odd为奇数，则把“剩下的”乘起来  

      }  

      n *= n; // 主体乘方  

      p /= 2; // 指数除以2  

     }  

    return n * odd; // 最后把主体和“剩下的”乘起来作为结果  

}  

========================================================

5."蒙格马利”快速幂模算法

后面我们会用到这样一种运算：(X^Y)%Z。但问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？

考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。

下面的讲解是给尚没有做出这样动作的同学们准备的：

X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……

如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法：

[cpp] view plain copy

unsigned Montgomery(unsigned n, unsigned p, unsigned m)  

{ // 快速计算 (n ^ e) % m 的值，与power算法极类似  

    unsigned r = n % m; // 这里的r可不能省  

    unsigned k = 1;  

    while (p > 1)  

    {  

        if ((p & 1)!=0)  

        {  

            k = (k * r) % m; // 直接取模  

        }  

        r = (r * r) % m; // 同上  

        p /= 2;  

    }  

    return (r * k) % m; // 还是同上  

}  

上面的代码还可以优化。下面是蒙格马利极速版：

[cpp]

unsigned Montgomery(unsigned n,unsigned p,unsigned m)  

{ //快速计算(n^p)%m的值  

      unsignedk=1;  

      n%=m;  

     while(p!=1)  

     {  

         if(0!=(p&1))k=(k*n)%m;  

         n=(n*n)%m;  

         p>>=1;  

    }  

    return(n*k)%m;  

}  

=====================================================

6.怎么判断一个数是否为素数？

1)笨蛋的作法：

[cpp]

bool IsPrime(unsigned n)  

{  

    if (n<2)  

    {     

        //小于2的数即不是合数也不是素数  

        throw 0;  

    }  

    for (unsigned i=2;i<n;++i)  

    {   

        //和比它小的所有的数相除，如果都除不尽，证明素数  

        if (n%i==0)  

        {  

            //除尽了，则是合数  

            return false;  

        }  

    }  

    return true;  

}  

一个数去除以比它的一半还要大的数，一定除不尽，所以还用判断吗？？

2)下面是小学生的做法：

[cpp]

bool IsPrime(unsigned n)  

{  

    if (n<2)  

    {  

        //小于2的数即不是合数也不是素数  

        throw 0;  

    }  

    for(unsigned i=2;i<n/2+1;++i)  

    {   

        // 和比它的一半小数相除，如果都除不尽，证明素数  

        if ( 0 == n % i )  

        {   

            // 除尽了，合数  

            return false;  

        }  

    }  

    return true; // 都没除尽，素数  

}  

一个合数必然可以由两个或多个质数相乘而得到。那么如果一个数不能被比它的一半小的所有的质数整除，那么比它一半小的所有的合数也一样不可能整除它。建立一个素数表是很有用的。

3)下面是中学生的做法：

[cpp]

bool IsPrime2(unsigned n)  

{  

    if ( n < 2 )  

    { // 小于2的数即不是合数也不是素数  

        throw 0;  

    }  

    static unsigned aPrimeList[] = { // 素数表  

        1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,  

        43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 113,   

        193, 241, 257, 337, 353, 401, 433, 449, 577, 593, 641,   

        673, 769, 881, 929, 977, 1009, 1153, 1201, 1217, 1249,   

        1297,1361, 1409, 1489, 1553, 1601, 1697, 1777, 1873,   

        1889, 2017, 2081, 2113, 2129, 2161, 2273, 2417, 2593,   

        2609, 2657, 2689, 2753, 2801, 2833, 2897, 3041, 3089,   

        3121, 3137, 3169, 3217, 3313, 3329, 3361, 3457, 3617,   

        3697, 3761, 3793, 3889, 4001, 4049, 4129, 4177, 4241,   

        4273, 4289, 4337, 4481, 4513, 4561, 4657, 4673, 4721,   

        4801, 4817, 4993, 5009, 5153, 5233, 5281, 5297, 5393,   

        5441, 5521, 5569, 5857, 5953, 6113, 6257, 6337, 6353,   

        6449, 6481, 6529, 6577, 6673, 6689, 6737, 6833, 6961,   

        6977, 7057, 7121, 7297, 7393, 7457, 7489, 7537, 7649,   

        7681, 7793, 7841, 7873, 7937, 8017, 8081, 8161, 8209,   

        8273, 8353, 8369, 8513, 8609, 8641, 8689, 8737, 8753,   

        8849, 8929, 9041, 9137, 9281, 9377, 9473, 9521, 9601,   

        9649, 9697, 9857   

    };  

      

    const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//计算素数表里元素的个数  

    for (unsigned i=2;i<nListNum;++i )  

    {   

        if(n/2+1<aPrimeList[i])  

        {  

            return true;  

        }  

        if(0==n%aPrimeList[i])  

        {  

            return false;  

        }  

    }  

    /*由于素数表中元素个数是有限的，那么对于用素数表判断不到的数，就只有用笨蛋办法了*/  

    for (unsigned i=aPrimeList[nListNum-1];i<n/2+1;i++ )  

    {   

        if (0==n%i)  

        {   

            // 除尽了，合数   

            return false;  

        }  

    }  

    return true;   

}  

        还是太糟了，我们现在要做的对于大型素数的判断，那个素数表倒顶个P用！当然，我们可以利用动态的素数表来进行优化，这就是大学生的做法了。但是动态生成素数表的策略又复杂又没有效率，所以我们还是直接跳跃到专家的做法吧：

        根据上面讲到的费马小定理，对于两个互质的素数N和P，必有：N^(P-1)%P=1 ,那么我们通过这个性质来判断素数吧，当然，你会担心当P很大的时候乘方会很麻烦。不用担心！我们上面不是有个快速的幂模算法么？好好的利用蒙格马利这位大数学家为我们带来

的快乐吧！

算法思路是这样的： 

        对于N，从素数表中取出任意的素数对其进行费马测试，如果取了很多个素数，N仍未测试失败，那么则认为N是素数。当然，测试次数越多越准确，但一般来讲50次就足够了。另外，预先用“小学生”的算法构造一个包括500个素数的数组，先对Q进行整除测试，

将会大大提高通过率，方法如下:

6)下面是专家的做法：

[cpp]

bool IsPrime3(unsigned n)  

{  

    if ( n < 2 )  

    {   

        // 小于2的数即不是合数也不是素数  

        throw 0;  

    }  

    static unsigned aPrimeList[] = {  

        2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,  

        43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97  

    };  

    const int nListNum = sizeof(aPrimeList) / sizeof(unsigned);  

    for (int i=0;i<nListNum;++i)  

    {   

        // 按照素数表中的数对当前素数进行判断  

        if (1!=Montgomery(aPrimeList[i],n-1,n)) // 蒙格马利算法  

        {  

            return false;  

        }  

    }  

    return true;  

}  

        OK，这就专家的作法了。

        等等，什么？好像有点怪，看一下这个数29341，它等于13 * 37 * 61，显然是一个合数，但是竟通过了测试！！哦，抱歉，我忘了在素数表中加入13，37，61这三个数，我其实是故意的，我只是想说明并费马测试并不完全可靠。

        现在我们发现了重要的一点，费马定理是素数的必要条件而非充分条件。这种不是素数，但又能通过费马测试的数字还有不少，数学上把它们称为卡尔麦克数，现在数学家们已经找到所有10 ^ 16以内的卡尔麦克数，最大的一个是9585921133193329。我们必须寻

找更为有效的测试方法。数学家们通过对费马小定理的研究，并加以扩展，总结出了多种快速有效的素数测试方法，目前最快的算法是拉宾米勒测试算法，下面介绍拉宾米勒测试。

================================================================

7.拉宾米勒测试

        拉宾米勒测试是一个不确定的算法，只能从概率意义上判定一个数可能是素数，但并不能确保。算法流程如下:

       1.选择T个随机数A，并且有A<N成立。

       2.找到R和M，使得N=2*R*M+1成立。

       快速得到R和M的方式：N用二进制数B来表示，令C=B-1。因为N为奇数（素数都是奇数），所以C的最低位为0，从C的最低位的0开始向高位统计，一直到遇到第一个1。这时0的个数即为R，M为B右移R位的值。

       3.如果A^M%N=1，则通过A对于N的测试，然后进行下一个A的测试

       4.如果A^M%N!=1，那么令i由0迭代至R，进行下面的测试

       5.如果A^((2^i)*M)%N=N-1则通过A对于N的测试，否则进行下一个i的测试 

       6.如果i=r，且尚未通过测试，则此A对于N的测试失败，说明N为合数。

       7.进行下一个A对N的测试，直到测试完指定个数的A

       通过验证得知，当T为素数，并且A是平均分布的随机数，那么测试有效率为1 / ( 4 ^ T )。如果T > 8那么测试失误的机率就会小于10^(-5)，这对于一般的应用是足够了。如果需要求的素数极大，或着要求更高的保障度，可以适当调高T的值。

下面是代码：

[cpp]

bool RabbinMillerTest( unsigned n )  

{  

    if (n<2)  

    {   

         // 小于2的数即不是合数也不是素数  

        throw 0;  

    }  

    const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数  

    for(int i=0;i<nPrimeListSize;++i)  

    {  

        // 按照素数表中的数对当前素数进行判断  

        if (n/2+1<=g_aPrimeList[i])  

        {  

            // 如果已经小于当前素数表的数，则一定是素数  

            return true;  

        }  

        if (0==n%g_aPrimeList[i])  

        {  

            // 余数为0则说明一定不是素数  

            return false;  

        }  

    }  

    // 找到r和m，使得n = 2^r * m + 1;  

    int r = 0, m = n - 1; // ( n - 1 ) 一定是合数  

    while ( 0 == ( m & 1 ) )  

    {  

        m >>= 1; // 右移一位  

        r++; // 统计右移的次数  

    }  

    const unsigned nTestCnt = 8; // 表示进行测试的次数  

    for ( unsigned i = 0; i < nTestCnt; ++i )  

    {   

        // 利用随机数进行测试，  

        int a = g_aPrimeList[ rand() % nPrimeListSize ];  

        if ( 1 != Montgomery( a, m, n ) )  

        {  

            int j = 0;  

            int e = 1;  

            for ( ; j < r; ++j )  

            {  

                if ( n - 1 == Montgomery( a, m * e, n ) )   

                {  

                    break;  

                }  

                e <<= 1;  

            }  

            if (j == r)  

            {  

                return false;  

            }  

        }  

    }  

    return true;  

}  

========