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;
}