diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2014
+
+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/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/accelerate-fourier.cabal b/accelerate-fourier.cabal
new file mode 100644
--- /dev/null
+++ b/accelerate-fourier.cabal
@@ -0,0 +1,96 @@
+Name:             accelerate-fourier
+Version:          0.0
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         http://code.haskell.org/~thielema/accelerate-fourier/
+Category:         Math
+Synopsis:         Fast Fourier transform and convolution using the Accelerate framework
+Description:
+  Fast Fourier transform and convolution via the @accelerate@ package.
+  It uses many of the common tricks to achieve optimal speed
+  for every size of the data set,
+  but it does not contain low-level optimizations.
+  .
+  * For CUDA specific optimizations please see
+    @cufft@ and @accelerate-fft@ and @accelerate-cufft@ packages.
+  .
+  * For CPU specific optimizations cf. @accelerate-fftw@ package.
+Tested-With:      GHC==7.8.3
+Cabal-Version:    >=1.14
+Build-Type:       Simple
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/accelerate-fourier/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/accelerate-fourier/
+
+Library
+  Build-Depends:
+    accelerate-arithmetic >=0.0.1 && <0.1,
+    accelerate-utility >=0.1 && <0.2,
+    accelerate >=0.15 && <0.16,
+    containers >=0.5 && <0.6,
+    transformers >=0.3 && <0.5,
+    utility-ht >=0.0.8 && <0.1,
+    QuickCheck >=2.4 && <2.8,
+    base >=4.5 && <4.8
+
+  GHC-Options:      -Wall -fwarn-missing-import-lists
+  Hs-Source-Dirs:   src
+  Default-Language: Haskell98
+  Exposed-Modules:
+    Data.Array.Accelerate.Fourier.Planned
+    Data.Array.Accelerate.Fourier.Real
+    Data.Array.Accelerate.Fourier.Adhoc
+    Data.Array.Accelerate.Fourier.Preprocessed
+    Data.Array.Accelerate.Fourier.Utility
+    Data.Array.Accelerate.Convolution.Adhoc
+    Data.Array.Accelerate.Convolution.Preprocessed
+    Data.Array.Accelerate.Cyclic
+  Other-Modules:
+    Data.Array.Accelerate.Convolution.Small
+    Data.Array.Accelerate.Convolution.Private
+    Data.Array.Accelerate.Permutation
+    Data.Array.Accelerate.NumberTheory
+    Data.Array.Accelerate.Fourier.Sign
+    Data.Array.Accelerate.Fourier.Private
+
+Test-Suite test
+  Type: exitcode-stdio-1.0
+  Main-Is:          Test.hs
+  GHC-Options:      -Wall -fwarn-missing-import-lists
+  Hs-Source-Dirs:   test
+  Default-Language: Haskell98
+  Build-Depends:
+    accelerate-fourier,
+    accelerate-arithmetic,
+    accelerate-utility,
+    accelerate,
+    utility-ht,
+    QuickCheck,
+    base
+  Other-Modules:
+    Test.Data.Array.Accelerate.Fourier
+
+Benchmark benchmark
+  Type: exitcode-stdio-1.0
+  Main-Is:          Main.hs
+  GHC-Options:      -Wall -fwarn-missing-import-lists -threaded
+  GHC-Prof-Options: -fprof-auto -rtsopts
+  Hs-Source-Dirs:   benchmark
+  Default-Language: Haskell98
+  Build-Depends:
+    criterion >=1.0 && <1.1,
+    accelerate-fourier,
+    accelerate-arithmetic,
+    accelerate-utility,
+    accelerate-cuda >=0.15 && <0.16,
+    accelerate,
+    utility-ht,
+    base
diff --git a/benchmark/Main.hs b/benchmark/Main.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/Main.hs
@@ -0,0 +1,75 @@
+module Main where
+
+import Criterion.Main (Benchmark, defaultMain, bgroup, bench, whnf, )
+
+import qualified Data.Array.Accelerate.Fourier.Planned as Planned
+import qualified Data.Array.Accelerate.Fourier.Preprocessed as Prep
+import qualified Data.Array.Accelerate.Fourier.Adhoc as Adhoc
+
+import qualified Data.Array.Accelerate.CUDA as CUDA
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate (Array, DIM2, Z(Z), (:.)((:.)), )
+import Data.Complex (Complex, )
+
+
+powersOfTwo ::
+   (Int ->
+    A.Acc (Array DIM2 (Complex Float)) ->
+    A.Acc (Array DIM2 (Complex Float))) ->
+   [Benchmark]
+powersOfTwo f =
+   take 6 $ flip map (iterate (2*) 1024) $ \len ->
+      bench (show len) $ whnf (CUDA.run1 (f len)) $
+      A.fromList (Z:.16:.len) $ repeat 0
+
+powersOfTwos :: [Benchmark]
+powersOfTwos =
+   bgroup "split-radix adhoc"
+      (powersOfTwo (const $ Adhoc.ditSplitRadix Adhoc.forward)) :
+   bgroup "split-radix preprocessed"
+      (powersOfTwo (Prep.ditSplitRadix Prep.forward)) :
+   bgroup "dit2 adhoc"
+      (powersOfTwo (const $ Adhoc.dit2 Adhoc.forward)) :
+   bgroup "dit2 preprocessed"
+      (powersOfTwo (Prep.dit2 Prep.forward)) :
+   bgroup "dif2 preprocessed"
+      (powersOfTwo (Prep.dif2 Prep.forward)) :
+   bgroup "dit2 planned"
+      (powersOfTwo (Planned.transform Planned.forward)) :
+   []
+
+
+arbitrary ::
+   (Int ->
+    A.Acc (Array DIM2 (Complex Float)) ->
+    A.Acc (Array DIM2 (Complex Float))) ->
+   [Benchmark]
+arbitrary f =
+   take 128 $ flip map (iterate (1+) 1) $ \len ->
+      bench (show len) $ whnf (CUDA.run1 (f len)) $
+      A.fromList (Z:.4096:.len) $ repeat 0
+
+arbitraryLengths :: [Benchmark]
+arbitraryLengths =
+   bgroup "auto planned"
+      (arbitrary (Planned.transform Planned.forward)) :
+   bgroup "decompose planned"
+      (arbitrary (Planned.transformDecompose Planned.forward)) :
+   bgroup "chirp235 planned"
+      (arbitrary (Planned.transformChirp235 Planned.forward)) :
+   bgroup "chirp2 planned"
+      (arbitrary (Planned.transformChirp2 Planned.forward)) :
+   bgroup "auto adhoc"
+      (arbitrary (const $ Adhoc.transform Adhoc.forward)) :
+   bgroup "chirp235 adhoc"
+      (arbitrary (const $ Adhoc.transformChirp235 Adhoc.forward)) :
+   bgroup "chirp2 adhoc"
+      (arbitrary (const $ Adhoc.transformChirp2 Adhoc.forward)) :
+   []
+
+
+main :: IO ()
+main = defaultMain $
+   bgroup "2^n" powersOfTwos :
+   bgroup "any" arbitraryLengths :
+   []
diff --git a/src/Data/Array/Accelerate/Convolution/Adhoc.hs b/src/Data/Array/Accelerate/Convolution/Adhoc.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Convolution/Adhoc.hs
@@ -0,0 +1,138 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.Convolution.Adhoc (
+   Transform2,
+   karatsuba,
+   cyclic,
+   complex,
+   ) where
+
+import Data.Array.Accelerate.Convolution.Private (Transform2, indexPad, )
+import Data.Array.Accelerate.Fourier.Private (Transform, )
+
+import qualified Data.Array.Accelerate.Utility.Sliced1 as Sliced1
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+import Data.Array.Accelerate.Utility.Lift.Acc (acc)
+
+import qualified Data.Array.Accelerate.Data.Complex as Complex
+import Data.Array.Accelerate.Data.Complex (Complex((:+)), )
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Exp, Acc, Array, IsNum, Elt,
+           Z(Z), (:.)((:.)), Any(Any), All(All), Slice, Shape,
+           (!), (>*), )
+
+
+{- |
+Both arrays must have the same size.
+-}
+karatsuba ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Transform2 (sh :. Int) a
+karatsuba x y =
+   flip A.slice (A.lift $ Any :. (0::Int) :. All)
+   .
+   A.afst
+   .
+   A.awhile
+      (\arrs -> A.unit $ (Sliced.length $ A.asnd arrs) >* 1)
+      (Acc.modify (acc, acc) $
+       \(z, lens) ->
+          (karatsubaGo (lens ! A.index1 0) (2*(lens ! A.index1 1)-1) z,
+           Sliced.tail lens))
+   .
+   (Acc.modify ((acc, acc), acc) $
+    \((x0,y0), lens) -> (A.zipWith (*) x0 y0, lens))
+   .
+   A.awhile
+      (\arrs -> A.unit $ (Sliced.length $ A.afst $ A.afst arrs) >* 1)
+      (Acc.modify ((acc, acc), acc) $
+       \((x0,y0), lens) ->
+          let (x1,y1) = karatsubaReorder (x0,y0)
+          in  ((x1,y1), Sliced.consExp (Sliced.length x1) lens))
+   $
+   A.lift
+      ((A.replicate (A.lift $ Any :. (1::Int) :. All) x,
+        A.replicate (A.lift $ Any :. (1::Int) :. All) y),
+       A.fill (A.constant $ Z:.1) (Sliced.length x))
+
+karatsubaReorder ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   (Acc (Array (sh :. Int :. Int) a),
+    Acc (Array (sh :. Int :. Int) a)) ->
+   (Acc (Array (sh :. Int :. Int) a),
+    Acc (Array (sh :. Int :. Int) a))
+karatsubaReorder (x,y) =
+   let len2 = - div (- Sliced.length x) 2
+       xl = Sliced.take len2 x
+       yl = Sliced.take len2 y
+       xr = Sliced.pad 0 len2 $ Sliced.drop len2 x
+       yr = Sliced.pad 0 len2 $ Sliced.drop len2 y
+   in  (Sliced1.append3 xl (A.zipWith (+) xl xr) xr,
+        Sliced1.append3 yl (A.zipWith (+) yl yr) yr)
+
+karatsubaGo ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Exp Int ->
+   Exp Int ->
+   Transform (sh :. Int :. Int) a
+karatsubaGo xlen zlen zmerged =
+   let (sh:.n:._m) = Exp.unlift (expr:.expr:.expr) $ A.shape zmerged
+       n3 = div n 3
+       zl = Sliced1.take n3 zmerged
+       zm = Sliced1.drop n3 zmerged
+       zr = Sliced1.drop (2*n3) zmerged
+       zc = A.zipWith (-) zm $ A.zipWith (+) zl zr
+   in  A.generate (A.lift $ sh :. n3 :. zlen) $
+       Exp.modify (expr:.expr) $
+       \(ix:.k) ->
+          indexPad (ix:.k)        zl +
+          indexPad (ix:.k-xlen)   zc +
+          indexPad (ix:.k-xlen*2) zr
+
+
+{- |
+Turn an ordinary convolution into a cyclic convolution of the same length.
+-}
+cyclic ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Transform2 (sh :. Int) a ->
+   Transform2 (sh :. Int) a
+cyclic conv x y =
+   let z = conv x y
+       len = Sliced.length x
+   in  A.zipWith (+) z $ Sliced.pad 0 len $ Sliced.drop len z
+
+
+{- |
+Turn a real-valued convolution into a complex-valued convolution.
+Can be removed when we get @instance IsNum (Complex a)@.
+-}
+complex, _complex ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Transform2 (sh :. Int) a ->
+   Transform2 (sh :. Int) (Complex a)
+complex conv x y =
+   let xr = A.map Complex.real x; xi = A.map Complex.imag x
+       yr = A.map Complex.real y; yi = A.map Complex.imag y
+       xm = A.zipWith (+) xr xi
+       ym = A.zipWith (+) yr yi
+       xryr = conv xr yr
+       xiyi = conv xi yi
+       xmym = conv xm ym
+   in  A.zipWith
+          (Exp.modify2 expr expr (:+))
+          (A.zipWith (-) xryr xiyi)
+          (A.zipWith (-) xmym $ A.zipWith (+) xryr xiyi)
+
+_complex conv x y =
+   let xr = A.map Complex.real x; xi = A.map Complex.imag x
+       yr = A.map Complex.real y; yi = A.map Complex.imag y
+   in  A.zipWith
+          (Exp.modify2 expr expr (:+))
+          (A.zipWith (-) (conv xr yr) (conv xi yi))
+          (A.zipWith (+) (conv xr yi) (conv xi yr))
diff --git a/src/Data/Array/Accelerate/Convolution/Preprocessed.hs b/src/Data/Array/Accelerate/Convolution/Preprocessed.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Convolution/Preprocessed.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.Convolution.Preprocessed (
+   Transform2,
+   karatsuba,
+   ) where
+
+import Data.Array.Accelerate.Convolution.Private (Transform2, indexPad, )
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (IsNum, Elt,
+           (:.)((:.)), Any(Any), All(All), Slice, Shape, )
+
+
+{- |
+Both arrays must have the same size.
+
+There is not much to preprocess,
+thus you should prefer 'Data.Array.Accelerate.Convolution.Adhoc.karatsuba'.
+-}
+karatsuba ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Int -> Transform2 (sh :. Int) a
+karatsuba len x y =
+   if len <= 1
+     then A.zipWith (*) x y
+     else
+        let len2 = - div (-len) 2
+            elen2 = A.constant len2
+            xl = Sliced.take elen2 x
+            yl = Sliced.take elen2 y
+            xr = Sliced.pad 0 elen2 $ Sliced.drop elen2 x
+            yr = Sliced.pad 0 elen2 $ Sliced.drop elen2 y
+            zmerged =
+               karatsuba len2
+                  (Sliced.stack3 xl (A.zipWith (+) xl xr) xr)
+                  (Sliced.stack3 yl (A.zipWith (+) yl yr) yr)
+            zl = A.slice zmerged $ A.lift $ Any :. (0::Int) :. All
+            zm = A.slice zmerged $ A.lift $ Any :. (1::Int) :. All
+            zr = A.slice zmerged $ A.lift $ Any :. (2::Int) :. All
+            zc = A.zipWith (-) zm $ A.zipWith (+) zl zr
+            sh = A.indexTail $ A.shape zc
+        in  A.generate (A.lift $ sh :. 2*len-1) $
+            Exp.modify (expr:.expr) $
+            \(ix:.k) ->
+               indexPad (ix:.k)         zl +
+               indexPad (ix:.k-elen2)   zc +
+               indexPad (ix:.k-elen2*2) zr
diff --git a/src/Data/Array/Accelerate/Convolution/Private.hs b/src/Data/Array/Accelerate/Convolution/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Convolution/Private.hs
@@ -0,0 +1,25 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.Convolution.Private where
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Exp, Acc, Array, IsNum, Elt,
+           (:.)((:.)), Slice, Shape,
+           (!), (?), (&&*), (<*), (<=*), )
+
+
+type Transform2 sh a =
+   Acc (Array sh a) ->
+   Acc (Array sh a) ->
+   Acc (Array sh a)
+
+
+indexPad ::
+   (Shape sh, Slice sh, Elt a, IsNum a) =>
+   Exp sh :. Exp Int ->
+   Acc (Array (sh:.Int) a) -> Exp a
+indexPad (ix:.k) xs =
+   0 <=* k &&* k <* Sliced.length xs ? (xs ! A.lift (ix:.k), 0)
diff --git a/src/Data/Array/Accelerate/Convolution/Small.hs b/src/Data/Array/Accelerate/Convolution/Small.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Convolution/Small.hs
@@ -0,0 +1,85 @@
+-- cf. synthesizer-core:Synthesizer.Generic.Cyclic
+{- |
+Some small size convolutions using the Karatsuba trick.
+We do not use Toom-3 multiplication,
+because this requires division by 2 and 6,
+and thus 'Fractional' constraints.
+-}
+module Data.Array.Accelerate.Convolution.Small where
+
+
+type Pair a = (a,a)
+
+convolvePair ::
+   (Num a) =>
+   Pair a -> Pair a -> Pair a
+convolvePair a b =
+   snd $ sumAndConvolvePair a b
+
+sumAndConvolvePair ::
+   (Num a) =>
+   Pair a -> Pair a -> ((a,a), Pair a)
+sumAndConvolvePair (a0,a1) (b0,b1) =
+   let sa01 = a0+a1
+       sb01 = b0+b1
+       ab0ab1 = a0*b0+a1*b1
+   in  ((sa01, sb01), (ab0ab1, sa01*sb01-ab0ab1))
+
+
+type Triple a = (a,a,a)
+
+convolveTriple ::
+   (Num a) =>
+   Triple a -> Triple a -> Triple a
+convolveTriple a b =
+   snd $ sumAndConvolveTriple a b
+
+sumAndConvolveTriple ::
+   (Num a) =>
+   Triple a -> Triple a -> ((a,a), Triple a)
+sumAndConvolveTriple (a0,a1,a2) (b0,b1,b2) =
+   let ab0 = a0*b0
+       dab12 = a1*b1 - a2*b2
+       sa01 = a0+a1; sb01 = b0+b1; tab01 = sa01*sb01 - ab0
+       sa02 = a0+a2; sb02 = b0+b2; tab02 = sa02*sb02 - ab0
+       sa012 = sa01+a2
+       sb012 = sb01+b2
+
+       d0 = sa012*sb012 - tab01 - tab02
+       d1 = tab01 - dab12
+       d2 = tab02 + dab12
+   in  ((sa012, sb012), (d0, d1, d2))
+
+
+type Quadruple a = (a,a,a,a)
+
+convolveQuadruple ::
+   (Num a) =>
+   Quadruple a -> Quadruple a -> Quadruple a
+convolveQuadruple a b =
+   snd $ sumAndConvolveQuadruple a b
+
+sumAndConvolveQuadruple ::
+   (Num a) =>
+   Quadruple a -> Quadruple a -> ((a,a), Quadruple a)
+sumAndConvolveQuadruple (a0,a1,a2,a3) (b0,b1,b2,b3) =
+   let ab0 = a0*b0
+       ab1 = a1*b1
+       sa01 = a0+a1; sb01 = b0+b1
+       ab01 = sa01*sb01 - (ab0+ab1)
+       ab2 = a2*b2
+       ab3 = a3*b3
+       sa23 = a2+a3; sb23 = b2+b3
+       ab23 = sa23*sb23 - (ab2+ab3)
+       c0 = ab0  + ab2 - (ab1 + ab3)
+       c1 = ab01 + ab23
+       ab02 = (a0+a2)*(b0+b2)
+       ab13 = (a1+a3)*(b1+b3)
+       sa0123 = sa01+sa23
+       sb0123 = sb01+sb23
+       ab0123 = sa0123*sb0123 - (ab02+ab13)
+       d0 = ab13   + c0
+       d1 =          c1
+       d2 = ab02   - c0
+       d3 = ab0123 - c1
+   in  ((sa0123, sb0123), (d0, d1, d2, d3))
diff --git a/src/Data/Array/Accelerate/Cyclic.hs b/src/Data/Array/Accelerate/Cyclic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Cyclic.hs
@@ -0,0 +1,39 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.Cyclic (
+   Transform,
+   reverse,
+   reverse2d,
+   ) where
+
+import Data.Array.Accelerate.Fourier.Private (Transform)
+
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate ((:.)((:.)), )
+
+import Prelude (Int, mod, ($), )
+
+
+reverse ::
+   (A.Shape sh, A.Slice sh, A.Elt a) =>
+   Transform (sh :. Int) a
+reverse arr =
+   let sh = A.shape arr
+   in  A.backpermute sh
+          (Exp.modify (expr:.expr) $
+           \(ix:.k) -> ix :. mod (-k) (A.indexHead sh))
+          arr
+
+reverse2d ::
+   (A.Shape sh, A.Slice sh, A.Elt a) =>
+   Transform (sh :. Int :. Int) a
+reverse2d arr =
+   let sh = A.shape arr
+       (_z:.height:.width) = Exp.unlift (expr:.expr:.expr) sh
+   in  A.backpermute sh
+          (Exp.modify (expr:.expr:.expr) $
+           \(ix:.y:.x) -> ix :. mod (-y) height :. mod (-x) width)
+          arr
diff --git a/src/Data/Array/Accelerate/Fourier/Adhoc.hs b/src/Data/Array/Accelerate/Fourier/Adhoc.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Adhoc.hs
@@ -0,0 +1,380 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{- |
+The implementations in this module work entirely in the 'A.Acc' domain.
+This means that they can be applied to any array
+without knowing their extent on the Haskell side.
+The downside is, that they cannot share any preprocessing.
+-}
+module Data.Array.Accelerate.Fourier.Adhoc (
+   Transform,
+   transform,
+
+   ditSplitRadix,
+   dit2,
+   dit235,
+
+   ceiling5Smooth,
+   transformChirp2,
+   transformChirp235,
+
+   Sign,
+   forward,
+   inverse,
+
+   transform2d,
+   transform3d,
+
+   SubTransform(SubTransform),
+   ) where
+
+import qualified Data.Array.Accelerate.Fourier.Private as Fourier
+import qualified Data.Array.Accelerate.Fourier.Sign as Sign
+import Data.Array.Accelerate.Fourier.Private
+          (SubTransform(SubTransform), SubTransformPair(SubTransformPair),
+           Transform, twist, )
+import Data.Array.Accelerate.Fourier.Utility (scaleDown, )
+import Data.Array.Accelerate.Fourier.Sign (Sign, )
+
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import Data.Array.Accelerate.Arithmetic.LinearAlgebra
+          (zipExtrudedVectorWith, zipExtrudedMatrixWith, )
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Sliced1 as Sliced1
+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+import Data.Array.Accelerate.Utility.Lift.Acc (acc)
+import Data.Array.Accelerate.Utility.Ord (argminimum)
+
+import qualified Data.Array.Accelerate.Data.Complex as Complex
+import Data.Array.Accelerate.Data.Complex (Complex, )
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Slice, Shape, DIM1, Z(Z), (:.)((:.)),
+           Exp, Elt, IsFloating, (<=*), (>*), (==*), )
+
+import Data.Bits ((.&.))
+
+
+forward, inverse :: (Elt a, A.IsNum a) => Exp (Sign a)
+forward = Sign.forwardExp
+inverse = Sign.inverseExp
+
+
+{- |
+Automatically choose transformation algorithms
+according to the size of the array.
+However, they are not as sophisticated as the algorithms in
+"Data.Array.Accelerate.Fourier.Planned".
+-}
+transform ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int) (Complex a)
+transform sign arr =
+   let len = Sliced.length arr
+   in  A.acond (len <=* 1) arr $
+          let (pow2, smooth5) = is2or5smooth len
+          in  A.acond pow2 (ditSplitRadixLoop sign arr) $
+              A.acond smooth5 (dit235 sign arr) $
+              transformChirp235 sign arr
+
+is2or5smooth :: Exp Int -> (Exp Bool, Exp Bool)
+is2or5smooth len =
+   let maxPowerOfTwo = len .&. negate len
+       lenOdd = div len maxPowerOfTwo
+   in  (lenOdd ==* 1,
+        (divideMaxPower 5 $ divideMaxPower 3 lenOdd) ==* 1)
+
+divideMaxPower :: Exp Int -> Exp Int -> Exp Int
+divideMaxPower fac =
+   A.while (\n -> mod n fac ==* 0) (flip div fac)
+
+
+{- |
+Split-Radix for power-of-two sizes.
+-}
+ditSplitRadix ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int) (Complex a)
+ditSplitRadix sign arr =
+   A.acond
+      (Sliced.length arr <=* 1)
+      arr (ditSplitRadixLoop sign arr)
+
+ditSplitRadixLoop ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int) (Complex a)
+ditSplitRadixLoop sign =
+   Fourier.finishSplitRadix . A.afst
+   .
+   A.awhile
+      (\x -> A.unit $ (Sliced1.length $ A.asnd x) >* 0)
+      (Acc.modify (acc, acc) $
+       \(arr2, arr1) ->
+         Fourier.ditSplitRadixStep
+            (Fourier.imagSplitRadix sign)
+            (Fourier.twiddleFactorsSRPair sign $ Sliced.length arr1)
+            (arr2, arr1))
+   .
+   Acc.modify (acc, acc) Fourier.ditSplitRadixBase
+   .
+   A.awhile
+      (\x -> A.unit $ (Sliced.length $ A.asnd x) >* 1)
+      (Acc.modify (acc, acc) Fourier.ditSplitRadixReorder)
+   .
+   A.lift . Fourier.initSplitRadix
+
+
+{- |
+Decimation in time for power-of-two sizes.
+-}
+dit2 ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int) (Complex a)
+dit2 sign =
+   flip A.slice (A.lift $ A.Any :. (0::Int) :. A.All)
+   .
+   A.awhile
+      (\x -> A.unit $ Sliced1.length x >* 1)
+      (ditStep sign)
+   .
+   A.awhile
+      (\x -> A.unit $ Sliced.length x >* 1)
+      (twist 2)
+   .
+   A.replicate (A.lift $ A.Any :. (1::Int) :. A.All)
+
+ditStep ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int:.Int) (Complex a)
+ditStep sign x =
+   let twiddles = Fourier.twiddleFactors2 sign $ Sliced.length x
+       evens = Sliced1.sieve 2 0 x
+       odds = zipExtrudedVectorWith (*) twiddles $ Sliced1.sieve 2 1 x
+   in  A.zipWith (+) evens odds  A.++  A.zipWith (-) evens odds
+
+
+{- |
+Decimation in time for sizes that are composites of the factors 2, 3 and 5.
+These sizes are known as 5-smooth numbers or the Hamming sequence.
+<http://oeis.org/A051037>.
+-}
+dit235 ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Transform (sh:.Int) (Complex a)
+dit235 sign =
+   flip A.slice (A.lift $ A.Any :. (0::Int) :. A.All)
+   .
+   A.afst
+   .
+   A.awhile
+      (\x -> A.unit $ (A.length $ A.asnd x) >* 0)
+      (Acc.modify (acc,acc) $
+       \(arr,factors) ->
+         let fac = factors A.! A.index1 0
+         in  (dit235Step sign fac arr, Sliced.tail factors))
+   .
+   A.awhile
+      (\x -> A.unit $ (Sliced.length $ A.afst x) >* 1)
+      (Acc.modify (acc,acc) $
+       \(arr,factors) ->
+         let divides k n = mod n k ==* 0
+             caseFactor k = (divides k, k)
+             len = Sliced.length arr
+             factor =
+                flip (A.caseof len) 2 $
+                   caseFactor 3 :
+                   caseFactor 4 :
+                   caseFactor 5 :
+                   []
+         in  (twist factor arr, Sliced.consExp factor factors))
+   .
+   A.lift . flip (,) (A.fill (A.index1 0) 0)
+   .
+   A.replicate (A.lift $ A.Any :. (1::Int) :. A.All)
+
+dit235Step ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Exp Int ->
+   Transform (sh:.Int:.Int) (Complex a)
+dit235Step sign fac x =
+   let (sh:.count:.len) = Exp.unlift (expr:.expr:.expr) $ A.shape x
+       twiddled =
+          LinAlg.transpose .
+          zipExtrudedMatrixWith (*) (Fourier.twiddleFactors sign fac len) .
+          A.reshape (A.lift (sh :. div count fac :. fac :. len))
+          $
+          x
+   in  Fourier.merge $
+       A.acond (fac ==* 5) (Fourier.transform5 (Fourier.cache5 sign) twiddled) $
+       A.acond (fac ==* 4) (Fourier.transform4 (Fourier.cache4 sign) twiddled) $
+       A.acond (fac ==* 3) (Fourier.transform3 (Fourier.cache3 sign) twiddled) $
+       Fourier.transform2 (Fourier.cache2 sign) twiddled
+
+
+{- |
+Next greater or equal 5-smooth number as needed by 'dit235'.
+-}
+ceiling5Smooth :: Exp Int -> Exp Int
+ceiling5Smooth n =
+   Exp.modify (expr,expr,expr)
+      (\(e2, e3, e5) -> pow e2 2 * pow e3 3 * pow e5 5) $
+   A.snd $ ceiling5SmoothFloat $
+   (A.fromIntegral n :: Exp Double)
+
+ceiling5SmoothFloat ::
+   (Elt a, IsFloating a) =>
+   Exp a -> Exp (a, (Int, Int, Int))
+ceiling5SmoothFloat n =
+   let d3 = ceilingLogBase 3 n
+       d5 = ceilingLogBase 5 n
+   in  A.the $ argminimum $
+       A.generate (A.lift $ Z:.d5:.d3) $
+       Exp.modify (expr:.expr:.expr) $
+       \(_z:.e5:.e3) ->
+          let p53 = 5 ** A.fromIntegral e5 * 3 ** A.fromIntegral e3
+              e2 = max 0 $ ceilingLogBase 2 $ n/p53
+          in  (p53 * 2 ** A.fromIntegral e2, (e2, e3, e5))
+
+{-
+Should be more efficient than ceiling5SmoothFloat,
+but sometimes misses optimal results due to rounding errors.
+-}
+_ceiling5SmoothLog ::
+   (Elt a, IsFloating a) =>
+   Exp a -> Exp (a, (Int, Int, Int))
+_ceiling5SmoothLog n =
+   let log3 = logBase 2 3
+       log5 = logBase 2 5
+       logN = logBase 2 n
+       d3 = max 1 $ A.ceiling $ logN / log3
+       d5 = max 1 $ A.ceiling $ logN / log5
+   in  A.the $ argminimum $
+       A.generate (A.lift $ Z:.d5:.d3) $
+       Exp.modify (expr:.expr:.expr) $
+       \(_z:.e5:.e3) ->
+          let logP53 = log5 * A.fromIntegral e5 + log3 * A.fromIntegral e3
+              e2 = max 0 $ A.ceiling $ logN-logP53
+          in  (logP53 + A.fromIntegral e2, (e2, e3, e5))
+
+_ceiling5SmoothFloat ::
+   (Elt a, IsFloating a) =>
+   Exp a -> Exp (a, (Int, Int, Int))
+_ceiling5SmoothFloat n =
+   let powers base =
+          A.scanl (*) 1 $
+          A.fill (A.index1 $ ceilingLogBase base n) $
+          A.fromIntegral base
+   in  A.the $ argminimum $
+       Arrange.mapWithIndex
+          (Exp.modify2 (expr:.expr:.expr) expr $
+           \(_z:.e5:.e3) p53 ->
+              let e2 = max 0 $ ceilingLogBase 2 $ n/p53
+              in  (p53 * 2 ** A.fromIntegral e2, (e2, e3, e5))) $
+       LinAlg.outer (powers 5) (powers 3)
+
+ceilingLogBase ::
+   (Elt a, IsFloating a) =>
+   Exp Int -> Exp a -> Exp Int
+ceilingLogBase base x =
+   A.ceiling $ logBase (A.fromIntegral base) x
+
+pow :: Exp Int -> Exp Int -> Exp Int
+pow e n = A.the $ A.product $ A.fill (A.index1 e) n
+
+
+_transformChirp ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Exp Int ->
+   (Transform DIM1 (Complex a),
+    Transform (sh:.Int) (Complex a),
+    Transform (sh:.Int) (Complex a)) ->
+   Transform (sh:.Int) (Complex a)
+_transformChirp sign padLen (analysis1,analysis,synthesis) arr =
+   let len = Sliced.length arr
+       chirp = Fourier.chirp sign padLen $ A.fromIntegral len
+   in  A.acond (len<=*1) arr $
+       Sliced.take len $ scaleDown $
+       LinAlg.zipExtrudedVectorWith (*) chirp $ synthesis $
+       LinAlg.zipExtrudedVectorWith (*)
+          (analysis1 $ A.map Complex.conjugate chirp)
+          (analysis $ Sliced.pad 0 padLen $
+           LinAlg.zipExtrudedVectorWith (*) chirp arr)
+
+transformChirp ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Exp Int ->
+   SubTransformPair (Complex a) ->
+   Transform (sh:.Int) (Complex a)
+transformChirp sign padLen (SubTransformPair analysis synthesis) arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+       chirp = Fourier.chirp sign padLen $ A.fromIntegral len
+       spectrum =
+          analysis $
+          Sliced1.cons
+             (A.map Complex.conjugate chirp)
+             (A.reshape (A.lift $ A.index1 (A.shapeSize sh) :. padLen) $
+              Sliced.pad 0 padLen $
+              LinAlg.zipExtrudedVectorWith (*) chirp arr)
+   in  A.acond (len<=*1) arr $
+       Sliced.take len $ scaleDown $
+       LinAlg.zipExtrudedVectorWith (*) chirp $ synthesis $
+       LinAlg.zipExtrudedVectorWith (*)
+          (Sliced1.head spectrum)
+          (A.reshape (A.lift $ sh:.padLen) $ Sliced1.tail spectrum)
+
+{- |
+Transformation of arbitrary length based on Bluestein on a power-of-two size.
+-}
+transformChirp2 ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   Exp (Sign a) -> Transform (sh:.Int) (Complex a)
+transformChirp2 sign arr =
+   transformChirp sign
+      (let n = Sliced.length arr
+       in  pow (ceilingLogBase 2 (A.fromIntegral (2*n-1) :: Exp Double)) 2)
+      (SubTransformPair (ditSplitRadix forward) (ditSplitRadix inverse))
+      arr
+
+{- |
+Transformation of arbitrary length based on Bluestein on a 5-smooth size.
+-}
+transformChirp235 ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   Exp (Sign a) -> Transform (sh:.Int) (Complex a)
+transformChirp235 sign arr =
+   transformChirp sign
+      (ceiling5Smooth (2 * Sliced.length arr))
+      (SubTransformPair (dit235 forward) (dit235 inverse))
+      arr
+
+
+transform2d ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   SubTransform (Complex a) ->
+   Transform (sh:.Int:.Int) (Complex a)
+transform2d (SubTransform trans) =
+   LinAlg.transpose . trans .
+   LinAlg.transpose . trans
+
+transform3d ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   SubTransform (Complex a) ->
+   Transform (sh:.Int:.Int:.Int) (Complex a)
+transform3d (SubTransform trans) =
+   Fourier.cycleDim3 . trans .
+   Fourier.cycleDim3 . trans .
+   Fourier.cycleDim3 . trans
diff --git a/src/Data/Array/Accelerate/Fourier/Planned.hs b/src/Data/Array/Accelerate/Fourier/Planned.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Planned.hs
@@ -0,0 +1,823 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{- |
+Like "Data.Array.Accelerate.Fourier.Preprocessed"
+this module allows to factor out some preprocessing.
+Additionally it gives you concrete objects (plans and caches)
+for sharing preprocessed data between transforms.
+You cannot only share the preprocessing between transforms of the same size,
+but across all array sizes.
+This implementation also has the largest collection of algorithms
+and thus should be generally fastest among all implementations in this package.
+-}
+module Data.Array.Accelerate.Fourier.Planned (
+   -- * Transforms
+   Transform,
+
+   transform,
+
+   transformDecompose,
+   transformChirp2,
+   transformChirp235,
+
+   convolveCyclic,
+
+   -- * Planning
+   Plan,
+   plan,
+   transformWithPlanner,
+
+   PlanMap,
+   planWithMapUpdate,
+   planDecomposeWithMapUpdate,
+   planChirpWithMapUpdate,
+   smallPlanMap,
+
+   -- * Caching
+   Cache,
+   cache,
+   cacheDuplex,
+   transformWithCache,
+
+   CacheMap,
+   Direction(..),
+   cacheFromPlanWithMapUpdate, directionMode,
+   cacheFromPlanWithMapUpdate2, directionModes,
+
+   -- * Miscellaneous
+   Sign.Sign,
+   Sign.forward,
+   Sign.inverse,
+   ) where
+
+import qualified Data.Array.Accelerate.Convolution.Adhoc as Convolution
+import qualified Data.Array.Accelerate.Permutation as Permutation
+import qualified Data.Array.Accelerate.NumberTheory as NumberTheory
+import qualified Data.Array.Accelerate.Fourier.Private as Fourier
+import qualified Data.Array.Accelerate.Fourier.Sign as Sign
+import Data.Array.Accelerate.Fourier.Private
+          (SubTransform(SubTransform), SubTransformPair(SubTransformPair),
+           SubPairTransform(SubPairTransform),
+           PairTransform, Transform, )
+import Data.Array.Accelerate.Fourier.Utility (scaleDown, )
+import Data.Array.Accelerate.Fourier.Sign (Sign(Sign))
+
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import Data.Array.Accelerate.Arithmetic.LinearAlgebra
+          (zipExtrudedVectorWith, zipExtrudedMatrixWith, )
+
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Data.Complex (Complex, conjugate, )
+import Data.Array.Accelerate
+          (Exp, Acc, Array, DIM1, DIM2, IsFloating, Elt,
+           (:.)((:.)), Slice, Shape, )
+
+import qualified Control.Monad.Trans.State as State
+import Control.Monad (liftM2, )
+import Control.Applicative ((<$>), )
+import Data.Traversable (for, )
+
+import qualified Data.Map as Map
+import Data.Tuple.HT (mapPair, )
+
+
+{- |
+Fourier transform of arbitrary size.
+Sign can be
+
+* @forward@: from time domain to frequency spectrum
+
+* @inverse@: from frequency spectrum to time domain
+
+You may share @transform sign n@ between several calls
+in order to run some preprocessing only once.
+You must make sure that the @length@
+is equal to the extent of the inner dimension of every transformed array.
+-}
+transform ::
+   (Slice sh, Shape sh, RealFloat a, Elt a, IsFloating a) =>
+   Sign a -> Int -> Transform (sh:.Int) (Complex a)
+transform sign len = transformWithCache $ cache sign len
+
+
+{- |
+Transform using only Cooley-Tukey, Good-Thomas, Rader, Split-Radix,
+but no Bluestein.
+This is more for testing and benchmarking than for real use.
+-}
+transformDecompose ::
+   (Slice sh, Shape sh, Elt a, IsFloating a, RealFloat a) =>
+   Sign a -> Int ->
+   Transform (sh :. Int) (Complex a)
+transformDecompose =
+   transformWithPlanner planDecomposeWithMapUpdate
+
+transformWithPlanner ::
+   (Slice sh, Shape sh, Elt a, IsFloating a, RealFloat a) =>
+   (Integer -> State.State PlanMap Plan) ->
+   Sign a -> Int ->
+   Transform (sh :. Int) (Complex a)
+transformWithPlanner planner sign len =
+   transformWithCache $
+   cacheFromPlan
+      (flip State.evalState smallPlanMap $ planner $ fromIntegral len) $
+   directionMode sign len
+
+
+{- |
+The size and type of the signal must match the parameters,
+that the cache was generated for.
+-}
+transformWithCache ::
+   (Slice sh, Shape sh, Elt a, IsFloating a) =>
+   Cache (Complex a) -> Transform (sh:.Int) (Complex a)
+transformWithCache ch =
+   case ch of
+      CacheIdentity -> id
+      CacheSmall size ->
+         case size of
+            LevelCache2 zs -> Fourier.transform2 zs
+            LevelCache3 zs -> Fourier.transform3 zs
+            LevelCache4 zs -> Fourier.transform4 zs
+            LevelCache5 zs -> Fourier.transform5 zs
+      CacheRadix2 level subCache ->
+         transformRadix2InterleavedTime level $
+         subTransformWithCache subCache
+      CacheSplitRadix chain ->
+         Fourier.finishSplitRadix . fst .
+         transformSplitRadixInterleavedTimeChain chain .
+         Fourier.initSplitRadix
+      CachePrime level subCaches ->
+         transformPrime level $
+         fmap subTransformPairWithCache subCaches
+      CacheCoprime level subCaches ->
+         transformCoprime level $
+         subTransformPairWithCache subCaches
+      CacheComposite level subCaches ->
+         transformComposite level $
+         subTransformPairWithCache subCaches
+      CacheChirp level subCaches ->
+         transformChirp level $
+         subTransformPairWithCache subCaches
+
+subTransformWithCache ::
+   (Elt a, IsFloating a) =>
+   Cache (Complex a) -> SubTransform (Complex a)
+subTransformWithCache ch = SubTransform (transformWithCache ch)
+
+subTransformPairWithCache ::
+   (Elt a, IsFloating a) =>
+   (Cache (Complex a), Cache (Complex a)) -> SubTransformPair (Complex a)
+subTransformPairWithCache (ch0,ch1) =
+   SubTransformPair (transformWithCache ch0) (transformWithCache ch1)
+
+
+{- |
+Memorize factorizations of the data size and permutation vectors.
+-}
+data Plan = Plan Integer PlanStructure
+   deriving (Show)
+
+data PlanStructure =
+     PlanIdentity
+   | PlanSmall LevelSmall
+   | PlanRadix2 Plan
+   | PlanSplitRadix Plan
+   | PlanPrime (Maybe Plan)
+   | PlanCoprime (Plan, Plan)
+   | PlanComposite (Plan, Plan)
+   | PlanChirp Plan
+   deriving (Show)
+
+
+{- |
+Plan transform algorithms for a certain array size.
+-}
+plan :: Integer -> Plan
+plan n =
+   State.evalState (planWithMapUpdate n) smallPlanMap
+
+{- |
+Too many nested Rader transformations slow down the transform,
+up to quadratic time in the worst case.
+As a heuristic we allow at most nesting depth two,
+and switch to Bluestein transformation otherwise.
+We could compute more precise operation counts
+and base our decision on these,
+but we found that the actual execution time
+differs considerably from the operation counts.
+-}
+planWithMapUpdate :: Integer -> State.State PlanMap Plan
+planWithMapUpdate n = do
+   p <- planDecomposeWithMapUpdate n
+   if planCountPrimes p < 3
+     then return p
+     else planChirpWithMapUpdate NumberTheory.ceiling5Smooth n
+
+planCountPrimes :: Plan -> Int
+planCountPrimes (Plan _ struct) =
+   case struct of
+      PlanIdentity -> 0
+      PlanSmall _ -> 0
+      PlanRadix2 p -> planCountPrimes p
+      PlanSplitRadix p -> planCountPrimes p
+      PlanPrime mp -> 1 + maybe 0 planCountPrimes mp
+      PlanCoprime (m, n) -> max (planCountPrimes m) (planCountPrimes n)
+      PlanComposite (m, n) -> max (planCountPrimes m) (planCountPrimes n)
+      PlanChirp p -> planCountPrimes p
+
+type PlanMap = Map.Map Integer Plan
+
+{- |
+Map of primitive transforms.
+You should use this as the initial map
+when evaluating a planning sequence using 'State.evalState'.
+-}
+smallPlanMap :: PlanMap
+smallPlanMap =
+   Map.fromAscList $ zipWith (\n struct -> (n, Plan n struct)) [0..] $
+   PlanIdentity :
+   PlanIdentity :
+   PlanSmall Level2 :
+   PlanSmall Level3 :
+   PlanSmall Level4 :
+   PlanSmall Level5 :
+   []
+
+{- |
+Detect and re-use common sub-plans.
+-}
+planDecomposeWithMap :: Integer -> State.State PlanMap Plan
+planDecomposeWithMap n =
+   fmap (Plan n) $
+   case divMod n 2 of
+      (n2,0) ->
+         case (0::Int) of
+            0 -> PlanSplitRadix <$> planDecomposeWithMapUpdate n2
+            1 -> PlanComposite <$> planDecomposeWithMapUpdate2 (2,n2)
+            _ -> PlanRadix2 <$> planDecomposeWithMapUpdate n2
+      _ ->
+         let facs = NumberTheory.fermatFactors n
+         in  -- find unitary divisors
+             case filter (\(a,b) -> a>1 && gcd a b == 1) facs of
+                q2 : _ -> PlanCoprime <$> planDecomposeWithMapUpdate2 q2
+                _ ->
+                   let (q2 : _) = facs
+                   in  if fst q2 == 1
+                         then
+                            PlanPrime <$>
+                            if False
+                              then return Nothing
+                              else Just <$> planDecomposeWithMapUpdate (n-1)
+                         else PlanComposite <$> planDecomposeWithMapUpdate2 q2
+
+planDecomposeWithMapUpdate :: Integer -> State.State PlanMap Plan
+planDecomposeWithMapUpdate n = do
+   item <- State.gets (Map.lookup n)
+   case item of
+      Just p -> return p
+      Nothing -> do
+         m <- planDecomposeWithMap n
+         State.modify (Map.insert n m)
+         return m
+
+planDecomposeWithMapUpdate2 ::
+   (Integer, Integer) -> State.State PlanMap (Plan, Plan)
+planDecomposeWithMapUpdate2 =
+   uncurry (liftM2 (,)) .
+   mapPair (planDecomposeWithMapUpdate,planDecomposeWithMapUpdate)
+
+
+{- |
+Cache arrays of twiddle factors,
+i.e. powers of the primitive root of unity.
+-}
+data Cache a =
+     CacheIdentity
+   | CacheSmall (LevelCacheSmall (Exp a))
+   | CacheRadix2 (LevelCacheRadix2 a) (Cache a)
+   | CacheSplitRadix (CacheSplitRadixChain a)
+   | CachePrime (LevelCachePrime a) (Maybe (Cache a, Cache a))
+   | CacheCoprime LevelCacheCoprime (Cache a, Cache a)
+   | CacheComposite (LevelCacheComposite a) (Cache a, Cache a)
+   | CacheChirp (LevelCacheChirp a) (Cache a, Cache a)
+   deriving (Show)
+
+data CacheSplitRadixChain a =
+     CacheSplitRadixCons (LevelCacheSplitRadix a) (CacheSplitRadixChain a)
+   | CacheSplitRadixEnd (Cache a) (Cache a)
+   deriving (Show)
+
+{- |
+The expression @cache sign len@
+precomputes all data that is needed for Fourier transforms
+for signals of length @len@.
+You can use this cache in 'transformWithCache'.
+-}
+cache ::
+   (RealFloat a, Elt a, IsFloating a) =>
+   Sign a -> Int -> Cache (Complex a)
+cache sign len =
+   cacheFromPlan
+      (plan $ fromIntegral len)
+      (directionMode sign len)
+
+
+{- |
+It is @(cache inverse x, cache forward x) = cacheDuplex x@
+but 'cacheDuplex' shares common data of both caches.
+-}
+cacheDuplex ::
+   (a ~ Complex b, RealFloat b, Elt b, IsFloating b) =>
+   Int -> (Cache a, Cache a)
+cacheDuplex len =
+   let p = plan $ fromIntegral len
+   in  flip State.evalState Map.empty $
+          cacheFromPlanWithMapUpdate2 (p,p) (directionModes len)
+
+
+data Direction = Forward | Inverse
+   deriving (Show, Eq, Ord)
+
+
+type CacheMap a = Map.Map (Integer,Direction) (Cache a)
+
+cacheFromPlan ::
+   (a ~ Complex b, RealFloat b, Elt b, IsFloating b) =>
+   Plan -> (Direction, Sign b) -> Cache a
+cacheFromPlan p z =
+   State.evalState (cacheFromPlanWithMapUpdate p z) Map.empty
+
+
+{- |
+Detect and re-use common sub-caches.
+-}
+cacheFromPlanWithMap ::
+   (a ~ Complex b, RealFloat b, Elt b, IsFloating b) =>
+   Plan -> (Direction, Sign b) ->
+   State.State (CacheMap a) (Cache a)
+cacheFromPlanWithMap (Plan len struct) dsign@(_d,sign) =
+   case struct of
+      PlanIdentity -> return $ CacheIdentity
+      PlanSmall size -> return $ CacheSmall $
+         case size of
+            Level2 -> LevelCache2 $ Fourier.cache2 $ A.constant sign
+            Level3 -> LevelCache3 $ Fourier.cache3 $ A.constant sign
+            Level4 -> LevelCache4 $ Fourier.cache4 $ A.constant sign
+            Level5 -> LevelCache5 $ Fourier.cache5 $ A.constant sign
+      PlanRadix2 subPlan@(Plan len2 _) ->
+         CacheRadix2 (levelCacheRadix2 len2 sign) <$>
+         cacheFromPlanWithMapUpdate subPlan dsign
+      PlanSplitRadix subPlan@(Plan len2 subStruct) -> do
+         subCache <- cacheFromPlanWithMapUpdate subPlan dsign
+         case subCache of
+            CacheSplitRadix chain ->
+               return $
+                  CacheSplitRadix $
+                  CacheSplitRadixCons (levelCacheSplitRadix len2 sign) chain
+            _ ->
+               case subStruct of
+                  PlanSplitRadix subsubPlan -> do
+                     subsubCache <- cacheFromPlanWithMapUpdate subsubPlan dsign
+                     return $ CacheSplitRadix $
+                        CacheSplitRadixCons (levelCacheSplitRadix len2 sign) $
+                        CacheSplitRadixEnd subCache subsubCache
+                  _ ->
+                     return $ CacheRadix2 (levelCacheRadix2 len2 sign) subCache
+      PlanPrime maybeSubPlan ->
+         (\maybeSubCaches ->
+            CachePrime
+               (levelCachePrime len
+                  (fmap (subTransformWithCache . fst) maybeSubCaches) sign)
+               maybeSubCaches)
+         <$>
+         for maybeSubPlan
+            (\subPlan ->
+               cacheFromPlanWithMapUpdate2 (subPlan,subPlan)
+                  (directionModes $ fromInteger $ len-1))
+      PlanCoprime subPlans@(Plan n _, Plan m _) ->
+         CacheCoprime (levelCacheCoprime (n,m)) <$>
+         cacheFromPlanWithMapUpdate2 subPlans (dsign, dsign)
+      PlanComposite subPlans@(Plan n _, Plan m _) ->
+         CacheComposite (levelCacheComposite (n,m) sign)
+         <$>
+         cacheFromPlanWithMapUpdate2 subPlans (dsign, dsign)
+      PlanChirp subPlan@(Plan padlen _) ->
+         (\subCaches ->
+            CacheChirp
+               (levelCacheChirp len padlen
+                  (subTransformWithCache (fst subCaches)) sign)
+               subCaches)
+         <$>
+         cacheFromPlanWithMapUpdate2 (subPlan,subPlan)
+            (directionModes $ fromInteger padlen)
+
+cacheFromPlanWithMapUpdate ::
+   (a ~ Complex b, RealFloat b, Elt b, IsFloating b) =>
+   Plan -> (Direction, Sign b) ->
+   State.State (CacheMap a) (Cache a)
+cacheFromPlanWithMapUpdate p@(Plan len _) z = do
+   let key = (len, fst z)
+   item <- State.gets (Map.lookup key)
+   case item of
+      Just c -> return c
+      Nothing -> do
+         m <- cacheFromPlanWithMap p z
+         State.modify (Map.insert key m)
+         return m
+
+cacheFromPlanWithMapUpdate2 ::
+   (a ~ Complex b, RealFloat b, Elt b, IsFloating b) =>
+   (Plan, Plan) -> ((Direction, Sign b), (Direction, Sign b)) ->
+   State.State (CacheMap a) (Cache a, Cache a)
+cacheFromPlanWithMapUpdate2 (p0,p1) (dm0,dm1) =
+   liftM2 (,)
+      (cacheFromPlanWithMapUpdate p0 dm0)
+      (cacheFromPlanWithMapUpdate p1 dm1)
+
+
+directionMode ::
+   (Num a, Ord a) =>
+   Sign a -> Int -> (Direction, Sign a)
+directionMode (Sign sign) len =
+   (if sign>0 then fst else snd) $ directionModes len
+
+directionModes ::
+   (Num a) =>
+   Int -> ((Direction, Sign a), (Direction, Sign a))
+directionModes _len =
+   ((Inverse, Sign.inverse), (Forward, Sign.forward))
+
+
+
+data LevelSmall = Level2 | Level3 | Level4 | Level5
+   deriving (Show, Eq, Ord, Enum)
+
+data LevelCacheSmall a =
+     LevelCache2 a
+   | LevelCache3 (a,a)
+   | LevelCache4 (a,a,a)
+   | LevelCache5 (a,a,a,a)
+   deriving (Show)
+
+
+
+data LevelCacheRadix2 a = LevelCacheRadix2 (Acc (Array DIM1 a))
+   deriving (Show)
+
+levelCacheRadix2 ::
+   (Elt a, IsFloating a) =>
+   Integer -> Sign a -> LevelCacheRadix2 (Complex a)
+levelCacheRadix2 n2 sign =
+   LevelCacheRadix2 $
+   Fourier.twiddleFactors2 (A.constant sign) (expInteger n2)
+
+transformRadix2InterleavedTime ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCacheRadix2 a ->
+   SubTransform a ->
+   Transform (sh:.Int) a
+transformRadix2InterleavedTime
+      (LevelCacheRadix2 twiddles) (SubTransform subTrans) =
+   Fourier.transformRadix2InterleavedTime twiddles subTrans
+
+
+data LevelCacheSplitRadix a =
+   LevelCacheSplitRadix a (Acc (Array DIM1 a), Acc (Array DIM1 a))
+   deriving (Show)
+
+levelCacheSplitRadix ::
+   (Elt a, IsFloating a) =>
+   Integer -> Sign a -> LevelCacheSplitRadix (Complex a)
+levelCacheSplitRadix n2 sign =
+   LevelCacheSplitRadix (Fourier.imagSplitRadixPlain sign) $
+   Fourier.twiddleFactorsSRPair (A.constant sign) (expInteger (div n2 2))
+
+transformSplitRadixInterleavedTime ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCacheSplitRadix a ->
+   SubPairTransform a ->
+   PairTransform (sh:.Int:.Int) a
+transformSplitRadixInterleavedTime
+      (LevelCacheSplitRadix imag twiddles) (SubPairTransform subTrans) =
+   Fourier.ditSplitRadixStep (A.constant imag) twiddles .
+   subTrans .
+   Fourier.ditSplitRadixReorder
+
+transformSplitRadixInterleavedTimeChain ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   CacheSplitRadixChain a ->
+   PairTransform (sh:.Int:.Int) a
+transformSplitRadixInterleavedTimeChain chain =
+   case chain of
+      CacheSplitRadixCons level remChain ->
+         transformSplitRadixInterleavedTime level $
+         SubPairTransform (transformSplitRadixInterleavedTimeChain remChain)
+      CacheSplitRadixEnd subCache2 subCache1 ->
+         mapPair (transformWithCache subCache2, transformWithCache subCache1)
+
+
+newtype LevelCacheComposite a =
+   LevelCacheComposite (Acc (Array DIM2 a))
+   deriving (Show)
+
+levelCacheComposite ::
+   (Elt a, IsFloating a) =>
+   (Integer, Integer) -> Sign a -> LevelCacheComposite (Complex a)
+levelCacheComposite (n,m) sign =
+   LevelCacheComposite $
+   Fourier.twiddleFactors (A.constant sign) (expInteger n) (expInteger m)
+
+
+{- |
+For @transformComposite z (n,m) sig@,
+it must hold @n*m == length sig@ and @z ^ length sig == 1@.
+
+Cooley-Tukey-algorithm
+-}
+transformComposite ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCacheComposite a ->
+   SubTransformPair a ->
+   Transform (sh:.Int) a
+transformComposite
+      (LevelCacheComposite twiddles)
+      (SubTransformPair subTransN subTransM) =
+
+   Fourier.merge .
+   subTransN .
+   LinAlg.transpose .
+   zipExtrudedMatrixWith (*) twiddles .
+   subTransM .
+   Sliced.sliceHorizontal (A.shape twiddles)
+
+
+newtype LevelCacheCoprime = LevelCacheCoprime (Integer, Integer)
+   deriving (Show)
+
+{-
+Fourier exponent matrix of a signal of size 6.
+
+0 0 0 0 0 0     0               0   0   0       0     0
+0 1 2 3 4 5               0       2   0   4       3     0
+0 2 4 0 2 4  =          0    *  0   2   4    *      0     0
+0 3 0 3 0 3           0           0   0   0     0     3
+0 4 2 0 4 2         0           0   4   2         0     0
+0 5 4 3 2 1       0               4   0   2         0     3
+-}
+levelCacheCoprime :: (Integer, Integer) -> LevelCacheCoprime
+levelCacheCoprime = LevelCacheCoprime
+
+
+{- |
+For @transformCoprime z (n,m) sig@,
+the parameters @n@ and @m@ must be relatively prime
+and @n*m == length sig@ and @z ^ length sig == 1@.
+
+Good-Thomas algorithm
+-}
+transformCoprime ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCacheCoprime ->
+   SubTransformPair a ->
+   Transform (sh:.Int) a
+transformCoprime
+      (LevelCacheCoprime (n,m)) (SubTransformPair subTransN subTransM) =
+   permuteSkewGridInv .
+   subTransM .
+   LinAlg.transpose .
+   subTransN .
+   permuteSkewGrid (expInteger m) (expInteger n)
+
+permuteSkewGrid ::
+   (Slice sh, Shape sh, Elt a) =>
+   Exp Int -> Exp Int -> LinAlg.Vector sh a -> LinAlg.Matrix sh a
+permuteSkewGrid m n arr =
+   let (sh:.nm) = Exp.unlift (expr:.expr) $ A.shape arr
+   in  A.backpermute
+          (A.lift (sh :. m :. n))
+          (Exp.modify (expr:.expr:.expr) $
+           \(ix:.k:.j) -> ix :. mod (n*k + m*j) nm)
+          arr
+
+permuteSkewGridInv ::
+   (Slice sh, Shape sh, Elt a) =>
+   LinAlg.Matrix sh a -> LinAlg.Vector sh a
+permuteSkewGridInv arr =
+   let (sh:.m:.n) = Exp.unlift (expr:.expr:.expr) $ A.shape arr
+   in  A.backpermute
+          (A.lift (sh :. n*m))
+          (Exp.modify (expr:.expr) $
+           \(ix:.k) -> ix :. mod k m :. mod k n)
+          arr
+
+
+
+{-
+Fourier exponent matrix of a signal of size 7.
+
+0 0 0 0 0 0 0
+0 1 2 3 4 5 6
+0 2 4 6 1 3 5
+0 3 6 2 5 1 4
+0 4 1 5 2 6 3
+0 5 3 1 6 4 2
+0 6 5 4 3 2 1
+
+multiplicative generator in Z7: 3
+permutation of rows and columns by powers of 3: 1 3 2 6 4 5
+
+0 0 0 0 0 0 0
+0 1 3 2 6 4 5
+0 3 2 6 4 5 1
+0 2 6 4 5 1 3
+0 6 4 5 1 3 2
+0 4 5 1 3 2 6
+0 5 1 3 2 6 4
+
+Inverse permutation: 1 3 2 5 6 4
+The inverse permutations seems not to be generated by a multiplication.
+-}
+data LevelCachePrime a =
+   LevelCachePrime (Permutation.T, Permutation.T) (Acc (Array DIM1 a))
+      deriving (Show)
+
+levelCachePrime ::
+   (RealFloat a, IsFloating a, Elt a) =>
+   Integer ->
+   Maybe (SubTransform (Complex a)) ->
+   Sign a -> LevelCachePrime (Complex a)
+levelCachePrime n maybeSubTrans sign =
+   let len = fromInteger n
+       perm = A.use $ Permutation.multiplicative len
+       kernel =
+          A.map (Sign.cisRat (A.constant sign) (A.constant len) . (1+)) perm
+   in  LevelCachePrime
+          (Permutation.reverse perm, Permutation.inverse perm)
+          (maybe id
+             (\(SubTransform subTrans) -> scaleDown . subTrans)
+             maybeSubTrans kernel)
+
+{- |
+Rader's algorithm for prime length signals.
+-}
+transformPrime ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCachePrime a ->
+   Maybe (SubTransformPair a) ->
+   Transform (sh:.Int) a
+transformPrime (LevelCachePrime (rev, inv) zs) maybeSubTranss =
+   let conv =
+          case maybeSubTranss of
+             Nothing ->
+                \xs ->
+                   Convolution.complex
+                      (Convolution.cyclic Convolution.karatsuba)
+                      (LinAlg.extrudeVector (A.indexTail $ A.shape xs) zs)
+                      xs
+             Just subTranss ->
+                convolveSingleSpectrumCyclicCache subTranss zs
+   in  \arr ->
+          let x0 = Sliced.head arr
+              res = Sliced.tail arr
+          in  LinAlg.zipScalarVectorWith (+) x0 $
+              Sliced.cons (A.fold1 (+) res) $
+              Permutation.apply inv $
+              conv $
+              Permutation.apply rev res
+
+
+
+{- |
+Fourier transform for arbitrary lengths
+based on the Bluestein transform or chirp z-transform
+on an array with power-of-two size.
+It may be faster than 'transform' for certain prime factors.
+Find bad factors e.g. in <http://oeis.org/A061092> and <http://oeis.org/A059411>
+and nicer factors in <http://oeis.org/A061303>.
+-}
+transformChirp2 ::
+   (Slice sh, Shape sh, Elt a, IsFloating a, RealFloat a) =>
+   Sign a -> Int ->
+   Transform (sh :. Int) (Complex a)
+transformChirp2 = transformChirpComplete NumberTheory.ceilingPowerOfTwo
+
+{- |
+Fourier transform for arbitrary lengths
+based on the Bluestein transform
+on an array with 5-smooth size.
+(5-smooth = all prime factors are at most 5)
+-}
+transformChirp235 ::
+   (Slice sh, Shape sh, Elt a, IsFloating a, RealFloat a) =>
+   Sign a -> Int ->
+   Transform (sh :. Int) (Complex a)
+transformChirp235 = transformChirpComplete NumberTheory.ceiling5Smooth
+
+
+transformChirpComplete ::
+   (Slice sh, Shape sh, Elt a, IsFloating a, RealFloat a) =>
+   (Integer -> Integer) ->
+   Sign a -> Int ->
+   Transform (sh :. Int) (Complex a)
+transformChirpComplete padLength =
+   transformWithPlanner (planChirpWithMapUpdate padLength)
+
+planChirpWithMapUpdate ::
+   (Integer -> Integer) -> Integer -> State.State PlanMap Plan
+planChirpWithMapUpdate padLength len =
+   Plan len <$>
+   if len<2
+     then return PlanIdentity
+     else PlanChirp <$> planDecomposeWithMapUpdate (padLength (2*len-1))
+
+
+data LevelCacheChirp a =
+   LevelCacheChirp (Acc (Array DIM1 a)) (Acc (Array DIM1 a))
+      deriving (Show)
+
+levelCacheChirp ::
+   (RealFloat a, IsFloating a, Elt a) =>
+   Integer -> Integer ->
+   SubTransform (Complex a) ->
+   Sign a -> LevelCacheChirp (Complex a)
+levelCacheChirp len padlen (SubTransform subTrans) sign =
+   let chirp =
+          Fourier.chirp (A.constant sign) (expInteger padlen) (expInteger len)
+   in  LevelCacheChirp
+          (A.take (expInteger len) chirp)
+          (scaleDown $ subTrans $ A.map conjugate chirp)
+
+{- |
+Bluestein's algorithm for signals of arbitrary length
+and possibly slightly generalised basis vectors.
+-}
+transformChirp ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   LevelCacheChirp a ->
+   SubTransformPair a ->
+   Transform (sh:.Int) a
+transformChirp (LevelCacheChirp chirp chirpSpec) subTranss =
+   let conv = convolveSingleSpectrumCyclicCache subTranss chirpSpec
+       twistChirp = zipExtrudedVectorWith (*) chirp
+   in  \arr ->
+          twistChirp $
+          Sliced.take (Sliced.length arr) $
+          conv $
+          Sliced.pad 0 (A.length chirpSpec) $
+          twistChirp arr
+
+
+{- |
+Signals must have equal size and must not be empty.
+-}
+convolveCyclic ::
+   (Shape sh, Slice sh, a ~ Complex b, Elt b, IsFloating b, RealFloat b) =>
+   Int ->
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int) a)
+convolveCyclic leni =
+   let len = fromIntegral leni
+       (z,zInv) = directionModes leni
+   in  convolveCyclicCache $
+       subTransformPairWithCache
+          (cacheFromPlan (plan len) z,
+           cacheFromPlan (plan len) zInv)
+
+convolveCyclicCache ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   SubTransformPair a ->
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int) a)
+convolveCyclicCache transs@(SubTransformPair trans _) x =
+   convolveSpectrumCyclicCache transs $ scaleDown $ trans x
+
+convolveSingleSpectrumCyclicCache ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   SubTransformPair a ->
+   Acc (Array DIM1 a) -> Transform (sh:.Int) a
+convolveSingleSpectrumCyclicCache caches x y =
+   convolveSpectrumCyclicCache caches
+      (LinAlg.extrudeVector (A.indexTail $ A.shape y) x) y
+
+{- |
+This function does not apply scaling.
+That is you have to scale the spectrum by @recip (length x)@
+if you want a plain convolution.
+-}
+convolveSpectrumCyclicCache ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   SubTransformPair a ->
+   Acc (Array (sh:.Int) a) -> Transform (sh:.Int) a
+convolveSpectrumCyclicCache (SubTransformPair trans transInv) x y =
+   transInv $ A.zipWith (*) x (trans y)
+
+
+expInteger :: (Elt a, Num a) => Integer -> Exp a
+expInteger = A.constant . fromInteger
diff --git a/src/Data/Array/Accelerate/Fourier/Preprocessed.hs b/src/Data/Array/Accelerate/Fourier/Preprocessed.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Preprocessed.hs
@@ -0,0 +1,162 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{- |
+The implementations in this module require
+that you know the transformation data set size on the Haskell side.
+This knowledge is baked into the Accelerate code.
+The advantage is,
+that you can share preprocessing between calls to the Fourier transforms,
+like in:
+
+> let transform = dit2 1024
+> in  transform x ... transform y
+-}
+module Data.Array.Accelerate.Fourier.Preprocessed (
+   Transform,
+   ditSplitRadix,
+   dit2,
+   dif2,
+
+   Sign.Sign,
+   Sign.forward,
+   Sign.inverse,
+
+   transform2d,
+   transform3d,
+
+   SubTransformPair(SubTransformPair),
+   SubTransformTriple(SubTransformTriple),
+   ) where
+
+import qualified Data.Array.Accelerate.Fourier.Private as Fourier
+import qualified Data.Array.Accelerate.Fourier.Sign as Sign
+import Data.Array.Accelerate.Fourier.Sign (Sign, )
+import Data.Array.Accelerate.Fourier.Private
+          (SubTransformPair(SubTransformPair),
+           SubTransformTriple(SubTransformTriple),
+           Transform, PairTransform, )
+
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import Data.Array.Accelerate.Arithmetic.LinearAlgebra (zipExtrudedVectorWith, )
+import Data.Array.Accelerate.Data.Complex (Complex, )
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Slice, Shape, (:.), Exp, Elt, IsFloating, )
+
+
+{- |
+Decimation in time for power-of-two using the split-radix algorithm.
+Should be faster than 'dit2'.
+-}
+ditSplitRadix ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Sign a ->
+   Int ->
+   Transform (sh:.Int) (Complex a)
+ditSplitRadix mode len =
+   if len<2
+     then id
+     else
+        Fourier.finishSplitRadix . fst .
+        ditSplitRadixGo (A.constant mode) (div len 2) .
+        Fourier.initSplitRadix
+
+{- |
+Compute the Fourier transforms
+of a collection of 2N length signals
+and a collection of N length signals
+and share some computations between them.
+The global extent of @sh@ of all arrays must be equal.
+First array must have extent @sh:.count2:.2*len@
+and second array must have extent @sh:.count1:.len@.
+If this is a restriction for you,
+you may use 'Fourier.finishSplitRadixFlat' and 'Fourier.initSplitRadixFlat'
+which merge the global shape with our auxiliary dimension
+and then work with @sh = Z@.
+-}
+ditSplitRadixGo ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Exp (Sign a) ->
+   Int ->
+   PairTransform (sh:.Int:.Int) (Complex a)
+ditSplitRadixGo mode len =
+   if len<=1
+     then Fourier.ditSplitRadixBase
+     else
+        let len2 = div len 2
+            twiddles = Fourier.twiddleFactorsSRPair mode (A.constant len2)
+            imag = Fourier.imagSplitRadix mode
+        in  Fourier.ditSplitRadixStep imag twiddles .
+            ditSplitRadixGo mode len2 .
+            Fourier.ditSplitRadixReorder
+
+
+{- |
+Decimation in time for power-of-two sizes.
+-}
+dit2 ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Sign a ->
+   Int ->
+   Transform (sh:.Int) (Complex a)
+dit2 mode len =
+   if len<=1
+     then id
+     else
+        let len2 = div len 2
+        in  Fourier.transformRadix2InterleavedTime
+               (Fourier.twiddleFactors2 (A.constant mode) (A.constant len2))
+               (dit2 mode len2)
+
+
+{- |
+Decimation in frequency for power-of-two sizes.
+-}
+dif2 ::
+   (Slice sh, Shape sh, IsFloating a, Elt a) =>
+   Sign a ->
+   Int ->
+   Transform (sh:.Int) (Complex a)
+dif2 mode len =
+   if len<=1
+     then id
+     else
+        let len2 = div len 2
+            twiddles = Fourier.twiddleFactors2 (A.constant mode) (A.constant len2)
+        in  \arr ->
+              let part0 = Sliced.take (A.constant len2) arr
+                  part1 = Sliced.drop (A.constant len2) arr
+                  evens = A.zipWith (+) part0 part1
+                  odds =
+                     zipExtrudedVectorWith (*) twiddles $
+                     A.zipWith (-) part0 part1
+              in  Fourier.merge $ dif2 mode len2 $ Fourier.stack evens odds
+
+
+{- |
+Transforms in 'SubTransformPair'
+are ordered from least-significant to most-significant dimension.
+-}
+transform2d ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   SubTransformPair (Complex a) ->
+   Transform (sh:.Int:.Int) (Complex a)
+transform2d (SubTransformPair transform0 transform1) =
+   LinAlg.transpose . transform1 .
+   LinAlg.transpose . transform0
+
+{- |
+Transforms in 'SubTransformTriple'
+are ordered from least-significant to most-significant dimension.
+-}
+transform3d ::
+   (Shape sh, Slice sh, IsFloating a, Elt a) =>
+   SubTransformTriple (Complex a) ->
+   Transform (sh:.Int:.Int:.Int) (Complex a)
+transform3d (SubTransformTriple transform0 transform1 transform2) =
+   Fourier.cycleDim3 . transform2 .
+   Fourier.cycleDim3 . transform1 .
+   Fourier.cycleDim3 . transform0
diff --git a/src/Data/Array/Accelerate/Fourier/Private.hs b/src/Data/Array/Accelerate/Fourier/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Private.hs
@@ -0,0 +1,431 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE Rank2Types #-}
+module Data.Array.Accelerate.Fourier.Private where
+
+import qualified Data.Array.Accelerate.Fourier.Sign as Sign
+import qualified Data.Array.Accelerate.Convolution.Small as Cyclic
+import Data.Array.Accelerate.Fourier.Sign (Sign, )
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Sliced1 as Sliced1
+
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import Data.Array.Accelerate.Arithmetic.LinearAlgebra (zipExtrudedVectorWith, )
+
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Data.Complex (Complex((:+)), )
+import Data.Array.Accelerate
+          (Exp, Acc, Array, DIM1, DIM2, IsNum, IsFloating, Elt,
+           Z(Z), (:.)((:.)), Slice, Shape, (!), (?), (==*), (<*), )
+
+
+type Transform sh a = Acc (Array sh a) -> Acc (Array sh a)
+
+data SubTransform a =
+        SubTransform
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+
+data SubTransformPair a =
+        SubTransformPair
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+
+data SubTransformTriple a =
+        SubTransformTriple
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+           (forall sh. (Shape sh, Slice sh) => Transform (sh:.Int) a)
+
+
+type PairTransform sh a =
+        (Acc (Array sh a), Acc (Array sh a)) ->
+        (Acc (Array sh a), Acc (Array sh a))
+
+data SubPairTransform a =
+        SubPairTransform
+           (forall sh. (Shape sh, Slice sh) => PairTransform (sh:.Int:.Int) a)
+
+
+cache2 :: (sign ~ Exp (Sign b), a ~ Exp (Complex b), Elt b, IsFloating b) =>
+   sign -> a
+cache3 :: (sign ~ Exp (Sign b), a ~ Exp (Complex b), Elt b, IsFloating b) =>
+   sign -> (a,a)
+cache4 :: (sign ~ Exp (Sign b), a ~ Exp (Complex b), Elt b, IsFloating b) =>
+   sign -> (a,a,a)
+cache5 :: (sign ~ Exp (Sign b), a ~ Exp (Complex b), Elt b, IsFloating b) =>
+   sign -> (a,a,a,a)
+
+cache2 _sign = -1
+
+cache3 sign =
+   let sqrt3d2 = sqrt 3 / 2
+       mhalf = -1/2
+       s = Sign.toSign sign
+   in  (A.lift $ mhalf :+ s*sqrt3d2,
+        A.lift $ mhalf :+ (-s)*sqrt3d2)
+
+cache4 sign =
+   let s = Sign.toSign sign
+   in  (A.lift $ 0 :+ s, -1, A.lift $ 0 :+ (-s))
+
+cache5 sign =
+   let z = Sign.cisRat sign 5
+   in  (z 1, z 2, z 3, z 4)
+
+
+flatten2 ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array sh (a,a)) ->
+   Acc (Array (sh:.Int) a)
+flatten2 x =
+   A.generate
+      (Exp.indexCons (A.shape x) (A.constant 2))
+      (Exp.modify (expr :. expr) $
+       \(ix :. k)  ->  let xi = x ! ix in k ==* 0 ? (A.fst xi, A.snd xi))
+
+transform2 ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Exp a -> Transform (sh:.Int) a
+transform2 z arr =
+   flatten2 $
+   A.zipWith (\x0 x1 -> A.lift (x0+x1, x0+z*x1))
+      (A.slice arr (A.lift $ A.Any :. (0::Int)))
+      (A.slice arr (A.lift $ A.Any :. (1::Int)))
+
+
+flatten3 ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array sh (a,a,a)) ->
+   Acc (Array (sh:.Int) a)
+flatten3 x =
+   A.generate
+      (Exp.indexCons (A.shape x) (A.constant (3::Int)))
+      (Exp.modify (expr :. expr) $
+       \(ix :. k)  ->
+          let (x0,x1,x2) = A.unlift $ x ! ix
+          in  flip (A.caseof k) x0 $
+                 ((==*1), x1) :
+                 ((==*2), x2) :
+                 [])
+
+transform3 ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   (Exp a, Exp a) -> Transform (sh:.Int) a
+transform3 (z,z2) arr =
+   flatten3 $
+   A.zipWith3
+      (\x0 x1 x2 ->
+         let ((s,_), (zx1,zx2)) = Cyclic.sumAndConvolvePair (x1,x2) (z,z2)
+         in  A.lift (x0+s, x0+zx1, x0+zx2))
+      (A.slice arr (A.lift $ A.Any :. (0::Int)))
+      (A.slice arr (A.lift $ A.Any :. (1::Int)))
+      (A.slice arr (A.lift $ A.Any :. (2::Int)))
+
+
+flatten4 ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array sh (a,a,a,a)) ->
+   Acc (Array (sh:.Int) a)
+flatten4 x =
+   A.generate
+      (Exp.indexCons (A.shape x) (A.constant (4::Int)))
+      (Exp.modify (expr :. expr) $
+       \(ix :. k)  ->
+          let (x0,x1,x2,x3) = A.unlift $ x ! ix
+          in  flip (A.caseof k) x0 $
+                 ((==*1), x1) :
+                 ((==*2), x2) :
+                 ((==*3), x3) :
+                 [])
+
+transform4 ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   (Exp a, Exp a, Exp a) -> Transform (sh:.Int) a
+transform4 (z,z2,z3) arr =
+   flatten4 $
+   A.zipWith4
+      (\x0 x1 x2 x3 ->
+         let x02a = x0+x2; x02b = x0+z2*x2
+             x13a = x1+x3; x13b = x1+z2*x3
+         in  A.lift (x02a+   x13a, x02b+z *x13b,
+                     x02a+z2*x13a, x02b+z3*x13b))
+      (A.slice arr (A.lift $ A.Any :. (0::Int)))
+      (A.slice arr (A.lift $ A.Any :. (1::Int)))
+      (A.slice arr (A.lift $ A.Any :. (2::Int)))
+      (A.slice arr (A.lift $ A.Any :. (3::Int)))
+
+
+flatten5 ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array sh (a,a,a,a,a)) ->
+   Acc (Array (sh:.Int) a)
+flatten5 x =
+   A.generate
+      (Exp.indexCons (A.shape x) (A.constant (5::Int)))
+      (Exp.modify (expr :. expr) $
+       \(ix :. k)  ->
+          let (x0,x1,x2,x3,x4) = A.unlift $ x ! ix
+          in  flip (A.caseof k) x0 $
+                 ((==*1), x1) :
+                 ((==*2), x2) :
+                 ((==*3), x3) :
+                 ((==*4), x4) :
+                 [])
+
+
+{-
+Use Rader's trick for mapping the transform to a convolution
+and apply Karatsuba's trick at two levels (i.e. total three times)
+to that convolution.
+
+0 0 0 0 0
+0 1 2 3 4
+0 2 4 1 3
+0 3 1 4 2
+0 4 3 2 1
+
+Permutation.T: 0 1 2 4 3
+
+0 0 0 0 0
+0 1 2 4 3
+0 2 4 3 1
+0 4 3 1 2
+0 3 1 2 4
+-}
+transform5 ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   (Exp a, Exp a, Exp a, Exp a) -> Transform (sh:.Int) a
+transform5 (z1,z2,z3,z4) arr =
+   flatten5 $
+   A.zipWith5
+      (\x0 x1 x2 x3 x4 ->
+         let ((s,_), (d1,d2,d4,d3)) =
+                Cyclic.sumAndConvolveQuadruple (x1,x3,x4,x2) (z1,z2,z4,z3)
+         in  A.lift (x0+s, x0+d1, x0+d2, x0+d3, x0+d4))
+      (A.slice arr (A.lift $ A.Any :. (0::Int)))
+      (A.slice arr (A.lift $ A.Any :. (1::Int)))
+      (A.slice arr (A.lift $ A.Any :. (2::Int)))
+      (A.slice arr (A.lift $ A.Any :. (3::Int)))
+      (A.slice arr (A.lift $ A.Any :. (4::Int)))
+
+
+twist ::
+   (Shape sh, Slice sh, Elt a) =>
+   Exp Int -> Transform (sh:.Int:.Int) a
+twist fac x =
+   let sh :. m :. n = Exp.unlift (expr :. expr :. expr) $ A.shape x
+   in  A.backpermute
+          (A.lift $ sh :. fac*m :. div n fac)
+          (Exp.modify (expr :. expr :. expr) $
+           \(globalIx :. k :. j) -> globalIx :. div k fac :. fac*j + mod k fac)
+          x
+
+
+merge ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array (sh:.Int:.Int) a) ->
+   Acc (Array (sh:.Int) a)
+merge x =
+   let sh :. m :. n = Exp.unlift (expr :. expr :. expr) $ A.shape x
+   in  A.backpermute
+          (A.lift $ sh :. m*n)
+          (Exp.modify (expr :. expr) $
+           \(ix :. k) -> ix :. mod k m :. div k m)
+          x
+
+stack ::
+   (Shape sh, Slice sh, Elt a) =>
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int) a) ->
+   Acc (Array (sh:.Int:.Int) a)
+stack x y =
+   A.generate
+      (Exp.modify (expr :. expr)
+         (\(sh :. n) -> sh :. (2::Int) :. n)
+         (A.shape x))
+      (Exp.modify (expr :. expr :. expr) $
+       \(globalIx :. evenOdd :. k) ->
+          let ix = A.lift $ globalIx :. k
+          in  evenOdd ==* 0 ? (x ! ix, y ! ix))
+
+
+{- |
+twiddle factors for radix-2 Cooley-Tukey transforms
+-}
+twiddleFactors2 ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Acc (A.Vector (Complex a))
+twiddleFactors2 sign len2 =
+   A.generate (A.lift $ Z:.len2) $ twiddle2 sign len2 . A.indexHead
+
+twiddle2 ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Exp (Complex a)
+twiddle2 sign n2i ki =
+   let n2 = A.fromIntegral n2i
+       k  = A.fromIntegral ki
+   in  Sign.cis sign $ pi*k/n2
+
+
+twiddleFactors ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Acc (Array DIM2 (Complex a))
+twiddleFactors sign lenk lenj =
+   A.generate (A.lift $ Z:.lenk:.lenj) $
+   Exp.modify (expr :. expr :. expr) $
+   \(_z :. k :. j) -> twiddle sign (lenk*lenj) k j
+
+twiddle ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Exp Int -> Exp (Complex a)
+twiddle sign n k j =
+   Sign.cisRat sign n $ mod (k*j) n
+
+
+transformRadix2InterleavedTime ::
+   (Shape sh, Slice sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Acc (Array DIM1 a) ->
+   Transform (sh:.Int:.Int) a ->
+   Transform (sh:.Int) a
+transformRadix2InterleavedTime twiddles subTransform arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+       len2 = div len 2
+       subs =
+          subTransform $
+          if True
+            then Sliced.sliceHorizontal (A.lift $ Z:.(2::Int):.len2) arr
+            else
+               LinAlg.transpose $
+               A.reshape (A.lift $ sh:.len2:.(2::Int)) arr
+       evens = A.slice subs (A.lift $ A.Any :. (0::Int) :. A.All)
+       odds =
+          zipExtrudedVectorWith (*) twiddles $
+          A.slice subs (A.lift $ A.Any :. (1::Int) :. A.All)
+   in  A.zipWith (+) evens odds  A.++  A.zipWith (-) evens odds
+
+
+initSplitRadix ::
+   (Slice sh, Shape sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Acc (Array (sh:.Int) a) ->
+   (Acc (Array (sh:.Int:.Int) a), Acc (Array (sh:.Int:.Int) a))
+initSplitRadix arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+   in  (A.replicate (A.lift $ A.Any :. (1::Int) :. A.All) arr,
+        A.fill (A.lift $ sh:.(0::Int):.div len 2) 0)
+
+finishSplitRadix ::
+   (Slice sh, Shape sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Acc (Array (sh:.Int:.Int) a) -> Acc (Array (sh:.Int) a)
+finishSplitRadix =
+   flip A.slice (A.lift $ A.Any :. (0::Int) :. A.All)
+
+
+initSplitRadixFlat ::
+   (Slice sh, Shape sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Acc (Array (sh:.Int) a) ->
+   (Acc (Array DIM2 a), Acc (Array DIM2 a))
+initSplitRadixFlat arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+   in  (A.reshape (A.lift $ Z :. A.shapeSize sh :. len) arr,
+        A.fill (A.lift $ Z:.(0::Int):.div len 2) 0)
+
+finishSplitRadixFlat ::
+   (Slice sh, Shape sh, a ~ Complex b, IsFloating b, Elt b) =>
+   Exp (sh:.Int) -> Acc (Array DIM2 a) -> Acc (Array (sh:.Int) a)
+finishSplitRadixFlat = A.reshape
+
+
+imagSplitRadixPlain ::
+   (Elt a, IsNum a) =>
+   Sign a -> Complex a
+imagSplitRadixPlain sign = 0 :+ Sign.getSign sign
+
+imagSplitRadix ::
+   (Elt a, IsNum a) =>
+   Exp (Sign a) -> Exp (Complex a)
+imagSplitRadix sign =
+   A.lift (0 :+ Sign.toSign sign)
+
+ditSplitRadixReorder ::
+   (Slice sh, Shape sh, Elt a) =>
+   PairTransform (sh:.Int:.Int) a
+ditSplitRadixReorder (arr2, arr1) =
+   let evens = Sliced.sieve 2 0 arr2
+       odds  = Sliced.sieve 2 1 arr2
+   in  (Sliced1.append evens arr1, twist 2 odds)
+
+ditSplitRadixBase ::
+   (Slice sh, Shape sh, Elt a, IsFloating a) =>
+   PairTransform (sh:.Int:.Int) (Complex a)
+ditSplitRadixBase (arr2, arr1) = (transform2 (-1) arr2, arr1)
+
+ditSplitRadixStep ::
+   (Slice sh, Shape sh, a ~ Complex b, Elt b, IsFloating b) =>
+   Exp a ->
+   (Acc (Array DIM1 a), Acc (Array DIM1 a)) ->
+   PairTransform (sh:.Int:.Int) a
+ditSplitRadixStep imag (twiddles1, twiddles3) (u, zIntl) =
+   let twiddledZEven =
+          zipExtrudedVectorWith (*) twiddles1 $ Sliced1.sieve 2 0 zIntl
+       twiddledZOdd =
+          zipExtrudedVectorWith (*) twiddles3 $ Sliced1.sieve 2 1 zIntl
+       zSum = A.zipWith (+) twiddledZEven twiddledZOdd
+       zDiff = A.map (imag *) $ A.zipWith (-) twiddledZEven twiddledZOdd
+       zComplete = zSum A.++ zDiff
+   in  (A.zipWith (+) u zComplete
+        A.++
+        A.zipWith (-) u zComplete,
+           Sliced1.drop (Sliced1.length zComplete) u)
+
+
+twiddleSR ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Exp Int -> Exp (Complex a)
+twiddleSR sign n4i ki ji =
+   let n4 = A.fromIntegral n4i
+       k  = A.fromIntegral ki
+       j  = A.fromIntegral ji
+   in  Sign.cis sign $ pi*(k*j)/(2*n4)
+
+twiddleFactorsSR ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Acc (Array DIM1 (Complex a))
+twiddleFactorsSR sign len4 k =
+   A.generate (A.lift $ Z:.len4) $ twiddleSR sign len4 k . A.indexHead
+
+twiddleFactorsSRPair ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int ->
+   (Acc (Array DIM1 (Complex a)), Acc (Array DIM1 (Complex a)))
+twiddleFactorsSRPair sign len4 =
+   (twiddleFactorsSR sign len4 1,
+    twiddleFactorsSR sign len4 3)
+
+
+cycleDim3 ::
+   (Slice sh, Shape sh, Elt a) =>
+   Acc (Array (sh:.Int:.Int:.Int) a) ->
+   Acc (Array (sh:.Int:.Int:.Int) a)
+cycleDim3 arr =
+   A.backpermute
+      (Exp.modify (expr:.expr:.expr:.expr)
+         (\(sh:.k:.m:.n) -> (sh:.n:.k:.m)) $
+       A.shape arr)
+      (Exp.modify (expr:.expr:.expr:.expr)
+         (\(ix:.n:.k:.m) -> (ix:.k:.m:.n)))
+      arr
+
+
+chirp ::
+   (Elt a, IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp a -> A.Acc (A.Array DIM1 (Complex a))
+chirp sign padLen lenFloat =
+   A.generate (A.index1 padLen) $
+   \ix ->
+      let k = A.unindex1 ix
+          sk = A.fromIntegral (2*k <* padLen ? (k, k-padLen))
+      in  Sign.cis sign (pi*sk*sk/lenFloat)
diff --git a/src/Data/Array/Accelerate/Fourier/Real.hs b/src/Data/Array/Accelerate/Fourier/Real.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Real.hs
@@ -0,0 +1,239 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{- |
+Compute transforms on real data based on complex-valued transforms.
+-}
+module Data.Array.Accelerate.Fourier.Real (
+   toSpectrum,
+   fromSpectrum,
+
+   twoToSpectrum,
+   twoToSpectrum2d,
+   untangleSpectra,
+   untangleSpectra2d,
+   untangleCoefficient,
+
+   twoFromSpectrum,
+   twoFromSpectrum2d,
+   entangleSpectra,
+   entangleSpectra2d,
+   entangleCoefficient,
+   ) where
+
+import qualified Data.Array.Accelerate.Fourier.Sign as Sign
+import qualified Data.Array.Accelerate.Fourier.Private as Fourier
+import qualified Data.Array.Accelerate.Cyclic as Cyclic
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import Data.Array.Accelerate.Arithmetic.LinearAlgebra (zipExtrudedVectorWith, )
+
+import qualified Data.Array.Accelerate.Data.Complex as Complex
+import Data.Array.Accelerate.Data.Complex (Complex((:+)), )
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Acc, Array, Exp, Elt, IsFloating, Slice, Shape, (:.)((:.)),
+           (!), (?), (==*), )
+
+
+{- |
+Perform a real-to-complex transform
+using a complex-to-complex transform of half size.
+Input must have an even size.
+Result has the same size as the input, i.e. it is not halved.
+-}
+toSpectrum ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int) (Complex a) ->
+   Acc (Array (sh:.Int) a) -> Acc (Array (sh:.Int) (Complex a))
+toSpectrum subTrans arr =
+   let n2 = div (Sliced.length arr) 2
+       x = subTrans $ complexDeinterleave arr
+       xp = A.map (/2) $ A.zipWith (+) x (Cyclic.reverse x)
+       xm = A.map (/2) $ A.zipWith (-) x (Cyclic.reverse x)
+       twiddles =
+          A.map (imagUnit*) $
+          Fourier.twiddleFactors2 Sign.forwardExp n2
+       evens =
+          A.zipWith
+             (\xpk xmk -> A.lift $ Complex.real xpk :+ Complex.imag xmk) xp xm
+       odds  =
+          zipExtrudedVectorWith (*) twiddles $
+          A.zipWith
+             (\xpk xmk -> A.lift $ Complex.real xmk :+ Complex.imag xpk) xp xm
+   in  A.zipWith (-) evens odds
+       A.++
+       A.zipWith (+) evens odds
+
+complexDeinterleave ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int) a) -> Acc (Array (sh:.Int) (Complex a))
+complexDeinterleave arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+   in  A.generate
+          (A.lift $ sh :. div len 2)
+          (Exp.modify (expr:.expr) $
+           \(ix:.j) ->
+             arr ! A.lift (ix:.2*j)
+             :+
+             arr ! A.lift (ix:.2*j+1))
+
+
+{- |
+Perform a complex-to-real transform
+using a complex-to-complex of half size.
+Input must be self-adjoint and must have an even size.
+Result has the same size as the input, i.e. it is not doubled.
+-}
+fromSpectrum ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int) (Complex a) ->
+   Acc (Array (sh:.Int) (Complex a)) -> Acc (Array (sh:.Int) a)
+fromSpectrum subTrans spec =
+   let n2 = div (Sliced.length spec) 2
+       twiddles =
+          A.map (imagUnit*) $
+          Fourier.twiddleFactors2 Sign.inverseExp n2
+       part0 = Sliced.take n2 spec
+       part1 = Sliced.drop n2 spec
+       fe = A.zipWith (+) part0 part1
+       fo =
+          zipExtrudedVectorWith (*) twiddles $
+          A.zipWith (-) part0 part1
+   in  complexInterleave $ subTrans $ A.zipWith (+) fe fo
+
+complexInterleave ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int) (Complex a)) -> Acc (Array (sh:.Int) a)
+complexInterleave arr =
+   let (sh:.len) = Exp.unlift (expr:.expr) $ A.shape arr
+   in  A.generate
+          (A.lift $ sh :. 2*len)
+          (Exp.modify (expr:.expr) $
+           \(ix:.j) ->
+             let k = div j 2
+                 r = mod j 2
+                 x = arr ! A.lift (ix:.k)
+             in  r==*0 ? (Complex.real x, Complex.imag x))
+
+
+{- |
+Perform a real-to-complex transform of two real inputs
+using a complex-to-complex transform of the same size.
+Input can have arbitrary size.
+-}
+twoToSpectrum ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int) (Complex a) ->
+   Acc (Array (sh:.Int) (a,a)) ->
+   Acc (Array (sh:.Int) (Complex a, Complex a))
+twoToSpectrum subTrans =
+   untangleSpectra . subTrans .
+   A.map (Exp.modify (expr,expr) $ uncurry (:+))
+
+twoToSpectrum2d ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int:.Int) (Complex a) ->
+   Acc (Array (sh:.Int:.Int) (a,a)) ->
+   Acc (Array (sh:.Int:.Int) (Complex a, Complex a))
+twoToSpectrum2d subTrans =
+   untangleSpectra2d . subTrans .
+   A.map (Exp.modify (expr,expr) $ uncurry (:+))
+
+{- |
+You can transform two real data sets using one complex transform.
+This function can be used to untangle the resulting spectrum.
+-}
+{-
+Let f and g be two real valued images.
+The spectrum of f+i*g is spec f + i * spec g.
+Let 'flip' be the spectrum with negated indices modulo image size.
+It holds: flip (spec f) = conj (spec f).
+
+(a + conj b) / 2
+  = (spec (f+i*g) + conj (flip (spec (f+i*g)))) / 2
+  = (spec f + i*spec g + conj (flip (spec f)) + conj (flip (spec (i*g)))) / 2
+  = (2*spec f + i*spec g + conj (i*flip (spec g))) / 2
+  = (2*spec f + i*spec g - i * conj (flip (spec g))) / 2
+  = spec f
+
+(a - conj b) * (-i/2)
+  = (-i*a + conj (-i*b)) / 2
+  -> this swaps role of f and g in the proof above
+-}
+untangleSpectra ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int) (Complex a)) ->
+   Acc (Array (sh:.Int) (Complex a, Complex a))
+untangleSpectra spec =
+   A.zipWith untangleCoefficient spec (Cyclic.reverse spec)
+
+untangleSpectra2d ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int:.Int) (Complex a)) ->
+   Acc (Array (sh:.Int:.Int) (Complex a, Complex a))
+untangleSpectra2d spec =
+   A.zipWith untangleCoefficient spec (Cyclic.reverse2d spec)
+
+untangleCoefficient ::
+   (IsFloating a, Elt a) =>
+   Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a, Complex a)
+untangleCoefficient a b =
+   let bc = Complex.conjugate b
+   in  A.lift ((a + bc) / 2, (a - bc) * (-imagUnit / 2))
+
+
+twoFromSpectrum ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int) (Complex a) ->
+   Acc (Array (sh:.Int) (Complex a, Complex a)) ->
+   Acc (Array (sh:.Int) (a,a))
+twoFromSpectrum subTrans =
+   A.map (Exp.modify (expr:+expr) $ \(x:+y) -> (x,y)) .
+   subTrans . entangleSpectra
+
+twoFromSpectrum2d ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Fourier.Transform (sh:.Int:.Int) (Complex a) ->
+   Acc (Array (sh:.Int:.Int) (Complex a, Complex a)) ->
+   Acc (Array (sh:.Int:.Int) (a,a))
+twoFromSpectrum2d subTrans =
+   A.map (Exp.modify (expr:+expr) $ \(x:+y) -> (x,y)) .
+   subTrans . entangleSpectra2d
+
+entangleSpectra ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int) (Complex a, Complex a)) ->
+   Acc (Array (sh:.Int) (Complex a))
+entangleSpectra = entangleSpectraGen
+
+entangleSpectra2d ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array (sh:.Int:.Int) (Complex a, Complex a)) ->
+   Acc (Array (sh:.Int:.Int) (Complex a))
+entangleSpectra2d = entangleSpectraGen
+
+entangleSpectraGen ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Acc (Array sh (Complex a, Complex a)) ->
+   Acc (Array sh (Complex a))
+entangleSpectraGen = A.map (A.fst . A.uncurry entangleCoefficient)
+
+
+{-
+2 *c = a + bc     a  = c + i*d
+2i*d = a - bc     bc = c - i*d
+-}
+entangleCoefficient ::
+   (IsFloating a, Elt a) =>
+   Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a, Complex a)
+entangleCoefficient c d =
+   let di = d * imagUnit
+   in  A.lift (c + di, Complex.conjugate (c - di))
+
+
+imagUnit :: (A.Elt a, A.IsNum a) => Exp (Complex a)
+imagUnit = A.constant $ 0 :+ 1
diff --git a/src/Data/Array/Accelerate/Fourier/Sign.hs b/src/Data/Array/Accelerate/Fourier/Sign.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Sign.hs
@@ -0,0 +1,77 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+module Data.Array.Accelerate.Fourier.Sign where
+
+import Data.Array.Accelerate.Data.Complex (Complex((:+)), )
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate (Lift(lift), Unlift(unlift), Plain, )
+import Data.Array.Accelerate.Smart (Exp(Exp), PreExp(Tuple, Prj), )
+import Data.Array.Accelerate.Tuple
+          (IsTuple(TupleRepr, fromTuple, toTuple),
+           Tuple(NilTup, SnocTup), TupleIdx(ZeroTupIdx), )
+import Data.Array.Accelerate.Array.Sugar
+          (Elt(eltType, toElt, fromElt, eltType', toElt', fromElt'),
+           EltRepr, EltRepr', )
+
+import Data.Typeable (Typeable, )
+
+import qualified Test.QuickCheck as QC
+
+
+newtype Sign a = Sign {getSign :: a}
+   deriving (Eq, Show, Typeable)
+
+type instance EltRepr  (Sign a) = EltRepr  a
+type instance EltRepr' (Sign a) = EltRepr' a
+
+instance Elt a => Elt (Sign a) where
+   eltType = eltType . getSign
+   toElt   = Sign . toElt
+   fromElt = fromElt . getSign
+
+   eltType' = eltType' . getSign
+   toElt'   = Sign . toElt'
+   fromElt' = fromElt' . getSign
+
+instance IsTuple (Sign a) where
+   type TupleRepr (Sign a) = ((), a)
+   fromTuple (Sign a) = ((), a)
+   toTuple ((), a)    = Sign a
+
+instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Sign a) where
+   type Plain (Sign a) = Sign (Plain a)
+   lift (Sign a) = Exp $ Tuple (NilTup `SnocTup` lift a)
+
+instance Elt a => Unlift Exp (Sign (Exp a)) where
+   unlift e = Sign $ Exp $ ZeroTupIdx `Prj` e
+
+
+forward, inverse :: Num a => Sign a
+forward = Sign (-1)
+inverse = Sign 1
+
+forwardExp, inverseExp :: (Elt a, A.IsNum a) => Exp (Sign a)
+forwardExp = lift $ Sign $ A.fromIntegral (-1 :: Exp Int)
+inverseExp = lift $ Sign $ A.fromIntegral ( 1 :: Exp Int)
+
+toSign :: (Elt a) => Exp (Sign a) -> Exp a
+toSign = getSign . unlift
+
+cis ::
+   (Elt a, A.IsFloating a) =>
+   Exp (Sign a) -> Exp a -> Exp (Complex a)
+cis sign w  =  A.lift $ cos w :+ toSign sign * sin w
+
+cisRat ::
+   (Elt a, A.IsFloating a) =>
+   Exp (Sign a) -> Exp Int -> Exp Int -> Exp (Complex a)
+cisRat sign denom numer =
+   cis sign $ 2*pi * A.fromIntegral numer / A.fromIntegral denom
+
+
+instance (Num a) => QC.Arbitrary (Sign a) where
+   arbitrary = QC.elements [forward, inverse]
diff --git a/src/Data/Array/Accelerate/Fourier/Utility.hs b/src/Data/Array/Accelerate/Fourier/Utility.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Fourier/Utility.hs
@@ -0,0 +1,28 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.Fourier.Utility (
+   scaleDown,
+   ) where
+
+import Data.Array.Accelerate.Fourier.Private (Transform, )
+
+import Data.Array.Accelerate.Data.Complex (Complex((:+)), )
+
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Exp, Elt, IsNum, IsFloating, Slice, Shape, (:.), )
+
+
+scaleDown ::
+   (Shape sh, Slice sh, Elt a, IsFloating a) =>
+   Transform (sh:.Int) (Complex a)
+scaleDown zs =
+   A.map (cscale (recip $ A.fromIntegral $ Sliced.length zs)) zs
+
+cscale ::
+   (IsNum a, Elt a) =>
+   Exp a -> Exp (Complex a) -> Exp (Complex a)
+cscale x z =
+   case A.unlift z of
+      r :+ i -> A.lift (x*r :+ x*i)
diff --git a/src/Data/Array/Accelerate/NumberTheory.hs b/src/Data/Array/Accelerate/NumberTheory.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/NumberTheory.hs
@@ -0,0 +1,124 @@
+-- duplicate of synthesizer-core:NumberTheory
+module Data.Array.Accelerate.NumberTheory where
+
+import qualified Data.List.HT as ListHT
+import Data.Maybe.HT (toMaybe, )
+import Data.Bits ((.&.), (.|.), shiftR, )
+import Data.List (unfoldr, )
+
+
+{- |
+List all factorizations of an odd number
+where the first factor is at most the second factor
+and the first factors are in descending order.
+-}
+fermatFactors :: Integer -> [(Integer,Integer)]
+fermatFactors n =
+   let root = floor $ sqrt (fromInteger n :: Double)
+   in  map (\(a,b) -> (b-a,b+a)) $
+       mergeAndFilter
+          (zip (scanl (+) n [1,3..]) [0 .. div (n-1) 2])
+          (zip (scanl (+) (root*root) $ iterate (2+) (2*root+1)) [root..])
+
+mergeAndFilter :: (Ord a) => [(a,b)] -> [(a,c)] -> [(b,c)]
+mergeAndFilter ((a0,b):a0s) ((a1,c):a1s) =
+   case compare a0 a1 of
+      LT -> mergeAndFilter a0s ((a1,c):a1s)
+      GT -> mergeAndFilter ((a0,b):a0s) a1s
+      EQ -> (b,c) : mergeAndFilter a0s a1s
+mergeAndFilter _ _ = []
+
+
+multiplicativeGenerator :: Integer -> Integer
+multiplicativeGenerator p =
+   head $ primitiveRootsOfUnity p (p-1)
+
+primitiveRootsOfUnity :: Integer -> Integer -> [Integer]
+primitiveRootsOfUnity modu order =
+   let greatDivisors = map (div order) $ uniquePrimeFactors order
+   in  filter
+          (\n ->
+             let pow y = modularPower modu y n
+             in  coprime n modu
+                 &&
+                 pow order == 1
+                 &&
+                 all (\y -> pow y /= 1) greatDivisors) $
+       [1 .. modu-1]
+
+
+coprime :: Integer -> Integer -> Bool
+coprime x y  =  gcd x y == 1
+
+modularPower :: Integer -> Integer -> Integer -> Integer
+modularPower modu =
+   let go 0 _ = 1
+       go expo n =
+          case divMod expo 2 of
+             (expo2, r) ->
+                let n2 = mod (n*n) modu
+                in  if r==0
+                      then go expo2 n2
+                      else mod (go expo2 n2 * n) modu
+   in  go
+
+uniquePrimeFactors :: Integer -> [Integer]
+uniquePrimeFactors n =
+   let oddFactors =
+          foldr
+             (\p go m ->
+                let (q,r) = divMod m p
+                in  if r==0
+                      then p : go (divideByMaximumPower p q)
+                      else
+                        if q >= p
+                          then go m
+                          else if m==1 then [] else m : [])
+             (error "uniquePrimeFactors: end of infinite list")
+             (iterate (2+) 3)
+   in  case powerOfTwoFactors n of
+          (1,m) -> oddFactors m
+          (_,m) -> 2 : oddFactors m
+
+divideByMaximumPower :: Integer -> Integer -> Integer
+divideByMaximumPower b n =
+   last $
+   n : unfoldr (\m -> case divMod m b of (q,r) -> toMaybe (r==0) (q,q)) n
+
+powerOfTwoFactors :: Integer -> (Integer, Integer)
+powerOfTwoFactors n =
+   let powerOfTwo = n .&. (-n)
+   in  (powerOfTwo, div n powerOfTwo)
+
+
+ceilingPowerOfTwo :: Integer -> Integer
+ceilingPowerOfTwo 0 = 1
+ceilingPowerOfTwo n =
+   (1+) $ fst $ head $
+   dropWhile (uncurry (/=)) $
+   ListHT.mapAdjacent (,) $
+   scanl (\m d -> shiftR m d .|. m) (n-1) $
+   iterate (2*) 1
+
+{-
+For every reasonable pair of powers of 3 and 5
+it computes the least power of 2,
+such that their product is above @n@.
+-}
+ceiling5Smooth :: Integer -> Integer
+ceiling5Smooth n =
+   minimum $ map (minimum . ceilingSmooths 2 5 n) $
+   ceilingSmooths 2 3 n $ ceilingPowerOfTwo n
+
+{- |
+@ceilingSmooths a b n m@
+replaces successively @a@ factors in @m@ by @b@ factors
+while keeping the product above @n@.
+-}
+ceilingSmooths :: Integer -> Integer -> Integer -> Integer -> [Integer]
+ceilingSmooths a b n =
+   let divMany k =
+          case divMod k a of
+             (q,r) -> if r==0 && q>=n then divMany q else k
+       go m  =  m : if mod m a == 0 then go $ divMany $ m*b else []
+   in  go
diff --git a/src/Data/Array/Accelerate/Permutation.hs b/src/Data/Array/Accelerate/Permutation.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/Permutation.hs
@@ -0,0 +1,80 @@
+{- |
+Permutations of signals as needed for Fast Fourier transforms.
+Most functions are independent of the Signal framework.
+We could move them as well to Synthesizer.Basic.
+-}
+module Data.Array.Accelerate.Permutation where
+
+import qualified Data.Array.Accelerate.NumberTheory as NumberTheory
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Exp, Acc, Array, DIM1, Elt, Z(Z), (:.)((:.)),
+           Slice, Shape, (!), )
+
+
+type T = Acc Plain
+type Plain = Array DIM1 Int
+
+
+apply ::
+   (Slice sh, Shape sh, Elt a) =>
+   T -> LinAlg.Vector sh a -> LinAlg.Vector sh a
+apply p xs =
+   A.generate
+      (Exp.modify2 (expr:.expr) (expr:.expr)
+         (\(_z:.n) (sh:._) -> (sh:.n))
+         (A.shape p) (A.shape xs)) $
+   Exp.modify (expr:.expr) $
+   \(ix:.k) -> xs ! A.lift (ix :. p ! A.index1 k)
+
+
+plainSize :: Plain -> Int
+plainSize arr =
+   case A.arrayShape arr of
+      _ :. n -> n
+
+size :: T -> Exp Int
+size = A.length
+
+
+{- |
+Beware of 0-based indices stored in the result vector.
+-}
+multiplicative :: Int -> Plain
+multiplicative ni =
+   let n = fromIntegral ni
+       gen = NumberTheory.multiplicativeGenerator n
+   in  A.fromList (Z :. ni-1) $
+       map (fromInteger . subtract 1) $
+       iterate (\x -> mod (gen * x) n) 1
+
+
+{- |
+We only need to compute the inverse permutation explicitly,
+because not all signal structures support write to arbitrary indices,
+thus Generic.Write does not support it.
+For strict StorableVector it would be more efficient
+to build the vector directly.
+
+It holds:
+
+> inverse . inverse == id
+-}
+inverse :: T -> T
+inverse perm =
+   A.permute (+)
+      (A.fill (A.shape perm) 0)
+      (A.index1 . (perm!))
+      (A.generate (A.shape perm) A.unindex1)
+
+reverse :: T -> T
+reverse perm =
+   A.backpermute
+      (A.shape perm)
+      (\ix -> A.index1 $ mod (- A.unindex1 ix) (A.unindex1 $ A.shape perm))
+      perm
diff --git a/test/Test.hs b/test/Test.hs
new file mode 100644
--- /dev/null
+++ b/test/Test.hs
@@ -0,0 +1,15 @@
+module Main where
+
+import qualified Test.Data.Array.Accelerate.Fourier as Fourier
+
+import Data.Tuple.HT (mapFst, )
+
+
+prefix :: String -> [(String, IO ())] -> [(String, IO ())]
+prefix msg =
+   map (mapFst (\str -> msg ++ "." ++ str))
+
+main :: IO ()
+main =
+   mapM_ (\(msg,io) -> putStr (msg++": ") >> io) $
+   prefix "Fourier" Fourier.tests
diff --git a/test/Test/Data/Array/Accelerate/Fourier.hs b/test/Test/Data/Array/Accelerate/Fourier.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Data/Array/Accelerate/Fourier.hs
@@ -0,0 +1,583 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Test.Data.Array.Accelerate.Fourier where -- (tests) where
+
+import qualified Data.Array.Accelerate.Fourier.Real as FourierReal
+import qualified Data.Array.Accelerate.Fourier.Preprocessed as Prep
+import qualified Data.Array.Accelerate.Fourier.Adhoc as Adhoc
+import qualified Data.Array.Accelerate.Fourier.Planned as Planned
+import qualified Data.Array.Accelerate.Convolution.Preprocessed as ConvPrep
+import qualified Data.Array.Accelerate.Convolution.Adhoc as Convolution
+import qualified Data.Array.Accelerate.Cyclic as Cyclic
+import qualified Data.Array.Accelerate.Interpreter as AI
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Fourier.Planned (Transform, )
+import Data.Array.Accelerate.Fourier.Utility (scaleDown, )
+import Data.Array.Accelerate
+          (Acc, Exp, Array, DIM1, DIM2, DIM3, Z(Z), (:.)((:.)),
+           (<=*), (==*),(&&*), )
+
+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
+import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck (Arbitrary, arbitrary, quickCheck, )
+
+import Control.Monad (liftM2, liftM3, guard, )
+
+import qualified Data.Array.Accelerate.Data.Complex as Complex
+import Data.Complex (Complex((:+)), cis, )
+
+
+tolerance :: Double
+tolerance = 1e-10
+
+approxEqualAbs ::
+   (A.Elt a, A.IsFloating a) =>
+   Exp a -> Exp a -> Exp a -> Exp Bool
+approxEqualAbs eps x y =
+   abs (x-y) <=* eps
+
+approxEqualComplexAbs ::
+   (A.Elt a, A.IsFloating a) =>
+   Exp a -> Exp (Complex a) -> Exp (Complex a) -> Exp Bool
+approxEqualComplexAbs eps x y =
+   Complex.magnitude (x-y) <=* eps
+
+
+genComplex :: QC.Gen (Complex Double)
+genComplex = liftM2 (:+) (QC.choose (-1,1)) (QC.choose (-1,1))
+
+newtype Normed0 = Normed0 (Exp (Complex Double))
+   deriving (Show)
+
+instance Arbitrary Normed0 where
+   arbitrary = fmap (Normed0 . A.constant) genComplex
+
+
+data Normed1 = Normed1 Int (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1 where
+   arbitrary =
+      fmap
+         (\xs ->
+            let len = length xs
+            in  Normed1 len $ A.use $ A.fromList (Z :. len) xs) $
+      QC.listOf genComplex
+
+
+floorPowerOfTwo :: Int -> Int
+floorPowerOfTwo len =
+   2 ^ (floor (logBase 2 (fromIntegral len :: Double)) :: Int)
+
+data Normed1PowerTwo = Normed1PowerTwo Int (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1PowerTwo where
+   arbitrary =
+      fmap
+         (\xs ->
+            let len = floorPowerOfTwo $ length xs
+            in  Normed1PowerTwo len $ A.use $ A.fromList (Z :. len) xs) $
+      liftM2 (:) genComplex (QC.listOf genComplex)
+
+
+{-
+For every reasonable pair of powers of 3 and 5
+it computes the largest power of 2,
+such that their product is below @n@.
+-}
+floor5Smooth :: Int -> Int
+floor5Smooth n =
+   fromInteger $
+   maximum $ map (maximum . floorSmooths 2 5 (fromIntegral n)) $
+   floorSmooths 2 3 (fromIntegral n) $
+   fromIntegral $ floorPowerOfTwo n
+
+{- |
+@floorSmooths a b n m@
+replaces successively @a@ factors in @m@ by @b@ factors
+while keeping the product below @n@.
+-}
+floorSmooths :: Integer -> Integer -> Integer -> Integer -> [Integer]
+floorSmooths a b n =
+   let divMany k =
+          case divMod k a of
+             (q,r) -> guard (r==0) >> if q>n then divMany q else go q
+       go m  =  m : divMany (m*b)
+   in  go
+
+
+data Normed1Smooth5 = Normed1Smooth5 Int (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1Smooth5 where
+   arbitrary =
+      fmap
+         (\xs ->
+            let len = floor5Smooth $ length xs
+            in  Normed1Smooth5 len $ A.use $ A.fromList (Z :. len) xs) $
+      liftM2 (:) genComplex (QC.listOf genComplex)
+
+
+data Normed1Even = Normed1Even Int (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1Even where
+   arbitrary =
+      fmap
+         (\xs ->
+            let len = 2 * length xs
+            in  Normed1Even len $ A.use $ A.fromList (Z :. len) $
+                concatMap (\(x0,x1) -> [x0,x1]) xs) $
+      QC.listOf $ liftM2 (,) genComplex genComplex
+
+
+data
+   Normed1Pair =
+      Normed1Pair Int
+         (Acc (Array DIM1 (Complex Double)))
+         (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1Pair where
+   arbitrary =
+      fmap
+         (\xys ->
+            let len = length xys
+                (xs, ys) = unzip xys
+            in  Normed1Pair len
+                   (A.use $ A.fromList (Z :. len) xs)
+                   (A.use $ A.fromList (Z :. len) ys)) $
+      QC.listOf $ liftM2 (,) genComplex genComplex
+
+
+data
+   Normed1Triple =
+      Normed1Triple Int
+         (Acc (Array DIM1 (Complex Double)))
+         (Acc (Array DIM1 (Complex Double)))
+         (Acc (Array DIM1 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed1Triple where
+   arbitrary =
+      fmap
+         (\xyzs ->
+            let len = length xyzs
+                (xs, ys, zs) = unzip3 xyzs
+            in  Normed1Triple len
+                   (A.use $ A.fromList (Z :. len) xs)
+                   (A.use $ A.fromList (Z :. len) ys)
+                   (A.use $ A.fromList (Z :. len) zs)) $
+      QC.listOf $ liftM3 (,,) genComplex genComplex genComplex
+
+
+data Normed2 = Normed2 Int Int (Acc (Array DIM2 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed2 where
+   arbitrary = do
+      xs <- liftM2 (:) QC.arbitrary $ QC.listOf genComplex
+      let len = length xs
+      height <- QC.choose (1, round (sqrt (fromIntegral len :: Double)))
+      let width = div len height
+      return $ Normed2 width height $
+         A.use $ A.fromList (Z :. height :. width) xs
+
+
+data Normed3 = Normed3 Int Int Int (Acc (Array DIM3 (Complex Double)))
+   deriving (Show)
+
+instance Arbitrary Normed3 where
+   arbitrary = do
+      xs <- liftM2 (:) QC.arbitrary $ QC.listOf genComplex
+      let len = length xs
+          lenThd = round $ (fromIntegral len :: Double) ** recip 3
+      height <- QC.choose (1, lenThd)
+      width <- QC.choose (1, lenThd)
+      let depth = div len (width*height)
+      return $ Normed3 width height depth $
+         A.use $ A.fromList (Z :. depth :. height :. width) xs
+
+
+-- duplicate of Private.cycleDim3
+cycleDim3 :: (A.Elt a) => Transform DIM3 a
+cycleDim3 arr =
+   A.backpermute
+      (Exp.modify (expr:.expr:.expr:.expr)
+         (\(sh:.k:.m:.n) -> (sh:.n:.k:.m)) $
+       A.shape arr)
+      (Exp.modify (expr:.expr:.expr:.expr)
+         (\(ix:.n:.k:.m) -> (ix:.k:.m:.n)))
+      arr
+
+
+basisVector :: Int -> Int -> Array DIM1 (Complex Double)
+basisVector len freq =
+   A.fromList (Z:.len) $
+   map (\k -> cis $ 2*pi * fromIntegral (k*freq) / fromIntegral len) $
+   iterate (1+) 0
+
+
+norm2 ::
+   Acc (Array DIM1 (Complex Double)) ->
+   Acc (A.Scalar Double)
+norm2  =
+   A.map sqrt . A.sum . A.map (Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i)
+
+scalarProduct ::
+   Acc (Array DIM1 (Complex Double)) ->
+   Acc (Array DIM1 (Complex Double)) ->
+   Acc (A.Scalar (Complex Double))
+scalarProduct xs ys =
+   A.foldAll (+) 0 $ A.zipWith (*) xs (A.map Complex.conjugate ys)
+
+complexFromReal ::
+   Acc (Array DIM1 Double) ->
+   Acc (Array DIM1 (Complex Double))
+complexFromReal = A.map (A.lift . (:+0))
+
+toSelfAdjoint :: Transform DIM1 (Complex Double)
+toSelfAdjoint x =
+   A.zipWith (+) x $
+   A.map Complex.conjugate $ Cyclic.reverse x
+
+
+
+infixl 6 <+>
+
+(<+>) ::
+   Acc (Array DIM1 (Complex Double)) ->
+   Acc (Array DIM1 (Complex Double)) ->
+   Acc (Array DIM1 (Complex Double))
+(<+>) = A.zipWith (+)
+
+
+
+{-
+should be replaced by (==*) in future
+-}
+class (A.Shape sh, A.Slice sh) => EqShape sh where
+   eqShape :: Exp sh -> Exp sh -> Exp Bool
+
+instance EqShape Z where
+   eqShape _ _ = A.constant True
+
+instance (EqShape sh, i ~ Int) => EqShape (sh:.i) where
+   eqShape =
+      Exp.modify2 (expr:.expr) (expr:.expr) $
+      \(sh0:.n0) (sh1:.n1) ->
+         n0 ==* n1  &&*  eqShape sh0 sh1
+
+
+infix 4 =~=
+
+(=~=) ::
+   (EqShape sh) =>
+   Acc (Array sh (Complex Double)) ->
+   Acc (Array sh (Complex Double)) ->
+   Acc (A.Scalar Bool)
+(=~=) xs ys =
+   A.map (eqShape (A.shape xs) (A.shape ys)  &&*) $
+   A.and $ A.zipWith (approxEqualComplexAbs (A.constant tolerance)) xs ys
+
+
+run :: Acc (A.Scalar Bool) -> Bool
+run = Acc.the . AI.run
+
+
+tests :: [(String, IO ())]
+tests =
+   ("fourier generic vs. preprocessed dit2",
+      quickCheck $ \sign (Normed1PowerTwo len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Prep.dit2 sign len x) :
+   ("fourier generic vs. preprocessed dif2",
+      quickCheck $ \sign (Normed1PowerTwo len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Prep.dif2 sign len x) :
+   ("fourier generic vs. preprocessed ditSplitRadix",
+      quickCheck $ \sign (Normed1PowerTwo len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Prep.ditSplitRadix sign len x) :
+   ("fourier generic vs. adhoc dit2",
+      quickCheck $ \sign (Normed1PowerTwo len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Adhoc.dit2 (A.constant sign) x) :
+   ("fourier generic vs. adhoc ditSplitRadix",
+      quickCheck $ \sign (Normed1PowerTwo len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Adhoc.ditSplitRadix (A.constant sign) x) :
+   ("fourier generic vs. adhoc dit235",
+      quickCheck $ \sign (Normed1Smooth5 len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Adhoc.dit235 (A.constant sign) x) :
+   ("fourier adhoc chirp 2 vs. chirp 235",
+      quickCheck $ \sign (Normed1 _len x) -> run $
+             Adhoc.transformChirp2 (A.constant sign) x
+             =~=
+             Adhoc.transformChirp235 (A.constant sign) x) :
+   ("fourier generic vs. adhoc auto",
+      quickCheck $ \sign (Normed1 len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Adhoc.transform (A.constant sign) x) :
+   ("fourier generic vs. adhoc chirp 235",
+      quickCheck $ \sign (Normed1 len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Adhoc.transformChirp235 (A.constant sign) x) :
+   ("fourier generic vs. chirp2",
+      quickCheck $ \sign (Normed1 len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Planned.transformChirp2 sign len x) :
+   ("fourier generic vs. chirp235",
+      quickCheck $ \sign (Normed1 len x) -> run $
+             Planned.transform sign len x
+             =~=
+             Planned.transformChirp235 sign len x) :
+   ("homogeneity",
+      quickCheck $ \sign (Normed0 x) (Normed1 len y) -> run $
+         let transform = Planned.transform sign len
+         in  transform (A.map (x*) y)
+             =~=
+             A.map (x*) (transform y)) :
+   ("additivity",
+      quickCheck $ \sign (Normed1Pair len x y) -> run $
+         let transform = Planned.transform sign len
+         in  A.zipWith (+) (transform x) (transform y)
+             =~=
+             transform (A.zipWith (+) x y)) :
+   ("basis vector",
+      quickCheck $ \(Normed1 len _x) kp -> run $
+         let transform = Planned.transform Planned.inverse len
+             k = mod kp len
+             unit =
+                A.use $ A.fromList (Z:.len) $
+                replicate k 0 ++ 1 : repeat 0
+         in  transform unit
+             =~=
+             A.use (basisVector len k)) :
+   ("fourier inverse",
+      quickCheck $ \(Normed1 len x) -> run $
+             x =~=
+             (scaleDown $
+              Planned.transform Planned.forward len $
+              Planned.transform Planned.inverse len x)) :
+   ("double fourier = reverse",
+      quickCheck $ \sign (Normed1 len x) -> run $
+         let transform = Planned.transform sign len
+         in  x =~=
+             (Cyclic.reverse $
+              scaleDown $
+              transform $
+              transform x)) :
+   ("fourier of reverse",
+      quickCheck $ \sign (Normed1 len x) -> run $
+         let transform = Planned.transform sign len
+         in  Cyclic.reverse (transform x) =~=
+             transform (Cyclic.reverse x)) :
+   ("fourier of conjugate",
+      quickCheck $ \sign (Normed1 len x) -> run $
+         let transform = Planned.transform sign len
+         in  (A.map Complex.conjugate $ transform x)
+             =~=
+             (transform $
+              A.map Complex.conjugate $ Cyclic.reverse x)) :
+   ("isometry",
+      quickCheck $ \sign (Normed1 len x) -> run $
+         let transform = Planned.transform sign len
+         in  A.zipWith
+                (approxEqualAbs $ A.constant tolerance)
+                (norm2 $ transform x)
+                (A.map (A.constant (sqrt (fromIntegral len)) *) $ norm2 x)) :
+   ("unitarity",
+      quickCheck $ \sign (Normed1Pair len x y) -> run $
+         let transform = Planned.transform sign len
+         in  A.zipWith
+                (approxEqualComplexAbs $ A.constant tolerance)
+                (scalarProduct (transform x) (transform y))
+                (A.map (A.constant (fromIntegral len) *) $
+                 scalarProduct x y)) :
+   ("convolution commutativity",
+      quickCheck $ \(Normed1Pair len x y) -> run $
+         let (<*>) = Planned.convolveCyclic len
+         in  x <*> y
+             =~=
+             y <*> x) :
+   ("convolution associativity",
+      quickCheck $ \(Normed1Triple len x y z) -> run $
+         let (<*>) = Planned.convolveCyclic len
+         in  (x <*> y) <*> z
+             =~=
+             x <*> (y <*> z)) :
+   ("convolution distributivity",
+      quickCheck $ \(Normed1Triple len x y z) -> run $
+         let (<*>) = Planned.convolveCyclic len
+         in  x <*> (y <+> z)
+             =~=
+             (x <*> y) <+> (x <*> z)) :
+   ("convolution karatsuba rec vs. loop",
+      quickCheck $ \(Normed1 len xy) -> run $
+         let x = A.map Complex.real xy
+             y = A.map Complex.imag xy
+         in  complexFromReal (ConvPrep.karatsuba len x y)
+             =~=
+             complexFromReal (Convolution.karatsuba x y)) :
+{-
+   No instance for (A.IsNum (Complex Double))
+      arising from a use of 'ConvPrep.karatsuba'
+-}
+   ("convolution karatsuba",
+      quickCheck $ \(Normed1Pair len x y) -> run $
+         let resultLen = max 0 $ 2*len-1
+         in  Convolution.complex Convolution.karatsuba x y
+             =~=
+             Planned.convolveCyclic resultLen
+                (Sliced.pad 0 (A.constant resultLen) x)
+                (Sliced.pad 0 (A.constant resultLen) y)) :
+   ("convolution cyclic karatsuba",
+      quickCheck $ \(Normed1Pair len x y) -> run $
+             Convolution.complex
+                (Convolution.cyclic Convolution.karatsuba) x y
+             =~=
+             Planned.convolveCyclic len
+                (Sliced.pad 0 (A.constant len) x)
+                (Sliced.pad 0 (A.constant len) y)) :
+   ("real to spectrum",
+      quickCheck $ \(Normed1Even len x) -> run $
+         let xr = A.map Complex.real x
+         in  FourierReal.toSpectrum
+                (Planned.transform Planned.forward (div len 2)) xr
+             =~=
+             Planned.transform Planned.forward len (complexFromReal xr)) :
+   ("real from spectrum",
+      quickCheck $ \(Normed1Even len x) -> run $
+         let xSelfAdjoint = toSelfAdjoint x
+         in  (complexFromReal $
+              FourierReal.fromSpectrum
+                 (Planned.transform Planned.inverse (div len 2)) xSelfAdjoint)
+             =~=
+             Planned.transform Planned.inverse len xSelfAdjoint) :
+   ("real to and from spectrum",
+      quickCheck $ \(Normed1Even len x) -> run $
+         let xr = A.map Complex.real x
+             len2 = div len 2
+         in  (scaleDown $ complexFromReal $
+              FourierReal.fromSpectrum (Planned.transform Planned.inverse len2) $
+              FourierReal.toSpectrum (Planned.transform Planned.forward len2) xr)
+             =~=
+             complexFromReal xr) :
+   ("real from and to spectrum",
+      quickCheck $ \(Normed1Even len x) -> run $
+         let len2 = div len 2
+             xSelfAdjoint = toSelfAdjoint x
+         in  (scaleDown $
+              FourierReal.toSpectrum (Planned.transform Planned.forward len2) $
+              FourierReal.fromSpectrum (Planned.transform Planned.inverse len2) $
+              xSelfAdjoint)
+             =~=
+             xSelfAdjoint) :
+   ("double real to spectrum, even",
+      quickCheck $ \(Normed1Even len x) -> run $
+         let xr = A.map Complex.real x
+             xi = A.map Complex.imag x
+             transform = Planned.transform Planned.forward (div len 2)
+             (specr,speci) =
+                A.unzip $ FourierReal.untangleSpectra $
+                Planned.transform Planned.forward len x
+         in  A.zipWith (&&*)
+                (FourierReal.toSpectrum transform xr =~= specr)
+                (FourierReal.toSpectrum transform xi =~= speci)) :
+   ("double real to spectrum, arbitrary",
+      quickCheck $ \(Normed1 len x) -> run $
+         let xr = complexFromReal $ A.map Complex.real x
+             xi = complexFromReal $ A.map Complex.imag x
+             transform = Planned.transform Planned.forward len
+             (specr,speci) =
+                A.unzip $ FourierReal.untangleSpectra $ transform x
+         in  A.zipWith (&&*)
+                (transform xr =~= specr)
+                (transform xi =~= speci)) :
+   ("entangle and untangle spectrum of real data",
+      quickCheck $ \(Normed1Pair _len x y) -> run $
+         let (xt,yt) =
+                A.unzip $
+                A.map (A.uncurry FourierReal.untangleCoefficient) $
+                A.zipWith FourierReal.entangleCoefficient x y
+         in  A.zipWith (&&*) (x =~= xt) (y =~= yt)) :
+   ("double real from spectrum",
+      quickCheck $ \(Normed1 len x) -> run $
+         let imagUnit = A.constant (0:+1)
+             xSelfAdjoint = toSelfAdjoint x
+             ySelfAdjoint = toSelfAdjoint $ A.map (imagUnit*) x
+             transform = Planned.transform Planned.inverse len
+             (xSignal,ySignal) =
+                A.unzip $ FourierReal.twoFromSpectrum transform $
+                A.zip xSelfAdjoint ySelfAdjoint
+         in  A.zipWith (&&*)
+                (transform xSelfAdjoint =~= A.map (A.lift . (:+0)) xSignal)
+                (transform ySelfAdjoint =~= A.map (A.lift . (:+0)) ySignal)) :
+   ("transform2d vs. transposition, preprocessed",
+      quickCheck $ \sign (Normed2 width height x) -> run $
+         let transformH =
+                Prep.transform2d
+                   (Prep.SubTransformPair
+                      (Planned.transform sign width)
+                      (Planned.transform sign height))
+             transformV =
+                Prep.transform2d
+                   (Prep.SubTransformPair
+                      (Planned.transform sign height)
+                      (Planned.transform sign width))
+         in  LinAlg.transpose (transformH x)
+             =~=
+             transformV (LinAlg.transpose x)) :
+   ("transform2d vs. transposition, adhoc",
+      quickCheck $ \sign (Normed2 _width _height x) -> run $
+         let transform =
+                Adhoc.transform2d
+                   (Adhoc.SubTransform
+                      (Adhoc.transformChirp2 (A.constant sign)))
+         in  LinAlg.transpose (transform x)
+             =~=
+             transform (LinAlg.transpose x)) :
+   ("transform3d vs. transposition, preprocessed",
+      quickCheck $ \sign (Normed3 width height depth x) -> run $
+         let transformH =
+                Prep.transform3d
+                   (Prep.SubTransformTriple
+                      (Planned.transform sign width)
+                      (Planned.transform sign height)
+                      (Planned.transform sign depth))
+             transformV =
+                Prep.transform3d
+                   (Prep.SubTransformTriple
+                      (Planned.transform sign height)
+                      (Planned.transform sign depth)
+                      (Planned.transform sign width))
+         in  cycleDim3 (transformH x)
+             =~=
+             transformV (cycleDim3 x)) :
+   ("transform2d vs. transposition, adhoc",
+      quickCheck $ \sign (Normed3 _width _height _depth x) -> run $
+         let transform =
+                Adhoc.transform2d
+                   (Adhoc.SubTransform
+                      (Adhoc.transformChirp2 (A.constant sign)))
+         in  LinAlg.transpose (transform x)
+             =~=
+             transform (LinAlg.transpose x)) :
+   []