Функционально идиоматический БПФ

Question

Я написал этот FFT для radix-2 с целью сделать его функционально идиоматичным, не жертвуя при этом слишком высокой производительностью:

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 является самым большим вложенным хвостовым рекурсивным беспорядком, который я когда-либо писал. Я повторил алгоритм в статье в Википедии для алгоритма Кули-Тьюки , но единственные функциональные конструкции, которые я мог бы реализовать, используя функции высшего порядка, приводят к огромным ударам как по производительности, так и по эффективности памяти. Существуют ли другие решения, которые могли бы дать те же результаты, не приводя к значительным замедлениям, и при этом оставаться идиоматичными?

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

Это может быть не совсем верно (я сделал это просторефакторинг вашего кода), но я думаю, что это на самом деле более идиоматично, чем использование сложных вложенных хвостовых рекурсивных функций в этом случае.