diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,9 @@
 # Change log for the `patch-image` package
 
+## 0.3.1:
+
+ * Speed up computation by moving more stuff to Knead/LLVM.
+
 ## 0.3:
 
  * Allow to save program state, modify it manually
diff --git a/patch-image.cabal b/patch-image.cabal
--- a/patch-image.cabal
+++ b/patch-image.cabal
@@ -1,5 +1,5 @@
 Name:           patch-image
-Version:        0.3.0.1
+Version:        0.3.1
 License:        BSD3
 License-File:   LICENSE
 Author:         Henning Thielemann <haskell@henning-thielemann.de>
@@ -41,7 +41,7 @@
   README.md
 
 Source-Repository this
-  Tag:         0.3.0.1
+  Tag:         0.3.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/patch-image/
 
@@ -67,7 +67,9 @@
   Main-Is: Knead.hs
   Other-Modules:
     MatchImageBorders
-    KneadShape
+    Knead.CArray
+    Knead.Shape
+    Complex
     State
     LinearAlgebra
     Arithmetic
@@ -81,7 +83,7 @@
 
   If flag(llvm)
     Build-Depends:
-      knead >=0.2.1 && <0.3,
+      knead >=0.2.2 && <0.3,
       llvm-extra >=0.7 && <0.8,
       llvm-tf >=3.1 && <3.2,
       tfp >=1.0 && <1.1,
@@ -92,6 +94,7 @@
       enumset >=0.0.4 && <0.1,
       containers >=0.4.2 && <0.6,
       fft >=0.1.7 && <0.2,
+      storable-tuple >=0.0.3 && <0.1,
       carray >=0.1.5 && <0.2,
       array >=0.4 && <0.6,
       cassava >=0.4.5 && <0.5,
@@ -112,6 +115,7 @@
     State
     LinearAlgebra
     Arithmetic
+    Complex
     Degree
     Option.Utility
     Option
diff --git a/src/Accelerate.hs b/src/Accelerate.hs
--- a/src/Accelerate.hs
+++ b/src/Accelerate.hs
@@ -5,6 +5,7 @@
 import qualified State
 
 import qualified Arithmetic as Arith
+import qualified Complex as Komplex
 import qualified Degree
 import LinearAlgebra (
    absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,
@@ -724,11 +725,24 @@
 
 type GenDIM2 a = Z :. a :. a
 
-shrinkFactors :: (Integral a) => DIM2 -> GenDIM2 a -> GenDIM2 a -> GenDIM2 a
-shrinkFactors (Z:.heightPad:.widthPad)
+-- cf. Arithmetic.minimumOverlapAbsFromPortion
+minimumOverlapAbsFromPortion ::
+   (Num a, Ord i) => (a -> i, i -> a) -> a -> (i, i) -> i
+minimumOverlapAbsFromPortion
+   (afloor, fromInt) minOverlapPortion (width, height) =
+      afloor $ minOverlapPortion * fromInt (min width height)
+
+shrinkFactors ::
+   (Num a, Integral i) =>
+   (a -> i, i -> a) ->
+   DIM2 -> a -> GenDIM2 i -> GenDIM2 i -> GenDIM2 i
+shrinkFactors methods (Z:.heightPad:.widthPad) minOverlapPortion
    (Z :. heighta :. widtha) (Z :. heightb :. widthb) =
-      let yk = divUp (heighta+heightb) $ fromIntegral heightPad
-          xk = divUp (widtha +widthb)  $ fromIntegral widthPad
+      let minOverlap =
+            minimumOverlapAbsFromPortion methods minOverlapPortion
+               (min widtha widthb, min heighta heightb)
+          yk = divUp (heighta+heightb-minOverlap) $ fromIntegral heightPad
+          xk = divUp (widtha +widthb -minOverlap) $ fromIntegral widthPad
       in  Z :. yk :. xk
 
 {-
@@ -740,7 +754,8 @@
    let run =
           Run.with CUDA.run1 $ \minimumOverlap a b ->
              let factors@(_z:.yk:.xk) =
-                    shrinkFactors padExtent
+                    shrinkFactors (A.floor, A.fromIntegral) padExtent
+                       minimumOverlap
                        (A.unlift $ A.shape a) (A.unlift $ A.shape b)
                  scalePos =
                     Exp.modify (expr, (expr,expr)) $
@@ -792,7 +807,9 @@
           Run.with CUDA.run1 $ \minimumOverlap a b ->
              let shapeA = A.unlift $ A.shape a
                  shapeB = A.unlift $ A.shape b
-                 factors@(_z:.yk:.xk) = shrinkFactors padExtent shapeA shapeB
+                 factors@(_z:.yk:.xk) =
+                    shrinkFactors (A.floor, A.fromIntegral) padExtent
+                       minimumOverlap shapeA shapeB
                  coarsed@(coarsedx,coarsedy) =
                     mapPair ((xk*), (yk*)) $
                     Exp.unliftPair $ A.snd $ A.the $ argmaximum $
@@ -862,7 +879,8 @@
 
    in  \minimumOverlap mMaximumDiff a b ->
           let factors@(Z:.yk:.xk) =
-                 shrinkFactors padExtent (A.arrayShape a) (A.arrayShape b)
+                 shrinkFactors (floor, fromIntegral) padExtent
+                    minimumOverlap (A.arrayShape a) (A.arrayShape b)
 
               (_score, shrunkd@(shrunkdx, shrunkdy)) =
                  Acc.the $ overlapShrunk minimumOverlap factors a b
@@ -1597,7 +1615,7 @@
        map picColored pics,
        map
          (mapPair
-            (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac))
+            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))
          posRots)
 
 
diff --git a/src/Arithmetic.hs b/src/Arithmetic.hs
--- a/src/Arithmetic.hs
+++ b/src/Arithmetic.hs
@@ -1,5 +1,6 @@
 module Arithmetic where
 
