content: module Fourier where
  import Complex
  import Complex_Vectors
  import List (transpose)
  import Control.Parallel.Strategies (rpar, runEval)
  splitAtN :: Int -> [a] -> [[a]]
   
  splitAtN n [] = []
  splitAtN n xs
      = ys : (splitAtN n zs) where (ys, zs) = splitAt n xs
   
  ffth :: [ComplexF] -> [ComplexF] -> [ComplexF]
   
  ffth xs us
      | n > 1 =
          (ffmups (replikate fftEvn_2)
               (ffmus us (replikate fftOdd)))
      | n == 1 = xs
    where
        rvv = ffmus us (replikate fftOdd)
        rvv1 = replikate fftEvn_2
        fftEvn = ffth (evns xs) uEvns
        fftOdd = ffth (odds xs) uEvns
        uEvns = evns us
        evns = everyNth 2
        odds = (everyNth 2) . tail
        n = length xs
        fftEvn_2
            =
                runEval
                    (do fftEvn_2 <- rpar fftEvn
                        return fftEvn_2)
   
  ffmus :: [ComplexF] -> [ComplexF] -> [ComplexF]
   
  ffmus ns ms = zipWith (*) ns ms
   
  ffmups :: [ComplexF] -> [ComplexF] -> [ComplexF]
   
  ffmups ns ms = zipWith (+) ns ms
   
  dft :: [ComplexF] -> [ComplexF]
   
  dft xs
      = (map ((1 / (fromInt n)) *) (dfth fs xs us))
    where
        us = replikate (map conjugate (rootsOfUnity n))
        fs = factors n
        n = length xs
        fromInt = fromInteger . toInteger
   
  dftinv :: [ComplexF] -> [ComplexF]
   
  dftinv xs
      = dfth fs xs us
    where
        us = replikate (rootsOfUnity n)
        fs = factors n
        n = length xs
   
  dfth
      :: [Int] -> [ComplexF] -> [ComplexF] -> [ComplexF]
   
  dfth ((f : fs)) xs us
      | fs == [] = sfth f xs us
      | otherwise = (map sum (transpose pfts))
    where
        pfts
            =
                [(dftOfSection k) `times`
                     (us `toThe` k) | k <- [0 .. f - 1]]
        dftOfSection k
            = repl f (dfth fs (fsectionOfxs k) (us `toThe` f))
        fsectionOfxs k = everyNth f (drop k xs)
   
  sft :: [ComplexF] -> [ComplexF]
   
  sft xs
      = map ((1 / (fromInt n)) *) (sfth n xs us)
    where
        us = replikate (map conjugate (rootsOfUnity n))
        n = length xs
        fromInt = fromInteger . toInteger
   
  sftinv :: [ComplexF] -> [ComplexF]
   
  sftinv xs
      = sfth n xs us
    where
        us = replikate (rootsOfUnity n)
        n = length xs
   
  sfth :: Int -> [ComplexF] -> [ComplexF] -> [ComplexF]
   
  sfth n xs us
      =
          [sum (xs `times` upowers) | upowers <- [us `toThe`
                                                      k | k <- [0 .. n - 1]]]
   
  sct :: [Double] -> [Double]
   
  sct xs
      =
          concat
              (splitAtN 20
                   (map
                        (\ k ->
                             (xs_dot (map (cos . ((fromInt k) *)) (thetas n))))
                        [k | k <- [0 .. n - 1]]))
    where
        n = length xs
        xs_dot = sum . (zipWith (*) xs)
        fromInt = fromInteger . toInteger
   
  plus = zipWith (+)
   
  times = zipWith (*)
   
  replikate = cycle
   
  repl n = concat . ((take n) . repeat)
   
  everyNth n
      =
          (map head) .
              ((takeWhile (/= [])) . (iterate (drop n)))
   
  toThe :: [ComplexF] -> Int -> [ComplexF]
   
  toThe rootsOfUnity k = everyNth k rootsOfUnity
   
  factors :: Int -> [Int]
   
  factors n = factorsh n primes
   
  factorsh :: Int -> [Int] -> [Int]
   
  factorsh n ((p : ps))
      | n == 1 = []
      | (n `mod` p) == 0 =
          p : (factorsh (n `div` p) (p : ps))
      | otherwise = factorsh n ps
   
  primes :: [Int]
   
  primes = map head (iterate sieve [2 ..])
   
  sieve ((p : xs)) = [x | x <- xs, (x `mod` p) /= 0]
   