如何更改递归码迭代形式

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