+import qualified Complex as Komplex
 import qualified Data.Complex as Complex
 import Data.Complex (Complex, )
 
@@ -71,7 +72,7 @@
    let posRotPics =
          zipWith
             (\(angle,pic) (pos,rot) ->
-               (pos, (pairFromComplex (Complex.cis angle * rot), pic)))
+               (pos, (Komplex.toPair (Complex.cis angle * rot), pic)))
             picAngles floatPosRots
        bbox (rot, pic) =
          case extent pic of
@@ -307,16 +308,6 @@
 divUp :: (Integral a) => a -> a -> a
 divUp a b = - div (-a) b
 
-
-
-pairFromComplex :: (RealFloat a) => Complex a -> (a,a)
-pairFromComplex z = (Complex.realPart z, Complex.imagPart z)
-
-mapComplex :: (a -> b) -> Complex a -> Complex b
-mapComplex f (r Complex.:+ i)  =  f r Complex.:+ f i
-
-mulConj :: (RealFloat a) => Complex a -> Complex a -> Complex a
-mulConj x y = x * Complex.conjugate y
 
 
 -- ToDo: move to a new utility module
diff --git a/src/Complex.hs b/src/Complex.hs
new file mode 100644
--- /dev/null
+++ b/src/Complex.hs
@@ -0,0 +1,29 @@
+module Complex where
+
+import qualified Data.Complex as HComplex
+import Data.Complex (Complex((:+)))
+
+
+toPair :: (RealFloat a) => Complex a -> (a,a)
+toPair z = (HComplex.realPart z, HComplex.imagPart z)
+
+map :: (a -> b) -> Complex a -> Complex b
+map f (r :+ i)  =  f r :+ f i
+
+conjugate :: (Num a) => Complex a -> Complex a
+conjugate (r :+ i)  =  r :+ negate i
+
+add :: (Num a) => Complex a -> Complex a -> Complex a
+add (xr:+xi) (yr:+yi) = (xr+yr) :+ (xi+yi)
+
+sub :: (Num a) => Complex a -> Complex a -> Complex a
+sub (xr:+xi) (yr:+yi) = (xr-yr) :+ (xi-yi)
+
+mul :: (Num a) => Complex a -> Complex a -> Complex a
+mul (xr:+xi) (yr:+yi) = (xr*yr-xi*yi) :+ (xr*yi+xi*yr)
+
+mulConj :: (Num a) => Complex a -> Complex a -> Complex a
+mulConj x y = mul x $ conjugate y
+
+mulConj_ :: (RealFloat a) => Complex a -> Complex a -> Complex a
+mulConj_ x y = x * Complex.conjugate y
diff --git a/src/Knead.hs b/src/Knead.hs
--- a/src/Knead.hs
+++ b/src/Knead.hs
@@ -6,6 +6,8 @@
 
 import qualified MatchImageBorders
 import qualified Arithmetic as Arith
+import qualified Knead.CArray as KneadCArray
+import qualified Complex as Komplex
 import qualified Degree
 import MatchImageBorders (arrayCFromKnead, arrayKneadFromC)
 import Arithmetic (guardedPairs, maximum0)
@@ -13,7 +15,7 @@
    absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,
    layoutFromPairDisplacements, fixAtLeastOneAnglePosition,
    )
-import KneadShape
+import Knead.Shape
          (Size, Vec2(Vec2), Dim1, Dim2, Shape2, Index2, Ix2,
           verticalVal, horizontalVal)
 import Degree (Degree(Degree), getDegree)
@@ -31,7 +33,6 @@
 import Data.Array.Knead.Expression
          (Exp, (==*), (/=*), (<*), (<=*), (>=*), (&&*))
 
-import qualified Data.Array.CArray as CArray
 import Data.Array.IArray (amap)
 import Data.Array.CArray (CArray)
 import Data.Array.MArray (thaw)
@@ -45,7 +46,7 @@
 import qualified LLVM.Core as LLVM
 
 import qualified Data.Complex as Complex
-import Data.Complex (Complex((:+)), conjugate, realPart)
+import Data.Complex (Complex((:+)))
 
 import qualified Codec.Picture as Pic
 
@@ -61,7 +62,7 @@
 import Text.Printf (printf)
 
 import qualified Control.Monad.HT as MonadHT
-import Control.Monad (liftM2, when, join, foldM, (<=<))
+import Control.Monad (when, join, foldM, (<=<))
 import Control.Applicative (pure, (<$>), (<*>))
 
 import qualified Data.Foldable as Fold
@@ -71,12 +72,12 @@
 import Data.Monoid ((<>))
 import Data.Maybe.HT (toMaybe)
 import Data.Maybe (mapMaybe, isJust, isNothing)
+import Data.Bits (Bits)
 import Data.Traversable (forM)
 import Data.Foldable (forM_)
 import Data.Ord.HT (comparing)
 import Data.Tuple.HT
-         (mapPair, mapFst, mapSnd, mapTriple, swap,
-          mapThd3, fst3, thd3, uncurry3)
+         (mapPair, mapFst, mapSnd, mapTriple, swap, mapThd3, fst3, uncurry3)
 import Data.Word (Word8, Word32)
 
 
@@ -563,125 +564,63 @@
 
 
 -- counterpart to 'clip'
-pad ::
-   (MultiValue.C a) =>
-   Exp a -> Exp Dim2 -> SymbPlane a -> SymbPlane a
+pad :: (MultiValue.C a) => Exp a -> Exp Dim2 -> SymbPlane a -> SymbPlane a
 pad a sh img =
    let Vec2 height width = Expr.decompose atomDim2 $ Symb.shape img
    in  generate sh $ \p ->
          let Vec2 y x = Expr.decompose atomIx2 p
          in  Expr.ifThenElse (y<*height &&* x<*width) (img ! p) a
 
