求Q(x)模X^n + 1的余数多项式的FFT算法

贴上IFFT算法的C语言代码:https://blog.csdn.net/great978/article/details/84306827


 

Q\left ( x \right )X^{n}+1的余数多项式,最高次数不超过n,令其余数多项式的形式为f\left ( x \right )=a_{0}+a_{1}x+a_{2}x^{2}+......+a_{n-1}x^{n-1}

 

1、因为X^{n}+1有n个解,对应余数多项式有n个,分别是\left \{ f\left ( x_{j} \right ) \right \}其中,0<=j<=n-1,而x_{_{j}}则是X^{n}+1的n个解,而x_{_{j}}是容易求的,

x_{j} = e^{i\left ( \pi /n\right )\left ( 2j+1 \right )},但是如果每个都往\left \{ f\left ( x_{j} \right ) \right \}里面带入的话,会花费大量的计算量。

 

这里提一种基于FFT的简单运算,这里要求n=2^{k}。(这里对为什么是基于FFT的我也不是很明白,有一说是

\left ( a_{0} ,a_{1},...a_{n-1}\right )\rightarrow \left ( f\left ( x_{0} \right ) ,f\left ( x_{1} \right ),....,f\left ( x_{n-1} \right )\right )的多项式的离散傅里叶变换,用快速傅里叶变换FFT计算得来的,DFS的公式这儿就不列出了)

 

2、计算算法,因为n=2^{k},把f\left ( x \right )=a_{0}+a_{1}x+a_{2}x^{2}+......+a_{n-1}x^{n-1}按奇偶项分开,每次分开一半,直到最后只剩两项的时候进行计算,然后一直递归向上,这样就把大量的乘幂运算转化为了加法运算,下面给出C语言的代码

 


//

#include "stdafx.h"

#define _CRT_SECURE_NO_WARNINGS
#include "stdlib.h"
#include "math.h"
#include "stdio.h"


#define K 3

#define Pi  3.1415927   //定义圆周率Pi
#define LEN sizeof(struct Compx)  //定义复数结构体大小

//-----定义复数结构体-----------------------
struct Compx
{
	double real;
	double imag;
}Compx;

//-----复数乘法运算函数---------------------
struct Compx mult(struct Compx b1, struct Compx b2)
{
	struct Compx b3;
	b3.real = b1.real*b2.real - b1.imag*b2.imag;
	b3.imag = b1.real*b2.imag + b1.imag*b2.real;
	return(b3);
}

//-----复数加法运算函数---------------------
struct Compx add(struct Compx a, struct Compx b)
{
	struct Compx c;
	c.real = a.real + b.real;
	c.imag = a.imag + b.imag;
	return(c);
}

struct Compx FFT(struct Compx *t , int n , struct Compx root , struct Compx result ) ;

int main(  ) 
{
	int N,i;
	N =8;

    int x[8];
    for( i = 0 ; i < N; i++ )
	{
    	scanf("%d",&x[i]);
	}
	
    struct  Compx * Source = (struct Compx *)malloc(N*LEN);			//为结构体分配存储空间
    struct  Compx * Result = (struct Compx *)malloc(N*LEN);
    struct  Compx * Root = (struct Compx *)malloc(N*LEN);

    
    //初始化======================================= 
    printf("\nSource初始化:\n");
	for (i = 0; i<N; i++)
	{
		Source[i].real = x[i];
		Source[i].imag = 0;
		printf("%.4f ", Source[i].real);
		printf("+%.4fj  ", Source[i].imag);
		printf("\n");
	}
	printf("\nResult初始化:\n");
	for (i = 0; i<N; i++)
	{
		Result[i].real = 0;
		Result[i].imag = 0;
		printf("%.4f ", Result[i].real);
		printf("+%.4fj  ", Result[i].imag);
		printf("\n");
	}
	printf("\nRoot初始化:\n");
	for (i = 0; i<N; i++)
	{
		Root[i].real = cos( 2 * Pi * ( 2 * i + 1) / (2 * N));
		Root[i].imag = sin( 2 * Pi * ( 2 * i + 1) / (2 * N));
		printf("%.4f ", Root[i].real);
		printf("+%.4fj  ", Root[i].imag);
		printf("\n");
	}
	

    for( i = 0 ; i < N ; i++) 
    {
    	Result[i] = FFT(Source , N , Root[i] , Result[i]);
	}
	
    //结果表示
	printf("\nResult计算结果:\n");
	for (i = 0; i<N; i++)
	{
		printf("%.4f ", Result[i].real);
		printf("+%.4fj  ", Result[i].imag);
		printf("\n");
	}
    
	return 0;
}

struct Compx FFT(struct Compx *t , int n , struct Compx root , struct Compx result ) 
{
	int i,j;
	struct  Compx * even = (struct Compx *)malloc((n / 2) * LEN);	
	struct  Compx * odd = (struct Compx *)malloc((n / 2) * LEN);
	
	//划分奇偶项 
	for (i = 0 , j = 0 ; i < n; i += 2 , j++ )
	{
		even[j].real = t[i].real;
		even[j].imag = t[i].imag;
	}
	for (i = 1 , j = 0 ; i < n; i += 2 , j++)
	{
		odd[j].real = t[i].real;
		odd[j].imag = t[i].imag;
	}
	
	if(n == 2)
	{
		struct Compx s = add( even[0] , mult( root , odd[0]) );
		return add(result , s);
	}
	else
	{
		return add( FFT(even , n / 2 , mult(root , root) , result ) , mult( root , FFT(odd , n / 2 , mult(root , root) , result ) ) );
	}
}

2018年11月20日更新:代码中root根值计算有错误,上述代码中已修改 :

新的代码中是:

Root[i].real = cos( 2 * Pi * ( 2 * i + 1) / (2 * N));
Root[i].imag = sin( 2 * Pi * ( 2 * i + 1) / (2 * N));

之前代码是:

Root[i].real = cos( Pi * ( 2 * i + 1) / (2 * N));
Root[i].imag = sin( Pi * ( 2 * i + 1) / (2 * N));

欢迎指正!

posted @ 2018-11-15 10:22  great978  阅读(252)  评论(0编辑  收藏  举报