如何更改遞歸碼迭代形式

Question

我實現下面一個簡單的FFT（最終變被忽略）：

typedef complex<double> base; 
vector<base> w; 
int FFTN = 1024; 
void fft(vector<base> &fa){ 
    int n = fa.size(); 
    if (n==1) return; 
    int half = (n>>1); 
    vector<base> odd(half),even(half); 
    for(int i=0,j = 0;i<n;i+=2,j++) { 
     even[j] = fa[i]; 
     odd[j] = fa[i+1];  
    }  
    fft(odd); 
    fft(even);  
    int fact = FFTN/n;  
    for (int i=0;i<half;i++){   
     fa[i] = even[i] + odd[i] * w[i * fact]; 
     fa[i + half] = even[i] - odd[i] * w[i * fact]; 
    } 
} 

它運作良好。但我堅持將其轉換爲迭代形式。我已經嘗試到目前爲止：

int n = fa.size(); 
int fact = (FFTN>>1); 
int half = 1; 
while(half<n){ 
    for(int i=0;i<n/half;i+=2){ 
     base even = fa[i], odd = fa[i+1]; 
     fa[i] = even + odd * w[i*fact]; 
     fa[i+half] = even - odd*w[i*fact];        
    } 
    for(int j=0;j<n/half;j++) 
     fa[j] = fa[j+half]; 
    fact >>= 1; 
    half <<= 1; 
} 

有人可以幫助我的轉換技巧嗎？

Answer 1

我要做的第一件事就是讓你的函數「更加遞歸」。

void fft(base* fa, size_t stride, size_t n) { 
    if (n==1) return; 
    int half = (n>>1); 
    fft(fa+stride, stride*2, half); // odd 
    fft(fa, stride*2, half); // even 
    int fact = FFTN/n; 
    for (int i=0;i<half;i++){  
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact]; 
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact]; 
    } 
} 
void fft(std::vector<base>& fa){ fft(fa.data(), 1, fa.size()); } 

現在我們在緩衝區內就地執行我們的fft。

由於我們現在有兩個不同的fft實現，我們可以對它們進行測試。在這一點上建立一些單元測試，因此可以對已知的「良好」（或至少是穩定的）行爲進行進一步的更改測試。

接下來，我們可以檢查哪些元素在原始向量中組合的順序。檢查長度的緩衝器4

a b c d 

我們做遞歸奇數和偶數

a[e] b[o] a[e] d[o] 

然後做遞歸奇數和偶數

a[ee] b[oe] a[eo] d[oo] 

這些集尺寸1.他們離開單獨，然後我們結合奇數和偶數。

現在我們來看看8.經過兩年遞歸由裏的元素是 '擁有'：

0[ee] 1[oe] 2[eo] 3[oo] 4[ee] 5[oe] 6[eo] 7[oo] 

和後3：

0[eee] 1[oee] 2[eoe] 3[ooe] 4[eeo] 5[oeo] 6[eoo] 7[ooo] 

如果我們扭轉這些標籤，並呼籲e0和o1，我們得到：

0[000] 1[001] 2[010] 3[011] 4[100] 5[101] 6[110] 7[111] 

這是二進制計數。第一位被丟棄，並且現在相等的元素在第二次到最後一次遞歸調用中被合併。

然後丟棄前兩位，並組合具有匹配的最後一位的元素。

我們可以不看比特，看看每個組合的開始和步幅的長度。

第一個組合是步長等於陣列長度（每個1個元素）。

第二個是長度/ 2。第三是長度/ 4。

這繼續直到步幅1.

子陣列相結合的數量等於步幅長度。

所以

for(size_t stride = n; stride = stride/2; stride!=0) { 
    for (size_t offset = 0; offset != stride; ++offset) { 
    fft_process(array+offset, stride, n/stride); 
    } 
} 

其中fft_process是基於關閉的：

int fact = FFTN/n; 
    for (int i=0;i<half;i++){  
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact]; 
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact]; 
    } 

也許是這樣的：