-padCArray ::
-   (SV.Storable a) =>
-   a -> (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a
-padCArray a (height, width) img =
-   CArray.listArray ((0,0), (height-1, width-1)) (repeat a)
-   CArray.//
-   CArray.assocs img
-
-clipCArray ::
-   (SV.Storable a) => (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a
-clipCArray (height, width) =
-   CArray.ixmap ((0,0), (height-1, width-1)) id
-
-mapPairInt :: (Integral i, Integral j) => (i,i) -> (j,j)
-mapPairInt = mapPair (fromIntegral, fromIntegral)
-
-correlatePaddedSimpleCArray ::
-   (FFTWReal a) =>
-   (Int,Int) ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a
-correlatePaddedSimpleCArray sh =
-   let forward = FFT.dftRCN [0,1] . padCArray 0 sh
-       inverse = FFT.dftCRN [0,1]
-   in  \ a b ->
-         inverse $ CArray.liftArray2 Arith.mulConj (forward a) (forward b)
-
--- expects zero-based arrays
-cyclicReverse2d :: (SV.Storable a) => CArray (Int,Int) a -> CArray (Int,Int) a
+cyclicReverse2d :: (MultiValue.C a) => SymbPlane a -> SymbPlane a
 cyclicReverse2d spec =
-   let (height, width) = mapPair ((1+), (1+)) $ snd $ CArray.bounds spec
-   in  CArray.ixmap (CArray.bounds spec)
-         (\(y,x) -> (mod (-y) height, mod (-x) width)) spec
+   let (Vec2 height width) = Expr.decompose atomDim2 $ Symb.shape spec
+   in  Symb.backpermute (Symb.shape spec)
+         (Expr.modify atomIx2 $ \(Vec2 y x) ->
+            Vec2
+               (wrap height height (height-y))
+               (wrap width width (width-x)))
+         spec
 
-untangleCoefficient ::
-   (RealFloat a) => Complex a -> Complex a -> (Complex a, Complex a)
-untangleCoefficient a b =
-   let bc = conjugate b
-   in  ((a + bc) / 2, (a - bc) * (0 :+ (-1/2)))
+atomComplex :: Complex (Atom a)
+atomComplex = atom:+atom
 
--- ToDo: could be moved to fft package
 untangleSpectra2d ::
-   (RealFloat a, SV.Storable a) =>
-   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a, Complex a)
+   (MultiValue.C a, MultiValue.Field a,
+    MultiValue.Real a, MultiValue.RationalConstant a) =>
+   SymbPlane (Complex a) -> SymbPlane (Complex a, Complex a)
 untangleSpectra2d spec =
-   CArray.liftArray2 untangleCoefficient spec (cyclicReverse2d spec)
-
-{- |
-Equivalent to @amap (uncurry Arith.mulConj) . untangleSpectra2d@
-but much faster, since it avoids the slow @instance Storable (a,b)@
-based on @storable-tuple:storePair@.
--}
-mulConjUntangledSpectra2d ::
-   (RealFloat a, SV.Storable a) =>
-   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a)
-mulConjUntangledSpectra2d spec =
-   CArray.liftArray2
-      ((uncurry Arith.mulConj .) . untangleCoefficient)
+   Symb.zipWith
+      (Expr.modify2 atomComplex atomComplex KneadCArray.untangleCoefficient)
       spec (cyclicReverse2d spec)
 
-
-{-
-This is more efficient than 'correlatePaddedSimpleCArray'
-since it needs only one complex forward Fourier transform,
-where 'correlatePaddedSimpleCArray' needs two real transforms.
-Especially for odd sizes
-two real transforms are slower than a complex transform.
-For the analysis part,
-perform two real-valued Fourier transforms using one complex-valued transform.
-Afterwards we untangle the superposed spectra.
--}
-correlatePaddedComplexCArray ::
-   (FFTWReal a) =>
-   (Int,Int) ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a
-correlatePaddedComplexCArray sh a b =
-   amap realPart $ FFT.idftN [0,1] $
-   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $
-   CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b)
+correlatePadded ::
+   (FFTWReal a, MultiValue.Real a, MultiMem.C a,
+    MultiValue.Field a, MultiValue.RationalConstant a) =>
+   Dim2 -> IO (Plane a -> Plane a -> IO (Plane a))
+correlatePadded padExtent@(Vec2 height width) = do
+   let sh = Expr.cons padExtent
+   mergePlanes <-
+      RenderP.run $ \a b ->
+         Symb.zipWith Expr.consComplex (pad 0 sh a) (pad 0 sh b)
+   let (halfWidth,parity) = divMod width 2
+   let exprFromInt = Expr.cons . fromIntegral
+   mulSpecs <-
+      RenderP.run $
+         clip (0,0) (exprFromInt $ halfWidth+1, exprFromInt height) .
+         Symb.map
+            (Expr.modify (atomComplex, atomComplex) $ uncurry Komplex.mulConj) .
+         untangleSpectra2d
 
