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

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

我要做的第一件事是让你的函数"更递归"。

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]

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

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的两个早期版本）。

这是我

的实现：

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绘制小图片有助于理解它是正确的。

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

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