functional-programming - 功能惯用的 FFT

Question

我编写了这个 radix-2 FFT，其目标是使其在功能上符合习惯，而不牺牲太多性能:

let reverse x bits =
    let rec reverse' x bits y =
        match bits with
        | 0 -> y
        | _ -> ((y <<< 1) ||| (x &&& 1))
               |> reverse' (x >>> 1) (bits - 1) 
     reverse' x bits 0 

let radix2 (vector: Complex[]) (direction: int) =
    let z = vector.Length
    let depth = floor(Math.Log(double z, 2.0)) |> int
    if (1 <<< depth) <> z then failwith "Vector length is not a power of 2"

    // Complex roots of unity; "twiddle factors"
    let unity: Complex[] =
        let xpn = float direction * Math.PI / double z
        Array.Parallel.init<Complex> (z/2) (fun i ->
            Complex.FromPolarCoordinates(1.0, (float i) * xpn))

    // Permutes elements of input vector via bit-reversal permutation
    let pvec = Array.Parallel.init z (fun i -> vector.[reverse i depth])

    let outerLoop (vec: Complex[]) =
        let rec recLoop size =
            if size <= z then
                let mid, step = size / 2, z / size
                let rec inrecLoop i =
                    if i < z then
                        let rec bottomLoop idx k =
                            if idx < i + mid then
                                let temp = vec.[idx + mid] * unity.[k]
                                vec.[idx + mid] <- (vec.[idx] - temp)
                                vec.[idx] <- (vec.[idx] + temp)
                                bottomLoop (idx + 1) (k + step)
                        bottomLoop i 0
                        inrecLoop (i + size)
                inrecLoop 0
                recLoop (size * 2)
        recLoop 2
        vec

    outerLoop pvec

outerLoop段是我写过的最大的嵌套尾递归困惑。我复制了维基百科文章中的算法 Cooley-Tukey algorithm ，但我认为使用高阶函数实现的唯一函数结构会对性能和内存效率造成巨大影响。是否有其他解决方案可以产生相同的结果，而不会导致速度大幅下降，同时仍然是惯用的？

Answer 1

我不是该算法如何工作的专家，因此可能有一个很好的功能实现，但值得注意的是，使用局部突变在 F# 中是完全惯用的。

您的radix2函数从外部起作用 - 它接受vector数组作为输入，从不改变它，创建一个新数组pvec然后它初始化(沿途使用一些突变)然后返回它。这与 Array.map 等内置函数使用的模式类似(初始化一个新数组，对其进行变异，然后返回它)。这通常是一种明智的做法，因为某些算法最好使用变异来编写。

在这种情况下，使用局部可变变量和循环是完全合理的。与尾递归版本相比，这样做将使您的代码更具可读性。我还没有对此进行测试，但我对 outerLoop 函数的天真翻译只是使用三个嵌套循环 - 像这样:

let mutable size = 2
while size <= z do
    let mid, step = size / 2, z / size
    let mutable i = 0
    while i < z do
        for j in 0 .. mid - 1 do
            let idx, k = i + j, step * j
            let temp = pvec.[idx + mid] * unity.[k]
            pvec.[idx + mid] <- (pvec.[idx] - temp)
            pvec.[idx] <- (pvec.[idx] + temp)
        i <- i + size
    size <- size * 2

这可能不完全正确(我这样做只是为了重构你的代码)，但我认为它实际上比在这种情况下使用复杂的嵌套尾递归函数更惯用。