-{- |
-Should be yet a little bit more efficient than 'correlatePaddedComplexCArray'
-since it uses a real back transform.
--}
-correlatePaddedCArray ::
-   (FFTWReal a) =>
-   (Int,Int) ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a ->
-   CArray (Int,Int) a
-correlatePaddedCArray sh@(height,width) a b =
-   (case divMod width 2 of
-      (halfWidth,0) -> FFT.dftCRN [0,1] . clipCArray (height,halfWidth+1)
-      (halfWidth,_) -> FFT.dftCRON [0,1] . clipCArray (height,halfWidth+1)) $
-   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $
-   CArray.liftArray2 (:+) (padCArray 0 sh a) (padCArray 0 sh b)
+   return $ \ a b ->
+      liftCArray (if parity==0 then FFT.dftCRN [0,1] else FFT.dftCRON [0,1]) =<<
+      mulSpecs =<< liftCArray (FFT.dftN [0,1]) =<< mergePlanes a b
 
 
-liftCArray2 ::
-   (SV.Storable a) =>
-   (CArray (Int,Int) a -> CArray (Int,Int) a -> CArray (Int,Int) a) ->
-   Plane a -> Plane a -> IO (Plane a)
-liftCArray2 f a b =
-   arrayKneadFromC <$>
-   liftM2 f
-      (arrayCFromKnead a)
-      (arrayCFromKnead b)
+liftCArray ::
+   (SV.Storable a, SV.Storable b) =>
+   (CArray (Int,Int) a -> CArray (Int,Int) b) ->
+   Plane a -> IO (Plane b)
+liftCArray f a = arrayKneadFromC . f <$> arrayCFromKnead a
 
 
 fixArray :: Id (Symb.Array sh a)
@@ -776,17 +715,17 @@
 
 allOverlapsRun ::
    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Plane Word8))
-allOverlapsRun padExtent@(Vec2 height width) = do
+allOverlapsRun padExtent = do
    run <-
       RenderP.run $ \minOverlapPortion sha shb img ->
          imageByteFromFloat $
          Symb.map (0.0001*) $
          Symb.map Expr.fst $
          allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img
+   correlate <- correlatePadded padExtent
 
    return $ \overlap a b ->
-      run overlap (Phys.shape a) (Phys.shape b)
-         =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b
+      run overlap (Phys.shape a) (Phys.shape b) =<< correlate a b
 
 
 argmax ::
@@ -802,15 +741,15 @@
 
 optimalOverlap ::
    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))
-optimalOverlap padExtent@(Vec2 height width) = do
+optimalOverlap padExtent = do
    run <-
       RenderP.run $ \minOverlapPortion (sha, shb) img ->
          argmaximum $
          allOverlapsFromCorrelation padExtent minOverlapPortion sha shb img
+   correlate <- correlatePadded padExtent
 
    return $ \overlap a b ->
-      run overlap (Phys.shape a, Phys.shape b)
-         =<< liftCArray2 (correlatePaddedCArray $ mapPairInt (height, width)) a b
+      run overlap (Phys.shape a, Phys.shape b) =<< correlate a b
 
 
 shrink ::
@@ -826,14 +765,42 @@
       (Expr.modify (atomIx2, atomIx2) $
          \(Vec2 yi xi, Vec2 yj xj) -> Vec2 (yi*yk+yj) (xi*xk+xj))
 
-shrinkFactors :: (Integral a) => Dim2 -> Shape2 a -> Shape2 a -> Shape2 a
-shrinkFactors (Vec2 heightPad widthPad)
+{-
+The implementation accepts overlapping of at most minOverlapPortion
+of the two shrunken images.
+However, in practice this optimization is rarely effective.
+In most cases the shrink factors are the same
+independent from whether minOverlap is zero or not.
+-}
+shrinkFactors ::
+   (Integral a) => Dim2 -> Float -> Shape2 a -> Shape2 a -> Shape2 a
+shrinkFactors (Vec2 heightPad widthPad) minOverlapPortion
    (Vec2 heighta widtha) (Vec2 heightb widthb) =
-      Vec2
-         (Arith.divUp (heighta+heightb) $ fromIntegral heightPad)
-         (Arith.divUp (widtha +widthb)  $ fromIntegral widthPad)
+      let minOverlap =
+            Arith.minimumOverlapAbsFromPortion minOverlapPortion
+               (min widtha widthb, min heighta heightb)
+      in Vec2
+            (Arith.divUp (heighta+heightb-minOverlap) $ fromIntegral heightPad)
+            (Arith.divUp (widtha +widthb -minOverlap) $ fromIntegral widthPad)
 
+{-
+Should compute almost the same as shrinkFactors
+but is less optimized and more idiomatic.
+@correlationSize@ has a final @ceilingSmooth7@.
+This is not necessary here
+since we expect that the user chooses an FFT friendly target size.
+-}
+shrinkFactorsAlt ::
+   (Bits a, Integral a) => Float -> Dim2 -> Shape2 a -> Shape2 a -> Shape2 a
+shrinkFactorsAlt minOverlapPortion (Vec2 heightPad widthPad) a b =
+   let (widthc,heightc) =
+         Arith.correlationSize minOverlapPortion $
+         map (\(Vec2 height width) -> (width, height)) [a,b]
+   in Vec2
+         (Arith.divUp heightc $ fromIntegral heightPad)
+         (Arith.divUp widthc $ fromIntegral widthPad)
 
+
 optimalOverlapBig ::
    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))
 optimalOverlapBig padExtent = do
@@ -841,7 +808,7 @@
    optOverlap <- optimalOverlap padExtent
    return $ \minimumOverlap a b -> do
       let factors@(Vec2 yk xk) =
-            shrinkFactors padExtent (Phys.shape a) (Phys.shape b)
+            shrinkFactors padExtent minimumOverlap (Phys.shape a) (Phys.shape b)
       aSmall <- shrnk factors a
       bSmall <- shrnk factors b
       mapSnd (mapPair ((*xk), (*yk))) <$>
@@ -883,7 +850,7 @@
    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))
 optimalOverlapBigFine padExtent@(Vec2 heightPad widthPad) = do
    overlap <- optimalOverlap padExtent
-   -- optimalOverlap is rendered again here
+   -- optimalOverlap is compiled again here
    overlapBig <- optimalOverlapBig padExtent
    clp <- RenderP.run clip
    return $ \minimumOverlap a b -> do
@@ -936,7 +903,7 @@
 
    return $ \minimumOverlap mMaximumDiff a b -> do
       let factors@(Vec2 yk xk) =
