```-----------------------------------------------------------------------------
-- |
--
-- Stability   :  experimental
-- Portability :  portable
--
-- Rader's Algorithm for computing prime length FFT's
--
-----------------------------------------------------------------------------

import Data.List
import Data.Array
import Data.Complex

-------------------------------------------------------------------------------

-- Rader's Algorithm.  We define this two ways: using direct circular
-- convolution, and FFT circular convolution.  The algorithms and
-- implementations, are esentially the same, except for how hg is
-- computed.

-- | Rader's Algorithm using direct convolution

fft_rader1 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n]
-> a -- ^ N
-> Array a (Complex b) -- ^ X[k]

where h = listArray (0,n-2) [ f!(a ^* (n-(1+n'))) | n' <- [0..(n-2)] ]
g = listArray (0,n-2) [ w!(a ^* n') | n' <- [0..(n-2)] ]
hg = listArray (0,n-2) [ sum [ h!j * g!((i-j)`mod`(n-1)) | j <- [0..(n-2)] ] | i <- [0..(n-2)] ]
f' = array (0,n-1) ((0, sum [ f!i | i <- [0..(n-1)] ]) : [ (a ^* i, f!0 + hg!i) | i <- [0..(n-2)] ])
wn = cis (-2 * pi / fromIntegral n)
w = listArray (0,n-1) \$ iterate (* wn) 1
_ ^* 0 = 1
i ^* j = (i * (i ^* (j-1))) `mod` n
a = generator n

-- | Rader's Algorithm using FFT convolution

{-# specialize fft_rader2 :: Array Int (Complex Float) -> Int -> (Array Int (Complex Float) -> Array Int (Complex Float)) -> Array Int (Complex Float) #-}
{-# specialize fft_rader2 :: Array Int (Complex Double) -> Int -> (Array Int (Complex Double) -> Array Int (Complex Double)) -> Array Int (Complex Double) #-}

fft_rader2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n]
-> a -- ^ N
-> (Array a (Complex b) -> Array a (Complex b)) -- ^ FFT function
-> Array a (Complex b) -- ^ X[k]

fft_rader2 f n fft = f'
where h = listArray (0,n-2) [ f!(a ^* (n-(1+n'))) | n' <- [0..(n-2)] ]
g = listArray (0,n-2) [ w!(a ^* n') | n' <- [0..(n-2)] ]
h' = fft h
g' = fft g
hg' = listArray (0,n-2) [ h'!i * g'!i | i <- [0..(n-2)] ]
hg = ifft hg'
f' = array (0,n-1) ((0, sum [ f!i | i <- [0..(n-1)] ]) : [ (a ^* i, f!0 + hg!i) | i <- [0..(n-2)] ])
wn = cis (-2 * pi / fromIntegral n)
w = listArray (0,n-1) \$ iterate (* wn) 1
_ ^* 0 = 1
i ^* j = (i * (i ^* (j-1))) `mod` n
a = generator n
ifft b = fmap (/ fromIntegral (n-1)) \$ fmap swap \$ fft \$ fmap swap b
swap (x:+y) = (y:+x)

-- Haskell translation of find_generator from FFTW

{-# specialize generator :: Int -> Int #-}

generator :: (Integral a) => a -> a
generator p = findgen 1
where findgen 0 = error "rader: generator: no primitive root?"
findgen x | (period x x) == (p - 1) = x
| otherwise               = findgen ((x + 1) `mod` p)
period _ 1    = 1
period x prod = 1 + (period x (prod * x `mod` p))
```