super pi源码

经典super pi源码，c语言实现的FFT算法，用于测试cpu浮点计算性能。 Super PI是利用CPU的浮点运算能力来计算出π（圆周率），所以目前普遍被超频玩家用做测试系统稳定性和测试，CPU计算完后特定位数圆周率所需的时间。 make #如果是64位，请把 -march=i686 去掉； ./pi_css5 $((1<<20)) # 小数点后2^20次方的小数点位数

/*
Fast Fourier/Cosine/Sine Transform
    dimension   :one
    data length :power of 2
    decimation  :frequency
    radix       :split-radix
    data        :inplace
    table       :not use
functions
    cdft: Complex Discrete Fourier Transform
    rdft: Real Discrete Fourier Transform
    ddct: Discrete Cosine Transform
    ddst: Discrete Sine Transform
    dfct: Cosine Transform of RDFT （Real Symmetric DFT）
    dfst: Sine Transform of RDFT （Real Anti-symmetric DFT）
function prototypes
    void cdft（int int double *）;
    void rdft（int int double *）;
    void ddct（int int double *）;
    void ddst（int int double *）;
    void dfct（int double *）;
    void dfst（int double *）;


-------- Complex DFT （Discrete Fourier Transform） --------
    [definition]
        
            X[k] = sum_j=0^n-1 x[j]*exp（2*pi*i*j*k/n） 0<=k        
            X[k] = sum_j=0^n-1 x[j]*exp（-2*pi*i*j*k/n） 0<=k        （notes: sum_j=0^n-1 is a summation from j=0 to n-1）
    [usage]
        
            cdft（2*n 1 a）;
        
            cdft（2*n -1 a）;
    [parameters]
        2*n            :data length （int）
                        n >= 1 n = power of 2
        a[0...2*n-1]   :input/output data （double *）
                        input data
                            a[2*j] = Re（x[j]） 
                            a[2*j+1] = Im（x[j]） 0<=j                        output data
                            a[2*k] = Re（X[k]） 
                            a[2*k+1] = Im（X[k]） 0<=k    [remark]
        Inverse of 
            cdft（2*n -1 a）;
        is 
            cdft（2*n 1 a）;
            for （j = 0; j <= 2 * n - 1; j++） {
                a[j] *= 1.0 / n;
            }
        .


-------- Real DFT / Inverse of Real DFT --------
    [definition]
         RDFT
            R[k] = sum_j=0^n-1 a[j]*cos（2*pi*j*k/n） 0<=k<=n/2
            I[k] = sum_j=0^n-1 a[j]*sin（2*pi*j*k/n） 0         IRDFT （excluding scale）
            a[k] = （R[0] + R[n/2]*cos（pi*k））/2 + 
                   sum_j=1^n/2-1 R[j]*cos（2*pi*j*k/n） + 
                   sum_j=1^n/2-1 I[j]*sin（2*pi*j*k/n） 0<=k    [usage]
        
            rdft（n 1 a）;
        
            rdft（n -1 a）;
    [parameters]
        n              :data length （int）
                        n >= 2 n = power of 2
        a[0...n-1]     :input/output data （double *）
                        
                            output data
                                a[2*k] = R[k] 0<=k                                a[2*k+1] = I[k] 0                                a[1] = R[n/2]
                        
                            input data
                                a[2*j] = R[j] 0<=j                                a[2*j+1] = I[j] 0                                a[1] = R[n/2]
    [remark]
        Inverse of 
            rdft（n 1 a）;
        is 
            rdft（n -1 a）;
            for （j =