-            shrinkFactors padExtent (Phys.shape a) (Phys.shape b)
+            shrinkFactors padExtent minimumOverlap (Phys.shape a) (Phys.shape b)
       aSmall <- shrnk factors a
       bSmall <- shrnk factors b
 
@@ -1306,9 +1273,7 @@
     MultiValue.NativeFloating a ar) =>
    Exp (Geometry b) -> SymbPlane a -> SymbPlane Word8
 scaleDistanceMapGeom geom img =
-   let scale =
-         (4/) $ fromInt $
-         Expr.modify (atom,atom,(atom,atom)) (uncurry Expr.min . thd3) geom
+   let scale = (4/) $ fromInt $ Expr.uncurry Expr.min $ Expr.thd3 geom
    in  imageByteFromFloat $ Symb.map (scale*) img
 
 
@@ -1689,7 +1654,7 @@
        map picColored pics,
        map
          (mapPair
-            (mapPair (realToFrac, realToFrac), Arith.mapComplex realToFrac))
+            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))
          posRots)
 
 
@@ -1728,20 +1693,22 @@
    when False $ do
       notice "write fft"
       let pic0 : pic1 : _ = map snd rotated
-          size = (1024,768)
-      cpic0 <- arrayCFromKnead pic0
-      cpic1 <- arrayCFromKnead pic1
+          size = Vec2 1024 768
       makeByteImage <-
          RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . fixArray
+      runPad <- RenderP.run pad
       writeGrey (Option.quality opt) "/tmp/padded.jpeg" =<<
-         (makeByteImage 1 $ arrayKneadFromC $ padCArray 0 size cpic0)
+         (makeByteImage 1 =<< runPad 0 size pic0)
+      runMagnitude <-
+         RenderP.run $
+         Symb.map (Expr.modify atomComplex $ \(r:+i) -> Expr.sqrt$ r*r+i*i)
+            . fixArray
       writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" =<<
-         (makeByteImage 0.1 $ arrayKneadFromC $
-          CArray.liftArray Complex.magnitude $
-          FFT.dftRCN [0,1] $ padCArray 0 size cpic0)
+         (makeByteImage 0.1 =<< runMagnitude =<<
+          liftCArray (FFT.dftRCN [0,1]) =<< runPad 0 size pic0)
+      correlate <- correlatePadded size
       writeGrey (Option.quality opt) "/tmp/convolution.jpeg" =<<
-         (makeByteImage 0.1 $ arrayKneadFromC $
-          correlatePaddedCArray size cpic0 cpic1)
+         (makeByteImage 0.1 =<< correlate pic0 pic1)
 
    return $
       zipWith3
