diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2021
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Makefile b/Makefile
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,8 @@
+run-test:	update-test
+	runhaskell Setup configure --user --enable-tests
+	runhaskell Setup build
+	runhaskell Setup haddock
+	./dist/build/fftw-test/fftw-test
+
+update-test:
+	doctest-extract-0.1 -i src/ -o test/ --library-main=Main $$(cat test-module.list)
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/comfort-fftw.cabal b/comfort-fftw.cabal
new file mode 100644
--- /dev/null
+++ b/comfort-fftw.cabal
@@ -0,0 +1,82 @@
+Cabal-Version:    2.2
+Name:             comfort-fftw
+Version:          0.0
+License:          BSD-3-Clause
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         https://hub.darcs.net/thielema/comfort-fftw/
+Category:         Math
+Synopsis:         High-level interface to FFTW (Fast Fourier Transform) based on comfort-array
+Description:
+  FFTW claims to be the fastest Fourier Transform in the West.
+  This is a high-level interface to @libfftw@.
+  We re-use the type classes from @netlib-ffi@.
+  .
+  For rather simple examples
+  see the packages @align-audio@ and @morbus-meniere@.
+  .
+  See also package @fft@.
+Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.6.5
+Build-Type:       Simple
+Extra-Source-Files:
+  Makefile
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/comfort-fftw/
+
+Source-Repository head
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/comfort-fftw/
+
+Library
+  Build-Depends:
+    fftw-ffi >=0.0 && <0.2,
+    comfort-array >=0.5 && <0.6,
+    netlib-ffi >=0.0 && <0.2,
+    QuickCheck >=2 && <3,
+    deepseq >=1.3 && <1.5,
+    base >=4.5 && <5
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src
+  Default-Language: Haskell98
+  Exposed-Modules:
+    Numeric.FFTW.Rank1
+    Numeric.FFTW.Rank2
+    Numeric.FFTW.Rank3
+    Numeric.FFTW.RankN
+    Numeric.FFTW.Batch
+    Numeric.FFTW.Shape
+  Other-Modules:
+    Numeric.FFTW.Private
+
+Test-Suite fftw-test
+  Type: exitcode-stdio-1.0
+  Build-Depends:
+    comfort-fftw,
+    QuickCheck,
+    doctest-exitcode-stdio >=0.0 && <0.1,
+    doctest-lib >=0.1 && <0.2,
+    comfort-array,
+    netlib-ffi,
+    non-empty >=0.3.2 && <0.4,
+    storable-record >=0.0.6 && <0.1,
+    deepseq,
+    base
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   test
+  Default-Language: Haskell98
+  Main-Is:          Main.hs
+  Other-Modules:
+    Test.Main
+    Test.Numeric.FFTW.Rank1
+    Test.Numeric.FFTW.Rank2
+    Test.Numeric.FFTW.Rank3
+    Test.Numeric.FFTW.RankN
+    Test.Numeric.FFTW.Batch
+    Test.Numeric.FFTW.Shape
+    Test.Numeric.FFTW.Common
diff --git a/src/Numeric/FFTW/Batch.hs b/src/Numeric/FFTW/Batch.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Batch.hs
@@ -0,0 +1,184 @@
+module Numeric.FFTW.Batch (
+   fourier, Sign(..), flipSign,
+   fourierRC,
+   fourierCR,
+   ) where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+import Numeric.FFTW.Private
+         (Sign(..), flipSign, ffiSign, run, runCopiedArray, withDims)
+
+import Foreign.Marshal.Array (pokeArray)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (nullPtr)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import qualified Data.List as List
+import Data.Complex (Complex)
+import Data.Monoid ((<>))
+
+import Control.Monad (when)
+
+
+{- $setup
+>>> import Test.Numeric.FFTW.Common
+>>>    (approxReal, approxComplex, floatTol, doubleTol,
+>>>     genCyclicArray2, genCyclicArray3, immutable,
+>>>     arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+>>>     floatList, complexFloatList)
+>>>
+>>> import qualified Numeric.FFTW.Batch as Batch
+>>> import qualified Numeric.FFTW.Rank1 as Trafo1
+>>> import qualified Numeric.FFTW.Rank2 as Trafo2
+>>> import qualified Numeric.FFTW.Rank3 as Trafo3
+>>> import qualified Numeric.FFTW.Shape as Spectrum
+>>>
+>>> import qualified Numeric.Netlib.Class as Class
+>>>
+>>> import qualified Data.Array.Comfort.Boxed as BoxedArray
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import qualified Data.Array.Comfort.Shape as Shape
+>>> import qualified Data.Foldable as Fold
+>>> import Data.Array.Comfort.Storable (Array)
+>>> import Data.Complex (Complex)
+>>>
+>>> import qualified Test.QuickCheck as QC
+>>>
+>>> array1s :: (Shape.C sh0, Shape.C sh1) =>
+>>>    Array (sh0,sh1) a -> Array (sh0,(),sh1) a
+>>> array1s = Array.mapShape (\(sh0,sh1) -> (sh0,(),sh1))
+>>>
+>>> unarray1s :: (Shape.C sh0, Shape.C sh1) =>
+>>>    Array (sh0,(),sh1) a -> Array (sh0,sh1) a
+>>> unarray1s = Array.mapShape (\(sh0,(),sh1) -> (sh0,sh1))
+>>>
+>>> array2s :: (Shape.C sh0, Shape.C sh1, Shape.C sh2) =>
+>>>    Array (sh0,sh1,sh2) a -> Array (sh0,(sh1,sh2)) a
+>>> array2s = Array.mapShape (\(sh0,sh1,sh2) -> (sh0,(sh1,sh2)))
+>>>
+>>> allApproxReal ::
+>>>    (Shape.C sh0, Shape.C sh1, Class.Real a, Eq sh0, Eq sh1) =>
+>>>    a ->
+>>>    BoxedArray.Array sh0 (Array sh1 a) ->
+>>>    BoxedArray.Array sh0 (Array sh1 a) ->
+>>>    Bool
+>>> allApproxReal tol xs ys =
+>>>    Fold.and $ BoxedArray.zipWith (approxReal tol) xs ys
+>>>
+>>> allApproxComplex ::
+>>>    (Shape.C sh0, Shape.C sh1, Class.Real a, Eq sh0, Eq sh1) =>
+>>>    a ->
+>>>    BoxedArray.Array sh0 (Array sh1 (Complex a)) ->
+>>>    BoxedArray.Array sh0 (Array sh1 (Complex a)) ->
+>>>    Bool
+>>> allApproxComplex tol xs ys =
+>>>    Fold.and $ BoxedArray.zipWith (approxComplex tol) xs ys
+-}
+
+
+{- |
+>>> complexFloatList $ Batch.fourier Batch.Forward $ Array.fromList (Shape.ZeroBased (5::Int), Shape.Cyclic (5::Int)) [1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1]
+[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,(-0.81):+(-0.59),0.31:+0.95,0.31:+(-0.95),(-0.81):+0.59,1.00:+0.00,(-0.81):+0.59,0.31:+(-0.95),0.31:+0.95,(-0.81):+(-0.59),1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95)]
+>>> complexFloatList $ Batch.fourier Batch.Backward $ Array.fromList (Shape.ZeroBased (5::Int), Shape.Cyclic (5::Int)) [1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1]
+[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95),1.00:+0.00,(-0.81):+0.59,0.31:+(-0.95),0.31:+0.95,(-0.81):+(-0.59),1.00:+0.00,(-0.81):+(-0.59),0.31:+0.95,0.31:+(-0.95),(-0.81):+0.59,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95]
+
+prop> QC.forAll genCyclicArray2 $ \xs sign -> allApproxComplex floatTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo1.fourier sign) $ Array.toRowArray xs)
+prop> QC.forAll genCyclicArray2 $ \xs sign -> allApproxComplex doubleTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo1.fourier sign) $ Array.toRowArray xs)
+
+prop> QC.forAll (fmap array2s genCyclicArray3) $ \xs sign -> allApproxComplex floatTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo2.fourier sign) $ Array.toRowArray xs)
+prop> QC.forAll (fmap array2s genCyclicArray3) $ \xs sign -> allApproxComplex doubleTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo2.fourier sign) $ Array.toRowArray xs)
+
+prop> \sign -> QC.forAll genCyclicArray2 $ immutable (Batch.fourier sign) . arrayComplexFloat
+prop> \sign -> QC.forAll genCyclicArray2 $ immutable (Batch.fourier sign) . arrayComplexDouble
+
+prop> \sign -> QC.forAll genCyclicArray3 $ immutable (Batch.fourier sign) . array2s . arrayComplexFloat
+prop> \sign -> QC.forAll genCyclicArray3 $ immutable (Batch.fourier sign) . array2s . arrayComplexDouble
+-}
+fourier ::
+   (Shape.C loop, Spectrum.MultiCyclic sh, Class.Real a) =>
+   Sign -> Array (loop,sh) (Complex a) -> Array (loop,sh) (Complex a)
+fourier sign (Array (loop,sh) x) =
+   Array.unsafeCreateWithSize (loop,sh) $ \n yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $ run $
+   let sliceSize = fromIntegral $ Shape.size sh in
+   withDims (Spectrum.cyclicDimensions sh) $ \rank dimPtr ->
+      FFI.planManyDFT rank dimPtr (fromIntegral $ Shape.size loop)
+         xPtr nullPtr 1 sliceSize
+         yPtr nullPtr 1 sliceSize
+         (ffiSign sign) (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll genCyclicArray2 $ \xs -> allApproxComplex floatTol (Array.toRowArray $ unarray1s $ Batch.fourierRC $ array1s xs) (fmap Trafo1.fourierRC $ Array.toRowArray xs)
+prop> QC.forAll genCyclicArray2 $ \xs -> allApproxComplex doubleTol (Array.toRowArray $ unarray1s $ Batch.fourierRC $ array1s xs) (fmap Trafo1.fourierRC $ Array.toRowArray xs)
+
+prop> QC.forAll genCyclicArray3 $ \xs -> allApproxComplex floatTol (Array.toRowArray $ array2s $ Batch.fourierRC xs) (fmap Trafo2.fourierRC $ Array.toRowArray $ array2s xs)
+prop> QC.forAll genCyclicArray3 $ \xs -> allApproxComplex doubleTol (Array.toRowArray $ array2s $ Batch.fourierRC xs) (fmap Trafo2.fourierRC $ Array.toRowArray $ array2s xs)
+
+prop> QC.forAll genCyclicArray2 $ immutable Batch.fourierRC . array1s . arrayFloat
+prop> QC.forAll genCyclicArray2 $ immutable Batch.fourierRC . array1s . arrayDouble
+
+prop> QC.forAll genCyclicArray3 $ immutable Batch.fourierRC . arrayFloat
+prop> QC.forAll genCyclicArray3 $ immutable Batch.fourierRC . arrayDouble
+-}
+fourierRC ::
+   (Shape.C loop, Spectrum.MultiCyclic sh0, Integral n1, Class.Real a) =>
+   Array (loop, sh0, Shape.Cyclic n1) a ->
+   Array (loop, sh0, Spectrum.Half n1) (Complex a)
+fourierRC (Array (loop, sh0, sh1x@(Shape.Cyclic n1)) x) =
+   let sh1y = Spectrum.Half n1 in
+   Array.unsafeCreate (loop, sh0, sh1y) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   let dims = Spectrum.cyclicDimensions (sh0, sh1x) in
+   if any (<=0) dims
+      then pokeArray yPtr $ List.genericReplicate (Shape.size (loop,sh0)) 0
+      else run $
+         withDims dims $ \rank dimPtr ->
+         FFI.planManyDFTr2c rank dimPtr (fromIntegral $ Shape.size loop)
+            xPtr nullPtr 1 (fromIntegral $ Shape.size (sh0, sh1x))
+            yPtr nullPtr 1 (fromIntegral $ Shape.size (sh0, sh1y))
+            (FFI.estimate <> FFI.preserveInput)
+
+{- |
+>>> floatList $ Batch.fourierCR $ Array.fromList (Shape.ZeroBased (2::Int), (), Spectrum.Half (3::Int)) [1,0, 0,-1]
+[1.00,1.00,1.00,-2.00,1.00,1.00]
+
+>>> floatList $ Batch.fourierCR $ Array.fromList ((), Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0]
+[1.00,1.00,1.00]
+>>> floatList $ Batch.fourierCR $ Array.fromList ((), Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0]
+[1.00,1.00,1.00]
+
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> allApproxReal floatTol (Array.toRowArray $ unarray1s $ Batch.fourierCR $ array1s xs) (fmap Trafo1.fourierCR $ Array.toRowArray xs)
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> allApproxReal doubleTol (Array.toRowArray $ unarray1s $ Batch.fourierCR $ array1s xs) (fmap Trafo1.fourierCR $ Array.toRowArray xs)
+
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> allApproxReal floatTol (Array.toRowArray $ array2s $ Batch.fourierCR xs) (fmap Trafo2.fourierCR $ Array.toRowArray $ array2s xs)
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> allApproxReal doubleTol (Array.toRowArray $ array2s $ Batch.fourierCR xs) (fmap Trafo2.fourierCR $ Array.toRowArray $ array2s xs)
+
+
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Batch.fourierCR . array1s . arrayComplexFloat
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Batch.fourierCR . array1s . arrayComplexDouble
+
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Batch.fourierCR . arrayComplexFloat
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Batch.fourierCR . arrayComplexDouble
+-}
+fourierCR ::
+   (Shape.C loop, Spectrum.MultiCyclic sh0, Integral n1, Class.Real a) =>
+   Array (loop, sh0, Spectrum.Half n1) (Complex a) ->
+   Array (loop, sh0, Shape.Cyclic n1) a
+fourierCR arr@(Array (loop, sh0, sh1x@(Spectrum.Half n1)) _x) =
+   let sh1y = Shape.Cyclic n1 in
+   let sh = (sh0, sh1y) in
+   let dims = Spectrum.cyclicDimensions sh in
+   Array.unsafeCreate (loop, sh0, sh1y) $ \yPtr ->
+   when (all (>0) dims) $
+   runCopiedArray arr $ \xPtr ->
+   withDims dims $ \rank dimPtr ->
+      FFI.planManyDFTc2r rank dimPtr (fromIntegral $ Shape.size loop)
+         xPtr nullPtr 1 (fromIntegral $ Shape.size (sh0, sh1x))
+         yPtr nullPtr 1 (fromIntegral $ Shape.size (sh0, sh1y))
+         (FFI.estimate <> FFI.destroyInput)
diff --git a/src/Numeric/FFTW/Private.hs b/src/Numeric/FFTW/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Private.hs
@@ -0,0 +1,87 @@
+module Numeric.FFTW.Private where
+
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Foreign.C.Types as C
+import Foreign.Marshal.Array (copyArray, withArrayLen)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr, castPtr)
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import Control.Concurrent.MVar (MVar, newMVar, withMVar)
+import Control.Exception (bracket)
+
+import qualified Test.QuickCheck as QC
+
+
+{- |
+This lock must be taken during /planning/ of any transform.
+The FFTW library is not thread-safe in the planning phase.
+Thankfully, the lock is not needed during the execution phase.
+-}
+{-# NOINLINE lock #-}
+lock :: MVar ()
+lock = unsafePerformIO $ newMVar ()
+
+withLock :: IO a -> IO a
+withLock = withMVar lock . const
+
+
+run :: Class.Real a => IO (FFI.Plan a) -> IO ()
+run planner = bracket (withLock planner) FFI.destroyPlan FFI.execute
+
+runCopiedArray ::
+   (Shape.C sh, Class.Floating b, Class.Real a) =>
+   Array sh b -> (Ptr b -> IO (FFI.Plan a)) -> IO ()
+runCopiedArray (Array sh x) planner =
+   withForeignPtr x $ \ptr ->
+   let n = Shape.size sh in
+   allocaArray n $ \tmpPtr -> run $ do
+      plan <- planner tmpPtr
+      copyArray tmpPtr ptr n
+      return plan
+
+
+
+{- |
+Order is chosen such that the numeric sign is @(-1) ^ fromEnum sign@.
+-}
+data Sign = Backward | Forward
+   deriving (Eq, Ord, Enum, Show)
+
+instance QC.Arbitrary Sign where
+   arbitrary = QC.elements [Backward, Forward]
+
+flipSign :: Sign -> Sign
+flipSign Backward = Forward
+flipSign Forward = Backward
+
+ffiSign :: Sign -> FFI.Sign
+ffiSign Backward = FFI.backward
+ffiSign Forward = FFI.forward
+
+
+allocaArray :: (Class.Floating a) => Int -> (Ptr a -> IO b) -> IO b
+allocaArray n =
+   case mallocFree of
+      MallocFree alloc free -> bracket (alloc (fromIntegral n)) (free . castPtr)
+
+data MallocFree a = MallocFree (C.CSize -> IO (Ptr a)) (Ptr a -> IO ())
+
+mallocFree :: (Class.Floating a) => MallocFree a
+mallocFree =
+   Class.switchFloating
+      (MallocFree FFI.allocReal FFI.free)
+      (MallocFree FFI.allocReal FFI.free)
+      (MallocFree FFI.allocComplex FFI.freeComplex)
+      (MallocFree FFI.allocComplex FFI.freeComplex)
+
+
+withDims :: [C.CInt] -> (C.CInt -> Ptr C.CInt -> IO a) -> IO a
+withDims dims f =
+   withArrayLen dims $ \len dimPtr -> f (fromIntegral len) dimPtr
diff --git a/src/Numeric/FFTW/Rank1.hs b/src/Numeric/FFTW/Rank1.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Rank1.hs
@@ -0,0 +1,254 @@
+{-# LANGUAGE GADTs #-}
+module Numeric.FFTW.Rank1 (
+   fourier, Sign(..), flipSign,
+   fourierRC,
+   fourierCR,
+   cosine,
+   sine,
+   hartley,
+   ) where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+import Numeric.FFTW.Private (Sign(..), flipSign, ffiSign, run)
+
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Storable (poke)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import Data.Complex (Complex)
+import Data.Monoid ((<>))
+
+import Control.Monad (when)
+
+
+{- $setup
+>>> :set -XGADTs
+>>> import Test.Numeric.FFTW.Common
+>>>    (approxReal, approxComplex, floatTol, doubleTol, normInf,
+>>>     adjust, scalarProduct, split, genCyclicArray1, immutable,
+>>>     arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+>>>     floatList, complexFloatList)
+>>>
+>>> import qualified Numeric.FFTW.Rank1 as Trafo1
+>>> import qualified Numeric.FFTW.Shape as Spectrum
+>>> import Numeric.FFTW.Rank1 (flipSign)
+>>> import Numeric.FFTW.Shape
+>>>    (SymmetrySingleton(Even,Odd), ShiftSingleton(Exact,Halfway))
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import qualified Data.Array.Comfort.Shape as Shape
+>>> import qualified Data.Complex as Complex
+>>> import Data.Array.Comfort.Storable (Array)
+>>>
+>>> import qualified Numeric.Netlib.Class as Class
+>>>
+>>> import Foreign.Storable (Storable)
+>>>
+>>> import qualified Test.QuickCheck as QC
+>>>
+>>> genSymmetric ::
+>>>    (Spectrum.Symmetry symmetry) =>
+>>>    (Spectrum.Shift shiftTime) =>
+>>>    (Spectrum.Shift shiftSpectrum) =>
+>>>    (Spectrum.Symmetric symmetry shiftTime shiftSpectrum Int ~ sh) =>
+>>>    (QC.Arbitrary a, Storable a) =>
+>>>    Spectrum.SymmetrySingleton symmetry ->
+>>>    Spectrum.ShiftSingleton shiftTime ->
+>>>    Spectrum.ShiftSingleton shiftSpectrum ->
+>>>    QC.Gen (Array sh a)
+>>> genSymmetric symmetry shiftTime shiftSpectrum = do
+>>>    xs <- fmap (take 1000) $ QC.arbitrary
+>>>    return $ Array.fromList
+>>>       (Spectrum.Symmetric symmetry shiftTime shiftSpectrum $ length xs) xs
+>>>
+>>> adjustSymmetric ::
+>>>    (Spectrum.Symmetry symmetry) =>
+>>>    (Spectrum.Shift shiftTime) =>
+>>>    (Spectrum.Shift shiftSpectrum) =>
+>>>    (Spectrum.Symmetric symmetry shiftTime shiftSpectrum Int ~ sh) =>
+>>>    (Class.Floating a) =>
+>>>    Array sh a -> Array sh a
+>>> adjustSymmetric xs =
+>>>    let n = Spectrum.symmetricLogicalSize (Array.shape xs)
+>>>    in Array.map (fromIntegral n *) xs
+-}
+
+
+{- |
+>>> complexFloatList $ Trafo1.fourier Trafo1.Forward $ Array.fromList (Shape.Cyclic (5::Int)) [1,0,0,0,0]
+[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00]
+>>> complexFloatList $ Trafo1.fourier Trafo1.Forward $ Array.fromList (Shape.Cyclic (5::Int)) [0,1,0,0,0]
+ [1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95]
+>>> complexFloatList $ Trafo1.fourier Trafo1.Backward $ Array.fromList (Shape.Cyclic (5::Int)) [0,1,0,0,0]
+[1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95)]
+
+prop> QC.forAll genCyclicArray1 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo1.fourier sign (Trafo1.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo1.fourier sign (Trafo1.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray1 $ \xs sign -> approxComplex floatTol (Trafo1.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo1.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (Trafo1.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo1.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo1.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo1.fourier sign xs) (Trafo1.fourier sign ys))
+prop> QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo1.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo1.fourier sign xs) (Trafo1.fourier sign ys))
+
+prop> QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+prop> QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+
+prop> \sign -> QC.forAll genCyclicArray1 $ immutable (Trafo1.fourier sign) . arrayComplexFloat
+prop> \sign -> QC.forAll genCyclicArray1 $ immutable (Trafo1.fourier sign) . arrayComplexDouble
+-}
+fourier ::
+   (Integral n, Class.Real a) =>
+   Sign ->
+   Array (Shape.Cyclic n) (Complex a) ->
+   Array (Shape.Cyclic n) (Complex a)
+fourier sign (Array sh x) =
+   Array.unsafeCreateWithSize sh $ \n yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $
+   run $ FFI.planDFT1d (fromIntegral n) xPtr yPtr (ffiSign sign)
+            (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll genCyclicArray1 $ \xs -> approxReal floatTol (adjust xs) (Trafo1.fourierCR (Trafo1.fourierRC xs))
+prop> QC.forAll genCyclicArray1 $ \xs -> approxReal doubleTol (adjust xs) (Trafo1.fourierCR (Trafo1.fourierRC xs))
+
+prop> QC.forAll genCyclicArray1 $ immutable Trafo1.fourierRC . arrayFloat
+prop> QC.forAll genCyclicArray1 $ immutable Trafo1.fourierRC . arrayDouble
+-}
+fourierRC ::
+   (Integral n, Class.Real a) =>
+   Array (Shape.Cyclic n) a ->
+   Array (Spectrum.Half n) (Complex a)
+fourierRC (Array (Shape.Cyclic n) x) =
+   Array.unsafeCreate (Spectrum.Half n) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   if n<=0
+      then poke yPtr 0
+      else run $ FFI.planDFTr2c1d (fromIntegral n) xPtr yPtr
+                     (FFI.estimate <> FFI.preserveInput)
+
+{- |
+>>> floatList $ Trafo1.fourierCR $ Array.fromList (Spectrum.Half (5::Int)) [0,1,0]
+[2.00,0.62,-1.62,-1.62,0.62]
+
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable Trafo1.fourierCR . arrayComplexFloat
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable Trafo1.fourierCR . arrayComplexDouble
+-}
+fourierCR ::
+   (Integral n, Class.Real a) =>
+   Array (Spectrum.Half n) (Complex a) ->
+   Array (Shape.Cyclic n) a
+fourierCR (Array (Spectrum.Half n) x) =
+   Array.unsafeCreate (Shape.Cyclic n) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $
+   run $ FFI.planDFTc2r1d (fromIntegral n) xPtr yPtr
+            (FFI.estimate <> FFI.preserveInput)
+
+
+{- |
+@Symmetric Even Halfway Exact@ yields _the_ DCT,
+@Symmetric Even Exact Halfway@ yields _the_ inverse DCT.
+
+prop> QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Halfway Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Halfway Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Exact Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Exact Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Halfway Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+prop> QC.forAll (genSymmetric Even Halfway Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs))
+
+prop> QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ immutable Trafo1.cosine . arrayFloat
+prop> QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ immutable Trafo1.cosine . arrayDouble
+prop> QC.forAll (genSymmetric Even Halfway Exact) $ immutable Trafo1.cosine . arrayFloat
+prop> QC.forAll (genSymmetric Even Halfway Exact) $ immutable Trafo1.cosine . arrayDouble
+prop> QC.forAll (genSymmetric Even Exact Halfway) $ immutable Trafo1.cosine . arrayFloat
+prop> QC.forAll (genSymmetric Even Exact Halfway) $ immutable Trafo1.cosine . arrayDouble
+prop> QC.forAll (genSymmetric Even Halfway Halfway) $ immutable Trafo1.cosine . arrayFloat
+prop> QC.forAll (genSymmetric Even Halfway Halfway) $ immutable Trafo1.cosine . arrayDouble
+-}
+cosine ::
+   (Spectrum.Shift shiftTime, Spectrum.Shift shiftSpectrum,
+    Integral n, Class.Real a) =>
+   Array (Spectrum.Symmetric Spectrum.Even shiftTime shiftSpectrum n) a ->
+   Array (Spectrum.Symmetric Spectrum.Even shiftSpectrum shiftTime n) a
+cosine (Array (Spectrum.Symmetric symmetry shiftTime shiftSpectrum n) x) =
+   Array.unsafeCreate
+      (Spectrum.Symmetric symmetry shiftSpectrum shiftTime n) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $
+   let kind =
+         case (shiftTime,shiftSpectrum) of
+            (Spectrum.Halfway, Spectrum.Exact  ) -> FFI.redft10
+            (Spectrum.Exact,   Spectrum.Halfway) -> FFI.redft01
+            (Spectrum.Halfway, Spectrum.Halfway) -> FFI.redft11
+            (Spectrum.Exact,   Spectrum.Exact  ) ->
+               if n>1
+                  then FFI.redft00
+                  else error "DCT-1 must have at least two input data Exacts"
+   in run $ FFI.planR2r1d (fromIntegral n) xPtr yPtr kind
+               (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll (genSymmetric Odd Exact Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Exact Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Halfway Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Halfway Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Exact Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Exact Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Halfway Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+prop> QC.forAll (genSymmetric Odd Halfway Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs))
+
+prop> QC.forAll (genSymmetric Odd Exact Exact) $ immutable Trafo1.sine . arrayFloat
+prop> QC.forAll (genSymmetric Odd Exact Exact) $ immutable Trafo1.sine . arrayDouble
+prop> QC.forAll (genSymmetric Odd Halfway Exact) $ immutable Trafo1.sine . arrayFloat
+prop> QC.forAll (genSymmetric Odd Halfway Exact) $ immutable Trafo1.sine . arrayDouble
+prop> QC.forAll (genSymmetric Odd Exact Halfway) $ immutable Trafo1.sine . arrayFloat
+prop> QC.forAll (genSymmetric Odd Exact Halfway) $ immutable Trafo1.sine . arrayDouble
+prop> QC.forAll (genSymmetric Odd Halfway Halfway) $ immutable Trafo1.sine . arrayFloat
+prop> QC.forAll (genSymmetric Odd Halfway Halfway) $ immutable Trafo1.sine . arrayDouble
+-}
+sine ::
+   (Spectrum.Shift shiftTime, Spectrum.Shift shiftSpectrum,
+    Integral n, Class.Real a) =>
+   Array (Spectrum.Symmetric Spectrum.Odd shiftTime shiftSpectrum n) a ->
+   Array (Spectrum.Symmetric Spectrum.Odd shiftSpectrum shiftTime n) a
+sine (Array (Spectrum.Symmetric symmetry shiftTime shiftSpectrum n) x) =
+   Array.unsafeCreate
+      (Spectrum.Symmetric symmetry shiftSpectrum shiftTime n) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $
+   let kind =
+         case (shiftTime,shiftSpectrum) of
+            (Spectrum.Exact,   Spectrum.Exact  ) -> FFI.rodft00
+            (Spectrum.Halfway, Spectrum.Exact  ) -> FFI.rodft10
+            (Spectrum.Exact,   Spectrum.Halfway) -> FFI.rodft01
+            (Spectrum.Halfway, Spectrum.Halfway) -> FFI.rodft11
+   in run $ FFI.planR2r1d (fromIntegral n) xPtr yPtr kind
+               (FFI.estimate <> FFI.preserveInput)
+
+
+{- |
+prop> QC.forAll genCyclicArray1 $ \xs -> approxReal floatTol (adjust xs) (Trafo1.hartley (Trafo1.hartley xs))
+prop> QC.forAll genCyclicArray1 $ \xs -> approxReal doubleTol (adjust xs) (Trafo1.hartley (Trafo1.hartley xs))
+
+prop> QC.forAll genCyclicArray1 $ immutable Trafo1.hartley . arrayFloat
+prop> QC.forAll genCyclicArray1 $ immutable Trafo1.hartley . arrayDouble
+-}
+hartley ::
+   (Integral n, Class.Real a) =>
+   Array (Shape.Cyclic n) a ->
+   Array (Shape.Cyclic n) a
+hartley (Array sh@(Shape.Cyclic n) x) =
+   Array.unsafeCreate sh $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $
+   run $ FFI.planR2r1d (fromIntegral n) xPtr yPtr FFI.dht
+            (FFI.estimate <> FFI.preserveInput)
diff --git a/src/Numeric/FFTW/Rank2.hs b/src/Numeric/FFTW/Rank2.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Rank2.hs
@@ -0,0 +1,102 @@
+module Numeric.FFTW.Rank2 (
+   fourier, Sign(..), flipSign,
+   fourierRC,
+   fourierCR,
+   ) where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+import Numeric.FFTW.Private (Sign(..), flipSign, ffiSign, run, runCopiedArray)
+
+import Foreign.Marshal.Array (pokeArray)
+import Foreign.ForeignPtr (withForeignPtr)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import Data.Complex (Complex)
+import Data.Monoid ((<>))
+
+import Control.Monad (when)
+
+
+{- $setup
+>>> import Test.Numeric.FFTW.Common
+>>>    (approxReal, approxComplex, floatTol, doubleTol, normInf,
+>>>     adjust, scalarProduct, split, genCyclicArray2, immutable,
+>>>     arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble)
+>>>
+>>> import qualified Numeric.FFTW.Rank2 as Trafo2
+>>> import Numeric.FFTW.Rank1 (flipSign)
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import qualified Data.Complex as Complex
+>>>
+>>> import qualified Test.QuickCheck as QC
+-}
+
+
+{- |
+prop> QC.forAll genCyclicArray2 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo2.fourier sign (Trafo2.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo2.fourier sign (Trafo2.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray2 $ \xs sign -> approxComplex floatTol (Trafo2.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo2.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (Trafo2.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo2.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo2.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo2.fourier sign xs) (Trafo2.fourier sign ys))
+prop> QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo2.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo2.fourier sign xs) (Trafo2.fourier sign ys))
+
+prop> QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+prop> QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+
+prop> \sign -> QC.forAll genCyclicArray2 $ immutable (Trafo2.fourier sign) . arrayComplexFloat
+prop> \sign -> QC.forAll genCyclicArray2 $ immutable (Trafo2.fourier sign) . arrayComplexDouble
+-}
+fourier ::
+   (Integral n0, Integral n1, Class.Real a) =>
+   Sign ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1) (Complex a) ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1) (Complex a)
+fourier sign (Array sh@(Shape.Cyclic n0, Shape.Cyclic n1) x) =
+   Array.unsafeCreateWithSize sh $ \n yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $ run $
+      FFI.planDFT2d (fromIntegral n0) (fromIntegral n1)
+         xPtr yPtr (ffiSign sign) (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll genCyclicArray2 $ \xs -> approxReal floatTol (adjust xs) (Trafo2.fourierCR (Trafo2.fourierRC xs))
+prop> QC.forAll genCyclicArray2 $ \xs -> approxReal doubleTol (adjust xs) (Trafo2.fourierCR (Trafo2.fourierRC xs))
+
+prop> QC.forAll genCyclicArray2 $ immutable Trafo2.fourierRC . arrayFloat
+prop> QC.forAll genCyclicArray2 $ immutable Trafo2.fourierRC . arrayDouble
+-}
+fourierRC ::
+   (Integral n0, Integral n1, Class.Real a) =>
+   Array (Shape.Cyclic n0, Shape.Cyclic n1) a ->
+   Array (Shape.Cyclic n0, Spectrum.Half n1) (Complex a)
+fourierRC (Array (sh0@(Shape.Cyclic n0), Shape.Cyclic n1) x) =
+   Array.unsafeCreate (sh0, Spectrum.Half n1) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   if n0<=0 || n1<=0
+      then pokeArray yPtr $ replicate (fromIntegral n0) 0
+      else run $
+         FFI.planDFTr2c2d (fromIntegral n0) (fromIntegral n1)
+            xPtr yPtr (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Trafo2.fourierCR . arrayComplexFloat
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Trafo2.fourierCR . arrayComplexDouble
+-}
+fourierCR ::
+   (Integral n0, Integral n1, Class.Real a) =>
+   Array (Shape.Cyclic n0, Spectrum.Half n1) (Complex a) ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1) a
+fourierCR arr@(Array (sh0@(Shape.Cyclic n0), Spectrum.Half n1) _x) =
+   Array.unsafeCreate (sh0, Shape.Cyclic n1) $ \yPtr ->
+   when (n0>0 && n1>0) $
+   runCopiedArray arr $ \xPtr ->
+      FFI.planDFTc2r2d (fromIntegral n0) (fromIntegral n1)
+         xPtr yPtr (FFI.estimate <> FFI.destroyInput)
diff --git a/src/Numeric/FFTW/Rank3.hs b/src/Numeric/FFTW/Rank3.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Rank3.hs
@@ -0,0 +1,108 @@
+module Numeric.FFTW.Rank3 (
+   fourier, Sign(..), flipSign,
+   fourierRC,
+   fourierCR,
+   ) where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+import Numeric.FFTW.Private (Sign(..), flipSign, ffiSign, run, runCopiedArray)
+
+import Foreign.Marshal.Array (pokeArray)
+import Foreign.ForeignPtr (withForeignPtr)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import Data.Complex (Complex)
+import Data.Monoid ((<>))
+
+import Control.Monad (when)
+
+
+{- $setup
+>>> import Test.Numeric.FFTW.Common
+>>>    (approxReal, approxComplex, floatTol, doubleTol, normInf,
+>>>     adjust, scalarProduct, split, genCyclicArray3, immutable,
+>>>     arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble)
+>>>
+>>> import qualified Numeric.FFTW.Rank3 as Trafo3
+>>> import Numeric.FFTW.Rank1 (flipSign)
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import qualified Data.Complex as Complex
+>>>
+>>> import qualified Test.QuickCheck as QC
+-}
+
+
+{- |
+prop> QC.forAll genCyclicArray3 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo3.fourier sign (Trafo3.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo3.fourier sign (Trafo3.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray3 $ \xs sign -> approxComplex floatTol (Trafo3.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo3.fourier (flipSign sign) xs))
+prop> QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (Trafo3.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo3.fourier (flipSign sign) xs))
+
+prop> QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo3.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo3.fourier sign xs) (Trafo3.fourier sign ys))
+prop> QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo3.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo3.fourier sign xs) (Trafo3.fourier sign ys))
+
+prop> QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+prop> QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys
+
+prop> \sign -> QC.forAll genCyclicArray3 $ immutable (Trafo3.fourier sign) . arrayComplexFloat
+prop> \sign -> QC.forAll genCyclicArray3 $ immutable (Trafo3.fourier sign) . arrayComplexDouble
+-}
+fourier ::
+   (Integral n0, Integral n1, Integral n2, Class.Real a) =>
+   Sign ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2) (Complex a) ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2) (Complex a)
+fourier sign (Array sh@(Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2) x) =
+   Array.unsafeCreateWithSize sh $ \n yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $ run $
+      FFI.planDFT3d (fromIntegral n0) (fromIntegral n1) (fromIntegral n2)
+         xPtr yPtr (ffiSign sign) (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll genCyclicArray3 $ \xs -> approxReal floatTol (adjust xs) (Trafo3.fourierCR (Trafo3.fourierRC xs))
+prop> QC.forAll genCyclicArray3 $ \xs -> approxReal doubleTol (adjust xs) (Trafo3.fourierCR (Trafo3.fourierRC xs))
+
+prop> QC.forAll genCyclicArray3 $ immutable Trafo3.fourierRC . arrayFloat
+prop> QC.forAll genCyclicArray3 $ immutable Trafo3.fourierRC . arrayDouble
+-}
+fourierRC ::
+   (Integral n0, Integral n1, Integral n2, Class.Real a) =>
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2) a ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Spectrum.Half n2) (Complex a)
+fourierRC
+   (Array (sh0@(Shape.Cyclic n0), sh1@(Shape.Cyclic n1), Shape.Cyclic n2) x) =
+
+   Array.unsafeCreate (sh0, sh1, Spectrum.Half n2) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   if n0<=0 || n1<=0 || n2<=0
+      then pokeArray yPtr $ replicate (fromIntegral n0 * fromIntegral n1) 0
+      else run $
+         FFI.planDFTr2c3d (fromIntegral n0) (fromIntegral n1) (fromIntegral n2)
+            xPtr yPtr (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Trafo3.fourierCR . arrayComplexFloat
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Trafo3.fourierCR . arrayComplexDouble
+-}
+fourierCR ::
+   (Integral n0, Integral n1, Integral n2, Class.Real a) =>
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Spectrum.Half n2) (Complex a) ->
+   Array (Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2) a
+fourierCR
+   arr@(Array
+      (sh0@(Shape.Cyclic n0), sh1@(Shape.Cyclic n1), Spectrum.Half n2)
+      _x) =
+
+   Array.unsafeCreate (sh0, sh1, Shape.Cyclic n2) $ \yPtr ->
+   when (n0>0 && n1>0 && n2>0) $
+   runCopiedArray arr $ \xPtr ->
+      FFI.planDFTc2r3d (fromIntegral n0) (fromIntegral n1) (fromIntegral n2)
+         xPtr yPtr (FFI.estimate <> FFI.destroyInput)
diff --git a/src/Numeric/FFTW/RankN.hs b/src/Numeric/FFTW/RankN.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/RankN.hs
@@ -0,0 +1,155 @@
+module Numeric.FFTW.RankN (
+   fourier, Sign(..), flipSign,
+   fourierRC,
+   fourierCR,
+   ) where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.FFI as FFI
+import qualified Numeric.Netlib.Class as Class
+import Numeric.FFTW.Private
+         (Sign(..), flipSign, ffiSign, run, runCopiedArray, withDims)
+
+import Foreign.Marshal.Array (pokeArray)
+import Foreign.ForeignPtr (withForeignPtr)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import qualified Data.List as List
+import Data.Complex (Complex)
+import Data.Monoid ((<>))
+
+import Control.Monad (when)
+
+
+{- $setup
+>>> import Test.Numeric.FFTW.Common
+>>>    (approxReal, approxComplex, floatTol, doubleTol,
+>>>     genCyclicArray1, genCyclicArray2, genCyclicArray3, immutable,
+>>>     arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+>>>     floatList, complexFloatList)
+>>>
+>>> import qualified Numeric.FFTW.RankN as TrafoM
+>>> import qualified Numeric.FFTW.Rank1 as Trafo1
+>>> import qualified Numeric.FFTW.Rank2 as Trafo2
+>>> import qualified Numeric.FFTW.Rank3 as Trafo3
+>>> import qualified Numeric.FFTW.Shape as Spectrum
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import qualified Data.Array.Comfort.Shape as Shape
+>>> import Data.Array.Comfort.Storable (Array)
+>>>
+>>> import qualified Test.QuickCheck as QC
+>>>
+>>> array1 :: (Shape.C sh) => Array sh a -> Array ((), sh) a
+>>> array1 = Array.mapShape ((,) ())
+>>>
+>>> array3 :: (Shape.C sh0, Shape.C sh1, Shape.C sh2) =>
+>>>    Array (sh0,sh1,sh2) a -> Array ((sh0,sh1),sh2) a
+>>> array3 = Array.mapShape (\(sh0,sh1,sh2) -> ((sh0,sh1),sh2))
+-}
+
+
+{- |
+>>> complexFloatList $ TrafoM.fourier TrafoM.Forward $ Array.fromList (Shape.Cyclic (5::Int), Shape.Cyclic (5::Int)) [0,0,0,0,0, 0,1,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0]
+[1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,0.31:+(-0.95),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),0.31:+0.95,1.00:+(-0.00),0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59]
+
+prop> QC.forAll genCyclicArray1 $ \xs sign -> approxComplex floatTol (array1 $ Trafo1.fourier sign xs) (TrafoM.fourier sign $ array1 xs)
+prop> QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (array1 $ Trafo1.fourier sign xs) (TrafoM.fourier sign $ array1 xs)
+
+prop> QC.forAll genCyclicArray2 $ \xs sign -> approxComplex floatTol (Trafo2.fourier sign xs) (TrafoM.fourier sign xs)
+prop> QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (Trafo2.fourier sign xs) (TrafoM.fourier sign xs)
+
+prop> QC.forAll genCyclicArray3 $ \xs sign -> approxComplex floatTol (array3 $ Trafo3.fourier sign xs) (TrafoM.fourier sign $ array3 xs)
+prop> QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (array3 $ Trafo3.fourier sign xs) (TrafoM.fourier sign $ array3 xs)
+-}
+fourier ::
+   (Spectrum.MultiCyclic sh, Class.Real a) =>
+   Sign -> Array sh (Complex a) -> Array sh (Complex a)
+fourier sign (Array sh x) =
+   Array.unsafeCreateWithSize sh $ \n yPtr ->
+   withForeignPtr x $ \xPtr ->
+   when (n>0) $ run $
+   withDims (Spectrum.cyclicDimensions sh) $ \rank dimPtr ->
+      FFI.planDFT rank dimPtr xPtr yPtr (ffiSign sign)
+         (FFI.estimate <> FFI.preserveInput)
+
+{- |
+prop> QC.forAll genCyclicArray1 $ \xs -> approxComplex floatTol (array1 $ Trafo1.fourierRC xs) (TrafoM.fourierRC $ array1 xs)
+prop> QC.forAll genCyclicArray1 $ \xs -> approxComplex doubleTol (array1 $ Trafo1.fourierRC xs) (TrafoM.fourierRC $ array1 xs)
+
+prop> QC.forAll genCyclicArray2 $ \xs -> approxComplex floatTol (Trafo2.fourierRC xs) (TrafoM.fourierRC xs)
+prop> QC.forAll genCyclicArray2 $ \xs -> approxComplex doubleTol (Trafo2.fourierRC xs) (TrafoM.fourierRC xs)
+
+prop> QC.forAll genCyclicArray3 $ \xs -> approxComplex floatTol (array3 $ Trafo3.fourierRC xs) (TrafoM.fourierRC $ array3 xs)
+prop> QC.forAll genCyclicArray3 $ \xs -> approxComplex doubleTol (array3 $ Trafo3.fourierRC xs) (TrafoM.fourierRC $ array3 xs)
+
+
+prop> QC.forAll genCyclicArray1 $ immutable TrafoM.fourierRC . array1 . arrayFloat
+prop> QC.forAll genCyclicArray1 $ immutable TrafoM.fourierRC . array1 . arrayDouble
+
+prop> QC.forAll genCyclicArray2 $ immutable TrafoM.fourierRC . arrayFloat
+prop> QC.forAll genCyclicArray2 $ immutable TrafoM.fourierRC . arrayDouble
+
+prop> QC.forAll genCyclicArray3 $ immutable TrafoM.fourierRC . array3 . arrayFloat
+prop> QC.forAll genCyclicArray3 $ immutable TrafoM.fourierRC . array3 . arrayDouble
+-}
+fourierRC ::
+   (Spectrum.MultiCyclic sh0, Integral n1, Class.Real a) =>
+   Array (sh0, Shape.Cyclic n1) a ->
+   Array (sh0, Spectrum.Half n1) (Complex a)
+fourierRC (Array sh@(sh0, Shape.Cyclic n1) x) =
+   Array.unsafeCreate (sh0, Spectrum.Half n1) $ \yPtr ->
+   withForeignPtr x $ \xPtr ->
+   let dims = Spectrum.cyclicDimensions sh in
+   if any (<=0) dims
+      then pokeArray yPtr $ List.genericReplicate (Shape.size sh0) 0
+      else run $
+         withDims dims $ \rank dimPtr ->
+         FFI.planDFTr2c rank dimPtr xPtr yPtr
+            (FFI.estimate <> FFI.preserveInput)
+
+{- |
+>>> floatList $ Trafo1.fourierCR $ Array.fromList (Spectrum.Half (3::Int)) [1,0 ]
+[1.00,1.00,1.00]
+>>> floatList $ TrafoM.fourierCR $ Array.fromList ((), Spectrum.Half (3::Int)) [1,0]
+[1.00,1.00,1.00]
+
+>>> floatList $ Trafo2.fourierCR $ Array.fromList (Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0]
+[1.00,1.00,1.00]
+>>> floatList $ TrafoM.fourierCR $ Array.fromList (Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0]
+[1.00,1.00,1.00]
+
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ \xs -> approxReal floatTol (array1 $ Trafo1.fourierCR xs) (TrafoM.fourierCR $ array1 xs)
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ \xs -> approxReal doubleTol (array1 $ Trafo1.fourierCR xs) (TrafoM.fourierCR $ array1 xs)
+
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> approxReal floatTol (Trafo2.fourierCR xs) (TrafoM.fourierCR xs)
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> approxReal doubleTol (Trafo2.fourierCR xs) (TrafoM.fourierCR xs)
+
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> approxReal floatTol (array3 $ Trafo3.fourierCR xs) (TrafoM.fourierCR $ array3 xs)
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> approxReal doubleTol (array3 $ Trafo3.fourierCR xs) (TrafoM.fourierCR $ array3 xs)
+
+
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable TrafoM.fourierCR . array1 . arrayComplexFloat
+prop> QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable TrafoM.fourierCR . array1 . arrayComplexDouble
+
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable TrafoM.fourierCR . arrayComplexFloat
+prop> QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable TrafoM.fourierCR . arrayComplexDouble
+
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable TrafoM.fourierCR . array3 . arrayComplexFloat
+prop> QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable TrafoM.fourierCR . array3 . arrayComplexDouble
+-}
+fourierCR ::
+   (Spectrum.MultiCyclic sh0, Integral n1, Class.Real a) =>
+   Array (sh0, Spectrum.Half n1) (Complex a) ->
+   Array (sh0, Shape.Cyclic n1) a
+fourierCR arr@(Array (sh0, Spectrum.Half n1) _x) =
+   let sh = (sh0, Shape.Cyclic n1) in
+   let dims = Spectrum.cyclicDimensions sh in
+   Array.unsafeCreate sh $ \yPtr ->
+   when (all (>0) dims) $
+   runCopiedArray arr $ \xPtr ->
+   withDims dims $ \rank dimPtr ->
+      FFI.planDFTc2r rank dimPtr xPtr yPtr (FFI.estimate <> FFI.destroyInput)
diff --git a/src/Numeric/FFTW/Shape.hs b/src/Numeric/FFTW/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/FFTW/Shape.hs
@@ -0,0 +1,199 @@
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+module Numeric.FFTW.Shape (
+   Half(..),
+   MultiCyclic(cyclicDimensions),
+   Symmetric(..), symmetric,
+   symmetricLogicalSize,
+   Symmetry(..), SymmetrySingleton(..), Even, Odd,
+   Shift(..), ShiftSingleton(..), Exact, Halfway,
+   ) where
+
+import qualified Data.Array.Comfort.Shape as Shape
+
+import qualified Foreign.C.Types as C
+
+import Control.DeepSeq (NFData, rnf)
+
+
+newtype Half n = Half n
+   deriving (Eq, Show)
+
+instance (NFData n) => NFData (Half n) where
+   rnf (Half n) = rnf n
+
+instance (Integral n) => Shape.C (Half n) where
+   size (Half n) = halfSize $ fromIntegral n
+
+instance (Integral n) => Shape.Indexed (Half n) where
+   type Index (Half n) = n
+   indices (Half len) = Shape.indices $ Shape.ZeroBased $ halfSize len
+   unifiedOffset (Half len) =
+      Shape.unifiedOffset $ Shape.ZeroBased $ halfSize len
+   inBounds (Half len) ix = 0<=ix && ix<halfSize len
+
+instance (Integral n) => Shape.InvIndexed (Half n) where
+   unifiedIndexFromOffset (Half len) k0 = do
+      let k = fromIntegral k0
+      Shape.assertIndexFromOffset "Half" k0 $ 0<=k && k<halfSize len
+      return k
+
+halfSize :: Integral a => a -> a
+halfSize n = div n 2 + 1
+
+
+
+class (Shape.C sh) => MultiCyclic sh where
+   cyclicDimensions :: sh -> [C.CInt]
+
+instance (Integral n) => MultiCyclic (Shape.Cyclic n) where
+   cyclicDimensions (Shape.Cyclic n) = [fromIntegral n]
+
+instance MultiCyclic () where
+   cyclicDimensions () = []
+
+instance (MultiCyclic sh0, MultiCyclic sh1) => MultiCyclic (sh0,sh1) where
+   cyclicDimensions (sh0,sh1) = cyclicDimensions sh0 ++ cyclicDimensions sh1
+
+instance
+   (MultiCyclic sh0, MultiCyclic sh1, MultiCyclic sh2) =>
+      MultiCyclic (sh0,sh1,sh2) where
+   cyclicDimensions (sh0,sh1,sh2) =
+      cyclicDimensions sh0 ++ cyclicDimensions sh1 ++ cyclicDimensions sh2
+
+
+data Even
+data Odd
+data SymmetrySingleton symm where
+   Even :: SymmetrySingleton Even
+   Odd  :: SymmetrySingleton Odd
+
+class Symmetry symm where switchSymmetry :: f Even -> f Odd -> f symm
+instance Symmetry Even where switchSymmetry f _ = f
+instance Symmetry Odd  where switchSymmetry _ f = f
+
+autoSymmetry :: (Symmetry symm) => SymmetrySingleton symm
+autoSymmetry = switchSymmetry Even Odd
+
+instance Eq (SymmetrySingleton symm) where
+   x==y =
+      case (x,y) of
+         (Even, Even) -> True
+         (Odd, Odd) -> True
+
+instance Show (SymmetrySingleton symm) where
+   show Even = "Even"
+   show Odd = "Odd"
+
+instance NFData (SymmetrySingleton symm) where
+   rnf s = case s of Even -> (); Odd -> ()
+
+data Exact
+data Halfway
+data ShiftSingleton shift where
+   Exact   :: ShiftSingleton Exact
+   Halfway :: ShiftSingleton Halfway
+
+class Shift shift where switchShift :: f Exact -> f Halfway -> f shift
+instance Shift Exact   where switchShift f _ = f
+instance Shift Halfway where switchShift _ f = f
+
+autoShift :: (Shift shift) => ShiftSingleton shift
+autoShift = switchShift Exact Halfway
+
+instance Eq (ShiftSingleton shift) where
+   x==y =
+      case (x,y) of
+         (Exact, Exact) -> True
+         (Halfway, Halfway) -> True
+
+instance Show (ShiftSingleton shift) where
+   show Exact = "Exact"
+   show Halfway = "Halfway"
+
+instance NFData (ShiftSingleton shift) where
+   rnf s = case s of Exact -> (); Halfway -> ()
+
+
+{- |
+Shape for stored data of symmetric vectors.
+Even is for Cosine transform, Odd for Sine transform.
+@shiftTime@ refers to no or halfway shift of the data,
+@shiftSpectrum@ refers to no or halfway shift of the Cosine or Sine spectrum.
+
+0 means Exact, 1 means Halfway:
+
+* Even 0 0: even around 0 and even around n-1.
+
+* Even 1 0: even around -0.5 and even around n-0.5.
+
+* Even 0 1: even around 0 and odd around n.
+
+* Even 1 1: even around -0.5 and odd around n-0.5.
+
+* Odd  0 0: odd around -1 and odd around n.
+
+* Odd  1 0: odd around -0.5 and odd around n-0.5.
+
+* Odd  0 1: odd around -1 and even around n-1.
+
+* Odd  1 1: odd around -0.5 and even around n-0.5.
+
+
+We could pad data of Even symmetric vectors,
+but we cannot pad data of Odd symmetric vectors,
+because '!' would have to involve 'negate'.
+Thus we provide no padding, at all.
+-}
+data Symmetric symmetry shiftTime shiftSpectrum n =
+   Symmetric
+      (SymmetrySingleton symmetry)
+      (ShiftSingleton shiftTime)
+      (ShiftSingleton shiftSpectrum)
+      n
+   deriving (Eq, Show)
+
+symmetric ::
+   (Symmetry symmetry, Shift shiftTime, Shift shiftSpectrum) =>
+   n -> Symmetric symmetry shiftTime shiftSpectrum n
+symmetric = Symmetric autoSymmetry autoShift autoShift
+
+symmetricLogicalSize ::
+   (Num n) => Symmetric symmetry shiftTime shiftSpectrum n -> n
+symmetricLogicalSize (Symmetric symmetry shiftTime shiftSpectrum n) =
+   case (shiftTime, shiftSpectrum) of
+      (Exact, Exact) ->
+         case symmetry of
+            Even -> 2*n-2
+            Odd  -> 2*n+2
+      _ -> 2*n
+
+instance
+   (Symmetry symmetry, Shift shiftTime, Shift shiftSpectrum, NFData n) =>
+      NFData (Symmetric symmetry shiftTime shiftSpectrum n) where
+   rnf (Symmetric symmetry shiftTime shiftSpectrum n) =
+      rnf (symmetry, shiftTime, shiftSpectrum, n)
+
+instance
+   (Symmetry symmetry, Shift shiftTime, Shift shiftSpectrum, Integral n) =>
+      Shape.C (Symmetric symmetry shiftTime shiftSpectrum n) where
+   size (Symmetric _ _ _ n) = fromIntegral n
+
+instance
+   (Symmetry symmetry, Shift shiftTime, Shift shiftSpectrum, Integral n) =>
+      Shape.Indexed (Symmetric symmetry shiftTime shiftSpectrum n) where
+   type Index (Symmetric symmetry shiftTime shiftSpectrum n) = n
+   indices (Symmetric _ _ _ n) =
+      Shape.indices $ Shape.ZeroBased $ fromIntegral n
+   unifiedOffset (Symmetric _ _ _ n) =
+      Shape.unifiedOffset $ Shape.ZeroBased $ fromIntegral n
+   inBounds (Symmetric _ _ _ n) ix = 0<=ix && ix<fromIntegral n
+
+instance
+   (Symmetry symmetry, Shift shiftTime, Shift shiftSpectrum, Integral n) =>
+      Shape.InvIndexed (Symmetric symmetry shiftTime shiftSpectrum n) where
+   unifiedIndexFromOffset (Symmetric _ _ _ n) k0 = do
+      let k = fromIntegral k0
+      Shape.assertIndexFromOffset "Symmetric" k0 $ 0<=k && k<fromIntegral n
+      return k
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,15 @@
+module Main where
+
+import qualified Test.Numeric.FFTW.Shape as TestShape
+import qualified Test.Main as TestMain
+import Test.Numeric.FFTW.Common (prefix)
+
+import qualified Test.DocTest.Driver as DocTest
+
+
+main :: IO ()
+main = DocTest.run $ (>> TestMain.main) $
+   mapM_ (\(name,prop) ->
+            DocTest.printPrefix (name ++ ": ") >> DocTest.property prop) $
+   prefix "Shape" TestShape.tests ++
+   []
diff --git a/test/Test/Main.hs b/test/Test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Main.hs
@@ -0,0 +1,18 @@
+-- Do not edit! Automatically created with doctest-extract.
+module Test.Main where
+
+import qualified Test.Numeric.FFTW.Rank1
+import qualified Test.Numeric.FFTW.Rank2
+import qualified Test.Numeric.FFTW.Rank3
+import qualified Test.Numeric.FFTW.RankN
+import qualified Test.Numeric.FFTW.Batch
+
+import qualified Test.DocTest.Driver as DocTest
+
+main :: DocTest.T ()
+main = do
+    Test.Numeric.FFTW.Rank1.test
+    Test.Numeric.FFTW.Rank2.test
+    Test.Numeric.FFTW.Rank3.test
+    Test.Numeric.FFTW.RankN.test
+    Test.Numeric.FFTW.Batch.test
diff --git a/test/Test/Numeric/FFTW/Batch.hs b/test/Test/Numeric/FFTW/Batch.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Batch.hs
@@ -0,0 +1,214 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/FFTW/Batch.hs
+{-# LINE 28 "src/Numeric/FFTW/Batch.hs" #-}
+
+module Test.Numeric.FFTW.Batch where
+
+import Test.DocTest.Base
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 29 "src/Numeric/FFTW/Batch.hs" #-}
+import     Test.Numeric.FFTW.Common
+       (approxReal, approxComplex, floatTol, doubleTol,
+        genCyclicArray2, genCyclicArray3, immutable,
+        arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+        floatList, complexFloatList)
+
+import     qualified Numeric.FFTW.Batch as Batch
+import     qualified Numeric.FFTW.Rank1 as Trafo1
+import     qualified Numeric.FFTW.Rank2 as Trafo2
+import     qualified Numeric.FFTW.Rank3 as Trafo3
+import     qualified Numeric.FFTW.Shape as Spectrum
+
+import     qualified Numeric.Netlib.Class as Class
+
+import     qualified Data.Array.Comfort.Boxed as BoxedArray
+import     qualified Data.Array.Comfort.Storable as Array
+import     qualified Data.Array.Comfort.Shape as Shape
+import     qualified Data.Foldable as Fold
+import     Data.Array.Comfort.Storable (Array)
+import     Data.Complex (Complex)
+
+import     qualified Test.QuickCheck as QC
+
+array1s     :: (Shape.C sh0, Shape.C sh1) =>
+       Array (sh0,sh1) a -> Array (sh0,(),sh1) a
+array1s     = Array.mapShape (\(sh0,sh1) -> (sh0,(),sh1))
+
+unarray1s     :: (Shape.C sh0, Shape.C sh1) =>
+       Array (sh0,(),sh1) a -> Array (sh0,sh1) a
+unarray1s     = Array.mapShape (\(sh0,(),sh1) -> (sh0,sh1))
+
+array2s     :: (Shape.C sh0, Shape.C sh1, Shape.C sh2) =>
+       Array (sh0,sh1,sh2) a -> Array (sh0,(sh1,sh2)) a
+array2s     = Array.mapShape (\(sh0,sh1,sh2) -> (sh0,(sh1,sh2)))
+
+allApproxReal     ::
+       (Shape.C sh0, Shape.C sh1, Class.Real a, Eq sh0, Eq sh1) =>
+       a ->
+       BoxedArray.Array sh0 (Array sh1 a) ->
+       BoxedArray.Array sh0 (Array sh1 a) ->
+       Bool
+allApproxReal     tol xs ys =
+       Fold.and $ BoxedArray.zipWith (approxReal tol) xs ys
+
+allApproxComplex     ::
+       (Shape.C sh0, Shape.C sh1, Class.Real a, Eq sh0, Eq sh1) =>
+       a ->
+       BoxedArray.Array sh0 (Array sh1 (Complex a)) ->
+       BoxedArray.Array sh0 (Array sh1 (Complex a)) ->
+       Bool
+allApproxComplex     tol xs ys =
+       Fold.and $ BoxedArray.zipWith (approxComplex tol) xs ys
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.FFTW.Batch:90: "
+{-# LINE 90 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 90 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> allApproxComplex floatTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo1.fourier sign) $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:91: "
+{-# LINE 91 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 91 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> allApproxComplex doubleTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo1.fourier sign) $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:93: "
+{-# LINE 93 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 93 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap array2s genCyclicArray3) $ \xs sign -> allApproxComplex floatTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo2.fourier sign) $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:94: "
+{-# LINE 94 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 94 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap array2s genCyclicArray3) $ \xs sign -> allApproxComplex doubleTol (Array.toRowArray $ Batch.fourier sign xs) (fmap (Trafo2.fourier sign) $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:96: "
+{-# LINE 96 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 96 "src/Numeric/FFTW/Batch.hs" #-}
+     (\sign -> QC.forAll genCyclicArray2 $ immutable (Batch.fourier sign) . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:97: "
+{-# LINE 97 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 97 "src/Numeric/FFTW/Batch.hs" #-}
+     (\sign -> QC.forAll genCyclicArray2 $ immutable (Batch.fourier sign) . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:99: "
+{-# LINE 99 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 99 "src/Numeric/FFTW/Batch.hs" #-}
+     (\sign -> QC.forAll genCyclicArray3 $ immutable (Batch.fourier sign) . array2s . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:100: "
+{-# LINE 100 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 100 "src/Numeric/FFTW/Batch.hs" #-}
+     (\sign -> QC.forAll genCyclicArray3 $ immutable (Batch.fourier sign) . array2s . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:85: "
+{-# LINE 85 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.example
+{-# LINE 85 "src/Numeric/FFTW/Batch.hs" #-}
+   (complexFloatList $ Batch.fourier Batch.Forward $ Array.fromList (Shape.ZeroBased (5::Int), Shape.Cyclic (5::Int)) [1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1])
+  [ExpectedLine [LineChunk "[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,(-0.81):+(-0.59),0.31:+0.95,0.31:+(-0.95),(-0.81):+0.59,1.00:+0.00,(-0.81):+0.59,0.31:+(-0.95),0.31:+0.95,(-0.81):+(-0.59),1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95)]"]]
+ DocTest.printPrefix "Numeric.FFTW.Batch:87: "
+{-# LINE 87 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.example
+{-# LINE 87 "src/Numeric/FFTW/Batch.hs" #-}
+   (complexFloatList $ Batch.fourier Batch.Backward $ Array.fromList (Shape.ZeroBased (5::Int), Shape.Cyclic (5::Int)) [1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1])
+  [ExpectedLine [LineChunk "[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95),1.00:+0.00,(-0.81):+0.59,0.31:+(-0.95),0.31:+0.95,(-0.81):+(-0.59),1.00:+0.00,(-0.81):+(-0.59),0.31:+0.95,0.31:+(-0.95),(-0.81):+0.59,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95]"]]
+ DocTest.printPrefix "Numeric.FFTW.Batch:117: "
+{-# LINE 117 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 117 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> allApproxComplex floatTol (Array.toRowArray $ unarray1s $ Batch.fourierRC $ array1s xs) (fmap Trafo1.fourierRC $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:118: "
+{-# LINE 118 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 118 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> allApproxComplex doubleTol (Array.toRowArray $ unarray1s $ Batch.fourierRC $ array1s xs) (fmap Trafo1.fourierRC $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:120: "
+{-# LINE 120 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 120 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> allApproxComplex floatTol (Array.toRowArray $ array2s $ Batch.fourierRC xs) (fmap Trafo2.fourierRC $ Array.toRowArray $ array2s xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:121: "
+{-# LINE 121 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 121 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> allApproxComplex doubleTol (Array.toRowArray $ array2s $ Batch.fourierRC xs) (fmap Trafo2.fourierRC $ Array.toRowArray $ array2s xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:123: "
+{-# LINE 123 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 123 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable Batch.fourierRC . array1s . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:124: "
+{-# LINE 124 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 124 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable Batch.fourierRC . array1s . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:126: "
+{-# LINE 126 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 126 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable Batch.fourierRC . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:127: "
+{-# LINE 127 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 127 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable Batch.fourierRC . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:156: "
+{-# LINE 156 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 156 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> allApproxReal floatTol (Array.toRowArray $ unarray1s $ Batch.fourierCR $ array1s xs) (fmap Trafo1.fourierCR $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:157: "
+{-# LINE 157 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 157 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> allApproxReal doubleTol (Array.toRowArray $ unarray1s $ Batch.fourierCR $ array1s xs) (fmap Trafo1.fourierCR $ Array.toRowArray xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:159: "
+{-# LINE 159 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 159 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> allApproxReal floatTol (Array.toRowArray $ array2s $ Batch.fourierCR xs) (fmap Trafo2.fourierCR $ Array.toRowArray $ array2s xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:160: "
+{-# LINE 160 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 160 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> allApproxReal doubleTol (Array.toRowArray $ array2s $ Batch.fourierCR xs) (fmap Trafo2.fourierCR $ Array.toRowArray $ array2s xs))
+ DocTest.printPrefix "Numeric.FFTW.Batch:163: "
+{-# LINE 163 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 163 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Batch.fourierCR . array1s . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:164: "
+{-# LINE 164 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 164 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Batch.fourierCR . array1s . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:166: "
+{-# LINE 166 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 166 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Batch.fourierCR . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Batch:167: "
+{-# LINE 167 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.property
+{-# LINE 167 "src/Numeric/FFTW/Batch.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Batch.fourierCR . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Batch:148: "
+{-# LINE 148 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.example
+{-# LINE 148 "src/Numeric/FFTW/Batch.hs" #-}
+   (floatList $ Batch.fourierCR $ Array.fromList (Shape.ZeroBased (2::Int), (), Spectrum.Half (3::Int)) [1,0, 0,-1])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00,-2.00,1.00,1.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.Batch:151: "
+{-# LINE 151 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.example
+{-# LINE 151 "src/Numeric/FFTW/Batch.hs" #-}
+   (floatList $ Batch.fourierCR $ Array.fromList ((), Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.Batch:153: "
+{-# LINE 153 "src/Numeric/FFTW/Batch.hs" #-}
+ DocTest.example
+{-# LINE 153 "src/Numeric/FFTW/Batch.hs" #-}
+   (floatList $ Batch.fourierCR $ Array.fromList ((), Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
diff --git a/test/Test/Numeric/FFTW/Common.hs b/test/Test/Numeric/FFTW/Common.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Common.hs
@@ -0,0 +1,159 @@
+module Test.Numeric.FFTW.Common where
+
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.Complex as Complex
+import qualified Data.NonEmpty as NonEmpty
+import Data.Array.Comfort.Storable (Array)
+import Data.Complex (Complex((:+)))
+import Data.NonEmpty ((!:))
+
+import Foreign.Storable.Record.Tuple (Tuple(getTuple))
+import Foreign.Storable (Storable)
+
+import Control.DeepSeq (NFData, deepseq)
+import Control.Applicative ((<*>))
+
+import Text.Printf (printf)
+
+import qualified Test.QuickCheck as QC
+
+
+prefix :: String -> [(String, test)] -> [(String, test)]
+prefix msg =
+   map (\(str,test) -> (msg ++ "." ++ str, test))
+
+
+arrayFloat :: array Float -> array Float
+arrayFloat = id
+
+arrayDouble :: array Double -> array Double
+arrayDouble = id
+
+arrayComplexFloat :: array (Complex Float) -> array (Complex Float)
+arrayComplexFloat = id
+
+arrayComplexDouble :: array (Complex Double) -> array (Complex Double)
+arrayComplexDouble = id
+
+
+floatTol :: Float
+floatTol = 1e-4
+
+doubleTol :: Double
+doubleTol = 1e-10
+
+normInf ::
+   (Shape.C sh, Storable b, Class.Real a) => (b -> a) -> Array sh b -> a
+normInf mag xs =
+   NonEmpty.maximum $ 0 !: Array.toList (Array.map mag xs)
+
+approx ::
+   (Shape.C sh, Eq sh, Storable b, Num b, Class.Real a) =>
+   (b -> a) -> a -> Array sh b -> Array sh b -> Bool
+approx mag tol xs ys =
+   let absTol = tol * (normInf mag xs + normInf mag ys) in
+   all (<=absTol) $ Array.toList $ Array.map mag $ Array.zipWith (-) xs ys
+
+approxReal ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   a -> Array sh a -> Array sh a -> Bool
+approxReal = approx abs
+
+approxComplex ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   a -> Array sh (Complex a) -> Array sh (Complex a) -> Bool
+approxComplex = approx Complex.magnitude
+
+
+adjust :: (Shape.C sh, Class.Floating a) => Array sh a -> Array sh a
+adjust xs = Array.map (fromIntegral (Shape.size (Array.shape xs)) *) xs
+
+split ::
+   (Shape.C sh, Storable a, Storable b) =>
+   Array sh (Tuple (a, b)) -> (Array sh a, Array sh b)
+split xys = (Array.map (fst.getTuple) xys, Array.map (snd.getTuple) xys)
+
+scalarProduct ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   Array sh (Complex a) -> Array sh (Complex a) -> Complex a
+scalarProduct xs ys =
+   sum $ Array.toList $
+   Array.zipWith (\x y -> Complex.conjugate x * y) xs ys
+
+
+factors :: (Integral n) => n -> [(n,n)]
+factors n =
+   concatMap (\k ->
+      case divMod n k of
+         (q,0) -> if q==k then [(k,q)] else [(k,q),(q,k)]
+         _ -> []) $
+   takeWhile (\k -> k*k<=n) $ iterate (1+) 1
+
+genCyclicArray1 ::
+   (QC.Arbitrary a, Storable a) => QC.Gen (Array (Shape.Cyclic Int) a)
+genCyclicArray1 = do
+   xs <- fmap (take 1000) $ QC.arbitrary
+   return $ Array.fromList (Shape.Cyclic $ length xs) xs
+
+genCyclicArray2 ::
+   (QC.Arbitrary a, Storable a) =>
+   QC.Gen (Array (Shape.Cyclic Int, Shape.Cyclic Int) a)
+genCyclicArray2 = do
+   xys <- fmap (take 1000) $ QC.arbitrary
+   (n0,n1) <-
+      case xys of
+         _:_ -> QC.elements (factors $ length xys)
+         [] -> QC.elements [(,) 0, flip (,) 0] <*> QC.choose (0,10)
+   return $ Array.fromList (Shape.Cyclic n0, Shape.Cyclic n1) xys
+
+genCyclicArray3 ::
+   (QC.Arbitrary a, Storable a) =>
+   QC.Gen (Array (Shape.Cyclic Int, Shape.Cyclic Int, Shape.Cyclic Int) a)
+genCyclicArray3 = do
+   xyzs <- fmap (take 1000) $ QC.arbitrary
+   (n0,n1,n2) <-
+      case xyzs of
+         _:_ -> QC.elements $ do
+            (n0,m0) <- factors $ length xyzs
+            (n1,n2) <- factors m0
+            return (n0,n1,n2)
+         [] ->
+            QC.elements [(,,) 0, flip (,,) 0, \n0 n1 -> (n0,n1,0)]
+               <*> QC.choose (0,10)
+               <*> QC.choose (0,10)
+   let sh = (Shape.Cyclic n0, Shape.Cyclic n1, Shape.Cyclic n2)
+   return $ Array.fromList sh xyzs
+
+
+immutable ::
+   (Shape.C sh, Storable a, Eq sh, Eq a, NFData b) =>
+   (Array sh a -> b) -> (Array sh a -> Bool)
+immutable transform arr =
+   let copy = Array.map id arr
+   in deepseq (transform copy) (copy==arr)
+
+
+
+newtype FixedFloat = FixedFloat Float
+
+instance Show FixedFloat where
+   show (FixedFloat a) = printf "%.2f" a
+
+floatList :: (Shape.C sh) => Array sh Float -> [FixedFloat]
+floatList = map FixedFloat . Array.toList
+
+
+newtype FixedComplexFloat = FixedComplexFloat (Complex Float)
+
+instance Show FixedComplexFloat where
+   show (FixedComplexFloat (r:+i)) =
+      let str :: Float -> String
+          str x = if x>=0 then printf "%.2f" x else printf "(%.2f)" x
+      in printf "%s:+%s" (str r) (str i)
+
+complexFloatList ::
+   (Shape.C sh) => Array sh (Complex Float) -> [FixedComplexFloat]
+complexFloatList = map FixedComplexFloat . Array.toList
diff --git a/test/Test/Numeric/FFTW/Rank1.hs b/test/Test/Numeric/FFTW/Rank1.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Rank1.hs
@@ -0,0 +1,345 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/FFTW/Rank1.hs
+{-# LINE 29 "src/Numeric/FFTW/Rank1.hs" #-}
+
+{-# OPTIONS_GHC -XGADTs #-}
+module Test.Numeric.FFTW.Rank1 where
+
+import Test.DocTest.Base
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 31 "src/Numeric/FFTW/Rank1.hs" #-}
+import     Test.Numeric.FFTW.Common
+       (approxReal, approxComplex, floatTol, doubleTol, normInf,
+        adjust, scalarProduct, split, genCyclicArray1, immutable,
+        arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+        floatList, complexFloatList)
+
+import     qualified Numeric.FFTW.Rank1 as Trafo1
+import     qualified Numeric.FFTW.Shape as Spectrum
+import     Numeric.FFTW.Rank1 (flipSign)
+import     Numeric.FFTW.Shape
+       (SymmetrySingleton(Even,Odd), ShiftSingleton(Exact,Halfway))
+
+import     qualified Data.Array.Comfort.Storable as Array
+import     qualified Data.Array.Comfort.Shape as Shape
+import     qualified Data.Complex as Complex
+import     Data.Array.Comfort.Storable (Array)
+
+import     qualified Numeric.Netlib.Class as Class
+
+import     Foreign.Storable (Storable)
+
+import     qualified Test.QuickCheck as QC
+
+genSymmetric     ::
+       (Spectrum.Symmetry symmetry) =>
+       (Spectrum.Shift shiftTime) =>
+       (Spectrum.Shift shiftSpectrum) =>
+       (Spectrum.Symmetric symmetry shiftTime shiftSpectrum Int ~ sh) =>
+       (QC.Arbitrary a, Storable a) =>
+       Spectrum.SymmetrySingleton symmetry ->
+       Spectrum.ShiftSingleton shiftTime ->
+       Spectrum.ShiftSingleton shiftSpectrum ->
+       QC.Gen (Array sh a)
+genSymmetric     symmetry shiftTime shiftSpectrum = do
+       xs <- fmap (take 1000) $ QC.arbitrary
+       return $ Array.fromList
+          (Spectrum.Symmetric symmetry shiftTime shiftSpectrum $ length xs) xs
+
+adjustSymmetric     ::
+       (Spectrum.Symmetry symmetry) =>
+       (Spectrum.Shift shiftTime) =>
+       (Spectrum.Shift shiftSpectrum) =>
+       (Spectrum.Symmetric symmetry shiftTime shiftSpectrum Int ~ sh) =>
+       (Class.Floating a) =>
+       Array sh a -> Array sh a
+adjustSymmetric     xs =
+       let n = Spectrum.symmetricLogicalSize (Array.shape xs)
+       in Array.map (fromIntegral n *) xs
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.FFTW.Rank1:90: "
+{-# LINE 90 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 90 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo1.fourier sign (Trafo1.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:91: "
+{-# LINE 91 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 91 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo1.fourier sign (Trafo1.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:93: "
+{-# LINE 93 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 93 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign -> approxComplex floatTol (Trafo1.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo1.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:94: "
+{-# LINE 94 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 94 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (Trafo1.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo1.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:96: "
+{-# LINE 96 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 96 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo1.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:97: "
+{-# LINE 97 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 97 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo1.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:99: "
+{-# LINE 99 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 99 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:100: "
+{-# LINE 100 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 100 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo1.fourier sign xs) (Trafo1.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:102: "
+{-# LINE 102 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 102 "src/Numeric/FFTW/Rank1.hs" #-}
+     (\sign -> QC.forAll genCyclicArray1 $ immutable (Trafo1.fourier sign) . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:103: "
+{-# LINE 103 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 103 "src/Numeric/FFTW/Rank1.hs" #-}
+     (\sign -> QC.forAll genCyclicArray1 $ immutable (Trafo1.fourier sign) . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:83: "
+{-# LINE 83 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.example
+{-# LINE 83 "src/Numeric/FFTW/Rank1.hs" #-}
+   (complexFloatList $ Trafo1.fourier Trafo1.Forward $ Array.fromList (Shape.Cyclic (5::Int)) [1,0,0,0,0])
+  [ExpectedLine [LineChunk "[1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00,1.00:+0.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.Rank1:85: "
+{-# LINE 85 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.example
+{-# LINE 85 "src/Numeric/FFTW/Rank1.hs" #-}
+   (complexFloatList $ Trafo1.fourier Trafo1.Forward $ Array.fromList (Shape.Cyclic (5::Int)) [0,1,0,0,0])
+  [ExpectedLine [LineChunk "[1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95]"]]
+ DocTest.printPrefix "Numeric.FFTW.Rank1:87: "
+{-# LINE 87 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.example
+{-# LINE 87 "src/Numeric/FFTW/Rank1.hs" #-}
+   (complexFloatList $ Trafo1.fourier Trafo1.Backward $ Array.fromList (Shape.Cyclic (5::Int)) [0,1,0,0,0])
+  [ExpectedLine [LineChunk "[1.00:+0.00,0.31:+0.95,(-0.81):+0.59,(-0.81):+(-0.59),0.31:+(-0.95)]"]]
+ DocTest.printPrefix "Numeric.FFTW.Rank1:118: "
+{-# LINE 118 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 118 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxReal floatTol (adjust xs) (Trafo1.fourierCR (Trafo1.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:119: "
+{-# LINE 119 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 119 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxReal doubleTol (adjust xs) (Trafo1.fourierCR (Trafo1.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:121: "
+{-# LINE 121 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 121 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable Trafo1.fourierRC . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:122: "
+{-# LINE 122 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 122 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable Trafo1.fourierRC . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:140: "
+{-# LINE 140 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 140 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable Trafo1.fourierCR . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:141: "
+{-# LINE 141 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 141 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable Trafo1.fourierCR . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:137: "
+{-# LINE 137 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.example
+{-# LINE 137 "src/Numeric/FFTW/Rank1.hs" #-}
+   (floatList $ Trafo1.fourierCR $ Array.fromList (Spectrum.Half (5::Int)) [0,1,0])
+  [ExpectedLine [LineChunk "[2.00,0.62,-1.62,-1.62,0.62]"]]
+ DocTest.printPrefix "Numeric.FFTW.Rank1:159: "
+{-# LINE 159 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 159 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:160: "
+{-# LINE 160 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 160 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:161: "
+{-# LINE 161 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 161 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:162: "
+{-# LINE 162 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 162 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:163: "
+{-# LINE 163 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 163 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:164: "
+{-# LINE 164 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 164 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:165: "
+{-# LINE 165 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 165 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:166: "
+{-# LINE 166 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 166 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.cosine (Trafo1.cosine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:168: "
+{-# LINE 168 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 168 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ immutable Trafo1.cosine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:169: "
+{-# LINE 169 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 169 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Exact `QC.suchThat` ((>=2) . Shape.size . Array.shape)) $ immutable Trafo1.cosine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:170: "
+{-# LINE 170 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 170 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Exact) $ immutable Trafo1.cosine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:171: "
+{-# LINE 171 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 171 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Exact) $ immutable Trafo1.cosine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:172: "
+{-# LINE 172 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 172 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Halfway) $ immutable Trafo1.cosine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:173: "
+{-# LINE 173 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 173 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Exact Halfway) $ immutable Trafo1.cosine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:174: "
+{-# LINE 174 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 174 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Halfway) $ immutable Trafo1.cosine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:175: "
+{-# LINE 175 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 175 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Even Halfway Halfway) $ immutable Trafo1.cosine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:200: "
+{-# LINE 200 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 200 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:201: "
+{-# LINE 201 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 201 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:202: "
+{-# LINE 202 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 202 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Exact) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:203: "
+{-# LINE 203 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 203 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Exact) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:204: "
+{-# LINE 204 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 204 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:205: "
+{-# LINE 205 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 205 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:206: "
+{-# LINE 206 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 206 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Halfway) $ \xs -> approxReal floatTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:207: "
+{-# LINE 207 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 207 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Halfway) $ \xs -> approxReal doubleTol (adjustSymmetric xs) (Trafo1.sine (Trafo1.sine xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:209: "
+{-# LINE 209 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 209 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Exact) $ immutable Trafo1.sine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:210: "
+{-# LINE 210 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 210 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Exact) $ immutable Trafo1.sine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:211: "
+{-# LINE 211 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 211 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Exact) $ immutable Trafo1.sine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:212: "
+{-# LINE 212 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 212 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Exact) $ immutable Trafo1.sine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:213: "
+{-# LINE 213 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 213 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Halfway) $ immutable Trafo1.sine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:214: "
+{-# LINE 214 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 214 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Exact Halfway) $ immutable Trafo1.sine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:215: "
+{-# LINE 215 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 215 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Halfway) $ immutable Trafo1.sine . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:216: "
+{-# LINE 216 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 216 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll (genSymmetric Odd Halfway Halfway) $ immutable Trafo1.sine . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:239: "
+{-# LINE 239 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 239 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxReal floatTol (adjust xs) (Trafo1.hartley (Trafo1.hartley xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:240: "
+{-# LINE 240 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 240 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxReal doubleTol (adjust xs) (Trafo1.hartley (Trafo1.hartley xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank1:242: "
+{-# LINE 242 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 242 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable Trafo1.hartley . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank1:243: "
+{-# LINE 243 "src/Numeric/FFTW/Rank1.hs" #-}
+ DocTest.property
+{-# LINE 243 "src/Numeric/FFTW/Rank1.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable Trafo1.hartley . arrayDouble)
diff --git a/test/Test/Numeric/FFTW/Rank2.hs b/test/Test/Numeric/FFTW/Rank2.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Rank2.hs
@@ -0,0 +1,103 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/FFTW/Rank2.hs
+{-# LINE 25 "src/Numeric/FFTW/Rank2.hs" #-}
+
+module Test.Numeric.FFTW.Rank2 where
+
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 26 "src/Numeric/FFTW/Rank2.hs" #-}
+import     Test.Numeric.FFTW.Common
+       (approxReal, approxComplex, floatTol, doubleTol, normInf,
+        adjust, scalarProduct, split, genCyclicArray2, immutable,
+        arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble)
+
+import     qualified Numeric.FFTW.Rank2 as Trafo2
+import     Numeric.FFTW.Rank1 (flipSign)
+
+import     qualified Data.Array.Comfort.Storable as Array
+import     qualified Data.Complex as Complex
+
+import     qualified Test.QuickCheck as QC
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.FFTW.Rank2:42: "
+{-# LINE 42 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 42 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo2.fourier sign (Trafo2.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:43: "
+{-# LINE 43 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 43 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo2.fourier sign (Trafo2.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:45: "
+{-# LINE 45 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 45 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> approxComplex floatTol (Trafo2.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo2.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:46: "
+{-# LINE 46 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 46 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (Trafo2.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo2.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:48: "
+{-# LINE 48 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 48 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo2.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:49: "
+{-# LINE 49 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 49 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo2.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:51: "
+{-# LINE 51 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 51 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:52: "
+{-# LINE 52 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 52 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo2.fourier sign xs) (Trafo2.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:54: "
+{-# LINE 54 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 54 "src/Numeric/FFTW/Rank2.hs" #-}
+     (\sign -> QC.forAll genCyclicArray2 $ immutable (Trafo2.fourier sign) . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:55: "
+{-# LINE 55 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 55 "src/Numeric/FFTW/Rank2.hs" #-}
+     (\sign -> QC.forAll genCyclicArray2 $ immutable (Trafo2.fourier sign) . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:70: "
+{-# LINE 70 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 70 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> approxReal floatTol (adjust xs) (Trafo2.fourierCR (Trafo2.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:71: "
+{-# LINE 71 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 71 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> approxReal doubleTol (adjust xs) (Trafo2.fourierCR (Trafo2.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank2:73: "
+{-# LINE 73 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 73 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable Trafo2.fourierRC . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:74: "
+{-# LINE 74 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 74 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable Trafo2.fourierRC . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:90: "
+{-# LINE 90 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 90 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Trafo2.fourierCR . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank2:91: "
+{-# LINE 91 "src/Numeric/FFTW/Rank2.hs" #-}
+ DocTest.property
+{-# LINE 91 "src/Numeric/FFTW/Rank2.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable Trafo2.fourierCR . arrayComplexDouble)
diff --git a/test/Test/Numeric/FFTW/Rank3.hs b/test/Test/Numeric/FFTW/Rank3.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Rank3.hs
@@ -0,0 +1,103 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/FFTW/Rank3.hs
+{-# LINE 25 "src/Numeric/FFTW/Rank3.hs" #-}
+
+module Test.Numeric.FFTW.Rank3 where
+
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 26 "src/Numeric/FFTW/Rank3.hs" #-}
+import     Test.Numeric.FFTW.Common
+       (approxReal, approxComplex, floatTol, doubleTol, normInf,
+        adjust, scalarProduct, split, genCyclicArray3, immutable,
+        arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble)
+
+import     qualified Numeric.FFTW.Rank3 as Trafo3
+import     Numeric.FFTW.Rank1 (flipSign)
+
+import     qualified Data.Array.Comfort.Storable as Array
+import     qualified Data.Complex as Complex
+
+import     qualified Test.QuickCheck as QC
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.FFTW.Rank3:42: "
+{-# LINE 42 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 42 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign-> approxComplex floatTol (adjust xs) (Trafo3.fourier sign (Trafo3.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:43: "
+{-# LINE 43 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 43 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (adjust xs) (Trafo3.fourier sign (Trafo3.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:45: "
+{-# LINE 45 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 45 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign -> approxComplex floatTol (Trafo3.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo3.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:46: "
+{-# LINE 46 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 46 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (Trafo3.fourier sign (Array.map Complex.conjugate xs)) (Array.map Complex.conjugate (Trafo3.fourier (flipSign sign) xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:48: "
+{-# LINE 48 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 48 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in approxComplex floatTol (Trafo3.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:49: "
+{-# LINE 49 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 49 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in approxComplex doubleTol (Trafo3.fourier sign $ Array.zipWith (+) xs ys) (Array.zipWith (+) (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:51: "
+{-# LINE 51 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 51 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)) <= floatTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:52: "
+{-# LINE 52 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 52 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xys sign -> let (xs,ys) = split xys in Complex.magnitude (scalarProduct (adjust xs) ys - scalarProduct (Trafo3.fourier sign xs) (Trafo3.fourier sign ys)) <= doubleTol * normInf Complex.magnitude (adjust xs) * normInf Complex.magnitude ys)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:54: "
+{-# LINE 54 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 54 "src/Numeric/FFTW/Rank3.hs" #-}
+     (\sign -> QC.forAll genCyclicArray3 $ immutable (Trafo3.fourier sign) . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:55: "
+{-# LINE 55 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 55 "src/Numeric/FFTW/Rank3.hs" #-}
+     (\sign -> QC.forAll genCyclicArray3 $ immutable (Trafo3.fourier sign) . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:70: "
+{-# LINE 70 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 70 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> approxReal floatTol (adjust xs) (Trafo3.fourierCR (Trafo3.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:71: "
+{-# LINE 71 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 71 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> approxReal doubleTol (adjust xs) (Trafo3.fourierCR (Trafo3.fourierRC xs)))
+ DocTest.printPrefix "Numeric.FFTW.Rank3:73: "
+{-# LINE 73 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 73 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable Trafo3.fourierRC . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:74: "
+{-# LINE 74 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 74 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable Trafo3.fourierRC . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:92: "
+{-# LINE 92 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 92 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Trafo3.fourierCR . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.Rank3:93: "
+{-# LINE 93 "src/Numeric/FFTW/Rank3.hs" #-}
+ DocTest.property
+{-# LINE 93 "src/Numeric/FFTW/Rank3.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable Trafo3.fourierCR . arrayComplexDouble)
diff --git a/test/Test/Numeric/FFTW/RankN.hs b/test/Test/Numeric/FFTW/RankN.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/RankN.hs
@@ -0,0 +1,216 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/FFTW/RankN.hs
+{-# LINE 27 "src/Numeric/FFTW/RankN.hs" #-}
+
+module Test.Numeric.FFTW.RankN where
+
+import Test.DocTest.Base
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 28 "src/Numeric/FFTW/RankN.hs" #-}
+import     Test.Numeric.FFTW.Common
+       (approxReal, approxComplex, floatTol, doubleTol,
+        genCyclicArray1, genCyclicArray2, genCyclicArray3, immutable,
+        arrayFloat, arrayDouble, arrayComplexFloat, arrayComplexDouble,
+        floatList, complexFloatList)
+
+import     qualified Numeric.FFTW.RankN as TrafoM
+import     qualified Numeric.FFTW.Rank1 as Trafo1
+import     qualified Numeric.FFTW.Rank2 as Trafo2
+import     qualified Numeric.FFTW.Rank3 as Trafo3
+import     qualified Numeric.FFTW.Shape as Spectrum
+
+import     qualified Data.Array.Comfort.Storable as Array
+import     qualified Data.Array.Comfort.Shape as Shape
+import     Data.Array.Comfort.Storable (Array)
+
+import     qualified Test.QuickCheck as QC
+
+array1     :: (Shape.C sh) => Array sh a -> Array ((), sh) a
+array1     = Array.mapShape ((,) ())
+
+array3     :: (Shape.C sh0, Shape.C sh1, Shape.C sh2) =>
+       Array (sh0,sh1,sh2) a -> Array ((sh0,sh1),sh2) a
+array3     = Array.mapShape (\(sh0,sh1,sh2) -> ((sh0,sh1),sh2))
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.FFTW.RankN:59: "
+{-# LINE 59 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 59 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign -> approxComplex floatTol (array1 $ Trafo1.fourier sign xs) (TrafoM.fourier sign $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:60: "
+{-# LINE 60 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 60 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs sign -> approxComplex doubleTol (array1 $ Trafo1.fourier sign xs) (TrafoM.fourier sign $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:62: "
+{-# LINE 62 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 62 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> approxComplex floatTol (Trafo2.fourier sign xs) (TrafoM.fourier sign xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:63: "
+{-# LINE 63 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 63 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs sign -> approxComplex doubleTol (Trafo2.fourier sign xs) (TrafoM.fourier sign xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:65: "
+{-# LINE 65 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 65 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign -> approxComplex floatTol (array3 $ Trafo3.fourier sign xs) (TrafoM.fourier sign $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:66: "
+{-# LINE 66 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 66 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs sign -> approxComplex doubleTol (array3 $ Trafo3.fourier sign xs) (TrafoM.fourier sign $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:56: "
+{-# LINE 56 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.example
+{-# LINE 56 "src/Numeric/FFTW/RankN.hs" #-}
+   (complexFloatList $ TrafoM.fourier TrafoM.Forward $ Array.fromList (Shape.Cyclic (5::Int), Shape.Cyclic (5::Int)) [0,0,0,0,0, 0,1,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0])
+  [ExpectedLine [LineChunk "[1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,(-0.81):+(-0.59),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,0.31:+(-0.95),(-0.81):+0.59,0.31:+0.95,1.00:+0.00,0.31:+(-0.95),(-0.81):+(-0.59),0.31:+0.95,1.00:+(-0.00),0.31:+(-0.95),(-0.81):+(-0.59),(-0.81):+0.59]"]]
+ DocTest.printPrefix "Numeric.FFTW.RankN:80: "
+{-# LINE 80 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 80 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxComplex floatTol (array1 $ Trafo1.fourierRC xs) (TrafoM.fourierRC $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:81: "
+{-# LINE 81 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 81 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ \xs -> approxComplex doubleTol (array1 $ Trafo1.fourierRC xs) (TrafoM.fourierRC $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:83: "
+{-# LINE 83 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 83 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> approxComplex floatTol (Trafo2.fourierRC xs) (TrafoM.fourierRC xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:84: "
+{-# LINE 84 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 84 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ \xs -> approxComplex doubleTol (Trafo2.fourierRC xs) (TrafoM.fourierRC xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:86: "
+{-# LINE 86 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 86 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> approxComplex floatTol (array3 $ Trafo3.fourierRC xs) (TrafoM.fourierRC $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:87: "
+{-# LINE 87 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 87 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ \xs -> approxComplex doubleTol (array3 $ Trafo3.fourierRC xs) (TrafoM.fourierRC $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:90: "
+{-# LINE 90 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 90 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable TrafoM.fourierRC . array1 . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:91: "
+{-# LINE 91 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 91 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray1 $ immutable TrafoM.fourierRC . array1 . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:93: "
+{-# LINE 93 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 93 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable TrafoM.fourierRC . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:94: "
+{-# LINE 94 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 94 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray2 $ immutable TrafoM.fourierRC . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:96: "
+{-# LINE 96 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 96 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable TrafoM.fourierRC . array3 . arrayFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:97: "
+{-# LINE 97 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 97 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll genCyclicArray3 $ immutable TrafoM.fourierRC . array3 . arrayDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:125: "
+{-# LINE 125 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 125 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ \xs -> approxReal floatTol (array1 $ Trafo1.fourierCR xs) (TrafoM.fourierCR $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:126: "
+{-# LINE 126 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 126 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ \xs -> approxReal doubleTol (array1 $ Trafo1.fourierCR xs) (TrafoM.fourierCR $ array1 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:128: "
+{-# LINE 128 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 128 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> approxReal floatTol (Trafo2.fourierCR xs) (TrafoM.fourierCR xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:129: "
+{-# LINE 129 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 129 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ \xs -> approxReal doubleTol (Trafo2.fourierCR xs) (TrafoM.fourierCR xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:131: "
+{-# LINE 131 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 131 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> approxReal floatTol (array3 $ Trafo3.fourierCR xs) (TrafoM.fourierCR $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:132: "
+{-# LINE 132 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 132 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ \xs -> approxReal doubleTol (array3 $ Trafo3.fourierCR xs) (TrafoM.fourierCR $ array3 xs))
+ DocTest.printPrefix "Numeric.FFTW.RankN:135: "
+{-# LINE 135 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 135 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable TrafoM.fourierCR . array1 . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:136: "
+{-# LINE 136 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 136 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo1.fourierRC genCyclicArray1) $ immutable TrafoM.fourierCR . array1 . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:138: "
+{-# LINE 138 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 138 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable TrafoM.fourierCR . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:139: "
+{-# LINE 139 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 139 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo2.fourierRC genCyclicArray2) $ immutable TrafoM.fourierCR . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:141: "
+{-# LINE 141 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 141 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable TrafoM.fourierCR . array3 . arrayComplexFloat)
+ DocTest.printPrefix "Numeric.FFTW.RankN:142: "
+{-# LINE 142 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.property
+{-# LINE 142 "src/Numeric/FFTW/RankN.hs" #-}
+     (QC.forAll (fmap Trafo3.fourierRC genCyclicArray3) $ immutable TrafoM.fourierCR . array3 . arrayComplexDouble)
+ DocTest.printPrefix "Numeric.FFTW.RankN:115: "
+{-# LINE 115 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.example
+{-# LINE 115 "src/Numeric/FFTW/RankN.hs" #-}
+   (floatList $ Trafo1.fourierCR $ Array.fromList (Spectrum.Half (3::Int)) [1,0 ])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.RankN:117: "
+{-# LINE 117 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.example
+{-# LINE 117 "src/Numeric/FFTW/RankN.hs" #-}
+   (floatList $ TrafoM.fourierCR $ Array.fromList ((), Spectrum.Half (3::Int)) [1,0])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.RankN:120: "
+{-# LINE 120 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.example
+{-# LINE 120 "src/Numeric/FFTW/RankN.hs" #-}
+   (floatList $ Trafo2.fourierCR $ Array.fromList (Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
+ DocTest.printPrefix "Numeric.FFTW.RankN:122: "
+{-# LINE 122 "src/Numeric/FFTW/RankN.hs" #-}
+ DocTest.example
+{-# LINE 122 "src/Numeric/FFTW/RankN.hs" #-}
+   (floatList $ TrafoM.fourierCR $ Array.fromList (Shape.Cyclic (3::Int), Spectrum.Half (1::Int)) [1,0,0])
+  [ExpectedLine [LineChunk "[1.00,1.00,1.00]"]]
diff --git a/test/Test/Numeric/FFTW/Shape.hs b/test/Test/Numeric/FFTW/Shape.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Numeric/FFTW/Shape.hs
@@ -0,0 +1,49 @@
+module Test.Numeric.FFTW.Shape where
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import Numeric.FFTW.Shape
+         (SymmetrySingleton(Even,Odd), ShiftSingleton(Exact,Halfway))
+import Test.Numeric.FFTW.Common (prefix)
+
+import qualified Data.Array.Comfort.Shape.Test as ShapeTest
+
+import qualified Test.QuickCheck as QC
+
+import Control.Applicative ((<$>))
+
+
+genHalf :: Int -> QC.Gen (Spectrum.Half Int)
+genHalf n = Spectrum.Half <$> QC.choose (0,n)
+
+genSymmetric ::
+   (Spectrum.Symmetry symmetry) =>
+   (Spectrum.Shift shiftTime) =>
+   (Spectrum.Shift shiftSpectrum) =>
+   Spectrum.SymmetrySingleton symmetry ->
+   Spectrum.ShiftSingleton shiftTime ->
+   Spectrum.ShiftSingleton shiftSpectrum ->
+   Int -> QC.Gen (Spectrum.Symmetric symmetry shiftTime shiftSpectrum Int)
+genSymmetric symmetry shiftTime shiftSpectrum n =
+   Spectrum.Symmetric symmetry shiftTime shiftSpectrum <$> QC.choose (0,n)
+
+tests :: [(String, QC.Property)]
+tests =
+   prefix "Half"
+      (ShapeTest.tests $ genHalf 100) ++
+   prefix "Half Even Exact Exact"
+      (ShapeTest.tests $ genSymmetric Even Exact Exact 100) ++
+   prefix "Half Even Halfway Exact"
+      (ShapeTest.tests $ genSymmetric Even Halfway Exact 100) ++
+   prefix "Half Even Exact Halfway"
+      (ShapeTest.tests $ genSymmetric Even Exact Halfway 100) ++
+   prefix "Half Even Halfway Halfway"
+      (ShapeTest.tests $ genSymmetric Even Halfway Halfway 100) ++
+   prefix "Half Odd Exact Exact"
+      (ShapeTest.tests $ genSymmetric Odd Exact Exact 100) ++
+   prefix "Half Odd Halfway Exact"
+      (ShapeTest.tests $ genSymmetric Odd Halfway Exact 100) ++
+   prefix "Half Odd Exact Halfway"
+      (ShapeTest.tests $ genSymmetric Odd Exact Halfway 100) ++
+   prefix "Half Odd Halfway Halfway"
+      (ShapeTest.tests $ genSymmetric Odd Halfway Halfway 100) ++
+   []