【算法】数据结构与算法分析学习笔记——第三章习题选做快速傅里叶变换与多项式乘法

3.7编写一个函数将两个多项式相乘，用链表实现。必须保证输出的多项式按幂排列并且最多有一项为任意幂。

c.编写一个以O(MNlog(MN))时间执行乘法的程序。

这里介绍一种用快速傅里叶变化实现的方法，这篇真的写得很好极力推荐。可以参考以下链接：

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

现在讲讲是怎么用C语言去实现的

0.复数计算

其实库函数里面应该是有的，不过这么简单的功能，造造轮子比翻函数库要方便许多

// r = c1 + c2
static inline void complex_add(float *r, float *c1, float *c2)
{
    r[0] = c1[0] + c2[0];
    r[1] = c1[1] + c2[1];
}

// r = c1 - c2
static inline void complex_sub(float *r, float *c1, float *c2)
{
    r[0] = c1[0] - c2[0];
    r[1] = c1[1] - c2[1];
}

// r = c1 * c2
static inline void complex_mul(float *r, float *c1, float *c2)
{
    r[0] = c1[0] * c2[0] - c1[1] * c2[1];
    r[1] = c1[1] * c2[0] + c1[0] * c2[1];
}

1.离散傅里叶变换

对于一个自动化专业的学生，这东西还是见过的

参照上面推荐博客讲到的傅里叶变换，在多项式傅里叶变化这个例子中，我们做一个流程图方便理解

所以很容易能得出，这种计算方法的时间为O(N^2)，这是不符合题目要求的

对了，这里说一下为什么要用傅里叶变化是可行的，因为傅里叶变换对是一一对应的，而且有卷积定理

得出来的频域下的函数再进行傅里叶反变换就是多项式相乘的结果

这是C实现，用来检测是不是真的变快了以及结果的正确性

#if 1
static void dft(float *in, float *out, int n)
{
    int   i, j;
    float tmp[2];
    float *w = (float*)malloc(n*2*sizeof(float));
    for (i=0; i
    
    

2.快速傅里叶变换Cooley-Tukey算法

刚刚说到我们进行傅里叶变换的时间是O(N^2)，那么现在说的这个算法可以将时间降为O(NlogN)

至于原理上面推荐的那个博客已经讲得很清楚了

这个会把原理讲得更清楚一点：http://wenku.baidu.com/view/e9853c772e3f5727a5e962eb.html?from=search

大致就是利用对称性，进行计算得简化，也给出了递归与非递归的实现，参照作者的思路这里写一个非递归的C实现

注意：由于是不断地迭代，所以当n不是2的幂时需要补齐0到2的幂

#include 
     
     
      
      
#include 
      
      
       
       
#include 
       
       
        
        
#include 
        
        
         
         

void *fft_init   (int len);
void  fft_free   (void *c);
void  fft_execute(void *c, float *in, float *out);

typedef struct {
    int    N;
    float *W;
    float *T;
    int   *order;
} FFT_CONTEXT;

static int reverse_bits(int n)
{
    n = ((n & 0xAAAAAAAA) >> 1 ) | ((n & 0x55555555) << 1 );
    n = ((n & 0xCCCCCCCC) >> 2 ) | ((n & 0x33333333) << 2 );
    n = ((n & 0xF0F0F0F0) >> 4 ) | ((n & 0x0F0F0F0F) << 4 );
    n = ((n & 0xFF00FF00) >> 8 ) | ((n & 0x00FF00FF) << 8 );
    n = ((n & 0xFFFF0000) >> 16) | ((n & 0x0000FFFF) << 16);
    return n;
}

static void fft_execute_internal(FFT_CONTEXT *ctxt, float *data, int n, int w)
{
    int i;
    for (i=0; i
         
         
           W[i*w*2]); float T[2]; complex_add(C, A, B); complex_sub(T, A, B); complex_mul(D, T, W); } n /= 2; w *= 2; if (n > 1) { fft_execute_internal(ctxt, data + 0 , n, w); fft_execute_internal(ctxt, data + n * 2, n, w); } } void *fft_init(int n) { int i; FFT_CONTEXT *ctxt = (FFT_CONTEXT*)malloc(sizeof(FFT_CONTEXT)); if (!ctxt) return NULL; ctxt->N = n; ctxt->W = (float*)malloc(n * sizeof(float) * 1); ctxt->T = (float*)malloc(n * sizeof(float) * 2); ctxt->order = (int *)malloc(n * sizeof(int ) * 1); if (!ctxt->W || !ctxt->T || !ctxt->order) { fft_free(ctxt); return NULL; } for (i=0; i 
          
            N/2; i++) { ctxt->W[i * 2 + 0] = cos(2 * M_PI * i / ctxt->N); ctxt->W[i * 2 + 1] =-sin(2 * M_PI * i / ctxt->N); } int shift = (int)(32 - log(n)/log(2)); for (i=0; i 
           
             N; i++) { ctxt->order[i] = (unsigned)reverse_bits(i) >> shift; } return ctxt; } void fft_free(void *c) { FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c; if (!ctxt) return; if (ctxt->W ) free(ctxt->W ); if (ctxt->T ) free(ctxt->T ); if (ctxt->order) free(ctxt->order); free(ctxt); } void fft_execute(void *c, float *in, float *out) { int i; FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c; memcpy(ctxt->T, in, sizeof(float) * 2 * ctxt->N); fft_execute_internal(ctxt, ctxt->T, ctxt->N, 1); for (i=0; i 
            
              N; i++) { out[ctxt->order[i] * 2 + 0] = ctxt->T[i * 2 + 0]; out[ctxt->order[i] * 2 + 1] = ctxt->T[i * 2 + 1]; } } 
             
            
           
         
        
        
       
       
      
      
     
     

3.傅里叶逆变换

既然我们能从时域变频域，也要能从频域变时域

其他参考链接：

http://bbs.bccn.net/thread-465266-1-1.html