void fft_process(base* fa, size_t stride, size_t n) { 
    int fact = FFTN/n; // equals stride I think! Assuming outermost n is 1024. 
    for (int i=0;i<half;i++){  
    fa[i] = fa[stride*2*i] + fa[stride*2+i+stride] * w[i * fact]; 
    fa[i + half] = fa[stride*2*i] - fa[stride*2+i+stride] * w[i * fact]; 
    } 
} 

這一切都不是測試，但它給出了一個一步一步如何做到這一點的例子。您將希望在此迭代版本中釋放您之前編寫的單元測試（以測試fft的兩個早期版本）。

Answer 2

這是我的實現：

typedef complex<double> Data; 

const double PI = acos(-1); 

// Merges [low, (low + high)/2) with [(low + high)/2, high) parts. 
void merge(vector<Data>& b, int low, int high) { 
    int n = high - low; 
    Data cur(1), mul(cos(2. * PI/n), sin(2. * PI/n)); 
    for (int i = low; i < low + n/2; i++) { 
     Data temp = b[i + n/2] * cur; 
     b[i + n/2] = b[i] - temp; 
     b[i] = b[i] + temp; 
     cur = cur * mul; 
    } 
} 

// Computes FFT for the vector b. 
void do_fft(vector<Data>& b) { 
    int n = b.size(); 
    int hi = 0; 
    while ((1 << hi) < n) 
     hi++; 
    hi--; 
    // Permutes the input vector in a specific way. 
    vector<int> p(n); 
    for (int i = 0; i < n; i++) 
     for (int b = hi; b >= 0; b--) 
      if (i & (1 << b)) 
       p[i] |= (1 << (hi - b)); 
    vector<Data> buf(n); 
    for (int i = 0; i < n; i++) 
     buf[i] = b[p[i]]; 
    copy(buf.begin(), buf.end(), b.begin()); 
    for (int h = 2; h <= n; h *= 2) 
     for (int i = 0; i < n; i += h) 
      merge(b, i, i + h); 
} 

這個實現的想法是置換給定矢量中，我們需要相鄰子矢量在每個步驟合併（即，[0,0這樣的方式]與[1，1]，[2，2]與[3，3]等在第一步，[0，1]與[2,3]，[4,5]與[6,7]在第二步等等）。事實證明，這些元素應該按照以下方式進行置換：我們應該採用元素索引的二進制表示，將其逆轉，並將具有相反索引的元素放到當前位置。我無法證實這是嚴格的，但爲n = 8或n = 16繪製小圖可以幫助理解它是正確的。

Answer 3

這並不完全提供解決方案。但是可能會幫助一些解決類似問題的人將遞歸算法轉換爲迭代算法。遞歸在帶有堆棧的系統中實現。一種方法，每次遞歸調用推，以下信息到堆棧中：

1. 功能參數
2. 局部變量
3. 返回地址

如果程序員可以做上面有stack + while loop，我們可以實現迭代算法的遞歸算法。步驟將爲

1. 用於調用遞歸調用 調用的參數現在將推送到堆棧。
2. 然後，我們去了一個while循環（直到棧空）內，而從彈出堆棧 參數（LIFO），並調用核心邏輯
3. 繼續推動進一步的參數堆棧和重複（2）直到棧空。

使用上述方法進行迭代階乘計算的代碼示例。

int coreLogic(int current, int recursiveParameter) { 
    return current * recursiveParameter ; 
} 

int factorial(int n) { 

    std::stack<int> parameterStack ; 

    int tempFactorial = 1; 
    //parameters that would have been used to invoke the recursive call will now be pushed to stack 
    parameterStack.push(n); 

    while(!parameterStack.empty()) { 
     //popping arguments from stack 
     int current = parameterStack.top(); 
     parameterStack.pop(); 

     //and invoking core logic 
     tempFactorial = coreLogic(tempFactorial, current); 

     if(current > 1) { 
      //parameters that would have been used to invoke the recursive call will now be pushed to stack 
      parameterStack.push( current - 1); 
     } 

     /* 
     *if a divide and conquer algorithm like quick sort then again push right side args to stack 
     * - appers case in question 
     *if(condition) { 
     * parameterStack.push(args ); 
     *} 
     */ 
    } 
    return tempFactorial; 
}