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.傅里叶逆变换
既然我们能从时域变频域,也要能从频域变时域
其他参考链接: