HDU - 5958 New Signal Decomposition 原根下标变换+fft(未完待续)

题目链接点这里

过程看这里http://blog.csdn.net/u013368721/article/details/53001532

，，然而，，这位大佬的程序也gg了，，可能改过数据了，，不过思路是这样的，，

然而我与标程对拍，，答案总是差0.00几。。反正，，我写的fft自带低精度，大内存，高常数。。，，以后再看吧

#include<iostream>
#include<cstdio>
#include<math.h>
#include<algorithm>
#include<map>
#include<set>
#include<bitset>
#include<stack>
#include<queue>
#include<string.h>
#include<cstring>
#include<vector>
#include<time.h>
#include<stdlib.h>
using namespace std;
#define INF 0x3f3f3f3f
#define INFLL 0x3f3f3f3f3f3f3f3f
#define FIN freopen("input.txt","r",stdin)
#define mem(x,y) memset(x,y,sizeof(x))
typedef unsigned long long ULL;
typedef long long LL;
#define fuck(x) cout<<"x"<<endl;
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
typedef pair<pair<int,int>,int> PIII;
typedef pair<int,int> PII;
const double eps=1e-100;
const int MX=1111111;
int P;
int G;
double a[MX];
//fft
const double pi = acos(-1.0);
int len,res[MX],mx;//开大4倍
struct Complex
{
double r,i;
Complex(double r=0,double i=0):r(r),i(i) {};
Complex operator+(const Complex &rhs)
{
return Complex(r + rhs.r,i + rhs.i);
}
Complex operator-(const Complex &rhs)
{
return Complex(r - rhs.r,i - rhs.i);
}
Complex operator*(const Complex &rhs)
{
return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i);
}
} va[MX],vb[MX];
void rader(Complex F[],int len) //len = 2^M,reverse F[i] with F[j] j为i二进制反转
{
int j = len >> 1;
for(int i = 1; i < len - 1; ++i)
{
if(i < j) swap(F[i],F[j]); // reverse
int k = len>>1;
while(j>=k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}

void FFT(Complex F[],int len,int t)
{
rader(F,len);
for(int h=2; h<=len; h<<=1)
{
Complex wn(cos(-t*2*pi/h),sin(-t*2*pi/h));
for(int j=0; j<len; j+=h)
{
Complex E(1,0); //旋转因子
for(int k=j; k<j+h/2; ++k)
{
Complex u = F[k];
Complex v = E*F[k+h/2];
F[k] = u+v;
F[k+h/2] = u-v;
E=E*wn;
}
}
}
if(t==-1) //IDFT
for(int i=0; i<len; ++i)
F[i].r/=len;
}
/*
inline void FFT(Complex *a,int n,int r)
{
rader(a,len);
for (int i=1; i<n; i<<=1)
{
Complex wn(cos(pi/i),r*sin(pi/i));
for (int j=0; j<n; j+=(i<<1))
{
Complex w(1,0);
for (int k=0; k<i; k++,w=w*wn)
{
Complex x=a[j+k],y=w*a[j+k+i];
a[j+k]=x+y;
a[j+k+i]=x-y;
}
}
}
if (r==-1) for (int i=0; i<n; i++) a[i].r/=n;
}
*/
void Conv(Complex a[],Complex b[],int len) //求卷积
{
FFT(a,len,1);
FFT(b,len,1);
for(int i=0; i<len; ++i) a[i] = a[i]*b[i];
FFT(a,len,-1);
}
int change[MX];
double ans[MX];
double ss(int ret)
{
return pow(2,sin(2*pi*ret/P)*sin(2*pi*ret/P)*sin(2*pi*ret/P));
}
void gao()
{
mx=P-2+P-2;
len=1;
while(len<mx)len<<=1;//mx为卷积后最大下标
int ret=1;
va[P-2-0].r=a[1];
vb[0].r=ss(ret);
change[0]=1;
for(int i=1; i<P-1; i++)
{
ret=(ret*G)%P;
change[i]=ret;
va[P-2-i].r=a[ret];
va[P-2-i].i=0;
vb[i].r=ss(ret);
vb[i].i=0;
}
for(int i=P-1; i<len; i++)va[i].i=va[i].r=vb[i].i=vb[i].r=0;
Conv(va,vb,len);
ans[1]=va[P-2].r;
ans[0]=a[P-1];
for(int i=0; i<P-2; i++)
{
ans[change[i+1]]=va[i].r+va[i+P-1].r;
ans[0]+=a[i+1];
}
for(int i=0; i<=P-1; i++)
{
printf("%.3f ",ans[i]+a[0]);
}
puts("");
}
int main()
{
FIN;
freopen("output1.txt","w",stdout);
while(~scanf("%d",&P))
{
if(P==13)G=2;
else if(P==103)G=5;
else G=2;
for(int i=0; i<P; i++)scanf("%lf",&a[i]);
gao();
}
return 0;
}