packages feed

accelerate-fourier (empty) → 0.0

raw patch · 20 files changed

+3483/−0 lines, 20 filesdep +QuickCheckdep +acceleratedep +accelerate-arithmeticsetup-changed

Dependencies added: QuickCheck, accelerate, accelerate-arithmetic, accelerate-cuda, accelerate-fourier, accelerate-utility, base, containers, criterion, transformers, utility-ht

Files

+ LICENSE view
@@ -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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ accelerate-fourier.cabal view
@@ -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
+ benchmark/Main.hs view
@@ -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 :+   []
+ src/Data/Array/Accelerate/Convolution/Adhoc.hs view
@@ -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))
+ src/Data/Array/Accelerate/Convolution/Preprocessed.hs view
@@ -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
+ src/Data/Array/Accelerate/Convolution/Private.hs view
@@ -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)
+ src/Data/Array/Accelerate/Convolution/Small.hs view
@@ -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))
+ src/Data/Array/Accelerate/Cyclic.hs view
@@ -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
+ src/Data/Array/Accelerate/Fourier/Adhoc.hs view
@@ -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
+ src/Data/Array/Accelerate/Fourier/Planned.hs view
@@ -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
+ src/Data/Array/Accelerate/Fourier/Preprocessed.hs view
@@ -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
+ src/Data/Array/Accelerate/Fourier/Private.hs view
@@ -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)
+ src/Data/Array/Accelerate/Fourier/Real.hs view
@@ -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
+ src/Data/Array/Accelerate/Fourier/Sign.hs view
@@ -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]
+ src/Data/Array/Accelerate/Fourier/Utility.hs view
@@ -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)
+ src/Data/Array/Accelerate/NumberTheory.hs view
@@ -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
+ src/Data/Array/Accelerate/Permutation.hs view
@@ -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
+ test/Test.hs view
@@ -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
+ test/Test/Data/Array/Accelerate/Fourier.hs view
@@ -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)) :+   []