diff --git a/src/Knead/CArray.hs b/src/Knead/CArray.hs
new file mode 100644
--- /dev/null
+++ b/src/Knead/CArray.hs
@@ -0,0 +1,119 @@
+module Knead.CArray where
+
+import qualified Complex as Komplex
+
+import qualified Math.FFT as FFT
+import Math.FFT.Base (FFTWReal)
+
+import Foreign.Storable (Storable)
+import Foreign.Storable.Tuple ()
+
+import qualified Data.Array.CArray as CArray
+import Data.Array.IArray (amap)
+import Data.Array.CArray (CArray)
+
+import Data.Complex (Complex((:+)), realPart)
+
+import Data.Tuple.HT (mapPair)
+
+
+pad ::
+   (Storable a) => a -> (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a
+pad a (height, width) img =
+   CArray.listArray ((0,0), (height-1, width-1)) (repeat a)
+   CArray.//
+   CArray.assocs img
+
+clip :: (Storable a) => (Int,Int) -> CArray (Int,Int) a -> CArray (Int,Int) a
+clip (height, width) = CArray.ixmap ((0,0), (height-1, width-1)) id
+
+correlatePaddedSimple ::
+   (FFTWReal a) =>
+   (Int,Int) ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a
+correlatePaddedSimple sh =
+   let forward = FFT.dftRCN [0,1] . pad 0 sh
+       inverse = FFT.dftCRN [0,1]
+   in  \ a b ->
+         inverse $ CArray.liftArray2 Komplex.mulConj (forward a) (forward b)
+
+-- expects zero-based arrays
+cyclicReverse2d :: (Storable a) => CArray (Int,Int) a -> CArray (Int,Int) a
+cyclicReverse2d spec =
+   let (height, width) = mapPair ((1+), (1+)) $ snd $ CArray.bounds spec
+   in  CArray.ixmap (CArray.bounds spec)
+         (\(y,x) -> (mod (-y) height, mod (-x) width)) spec
+
+untangleCoefficient ::
+   (Fractional a) => Complex a -> Complex a -> (Complex a, Complex a)
+untangleCoefficient a b =
+   let bc = Komplex.conjugate b
+   in  (Komplex.mul (Komplex.add a bc) ((1/2) :+ 0),
+        Komplex.mul (Komplex.sub a bc) (0 :+ (-1/2)))
+
+untangleCoefficient_ ::
+   (RealFloat a) => Complex a -> Complex a -> (Complex a, Complex a)
+untangleCoefficient_ a b =
+   let bc = Komplex.conjugate b
+   in  ((a + bc) / 2, (a - bc) * (0 :+ (-1/2)))
+
+-- ToDo: could be moved to fft package
+untangleSpectra2d ::
+   (Fractional a, Storable a) =>
+   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a, Complex a)
+untangleSpectra2d spec =
+   CArray.liftArray2 untangleCoefficient spec (cyclicReverse2d spec)
+
+{- |
+Equivalent to @amap (uncurry Komplex.mulConj) . untangleSpectra2d@
+but much faster, since it avoids the slow @instance Storable (a,b)@
+based on @storable-tuple:storePair@.
+-}
+mulConjUntangledSpectra2d ::
+   (Fractional a, Storable a) =>
+   CArray (Int,Int) (Complex a) -> CArray (Int,Int) (Complex a)
+mulConjUntangledSpectra2d spec =
+   CArray.liftArray2
+      ((uncurry Komplex.mulConj .) . untangleCoefficient)
+      spec (cyclicReverse2d spec)
+
+
+{-
+This is more efficient than 'correlatePaddedSimpleCArray'
+since it needs only one complex forward Fourier transform,
+where 'correlatePaddedSimpleCArray' needs two real transforms.
+Especially for odd sizes
+two real transforms are slower than a complex transform.
+For the analysis part,
+perform two real-valued Fourier transforms using one complex-valued transform.
+Afterwards we untangle the superposed spectra.
+-}
+correlatePaddedComplex ::
+   (FFTWReal a) =>
+   (Int,Int) ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a
+correlatePaddedComplex sh a b =
+   amap realPart $ FFT.idftN [0,1] $
+   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $
+   CArray.liftArray2 (:+) (pad 0 sh a) (pad 0 sh b)
+
+{- |
+Should be yet a little bit more efficient than 'correlatePaddedComplexCArray'
+since it uses a real back transform.
+-}
+correlatePadded ::
+   (FFTWReal a) =>
+   (Int,Int) ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a ->
+   CArray (Int,Int) a
+correlatePadded sh@(height,width) a b =
+   (case divMod width 2 of
+      (halfWidth,0) -> FFT.dftCRN [0,1] . clip (height,halfWidth+1)
+      (halfWidth,_) -> FFT.dftCRON [0,1] . clip (height,halfWidth+1)) $
+   mulConjUntangledSpectra2d $ FFT.dftN [0,1] $
+   CArray.liftArray2 (:+) (pad 0 sh a) (pad 0 sh b)
diff --git a/src/Knead/Shape.hs b/src/Knead/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Knead/Shape.hs
@@ -0,0 +1,176 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE EmptyDataDecls #-}
+module Knead.Shape where
+
+import qualified Data.Array.Knead.Index.Nested.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+
+import qualified LLVM.Extra.Multi.Value.Memory as MultiMem
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Arithmetic as A
+import LLVM.Extra.Multi.Value (atom)
+
+import qualified LLVM.Core as LLVM
+
+import qualified Type.Data.Num.Decimal as TypeNum
+
+import Foreign.Storable
+         (Storable, sizeOf, alignment, poke, pokeElemOff, peek, peekElemOff)
+import Foreign.Ptr (Ptr, castPtr)
+
+import qualified Control.Monad.HT as Monad
+import Control.Monad (join)
+
+import Data.Int (Int64)
+
+
+{- |
+I choose a bit complicated Dim2 definition
+to make it distinct from size pairs with width and height swapped.
+Alternatives would be Index.Linear or intentionally complicated Shape types like:
+
+type Dim0 = ()
+type Dim1 = ((), Size)
+type Dim2 = ((), Size, Size)
+
+Problems with Index.Linear is that it is fixed to Word32 dimensions
+which causes trouble with negative coordinates
+that we encounter on rotations.
+
+The custom shape type requires lots of new definitions
+but it is certainly the cleanest solution.
+-}
+type Size = Int64
+type Dim0 = ()
+type Dim1 = Size
+type Dim2 = Shape2 Size
+type Ix2  = Index2 Size
+
+data Vec2 tag i = Vec2 {vertical, horizontal :: i}
+
+data ShapeTag
+data IndexTag
+
+type Shape2 = Vec2 ShapeTag
+type Index2 = Vec2 IndexTag
+
+
+
+squareShape :: n -> Vec2 tag n
+squareShape n = Vec2 n n
+
+castToElemPtr :: Ptr (Vec2 tag a) -> Ptr a
+castToElemPtr = castPtr
+
+instance (Storable n) => Storable (Vec2 tag n) where
+   -- cf. sample-frame:Frame.Stereo
+   sizeOf ~(Vec2 n m) =
+      sizeOf n + mod (- sizeOf n) (alignment m) + sizeOf m
+   alignment ~(Vec2 n _) = alignment n
+   poke p (Vec2 n m) =
+      let q = castToElemPtr p
+      in  poke q n >> pokeElemOff q 1 m
+   peek p =
+      let q = castToElemPtr p
+      in  Monad.lift2 Vec2 (peek q) (peekElemOff q 1)
+
+instance (MultiValue.C n) => MultiValue.C (Vec2 tag n) where
+   type Repr f (Vec2 tag n) = Vec2 tag (MultiValue.Repr f n)
+   cons (Vec2 n m) =
+      MultiValue.compose $ Vec2 (MultiValue.cons n) (MultiValue.cons m)
+   undef = MultiValue.compose $ squareShape MultiValue.undef
+   zero = MultiValue.compose $ squareShape MultiValue.zero
+   phis bb a =
+      case MultiValue.decompose (squareShape atom) a of
+         Vec2 a0 a1 ->
+            fmap MultiValue.compose $
+            Monad.lift2 Vec2 (MultiValue.phis bb a0) (MultiValue.phis bb a1)
+   addPhis bb a b =
+      case (MultiValue.decompose (squareShape atom) a,
+            MultiValue.decompose (squareShape atom) b) of
+         (Vec2 a0 a1, Vec2 b0 b1) ->
+            MultiValue.addPhis bb a0 b0 >>
+            MultiValue.addPhis bb a1 b1
+
+type instance
+   MultiValue.Decomposed f (Vec2 tag pat) =
+      Vec2 tag (MultiValue.Decomposed f pat)
+type instance
+   MultiValue.PatternTuple (Vec2 tag pat) =
+      Vec2 tag (MultiValue.PatternTuple pat)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Vec2 tag n) where
+   type Composed (Vec2 tag n) = Vec2 tag (MultiValue.Composed n)
+   compose (Vec2 n m) =
+      case (MultiValue.compose n, MultiValue.compose m) of
+         (MultiValue.Cons rn, MultiValue.Cons rm) ->
+            MultiValue.Cons (Vec2 rn rm)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Vec2 tag pn) where
+   decompose (Vec2 pn pm) (MultiValue.Cons (Vec2 n m)) =
+      Vec2
+         (MultiValue.decompose pn (MultiValue.Cons n))
+         (MultiValue.decompose pm (MultiValue.Cons m))
+
+instance (MultiMem.C i) => MultiMem.C (Vec2 tag i) where
+   type Struct (Vec2 tag i) =
+         LLVM.Struct (MultiMem.Struct i, (MultiMem.Struct i, ()))
+   decompose nm =
+      Monad.lift2 zipShape
+         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)
+         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)
+   compose nm =
+      case unzipShape nm of
+         Vec2 n m -> do
+            sn <- MultiMem.compose n
+            sm <- MultiMem.compose m
+            rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0
+            LLVM.insertvalue rn sm TypeNum.d1
+
+
+unzipShape :: MultiValue.T (Vec2 tag n) -> Vec2 tag (MultiValue.T n)
+unzipShape = MultiValue.decompose (squareShape atom)
+
+zipShape :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Vec2 tag n)
+zipShape y x = MultiValue.compose $ Vec2 y x
+
+instance (tag ~ ShapeTag, Shape.C i) => Shape.C (Vec2 tag i) where
+   type Index (Vec2 tag i) = Index2 (Shape.Index i)
+   intersectCode a b =
+      case (unzipShape a, unzipShape b) of
+         (Vec2 an am, Vec2 bn bm) ->
+            Monad.lift2 zipShape
+               (Shape.intersectCode an bn)
+               (Shape.intersectCode am bm)
+   sizeCode nm =
+      case unzipShape nm of
+         Vec2 n m ->
+            join $ Monad.lift2 A.mul (Shape.sizeCode n) (Shape.sizeCode m)
+   size (Vec2 n m) = Shape.size n * Shape.size m
+   flattenIndexRec nm ij =
+      case (unzipShape nm, unzipShape ij) of
+         (Vec2 n m, Vec2 i j) -> do
+            (ns, il) <- Shape.flattenIndexRec n i
+            (ms, jl) <- Shape.flattenIndexRec m j
+            Monad.lift2 (,)
+               (A.mul ns ms)
+               (A.add jl =<< A.mul ms il)
+   loop code nm =
+      case unzipShape nm of
+         Vec2 n m ->
+            Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n
+
+
+instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where
+   type Composed (Vec2 tag n) = Vec2 tag (Expr.Composed n)
+   compose (Vec2 n m) = Expr.lift2 zipShape (Expr.compose n) (Expr.compose m)
+
+instance (Expr.Decompose p) => Expr.Decompose (Vec2 tag p) where
+   decompose (Vec2 pn pm) vec =
+      Vec2
+         (Expr.decompose pn (verticalVal vec))
+         (Expr.decompose pm (horizontalVal vec))
+
+verticalVal, horizontalVal :: (Expr.Value val) => val (Vec2 tag n) -> val n
+verticalVal = Expr.lift1 (MultiValue.lift1 vertical)
+horizontalVal = Expr.lift1 (MultiValue.lift1 horizontal)
diff --git a/src/KneadShape.hs b/src/KneadShape.hs
deleted file mode 100644
--- a/src/KneadShape.hs
+++ /dev/null
@@ -1,176 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE EmptyDataDecls #-}
-module KneadShape where
-
-import qualified Data.Array.Knead.Index.Nested.Shape as Shape
-import qualified Data.Array.Knead.Expression as Expr
-
-import qualified LLVM.Extra.Multi.Value.Memory as MultiMem
-import qualified LLVM.Extra.Multi.Value as MultiValue
-import qualified LLVM.Extra.Arithmetic as A
-import LLVM.Extra.Multi.Value (atom)
-
-import qualified LLVM.Core as LLVM
-
-import qualified Type.Data.Num.Decimal as TypeNum
-
-import Foreign.Storable
-         (Storable, sizeOf, alignment, poke, pokeElemOff, peek, peekElemOff)
-import Foreign.Ptr (Ptr, castPtr)
-
-import qualified Control.Monad.HT as Monad
-import Control.Monad (join)
-
-import Data.Int (Int64)
-
-
-{- |
-I choose a bit complicated Dim2 definition
-to make it distinct from size pairs with width and height swapped.
-Alternatives would be Index.Linear or intentionally complicated Shape types like:
-
-type Dim0 = ()
-type Dim1 = ((), Size)
-type Dim2 = ((), Size, Size)
-
-Problems with Index.Linear is that it is fixed to Word32 dimensions
-which causes trouble with negative coordinates
-that we encounter on rotations.
-
-The custom shape type requires lots of new definitions
-but it is certainly the cleanest solution.
--}
-type Size = Int64
-type Dim0 = ()
-type Dim1 = Size
-type Dim2 = Shape2 Size
-type Ix2  = Index2 Size
-
-data Vec2 tag i = Vec2 {vertical, horizontal :: i}
-
-data ShapeTag
-data IndexTag
-
-type Shape2 = Vec2 ShapeTag
-type Index2 = Vec2 IndexTag
-
-
-
-squareShape :: n -> Vec2 tag n
-squareShape n = Vec2 n n
-
-castToElemPtr :: Ptr (Vec2 tag a) -> Ptr a
-castToElemPtr = castPtr
-
-instance (Storable n) => Storable (Vec2 tag n) where
-   -- cf. sample-frame:Frame.Stereo
-   sizeOf ~(Vec2 n m) =
-      sizeOf n + mod (- sizeOf n) (alignment m) + sizeOf m
-   alignment ~(Vec2 n _) = alignment n
-   poke p (Vec2 n m) =
-      let q = castToElemPtr p
-      in  poke q n >> pokeElemOff q 1 m
-   peek p =
-      let q = castToElemPtr p
-      in  Monad.lift2 Vec2 (peek q) (peekElemOff q 1)
-
-instance (MultiValue.C n) => MultiValue.C (Vec2 tag n) where
-   type Repr f (Vec2 tag n) = Vec2 tag (MultiValue.Repr f n)
-   cons (Vec2 n m) =
-      MultiValue.compose $ Vec2 (MultiValue.cons n) (MultiValue.cons m)
-   undef = MultiValue.compose $ squareShape MultiValue.undef
-   zero = MultiValue.compose $ squareShape MultiValue.zero
-   phis bb a =
-      case MultiValue.decompose (squareShape atom) a of
-         Vec2 a0 a1 ->
-            fmap MultiValue.compose $
-            Monad.lift2 Vec2 (MultiValue.phis bb a0) (MultiValue.phis bb a1)
-   addPhis bb a b =
-      case (MultiValue.decompose (squareShape atom) a,
-            MultiValue.decompose (squareShape atom) b) of
-         (Vec2 a0 a1, Vec2 b0 b1) ->
-            MultiValue.addPhis bb a0 b0 >>
-            MultiValue.addPhis bb a1 b1
-
-type instance
-   MultiValue.Decomposed f (Vec2 tag pat) =
-      Vec2 tag (MultiValue.Decomposed f pat)
-type instance
-   MultiValue.PatternTuple (Vec2 tag pat) =
-      Vec2 tag (MultiValue.PatternTuple pat)
-
-instance (MultiValue.Compose n) => MultiValue.Compose (Vec2 tag n) where
-   type Composed (Vec2 tag n) = Vec2 tag (MultiValue.Composed n)
-   compose (Vec2 n m) =
-      case (MultiValue.compose n, MultiValue.compose m) of
-         (MultiValue.Cons rn, MultiValue.Cons rm) ->
-            MultiValue.Cons (Vec2 rn rm)
-
-instance (MultiValue.Decompose pn) => MultiValue.Decompose (Vec2 tag pn) where
-   decompose (Vec2 pn pm) (MultiValue.Cons (Vec2 n m)) =
-      Vec2
-         (MultiValue.decompose pn (MultiValue.Cons n))
-         (MultiValue.decompose pm (MultiValue.Cons m))
-
-instance (MultiMem.C i) => MultiMem.C (Vec2 tag i) where
-   type Struct (Vec2 tag i) =
-         LLVM.Struct (MultiMem.Struct i, (MultiMem.Struct i, ()))
-   decompose nm =
-      Monad.lift2 zipShape
-         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d0)
-         (MultiMem.decompose =<< LLVM.extractvalue nm TypeNum.d1)
-   compose nm =
-      case unzipShape nm of
-         Vec2 n m -> do
-            sn <- MultiMem.compose n
-            sm <- MultiMem.compose m
-            rn <- LLVM.insertvalue (LLVM.value LLVM.undef) sn TypeNum.d0
-            LLVM.insertvalue rn sm TypeNum.d1
-
-
-unzipShape :: MultiValue.T (Vec2 tag n) -> Vec2 tag (MultiValue.T n)
-unzipShape = MultiValue.decompose (squareShape atom)
-
-zipShape :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Vec2 tag n)
-zipShape y x = MultiValue.compose $ Vec2 y x
-
-instance (tag ~ ShapeTag, Shape.C i) => Shape.C (Vec2 tag i) where
-   type Index (Vec2 tag i) = Index2 (Shape.Index i)
-   intersectCode a b =
-      case (unzipShape a, unzipShape b) of
-         (Vec2 an am, Vec2 bn bm) ->
-            Monad.lift2 zipShape
-               (Shape.intersectCode an bn)
-               (Shape.intersectCode am bm)
-   sizeCode nm =
-      case unzipShape nm of
-         Vec2 n m ->
-            join $ Monad.lift2 A.mul (Shape.sizeCode n) (Shape.sizeCode m)
-   size (Vec2 n m) = Shape.size n * Shape.size m
-   flattenIndexRec nm ij =
-      case (unzipShape nm, unzipShape ij) of
-         (Vec2 n m, Vec2 i j) -> do
-            (ns, il) <- Shape.flattenIndexRec n i
-            (ms, jl) <- Shape.flattenIndexRec m j
-            Monad.lift2 (,)
-               (A.mul ns ms)
-               (A.add jl =<< A.mul ms il)
-   loop code nm =
-      case unzipShape nm of
-         Vec2 n m ->
-            Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n
-
-
-instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where
-   type Composed (Vec2 tag n) = Vec2 tag (Expr.Composed n)
-   compose (Vec2 n m) = Expr.lift2 zipShape (Expr.compose n) (Expr.compose m)
-
-instance (Expr.Decompose p) => Expr.Decompose (Vec2 tag p) where
-   decompose (Vec2 pn pm) vec =
-      Vec2
-         (Expr.decompose pn (verticalVal vec))
-         (Expr.decompose pm (horizontalVal vec))
-
-verticalVal, horizontalVal :: (Expr.Value val) => val (Vec2 tag n) -> val n
-verticalVal = Expr.lift1 (MultiValue.lift1 vertical)
-horizontalVal = Expr.lift1 (MultiValue.lift1 horizontal)
diff --git a/src/MatchImageBorders.hs b/src/MatchImageBorders.hs
--- a/src/MatchImageBorders.hs
+++ b/src/MatchImageBorders.hs
@@ -10,7 +10,7 @@
 -}
 module MatchImageBorders where
 
-import KneadShape (Vec2(Vec2), Dim2)
+import Knead.Shape (Vec2(Vec2), Dim2)
 
 import qualified Data.Array.Knead.Simple.Physical as Phys