[求助]这里有一个pi算法，高手来看一下

George

用c++编的
#include <complex>

#define PI 3.141592653589793238

using namespace std;

template<class T>
void fft2(complex<T> *array,int size_log_2,bool positive)
{
        complex<T> swap;
        for (int i=0;i<(1<<size_log_2);i++)
        {
                int reverse=0;
                for (int j=0;j<size_log_2;j++) reverse^=(i&1<<j)>>j<<size_log_2-j-1;
                if (reverse>i)
                {
                        swap=array[reverse];
                        array[reverse]=array;
                        array=swap;
                }
        }
        complex<T> ca,cb;
        for (int s=0;s<size_log_2;s++)
        {
                complex<T> w=exp(complex<T>(0,(positive?-1:1)*PI/(1<<s)));
                for (int a=0;a<(1<<size_log_2-s-1);a++)
                {
                        complex<T> m=1;
                        for (int b=0;b<(1<<s);b++)
                        {
                                ca=array[a*(1<<s+1)+b];
                                cb=array[((a<<1)+1)*(1<<s)+b]*m;
                                array[a*(1<<s+1)+b]=ca+cb;
                                array[((a<<1)+1)*(1<<s)+b]=ca-cb;
                                m*=w;
                        }
                }
        }
        if (!positive) for (int c=0;c<(1<<size_log_2);c++) array[c]/=(1<<size_log_2);
}

template<class T>
int fft(complex<T> *array,int size,bool positive,int step=1)
{
        int N1,N2,m1,m2,k1,k2;
        if (step==1)
        {
                int inter=size;
                while (true)
                {
                        if (inter%2==0) inter=inter/2;
                        else if (inter%3==0) inter=inter/3;
                        else if (inter%5==0) inter=inter/5;
                        else if (inter%7==0) inter=inter/7;
                        else break;
                }
                if (inter!=1) return -1;
        }
        if (size%2==0) N1=2;
        else if (size%3==0) N1=3;
        else if (size%5==0) N1=5;
        else if (size%7==0) N1=7;
        else if (size==1) return 0;
        else return -1;
        N2=size/N1;
        complex<T> log_w=complex<T>(0,(positive?-1:1)*2*PI/size);
        complex<T> log_w1=complex<T>(0,(positive?-1:1)*2*PI/N1);
        if (N2>1)
                for (m1=0;m1<N1;m1++)
                {
                        fft<T>(array+m1*step,N2,positive,N1*step);
                        for (k2=0;k2<N2;k2++) array[(m1+k2*N1)*step]*=exp((T)(m1*k2)*log_w);
                }
        complex<T> *array_=new complex<T>[size];
        complex<T> *ar=new complex<T>[N1];
        for (k2=0;k2<N2;k2++)
        {
                for (k1=0;k1<N1;k1++) ar[k1]=0;
                for (k1=0;k1<N1;k1++)
                        for (m1=0;m1<N1;m1++) ar[k1]+=array[(k2*N1+m1)*step]*exp((T)(m1*k1)*log_w1);
                for (k1=0;k1<N1;k1++) array_[k2+k1*N2]=ar[k1];
        }
        for (int s=0;s<size;s++) array[s*step]=array_;
        delete []array_; delete []ar;
        if ((step==1)&&!positive) for (int c=0;c<size;c++) array[c]/=size;
        return 0;
}

详细讲解一下，谁有c++编译器编译一下传上来吧[em04][em06][em06][em06]