如何将递归代码更改为迭代形式

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

这是我的实现：

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 2

这并不能完全提供解决方案。但可能会帮助一些解决类似问题的人将递归算法转换为迭代算法。递归是在带有堆栈的系统中实现的。对方法的每次递归调用都会将以下信息压入堆栈：

函数参数
局部变量
Return地址

如果程序员可以使用 stack + while loop 执行上述操作，我们可以将递归算法实现为迭代算法。步骤将是

用于调用递归的参数现在调用将被推入堆栈。
然后我们进入一个 while 循环（直到堆栈为空），同时弹出来自堆栈（ LIFO ）的参数并调用核心逻辑
继续将更多参数压入堆栈并重复 (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;
}

Answer 3

我要做的第一件事就是让你的函数 "more recursive".

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 实现，我们可以相互测试它们。此时构建一些单元测试，以便可以针对已知的 "good"（或至少稳定的）行为集测试进一步的更改。

接下来，我们可以检查元素在原始向量中的组合顺序。检查长度为 4 的缓冲区。

a   b   c   d

我们递归做奇数和偶数

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

然后递归做奇数和偶数

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

这些集合的大小为 1。它们单独放置，然后我们将奇数和偶数合并。

现在我们看8，经过两次递归，元素为'owned' by:

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]

如果我们反转这些标签，并调用 e 0 和 o 1，我们得到：

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

这是二进制计数。第一位被丢弃，现在相等的元素在倒数第二个递归调用中合并。

然后舍弃前两位，合并最后一位匹配的元素

我们可以不看位，而是看每个组合的开始和步长。

第一个组合的步长等于数组长度（每个 1 个元素）。

第二个是length/2。第三个是length/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];
  }
}

none 已经过测试，但它给出了如何执行此操作的分步示例。您需要在此迭代版本上释放您之前编写的单元测试（以测试 fft 的两个早期版本）。

如何将递归代码更改为迭代形式

How to change a recursive code to iterative form

c++

iteration

algorithm

recursion

fft