{- |
Complete implementation for Fast Fourier Transform for any signal length.
Although defined for all kinds of signal storage,
we need fast access to arbitrary indices.
-}
{-
further thoughts
 - Test the algorithms using remainder polynomials
     with respect to [-1,0,...,0,1]
     Problem: We would need large polynomial degrees,
     namely LCM of the size of all sub-transforms.
     Those numbers are in the same magnitude
     as the integers we use for our integer residue class arithmetic.
 - a z-transform by convolving with a chirp would be nice,
     however we need a square of the primitive root of unity
     in order to compute cis((i/n)^2/2)
 - Can we write the Fourier transforms for lengths larger than the input signal length
     with implicit zero padding?
     This would be useful for Fourier based convolution.
     Our frequent use of 'rechunk' would be a problem, though.
     transformCoprime also needs explicit zero padding.
 - a type class could unify all Level generators
     and thus they would allow for a generic way to call a certain sub-transform
-}
{-# LANGUAGE NoImplicitPrelude #-}
module Synthesizer.Generic.Fourier (
   Element(..),
   -- * conversion between time and frequency domain (spectrum)
   transformForward,
   transformBackward,
   cacheForward,
   cacheBackward,
   cacheDuplex,
   transformWithCache,
   -- * convolution based on Fourier transform
   convolveCyclic,
   Window,
   window,
   convolveWithWindow,
   ) where

import qualified Synthesizer.Generic.Signal as SigG
import qualified Synthesizer.Generic.Cut as CutG
import qualified Synthesizer.Generic.Cyclic as Cyclic
import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG

import qualified Synthesizer.Generic.Permutation as Permutation
import qualified Synthesizer.Basic.NumberTheory as NumberTheory

import qualified Synthesizer.State.Analysis as Ana
import qualified Synthesizer.State.Signal as SigS

import qualified Control.Monad.Trans.State as State
import Control.Monad (liftM2, )
import Control.Applicative ((<$>), )

import qualified Data.Map as Map
import Data.Tuple.HT (mapPair, )

import qualified Algebra.Transcendental as Trans
import qualified Algebra.Ring as Ring
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.IntegralDomain as Integral

import qualified Number.ResidueClass.Check as RC
import Number.ResidueClass.Check ((/:), )

import qualified Number.Complex as Complex
import Number.Complex ((+:))

import NumericPrelude.Numeric
import NumericPrelude.Base hiding (head, )



class Ring.C y => Element y where
   recipInteger :: (SigG.Read sig y) => sig y -> y
   addId :: (SigG.Read sig y) => sig y -> y
   multId :: (SigG.Read sig y) => sig y -> y
   {- |
   It must hold:

   > uncurry (*) (conjugatePrimitiveRootsOfUnity n) = 1

   > mapPair ((^m), (^m)) (conjugatePrimitiveRootsOfUnity (n*m) y)
   >    == conjugatePrimitiveRootsOfUnity n y@

   since we need for caching that the cache is uniquely determined
   by singal length and transform direction.
   -}
   conjugatePrimitiveRootsOfUnity :: (SigG.Read sig y) => sig y -> (y,y)

instance Trans.C a => Element (Complex.T a) where
   recipInteger :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
recipInteger sig (T a)
sig = forall a. C a => a -> a
recip (forall a b. (C a, C b) => a -> b
fromIntegral (forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig)) forall a. a -> a -> T a
+: forall a. C a => a
zero
   addId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
addId sig (T a)
_sig = forall a. C a => a
zero
   multId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
multId sig (T a)
_sig = forall a. C a => a
one
   conjugatePrimitiveRootsOfUnity :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> (T a, T a)
conjugatePrimitiveRootsOfUnity sig (T a)
sig =
      (\T a
x -> (T a
x, forall a. C a => T a -> T a
Complex.conjugate T a
x)) forall a b. (a -> b) -> a -> b
$
      case forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig of
         Int
1 -> forall a. C a => a
one
         Int
2 -> forall a. C a => a -> a
negate forall a. C a => a
one
         Int
3 -> (forall a. C a => a -> a
negate forall a. C a => a
one forall a. a -> a -> T a
+: forall a. C a => a -> a
sqrt a
3) forall a. C a => a -> a -> a
/ T a
2
         Int
4 -> forall a. C a => a
zero forall a. a -> a -> T a
+: forall a. C a => a
one
         Int
5 ->
            let sqrt5 :: a
sqrt5 = forall a. C a => a -> a
sqrt a
5
            in  ((a
sqrt5 forall a. C a => a -> a -> a
- a
1) forall a. a -> a -> T a
+: forall a. C a => a -> a
sqrt a
2 forall a. C a => a -> a -> a
* forall a. C a => a -> a
sqrt(a
5 forall a. C a => a -> a -> a
+ a
sqrt5)) forall a. C a => a -> a -> a
/ T a
4
         Int
6 -> (forall a. C a => a
one forall a. a -> a -> T a
+: forall a. C a => a -> a
sqrt a
3) forall a. C a => a -> a -> a
/ T a
2
         Int
8 -> forall a. C a => a -> T a -> T a
Complex.scale (forall a. C a => a -> a
sqrt a
2 forall a. C a => a -> a -> a
/ a
2) (forall a. C a => a
one forall a. a -> a -> T a
+: forall a. C a => a
one)
         Int
12 -> (forall a. C a => a -> a
sqrt a
3 forall a. a -> a -> T a
+: forall a. C a => a
one) forall a. C a => a -> a -> a
/ T a
2
         Int
n -> forall a. C a => a -> T a
Complex.cis (a
2forall a. C a => a -> a -> a
*forall a. C a => a
pi forall a. C a => a -> a -> a
/ forall a b. (C a, C b) => a -> b
fromIntegral Int
n)

instance (NumberTheory.PrimitiveRoot a, PID.C a, Eq a) => Element (RC.T a) where
   recipInteger :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
recipInteger sig (T a)
sig =
      forall a. C a => a -> a
recip (forall a b. (C a, C b) => a -> b
fromIntegral (forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig) forall a. C a => a -> a -> T a
/: forall a. T a -> a
RC.modulus (forall (sig :: * -> *) y. Read sig y => sig y -> y
head sig (T a)
sig))
   addId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
addId sig (T a)
sig = forall a. C a => a
zero forall a. C a => a -> a -> T a
/: forall a. T a -> a
RC.modulus (forall (sig :: * -> *) y. Read sig y => sig y -> y
head sig (T a)
sig)
   multId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
multId sig (T a)
sig = forall a. C a => a
one forall a. C a => a -> a -> T a
/: forall a. T a -> a
RC.modulus (forall (sig :: * -> *) y. Read sig y => sig y -> y
head sig (T a)
sig)
   {-
   We cannot simply compute
     NumberTheory.primitiveRootsOfUnity modu (SigG.length sig)
   since we have to fulfill the laws.
   In order to fulfill them,
   we choose a root with maximum order,
   this will always be the same,
   and it is a root of all primitive roots
   of any possible order in that ring.
   -}
   conjugatePrimitiveRootsOfUnity :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> (T a, T a)
conjugatePrimitiveRootsOfUnity sig (T a)
sig =
      let modu :: a
modu = forall a. T a -> a
RC.modulus (forall (sig :: * -> *) y. Read sig y => sig y -> y
head sig (T a)
sig)
          order :: Order
order@(NumberTheory.Order Integer
expo) =
             forall a. PrimitiveRoot a => a -> Order
NumberTheory.maximumOrderOfPrimitiveRootsOfUnity a
modu
          a
r:[a]
_ = forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
NumberTheory.primitiveRootsOfUnity a
modu Order
order
          n :: Integer
n = forall a. (C a, C a) => a -> a -> a
Integral.divChecked Integer
expo (forall a b. (C a, C b) => a -> b
fromIntegral (forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig))
          z :: T a
z = (a
r forall a. C a => a -> a -> T a
/: a
modu) forall a. C a => a -> Integer -> a
^ Integer
n
      in  (T a
z, forall a. C a => a -> a
recip T a
z)


head :: (SigG.Read sig y) => sig y -> y
head :: forall (sig :: * -> *) y. Read sig y => sig y -> y
head =
   forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL (forall a. HasCallStack => [Char] -> a
error [Char]
"Generic.Signal.head: empty signal") forall a b. a -> b -> a
const forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState


directionPrimitiveRootsOfUnity ::
   (Element y, SigG.Read sig y) =>
   sig y -> ((Direction,y), (Direction,y))
directionPrimitiveRootsOfUnity :: forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
x =
   let (y
z,y
zInv) = forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> (y, y)
conjugatePrimitiveRootsOfUnity sig y
x
   in  ((Direction
Forward,y
z), (Direction
Backward,y
zInv))

transformForward ::
   (Element y, SigG.Transform sig y) =>
   sig y -> sig y
transformForward :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> sig y
transformForward sig y
xs =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> Cache sig y
cacheForward sig y
xs) sig y
xs

{- |
Shall we divide the result values by the length of the signal?
Our dimensional wrapper around the Fourier transform does not expect this.
-}
transformBackward ::
   (Element y, SigG.Transform sig y) =>
   sig y -> sig y
transformBackward :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> sig y
transformBackward sig y
xs =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> Cache sig y
cacheBackward sig y
xs) sig y
xs

{- |
The size of the signal must match the size, that the plan was generated for.
-}
_transformPlan ::
   (Element y, SigG.Transform sig y) =>
   Plan -> (Direction,y) -> sig y -> sig y
_transformPlan :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> sig y
_transformPlan Plan
p (Direction, y)
z sig y
xs =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan Plan
p (Direction, y)
z sig y
xs) sig y
xs

{- |
The size and type of the signal must match the parameters,
that the cache was generated for.
-}
transformWithCache ::
   (Element y, SigG.Transform sig y) =>
   Cache sig y -> sig y -> sig y
transformWithCache :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache sig y
xs =
   case Cache sig y
cache of
      Cache sig y
CacheIdentity -> sig y
xs
      CacheSmall LevelCacheSmall y
size ->
         case LevelCacheSmall y
size of
            LevelCache2 y
zs -> forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
transform2 y
zs sig y
xs
            LevelCache3 (y, y)
zs -> forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y) -> sig y -> sig y
transform3 (y, y)
zs sig y
xs
            LevelCache4 (y, y, y)
zs -> forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y, y) -> sig y -> sig y
transform4 (y, y, y)
zs sig y
xs
            LevelCache5 (y, y, y, y)
zs -> forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y, y, y) -> sig y -> sig y
transform5 (y, y, y, y)
zs sig y
xs
      CacheNaive LevelCacheNaive y
level ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheNaive y -> sig y -> sig y
transformNaive LevelCacheNaive y
level sig y
xs
      CacheRadix2 LevelCacheRadix2 sig y
level Cache sig y
subCache ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheRadix2 sig y -> Cache sig y -> sig y -> sig y
transformRadix2InterleavedFrequency LevelCacheRadix2 sig y
level Cache sig y
subCache sig y
xs
      CachePrime LevelCachePrime sig y
level (Cache sig y, Cache sig y)
subCaches ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCachePrime sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
transformPrime LevelCachePrime sig y
level (Cache sig y, Cache sig y)
subCaches sig y
xs
      CacheCoprime LevelCoprime
level (Cache sig y, Cache sig y)
subCaches ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCoprime -> (Cache sig y, Cache sig y) -> sig y -> sig y
transformCoprime LevelCoprime
level (Cache sig y, Cache sig y)
subCaches sig y
xs
      CacheComposite LevelCacheComposite sig y
level (Cache sig y, Cache sig y)
subCaches ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
transformComposite LevelCacheComposite sig y
level (Cache sig y, Cache sig y)
subCaches sig y
xs


{- |
Memorize factorizations of the data size and permutation vectors.
-}
data Plan =
     PlanIdentity
   | PlanSmall LevelSmall
   | PlanNaive  -- mainly for debugging
   | PlanRadix2 LevelRadix2 Plan
   | PlanPrime LevelPrime Plan
   | PlanCoprime LevelCoprime (Plan, Plan)
   | PlanComposite LevelComposite (Plan, Plan)
   deriving (Int -> Plan -> ShowS
[Plan] -> ShowS
Plan -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [Plan] -> ShowS
$cshowList :: [Plan] -> ShowS
show :: Plan -> [Char]
$cshow :: Plan -> [Char]
showsPrec :: Int -> Plan -> ShowS
$cshowsPrec :: Int -> Plan -> ShowS
Show)

{-
efficient swallow comparison
only correct for Plans generated by 'plan'.
-}
instance Eq Plan where
   Plan
p0 == :: Plan -> Plan -> Bool
== Plan
p1  =  forall a. Ord a => a -> a -> Ordering
compare Plan
p0 Plan
p1 forall a. Eq a => a -> a -> Bool
== Ordering
EQ

{-
Needed for keys in CacheMap
-}
instance Ord Plan where
   compare :: Plan -> Plan -> Ordering
compare Plan
p0 Plan
p1  =
      case (Plan
p0,Plan
p1) of
         (Plan
PlanIdentity, Plan
PlanIdentity) -> Ordering
EQ
         (Plan
PlanIdentity, Plan
_) -> Ordering
LT
         (Plan
_, Plan
PlanIdentity) -> Ordering
GT
         (PlanSmall LevelSmall
l0, PlanSmall LevelSmall
l1) -> forall a. Ord a => a -> a -> Ordering
compare LevelSmall
l0 LevelSmall
l1
         (PlanSmall LevelSmall
_, Plan
_) -> Ordering
LT
         (Plan
_, PlanSmall LevelSmall
_) -> Ordering
GT
         (Plan
PlanNaive, Plan
PlanNaive) -> Ordering
EQ
         (Plan
PlanNaive, Plan
_) -> Ordering
LT
         (Plan
_, Plan
PlanNaive) -> Ordering
GT
         (PlanRadix2 LevelRadix2
l0 Plan
_, PlanRadix2 LevelRadix2
l1 Plan
_) -> forall a. Ord a => a -> a -> Ordering
compare LevelRadix2
l0 LevelRadix2
l1
         (PlanRadix2 LevelRadix2
_ Plan
_, Plan
_) -> Ordering
LT
         (Plan
_, PlanRadix2 LevelRadix2
_ Plan
_) -> Ordering
GT
         (PlanPrime LevelPrime
l0 Plan
_, PlanPrime LevelPrime
l1 Plan
_) -> forall a. Ord a => a -> a -> Ordering
compare LevelPrime
l0 LevelPrime
l1
         (PlanPrime LevelPrime
_ Plan
_, Plan
_) -> Ordering
LT
         (Plan
_, PlanPrime LevelPrime
_ Plan
_) -> Ordering
GT
         (PlanCoprime LevelCoprime
l0 (Plan, Plan)
_, PlanCoprime LevelCoprime
l1 (Plan, Plan)
_) -> forall a. Ord a => a -> a -> Ordering
compare LevelCoprime
l0 LevelCoprime
l1
         (PlanCoprime LevelCoprime
_ (Plan, Plan)
_, Plan
_) -> Ordering
LT
         (Plan
_, PlanCoprime LevelCoprime
_ (Plan, Plan)
_) -> Ordering
GT
         (PlanComposite LevelComposite
l0 (Plan, Plan)
_, PlanComposite LevelComposite
l1 (Plan, Plan)
_) -> forall a. Ord a => a -> a -> Ordering
compare LevelComposite
l0 LevelComposite
l1


plan :: Integer -> Plan
plan :: Integer -> Plan
plan Integer
n =
   forall s a. State s a -> s -> a
State.evalState (Integer -> State PlanMap Plan
planWithMapUpdate Integer
n) PlanMap
smallPlanMap

type PlanMap = Map.Map Integer Plan

smallPlanMap :: PlanMap
smallPlanMap :: PlanMap
smallPlanMap =
   forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] forall a b. (a -> b) -> a -> b
$
   Plan
PlanIdentity forall a. a -> [a] -> [a]
:
   Plan
PlanIdentity forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level2 forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level3 forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level4 forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level5 forall a. a -> [a] -> [a]
:
   []

{- |
Detect and re-use common sub-plans.
-}
planWithMap :: Integer -> State.State PlanMap Plan
planWithMap :: Integer -> State PlanMap Plan
planWithMap Integer
n =
   case forall a. C a => a -> a -> (a, a)
divMod Integer
n Integer
2 of
      (Integer
n2,Integer
0) -> LevelRadix2 -> Plan -> Plan
PlanRadix2 (Integer -> LevelRadix2
levelRadix2 Integer
n2) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Integer -> State PlanMap Plan
planWithMapUpdate Integer
n2
      (Integer, Integer)
_ ->
         let facs :: [(Integer, Integer)]
facs = Integer -> [(Integer, Integer)]
NumberTheory.fermatFactors Integer
n
         in  -- find unitary divisors
             case forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
a,Integer
b) -> Integer
aforall a. Ord a => a -> a -> Bool
>Integer
1 Bool -> Bool -> Bool
&& forall a. C a => a -> a -> a
gcd Integer
a Integer
b forall a. Eq a => a -> a -> Bool
== Integer
1) [(Integer, Integer)]
facs of
                (Integer, Integer)
q2 : [(Integer, Integer)]
_ ->
                   LevelCoprime -> (Plan, Plan) -> Plan
PlanCoprime ((Integer, Integer) -> LevelCoprime
levelCoprime (Integer, Integer)
q2) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   (Integer, Integer) -> State PlanMap (Plan, Plan)
planWithMapUpdate2 (Integer, Integer)
q2
                [(Integer, Integer)]
_ ->
                   let ((Integer, Integer)
q2 : [(Integer, Integer)]
_) = [(Integer, Integer)]
facs
                   in  if forall a b. (a, b) -> a
fst (Integer, Integer)
q2 forall a. Eq a => a -> a -> Bool
== Integer
1
                         then LevelPrime -> Plan -> Plan
PlanPrime (Integer -> LevelPrime
levelPrime forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> b
snd (Integer, Integer)
q2) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                              Integer -> State PlanMap Plan
planWithMapUpdate (Integer
nforall a. C a => a -> a -> a
-Integer
1)
                         else LevelComposite -> (Plan, Plan) -> Plan
PlanComposite ((Integer, Integer) -> LevelComposite
levelComposite (Integer, Integer)
q2) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                              (Integer, Integer) -> State PlanMap (Plan, Plan)
planWithMapUpdate2 (Integer, Integer)
q2

planWithMapUpdate :: Integer -> State.State PlanMap Plan
planWithMapUpdate :: Integer -> State PlanMap Plan
planWithMapUpdate Integer
n = do
   Maybe Plan
item <- forall (m :: * -> *) s a. Monad m => (s -> a) -> StateT s m a
State.gets (forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
n)
   case Maybe Plan
item of
      Just Plan
p -> forall (m :: * -> *) a. Monad m => a -> m a
return Plan
p
      Maybe Plan
Nothing ->
         Integer -> State PlanMap Plan
planWithMap Integer
n forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Plan
m -> forall (m :: * -> *) s. Monad m => (s -> s) -> StateT s m ()
State.modify (forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert Integer
n Plan
m) forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall (m :: * -> *) a. Monad m => a -> m a
return Plan
m

planWithMapUpdate2 :: (Integer, Integer) -> State.State PlanMap (Plan, Plan)
planWithMapUpdate2 :: (Integer, Integer) -> State PlanMap (Plan, Plan)
planWithMapUpdate2 =
   forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (,)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair (Integer -> State PlanMap Plan
planWithMapUpdate,Integer -> State PlanMap Plan
planWithMapUpdate)


{- |
Cache powers of the primitive root of unity
in a storage compatible to the processed signal.
-}
data Cache sig y =
     CacheIdentity
   | CacheSmall (LevelCacheSmall y)
   | CacheNaive (LevelCacheNaive y)
   | CacheRadix2 (LevelCacheRadix2 sig y) (Cache sig y)
   | CachePrime (LevelCachePrime sig y) (Cache sig y, Cache sig y)
   | CacheCoprime LevelCoprime (Cache sig y, Cache sig y)
   | CacheComposite (LevelCacheComposite sig y) (Cache sig y, Cache sig y)
   deriving (Int -> Cache sig y -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Cache sig y -> ShowS
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Cache sig y] -> ShowS
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Cache sig y -> [Char]
showList :: [Cache sig y] -> ShowS
$cshowList :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Cache sig y] -> ShowS
show :: Cache sig y -> [Char]
$cshow :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Cache sig y -> [Char]
showsPrec :: Int -> Cache sig y -> ShowS
$cshowsPrec :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Cache sig y -> ShowS
Show)

{- |
The expression @cacheForward prototype@
precomputes all data that is needed for forward Fourier transforms
for signals of the type and length @prototype@.
You can use this cache in 'transformWithCache'.
-}
cacheForward ::
   (Element y, SigG.Transform sig y) =>
   sig y -> Cache sig y
cacheForward :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> Cache sig y
cacheForward sig y
xs =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan
      (Integer -> Plan
plan forall a b. (a -> b) -> a -> b
$ forall a b. (C a, C b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall sig. Read sig => sig -> Int
SigG.length sig y
xs)
      (forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
xs)
      sig y
xs

{- |
See 'cacheForward'.
-}
cacheBackward ::
   (Element y, SigG.Transform sig y) =>
   sig y -> Cache sig y
cacheBackward :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> Cache sig y
cacheBackward sig y
xs =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan
      (Integer -> Plan
plan forall a b. (a -> b) -> a -> b
$ forall a b. (C a, C b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall sig. Read sig => sig -> Int
SigG.length sig y
xs)
      (forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
xs)
      sig y
xs

{- |
It is @(cacheForward x, cacheBackward x) = cacheDuplex x@
but 'cacheDuplex' shared common data of both caches.
-}
cacheDuplex ::
   (Element y, SigG.Transform sig y) =>
   sig y -> (Cache sig y, Cache sig y)
cacheDuplex :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> (Cache sig y, Cache sig y)
cacheDuplex sig y
xs =
   let p :: Plan
p = Integer -> Plan
plan forall a b. (a -> b) -> a -> b
$ forall a b. (C a, C b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall sig. Read sig => sig -> Int
SigG.length sig y
xs
       ((Direction, y)
z,(Direction, y)
zInv) = forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
xs
   in  forall s a. State s a -> s -> a
State.evalState
          (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 (Plan
p,Plan
p) ((Direction, y)
z,(Direction, y)
zInv) (sig y
xs,sig y
xs)) forall a b. (a -> b) -> a -> b
$
       forall k a. Map k a
Map.empty


data Direction = Forward | Backward
   deriving (Int -> Direction -> ShowS
[Direction] -> ShowS
Direction -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [Direction] -> ShowS
$cshowList :: [Direction] -> ShowS
show :: Direction -> [Char]
$cshow :: Direction -> [Char]
showsPrec :: Int -> Direction -> ShowS
$cshowsPrec :: Int -> Direction -> ShowS
Show, Direction -> Direction -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Direction -> Direction -> Bool
$c/= :: Direction -> Direction -> Bool
== :: Direction -> Direction -> Bool
$c== :: Direction -> Direction -> Bool
Eq, Eq Direction
Direction -> Direction -> Bool
Direction -> Direction -> Ordering
Direction -> Direction -> Direction
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Direction -> Direction -> Direction
$cmin :: Direction -> Direction -> Direction
max :: Direction -> Direction -> Direction
$cmax :: Direction -> Direction -> Direction
>= :: Direction -> Direction -> Bool
$c>= :: Direction -> Direction -> Bool
> :: Direction -> Direction -> Bool
$c> :: Direction -> Direction -> Bool
<= :: Direction -> Direction -> Bool
$c<= :: Direction -> Direction -> Bool
< :: Direction -> Direction -> Bool
$c< :: Direction -> Direction -> Bool
compare :: Direction -> Direction -> Ordering
$ccompare :: Direction -> Direction -> Ordering
Ord)

type CacheMap sig y = Map.Map (Plan,Direction) (Cache sig y)

cacheFromPlan ::
   (Element y, SigG.Transform sig y) =>
   Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan Plan
p (Direction, y)
z sig y
xs =
   forall s a. State s a -> s -> a
State.evalState (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate Plan
p (Direction, y)
z sig y
xs) forall a b. (a -> b) -> a -> b
$
   forall k a. Map k a
Map.empty

{- |
Detect and re-use common sub-caches.
-}
cacheFromPlanWithMap ::
   (Element y, SigG.Transform sig y) =>
   Plan -> (Direction,y) -> sig y ->
   State.State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMap :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMap Plan
p (Direction
d,y
z) sig y
xs =
   case Plan
p of
      Plan
PlanIdentity -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall (sig :: * -> *) y. Cache sig y
CacheIdentity
      PlanSmall LevelSmall
size -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall (sig :: * -> *) y. LevelCacheSmall y -> Cache sig y
CacheSmall forall a b. (a -> b) -> a -> b
$
         case LevelSmall
size of
            LevelSmall
Level2 -> forall y. y -> LevelCacheSmall y
LevelCache2 forall a b. (a -> b) -> a -> b
$ forall y. C y => y -> y
cache2 y
z
            LevelSmall
Level3 -> forall y. (y, y) -> LevelCacheSmall y
LevelCache3 forall a b. (a -> b) -> a -> b
$ forall y. C y => y -> (y, y)
cache3 y
z
            LevelSmall
Level4 -> forall y. (y, y, y) -> LevelCacheSmall y
LevelCache4 forall a b. (a -> b) -> a -> b
$ forall y. C y => y -> (y, y, y)
cache4 y
z
            LevelSmall
Level5 -> forall y. (y, y, y, y) -> LevelCacheSmall y
LevelCache5 forall a b. (a -> b) -> a -> b
$ forall y. C y => y -> (y, y, y, y)
cache5 y
z
      Plan
PlanNaive ->
         forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall (sig :: * -> *) y. LevelCacheNaive y -> Cache sig y
CacheNaive forall a b. (a -> b) -> a -> b
$ forall y. y -> LevelCacheNaive y
LevelCacheNaive y
z
      PlanRadix2 level :: LevelRadix2
level@(LevelRadix2 Int
n2) Plan
subPlan ->
         let subxs :: sig y
subxs = forall sig. Transform sig => Int -> sig -> sig
CutG.take Int
n2 sig y
xs
         in  forall (sig :: * -> *) y.
LevelCacheRadix2 sig y -> Cache sig y -> Cache sig y
CacheRadix2 (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y
levelCacheRadix2 LevelRadix2
level y
z sig y
subxs) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
             forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate Plan
subPlan (Direction
d,y
zforall a. C a => a -> a -> a
*y
z) sig y
subxs
      PlanPrime level :: LevelPrime
level@(LevelPrime (T
perm,T
_,T
_)) Plan
subPlan ->
         (\(Cache sig y, Cache sig y)
subCaches ->
            forall (sig :: * -> *) y.
LevelCachePrime sig y -> (Cache sig y, Cache sig y) -> Cache sig y
CachePrime
               (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y
levelCachePrime LevelPrime
level (forall a b. (a, b) -> a
fst (Cache sig y, Cache sig y)
subCaches) y
z sig y
xs) (Cache sig y, Cache sig y)
subCaches)
         forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         let subxs :: sig y
subxs = forall sig. Transform sig => Int -> sig -> sig
CutG.take (T -> Int
Permutation.size T
perm) sig y
xs
         in  forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 (Plan
subPlan,Plan
subPlan)
                (forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
subxs)
                (sig y
subxs,sig y
subxs)
      PlanCoprime level :: LevelCoprime
level@(LevelCoprime (Integer
n,Integer
m) (T, T, T)
_) (Plan, Plan)
subPlans ->
         forall (sig :: * -> *) y.
LevelCoprime -> (Cache sig y, Cache sig y) -> Cache sig y
CacheCoprime LevelCoprime
level forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 (Plan, Plan)
subPlans ((Direction
d,y
zforall a. C a => a -> Integer -> a
^Integer
m), (Direction
d,y
zforall a. C a => a -> Integer -> a
^Integer
n))
            (forall sig. Transform sig => Int -> sig -> sig
CutG.take (forall a. C a => Integer -> a
fromInteger Integer
n) sig y
xs, forall sig. Transform sig => Int -> sig -> sig
CutG.take (forall a. C a => Integer -> a
fromInteger Integer
m) sig y
xs)
      PlanComposite level :: LevelComposite
level@(LevelComposite (Integer
n,Integer
m) (T, T)
_) (Plan, Plan)
subPlans ->
         forall (sig :: * -> *) y.
LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> Cache sig y
CacheComposite (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelComposite -> y -> sig y -> LevelCacheComposite sig y
levelCacheComposite LevelComposite
level y
z sig y
xs) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 (Plan, Plan)
subPlans ((Direction
d,y
zforall a. C a => a -> Integer -> a
^Integer
m), (Direction
d,y
zforall a. C a => a -> Integer -> a
^Integer
n))
            (forall sig. Transform sig => Int -> sig -> sig
CutG.take (forall a. C a => Integer -> a
fromInteger Integer
n) sig y
xs, forall sig. Transform sig => Int -> sig -> sig
CutG.take (forall a. C a => Integer -> a
fromInteger Integer
m) sig y
xs)

cacheFromPlanWithMapUpdate ::
   (Element y, SigG.Transform sig y) =>
   Plan -> (Direction,y) -> sig y ->
   State.State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate Plan
p (Direction, y)
z sig y
xs = do
   let key :: (Plan, Direction)
key = (Plan
p, forall a b. (a, b) -> a
fst (Direction, y)
z)
   Maybe (Cache sig y)
item <- forall (m :: * -> *) s a. Monad m => (s -> a) -> StateT s m a
State.gets (forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup (Plan, Direction)
key)
   case Maybe (Cache sig y)
item of
      Just Cache sig y
c -> forall (m :: * -> *) a. Monad m => a -> m a
return Cache sig y
c
      Maybe (Cache sig y)
Nothing ->
         forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMap Plan
p (Direction, y)
z sig y
xs forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Cache sig y
m ->
         forall (m :: * -> *) s. Monad m => (s -> s) -> StateT s m ()
State.modify (forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert (Plan, Direction)
key Cache sig y
m) forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>>
         forall (m :: * -> *) a. Monad m => a -> m a
return Cache sig y
m

cacheFromPlanWithMapUpdate2 ::
   (Element y, SigG.Transform sig y) =>
   (Plan, Plan) -> ((Direction,y),(Direction,y)) -> (sig y, sig y) ->
   State.State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
cacheFromPlanWithMapUpdate2 (Plan
p0,Plan
p1) ((Direction, y)
z0,(Direction, y)
z1) (sig y
xs0,sig y
xs1) =
   forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (,)
      (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate Plan
p0 (Direction, y)
z0 sig y
xs0)
      (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
cacheFromPlanWithMapUpdate Plan
p1 (Direction, y)
z1 sig y
xs1)


newtype LevelCacheNaive y =
      LevelCacheNaive y
   deriving (Int -> LevelCacheNaive y -> ShowS
forall y. Show y => Int -> LevelCacheNaive y -> ShowS
forall y. Show y => [LevelCacheNaive y] -> ShowS
forall y. Show y => LevelCacheNaive y -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelCacheNaive y] -> ShowS
$cshowList :: forall y. Show y => [LevelCacheNaive y] -> ShowS
show :: LevelCacheNaive y -> [Char]
$cshow :: forall y. Show y => LevelCacheNaive y -> [Char]
showsPrec :: Int -> LevelCacheNaive y -> ShowS
$cshowsPrec :: forall y. Show y => Int -> LevelCacheNaive y -> ShowS
Show)

transformNaive ::
   (Element y, SigG.Transform sig y) =>
   LevelCacheNaive y -> sig y -> sig y
transformNaive :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheNaive y -> sig y -> sig y
transformNaive (LevelCacheNaive y
z) sig y
sig =
   forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> T a -> T b
SigS.map
      (forall a. C a => T a -> T a -> a
scalarProduct1 (forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
sig) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> y -> T y
powers sig y
sig)
      (forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> y -> T y
powers sig y
sig y
z)

scalarProduct1 ::
   (Ring.C a) =>
   SigS.T a -> SigS.T a -> a
scalarProduct1 :: forall a. C a => T a -> T a -> a
scalarProduct1 T a
xs T a
ys =
   forall x. (x -> x -> x) -> T x -> x
SigS.foldL1 forall a. C a => a -> a -> a
(+) forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> T a -> T b -> T c
SigS.zipWith forall a. C a => a -> a -> a
(*) T a
xs T a
ys

_transformRing ::
   (Ring.C y, SigG.Transform sig y) =>
   y -> sig y -> sig y
_transformRing :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
_transformRing y
z sig y
sig =
   forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
   forall y. C y => y -> T y -> T y
Ana.chirpTransform y
z forall a b. (a -> b) -> a -> b
$ forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
sig

powers ::
   (Element y, SigG.Read sig y) =>
   sig y -> y -> SigS.T y
powers :: forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> y -> T y
powers sig y
sig y
c = forall a. (a -> a) -> a -> T a
SigS.iterate (y
cforall a. C a => a -> a -> a
*) forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
multId sig y
sig


data LevelSmall = Level2 | Level3 | Level4 | Level5
   deriving (Int -> LevelSmall -> ShowS
[LevelSmall] -> ShowS
LevelSmall -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelSmall] -> ShowS
$cshowList :: [LevelSmall] -> ShowS
show :: LevelSmall -> [Char]
$cshow :: LevelSmall -> [Char]
showsPrec :: Int -> LevelSmall -> ShowS
$cshowsPrec :: Int -> LevelSmall -> ShowS
Show, LevelSmall -> LevelSmall -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: LevelSmall -> LevelSmall -> Bool
$c/= :: LevelSmall -> LevelSmall -> Bool
== :: LevelSmall -> LevelSmall -> Bool
$c== :: LevelSmall -> LevelSmall -> Bool
Eq, Eq LevelSmall
LevelSmall -> LevelSmall -> Bool
LevelSmall -> LevelSmall -> Ordering
LevelSmall -> LevelSmall -> LevelSmall
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: LevelSmall -> LevelSmall -> LevelSmall
$cmin :: LevelSmall -> LevelSmall -> LevelSmall
max :: LevelSmall -> LevelSmall -> LevelSmall
$cmax :: LevelSmall -> LevelSmall -> LevelSmall
>= :: LevelSmall -> LevelSmall -> Bool
$c>= :: LevelSmall -> LevelSmall -> Bool
> :: LevelSmall -> LevelSmall -> Bool
$c> :: LevelSmall -> LevelSmall -> Bool
<= :: LevelSmall -> LevelSmall -> Bool
$c<= :: LevelSmall -> LevelSmall -> Bool
< :: LevelSmall -> LevelSmall -> Bool
$c< :: LevelSmall -> LevelSmall -> Bool
compare :: LevelSmall -> LevelSmall -> Ordering
$ccompare :: LevelSmall -> LevelSmall -> Ordering
Ord, Int -> LevelSmall
LevelSmall -> Int
LevelSmall -> [LevelSmall]
LevelSmall -> LevelSmall
LevelSmall -> LevelSmall -> [LevelSmall]
LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall]
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall]
$cenumFromThenTo :: LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall]
enumFromTo :: LevelSmall -> LevelSmall -> [LevelSmall]
$cenumFromTo :: LevelSmall -> LevelSmall -> [LevelSmall]
enumFromThen :: LevelSmall -> LevelSmall -> [LevelSmall]
$cenumFromThen :: LevelSmall -> LevelSmall -> [LevelSmall]
enumFrom :: LevelSmall -> [LevelSmall]
$cenumFrom :: LevelSmall -> [LevelSmall]
fromEnum :: LevelSmall -> Int
$cfromEnum :: LevelSmall -> Int
toEnum :: Int -> LevelSmall
$ctoEnum :: Int -> LevelSmall
pred :: LevelSmall -> LevelSmall
$cpred :: LevelSmall -> LevelSmall
succ :: LevelSmall -> LevelSmall
$csucc :: LevelSmall -> LevelSmall
Enum)

data LevelCacheSmall y =
     LevelCache2 y
   | LevelCache3 (y,y)
   | LevelCache4 (y,y,y)
   | LevelCache5 (y,y,y,y)
   deriving (Int -> LevelCacheSmall y -> ShowS
forall y. Show y => Int -> LevelCacheSmall y -> ShowS
forall y. Show y => [LevelCacheSmall y] -> ShowS
forall y. Show y => LevelCacheSmall y -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelCacheSmall y] -> ShowS
$cshowList :: forall y. Show y => [LevelCacheSmall y] -> ShowS
show :: LevelCacheSmall y -> [Char]
$cshow :: forall y. Show y => LevelCacheSmall y -> [Char]
showsPrec :: Int -> LevelCacheSmall y -> ShowS
$cshowsPrec :: forall y. Show y => Int -> LevelCacheSmall y -> ShowS
Show)

cache2 :: (Ring.C y) => y -> y
cache3 :: (Ring.C y) => y -> (y,y)
cache4 :: (Ring.C y) => y -> (y,y,y)
cache5 :: (Ring.C y) => y -> (y,y,y,y)

cache2 :: forall y. C y => y -> y
cache2 y
z = y
z
cache3 :: forall y. C y => y -> (y, y)
cache3 y
z = (y
z, y
zforall a. C a => a -> a -> a
*y
z)
cache4 :: forall y. C y => y -> (y, y, y)
cache4 y
z = let z2 :: y
z2=y
zforall a. C a => a -> a -> a
*y
z in (y
z,y
z2,y
zforall a. C a => a -> a -> a
*y
z2)
cache5 :: forall y. C y => y -> (y, y, y, y)
cache5 y
z = let z2 :: y
z2=y
zforall a. C a => a -> a -> a
*y
z in (y
z,y
z2,y
zforall a. C a => a -> a -> a
*y
z2,y
z2forall a. C a => a -> a -> a
*y
z2)


transform2 ::
   (Ring.C y, SigG.Transform sig y) =>
   y -> sig y -> sig y
transform2 :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
transform2 y
z sig y
sig =
   let y
x0:y
x1:[y]
_ = forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
   in  forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall y. [y] -> T y
SigS.fromList [y
x0forall a. C a => a -> a -> a
+y
x1, y
x0forall a. C a => a -> a -> a
+y
zforall a. C a => a -> a -> a
*y
x1]

transform3 ::
   (Ring.C y, SigG.Transform sig y) =>
   (y,y) -> sig y -> sig y
transform3 :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y) -> sig y -> sig y
transform3 (y
z,y
z2) sig y
sig =
   let y
x0:y
x1:y
x2:[y]
_ = forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
{- Rader's algorithm with convolution by 2-size-Fourier-transform
       xf1 = x1+x2
       xf2 = x1-x2
       zf1 = z+z2
       zf2 = z-z2
       xzf1 = xf1*zf1
       xzf2 = xf2*zf2
       xz1 = (xzf1+xzf2)/2
       xz2 = (xzf1-xzf2)/2
-}
{- naive
       [x0+x1+x2, x0+z*x1+z2*x2, x0+z2*x1+z*x2]
-}
       ((y
s,y
_), (y
zx1,y
zx2)) = forall y. C y => Pair y -> Pair y -> (Pair y, Pair y)
Cyclic.sumAndConvolvePair (y
x1,y
x2) (y
z,y
z2)
   in  forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall y. [y] -> T y
SigS.fromList [y
x0forall a. C a => a -> a -> a
+y
s, y
x0forall a. C a => a -> a -> a
+y
zx1, y
x0forall a. C a => a -> a -> a
+y
zx2]

transform4 ::
   (Ring.C y, SigG.Transform sig y) =>
   (y,y,y) -> sig y -> sig y
transform4 :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y, y) -> sig y -> sig y
transform4 (y
z,y
z2,y
z3) sig y
sig =
   let y
x0:y
x1:y
x2:y
x3:[y]
_ = forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
       x02a :: y
x02a = y
x0forall a. C a => a -> a -> a
+y
x2; x02b :: y
x02b = y
x0forall a. C a => a -> a -> a
+y
z2forall a. C a => a -> a -> a
*y
x2
       x13a :: y
x13a = y
x1forall a. C a => a -> a -> a
+y
x3; x13b :: y
x13b = y
x1forall a. C a => a -> a -> a
+y
z2forall a. C a => a -> a -> a
*y
x3
   in  forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall y. [y] -> T y
SigS.fromList [y
x02aforall a. C a => a -> a -> a
+   y
x13a, y
x02bforall a. C a => a -> a -> a
+y
z forall a. C a => a -> a -> a
*y
x13b,
                      y
x02aforall a. C a => a -> a -> a
+y
z2forall a. C a => a -> a -> a
*y
x13a, y
x02bforall a. C a => a -> a -> a
+y
z3forall a. C a => a -> a -> a
*y
x13b]
{-
This needs also five multiplications,
but in complex numbers it is z=i, and thus multiplications are cheap
and we should better make use of distributive law in order to save additions.

       x02a = x0+x2; x02b = x0+z2*x2
       x1_2 = z2*x1; x3_2 = z2*x3
   in  SigG.takeStateMatch sig $
       SigS.fromList [x02a + x1   + x3  , x02b+z*(x1   + x3_2),
                      x02a + x1_2 + x3_2, x02b+z*(x1_2 + x3  )]
-}

{-
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 ::
   (Ring.C y, SigG.Transform sig y) =>
   (y,y,y,y) -> sig y -> sig y
transform5 :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
(y, y, y, y) -> sig y -> sig y
transform5 (y
z1,y
z2,y
z3,y
z4) sig y
sig =
   let y
x0:y
x1:y
x2:y
x3:y
x4:[y]
_ = forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
       ((y
s,y
_), (y
d1,y
d2,y
d4,y
d3)) =
          forall y.
C y =>
Quadruple y -> Quadruple y -> ((y, y), Quadruple y)
Cyclic.sumAndConvolveQuadruple (y
x1,y
x3,y
x4,y
x2) (y
z1,y
z2,y
z4,y
z3)
   in  forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall y. [y] -> T y
SigS.fromList [y
x0forall a. C a => a -> a -> a
+y
s, y
x0forall a. C a => a -> a -> a
+y
d1, y
x0forall a. C a => a -> a -> a
+y
d2, y
x0forall a. C a => a -> a -> a
+y
d3, y
x0forall a. C a => a -> a -> a
+y
d4]

{-
transform7

Toom-3-multiplication at the highest level and Karatsuba below?
Toom-2.5-multiplication with manual addition of the missing parts?

Toom-3-multiplication with complex interpolation nodes?
Still requires division by 4 and then complex multiplication in the frequency domain.
A:=matrix(5,5,[1,0,0,0,0,1,1,1,1,1,1,-1,1,-1,1,1,I,-1,-I,1,0,0,0,0,1]);
A:=matrix(5,5,[1,0,0,0,0,1,1,1,1,1,1,-1,1,-1,1,1,I,-1,-I,1,1,-I,-1,I,1]);

Karatsuba at three levels for convolution of signal of size 8 with zero padding?

Modify the 3x3 Fourier matrix by multiplying a regular matrix
to make it more convenient to work with?
We will hardly get rid of the irrational numbers.
-}

newtype LevelRadix2 = LevelRadix2 Int
   deriving (Int -> LevelRadix2 -> ShowS
[LevelRadix2] -> ShowS
LevelRadix2 -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelRadix2] -> ShowS
$cshowList :: [LevelRadix2] -> ShowS
show :: LevelRadix2 -> [Char]
$cshow :: LevelRadix2 -> [Char]
showsPrec :: Int -> LevelRadix2 -> ShowS
$cshowsPrec :: Int -> LevelRadix2 -> ShowS
Show, LevelRadix2 -> LevelRadix2 -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: LevelRadix2 -> LevelRadix2 -> Bool
$c/= :: LevelRadix2 -> LevelRadix2 -> Bool
== :: LevelRadix2 -> LevelRadix2 -> Bool
$c== :: LevelRadix2 -> LevelRadix2 -> Bool
Eq, Eq LevelRadix2
LevelRadix2 -> LevelRadix2 -> Bool
LevelRadix2 -> LevelRadix2 -> Ordering
LevelRadix2 -> LevelRadix2 -> LevelRadix2
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
$cmin :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
max :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
$cmax :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
>= :: LevelRadix2 -> LevelRadix2 -> Bool
$c>= :: LevelRadix2 -> LevelRadix2 -> Bool
> :: LevelRadix2 -> LevelRadix2 -> Bool
$c> :: LevelRadix2 -> LevelRadix2 -> Bool
<= :: LevelRadix2 -> LevelRadix2 -> Bool
$c<= :: LevelRadix2 -> LevelRadix2 -> Bool
< :: LevelRadix2 -> LevelRadix2 -> Bool
$c< :: LevelRadix2 -> LevelRadix2 -> Bool
compare :: LevelRadix2 -> LevelRadix2 -> Ordering
$ccompare :: LevelRadix2 -> LevelRadix2 -> Ordering
Ord)

levelRadix2 :: Integer -> LevelRadix2
levelRadix2 :: Integer -> LevelRadix2
levelRadix2 =
   Int -> LevelRadix2
LevelRadix2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (C a, C b) => a -> b
fromIntegral


data LevelCacheRadix2 sig y =
   LevelCacheRadix2 Int (sig y)
   deriving (Int -> LevelCacheRadix2 sig y -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheRadix2 sig y -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheRadix2 sig y] -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheRadix2 sig y -> [Char]
showList :: [LevelCacheRadix2 sig y] -> ShowS
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheRadix2 sig y] -> ShowS
show :: LevelCacheRadix2 sig y -> [Char]
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheRadix2 sig y -> [Char]
showsPrec :: Int -> LevelCacheRadix2 sig y -> ShowS
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheRadix2 sig y -> ShowS
Show)

levelCacheRadix2 ::
   (Element y, SigG.Transform sig y) =>
   LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y
levelCacheRadix2 :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y
levelCacheRadix2 (LevelRadix2 Int
n2) y
z sig y
sig =
   forall (sig :: * -> *) y. Int -> sig y -> LevelCacheRadix2 sig y
LevelCacheRadix2 Int
n2
      (forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> y -> T y
powers sig y
sig y
z)


{- |
Cooley-Tukey specialised to one factor of the size being 2.

Size of the input signal must be even.
-}
transformRadix2InterleavedFrequency ::
   (Element y, SigG.Transform sig y) =>
   LevelCacheRadix2 sig y -> Cache sig y -> sig y -> sig y
transformRadix2InterleavedFrequency :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheRadix2 sig y -> Cache sig y -> sig y -> sig y
transformRadix2InterleavedFrequency
      (LevelCacheRadix2 Int
n2 sig y
twiddle) Cache sig y
subCache sig y
sig =
   let (sig y
xs0,sig y
xs1) = forall sig. Transform sig => Int -> sig -> (sig, sig)
SigG.splitAt Int
n2 sig y
sig
       fs0 :: sig y
fs0 = forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache forall a b. (a -> b) -> a -> b
$ forall (sig :: * -> *) a b c.
(Read sig a, Transform sig b, Transform sig c) =>
(a -> b -> c) -> sig a -> sig b -> sig c
SigG.zipWith forall a. C a => a -> a -> a
(+) sig y
xs0 sig y
xs1
       fs1 :: sig y
fs1 = forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache forall a b. (a -> b) -> a -> b
$
                forall (sig :: * -> *) a b c.
(Read sig a, Read sig b, Transform sig c) =>
(a -> b -> c -> c) -> sig a -> sig b -> sig c -> sig c
SigG.zipWith3
                   (\y
w y
x0 y
x1 -> y
wforall a. C a => a -> a -> a
*(y
x0forall a. C a => a -> a -> a
-y
x1))
                   sig y
twiddle sig y
xs0 sig y
xs1
   in  forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall y. T y -> T y -> T y
SigS.interleave (forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
fs0) (forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
fs1)


data LevelComposite =
   LevelComposite
      (Integer, Integer)
      (Permutation.T, Permutation.T)
   deriving (Int -> LevelComposite -> ShowS
[LevelComposite] -> ShowS
LevelComposite -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelComposite] -> ShowS
$cshowList :: [LevelComposite] -> ShowS
show :: LevelComposite -> [Char]
$cshow :: LevelComposite -> [Char]
showsPrec :: Int -> LevelComposite -> ShowS
$cshowsPrec :: Int -> LevelComposite -> ShowS
Show)

instance Eq LevelComposite where
   LevelComposite
a == :: LevelComposite -> LevelComposite -> Bool
== LevelComposite
b  =  forall a. Ord a => a -> a -> Ordering
compare LevelComposite
a LevelComposite
b forall a. Eq a => a -> a -> Bool
== Ordering
EQ

instance Ord LevelComposite where
   compare :: LevelComposite -> LevelComposite -> Ordering
compare (LevelComposite (Integer, Integer)
a (T, T)
_) (LevelComposite (Integer, Integer)
b (T, T)
_)  =
      forall a. Ord a => a -> a -> Ordering
compare (Integer, Integer)
a (Integer, Integer)
b

levelComposite :: (Integer, Integer) -> LevelComposite
levelComposite :: (Integer, Integer) -> LevelComposite
levelComposite (Integer
n,Integer
m) =
   let ni :: Int
ni = forall a. C a => Integer -> a
fromInteger Integer
n
       mi :: Int
mi = forall a. C a => Integer -> a
fromInteger Integer
m
   in  (Integer, Integer) -> (T, T) -> LevelComposite
LevelComposite (Integer
n,Integer
m)
          (Int -> Int -> T
Permutation.transposition Int
ni Int
mi,
           Int -> Int -> T
Permutation.transposition Int
mi Int
ni)


data LevelCacheComposite sig y =
   LevelCacheComposite
      (Integer, Integer)
      (Permutation.T, Permutation.T)
      (sig y)
   deriving (Int -> LevelCacheComposite sig y -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheComposite sig y -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheComposite sig y] -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheComposite sig y -> [Char]
showList :: [LevelCacheComposite sig y] -> ShowS
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheComposite sig y] -> ShowS
show :: LevelCacheComposite sig y -> [Char]
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheComposite sig y -> [Char]
showsPrec :: Int -> LevelCacheComposite sig y -> ShowS
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheComposite sig y -> ShowS
Show)

levelCacheComposite ::
   (Element y, SigG.Transform sig y) =>
   LevelComposite -> y -> sig y -> LevelCacheComposite sig y
levelCacheComposite :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelComposite -> y -> sig y -> LevelCacheComposite sig y
levelCacheComposite (LevelComposite (Integer
n,Integer
m) (T, T)
transpose) y
z sig y
sig =
   forall (sig :: * -> *) y.
(Integer, Integer) -> (T, T) -> sig y -> LevelCacheComposite sig y
LevelCacheComposite (Integer
n,Integer
m) (T, T)
transpose forall a b. (a -> b) -> a -> b
$
   forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
   forall a b c. (a -> b -> c) -> b -> a -> c
flip forall acc y. (acc -> (y, acc)) -> acc -> T y
SigS.generateInfinite (Integer
n, forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
multId sig y
sig, forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
multId sig y
sig) forall a b. (a -> b) -> a -> b
$ \(Integer
i,y
zi,y
zij) ->
   (y
zij,
    case forall a. Enum a => a -> a
pred Integer
i of
      Integer
0 -> (Integer
n, y
ziforall a. C a => a -> a -> a
*y
z, forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
multId sig y
sig)
      Integer
i1 -> (Integer
i1, y
zi, y
zijforall a. C a => a -> a -> a
*y
zi))
{-
   {-# SCC "levelCacheComposite:rechunk" #-}
   concatRechunk sig $
   {-# SCC "levelCacheComposite:subpowers" #-}
   SigS.map
      (SigG.takeStateMatch (SigG.take (fromIntegral n) sig) . powers sig)
      ({-# SCC "levelCacheComposite:powers" #-}
       powers sig z)
-}
{-
   SigS.map
      (SigG.takeStateMatch sig . SigS.take (fromIntegral n) . powers sig)
      ({-# SCC "levelCacheComposite:powers" #-}
       powers sig z)
-}
{- suffers from big inefficiency of repeated 'append'
   SigG.takeStateMatch sig $
   SigS.fold $
   SigS.map (SigS.take (fromIntegral n) . powers sig) $
   SigS.take (fromIntegral m) $ -- necessary for strict storable vectors
   powers sig z
-}

{- |
For @transformComposite z (n,m) sig@,
it must hold @n*m == length sig@ and @z ^ length sig == 1@.

Cooley-Tukey-algorithm
-}
transformComposite ::
   (Element y, SigG.Transform sig y) =>
   LevelCacheComposite sig y -> (Cache sig y, Cache sig y) -> sig y -> sig y
transformComposite :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
transformComposite
      (LevelCacheComposite (Integer
n,Integer
m) (T
transposeNM, T
transposeMN) sig y
twiddle)
      (Cache sig y
subCacheN,Cache sig y
subCacheM) sig y
sig =
   forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeMN forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall (sig :: * -> *) y.
Transform sig y =>
sig y -> T (sig y) -> sig y
concatRechunk sig y
sig forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall a b. (a -> b) -> T a -> T b
SigS.map (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCacheM) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (forall a. C a => Integer -> a
fromInteger Integer
m) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeNM forall b c a. (b -> c) -> (a -> b) -> a -> c
.
--       concatRechunk sig .
       forall (sig :: * -> *) a b c.
(Read sig a, Transform sig b, Transform sig c) =>
(a -> b -> c) -> sig a -> sig b -> sig c
SigG.zipWith forall a. C a => a -> a -> a
(*) sig y
twiddle forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall m. Monoid m => T m -> m
SigS.fold forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall a b. (a -> b) -> T a -> T b
SigS.map (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCacheN) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (forall a. C a => Integer -> a
fromInteger Integer
n) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeMN forall a b. (a -> b) -> a -> b
$
       sig y
sig


data LevelCoprime =
   LevelCoprime
      (Integer, Integer)
      (Permutation.T, Permutation.T, Permutation.T)
   deriving (Int -> LevelCoprime -> ShowS
[LevelCoprime] -> ShowS
LevelCoprime -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelCoprime] -> ShowS
$cshowList :: [LevelCoprime] -> ShowS
show :: LevelCoprime -> [Char]
$cshow :: LevelCoprime -> [Char]
showsPrec :: Int -> LevelCoprime -> ShowS
$cshowsPrec :: Int -> LevelCoprime -> ShowS
Show)

instance Eq LevelCoprime where
   LevelCoprime
a == :: LevelCoprime -> LevelCoprime -> Bool
== LevelCoprime
b  =  forall a. Ord a => a -> a -> Ordering
compare LevelCoprime
a LevelCoprime
b forall a. Eq a => a -> a -> Bool
== Ordering
EQ

instance Ord LevelCoprime where
   compare :: LevelCoprime -> LevelCoprime -> Ordering
compare (LevelCoprime (Integer, Integer)
a (T, T, T)
_) (LevelCoprime (Integer, Integer)
b (T, T, T)
_)  =
      forall a. Ord a => a -> a -> Ordering
compare (Integer, Integer)
a (Integer, Integer)
b

{-
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
-}
levelCoprime :: (Integer, Integer) -> LevelCoprime
levelCoprime :: (Integer, Integer) -> LevelCoprime
levelCoprime (Integer
n,Integer
m) =
   let ni :: Int
ni = forall a. C a => Integer -> a
fromInteger Integer
n
       mi :: Int
mi = forall a. C a => Integer -> a
fromInteger Integer
m
   in  (Integer, Integer) -> (T, T, T) -> LevelCoprime
LevelCoprime (Integer
n,Integer
m)
          (Int -> Int -> T
Permutation.skewGrid Int
mi Int
ni,
           Int -> Int -> T
Permutation.transposition Int
ni Int
mi,
           Int -> Int -> T
Permutation.skewGridCRTInv Int
ni Int
mi)


{- |
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
-}
{-
A very elegant way would be to divide the signal into chunks of size n,
define ring operations on these chunks
and perform one (length/n)-size-sub-transform in this chunk-ring.
This way we would also only have to plan the sub-transform once.
On StorableVectors the chunking could be performed in-place
in terms of a virtual reshape operation.
In the general case the performance can become very bad
if the chunks are very small, say 2 or 3 elements.
-}
transformCoprime ::
   (Element y, SigG.Transform sig y) =>
   LevelCoprime -> (Cache sig y, Cache sig y) -> sig y -> sig y
transformCoprime :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCoprime -> (Cache sig y, Cache sig y) -> sig y -> sig y
transformCoprime
      (LevelCoprime (Integer
n,Integer
m) (T
grid, T
transpose, T
gridInv)) (Cache sig y
subCacheN,Cache sig y
subCacheM) =
   let subTransform :: Cache sig y -> a -> sig y -> sig y
subTransform Cache sig y
cache a
j sig y
sig =
          forall (sig :: * -> *) y.
Transform sig y =>
sig y -> T (sig y) -> sig y
concatRechunk sig y
sig forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          forall a b. (a -> b) -> T a -> T b
SigS.map (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (forall a b. (C a, C b) => a -> b
fromIntegral a
j) forall a b. (a -> b) -> a -> b
$ sig y
sig
   in  forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
gridInv forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall {sig :: * -> *} {y} {a}.
(Transform sig y, Element y, C a) =>
Cache sig y -> a -> sig y -> sig y
subTransform Cache sig y
subCacheM Integer
m forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transpose forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall {sig :: * -> *} {y} {a}.
(Transform sig y, Element y, C a) =>
Cache sig y -> a -> sig y -> sig y
subTransform Cache sig y
subCacheN Integer
n forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
grid


-- concatenate and reorganize for faster indexing
concatRechunk ::
   (SigG.Transform sig y) =>
   sig y -> SigS.T (sig y) -> sig y
concatRechunk :: forall (sig :: * -> *) y.
Transform sig y =>
sig y -> T (sig y) -> sig y
concatRechunk sig y
pattern =
   forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
pattern forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall m. Monoid m => T m -> m
SigS.fold


data LevelPrime =
   LevelPrime (Permutation.T, Permutation.T, Permutation.T)
      deriving (Int -> LevelPrime -> ShowS
[LevelPrime] -> ShowS
LevelPrime -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [LevelPrime] -> ShowS
$cshowList :: [LevelPrime] -> ShowS
show :: LevelPrime -> [Char]
$cshow :: LevelPrime -> [Char]
showsPrec :: Int -> LevelPrime -> ShowS
$cshowsPrec :: Int -> LevelPrime -> ShowS
Show)

instance Eq LevelPrime where
   LevelPrime
a == :: LevelPrime -> LevelPrime -> Bool
== LevelPrime
b  =  forall a. Ord a => a -> a -> Ordering
compare LevelPrime
a LevelPrime
b forall a. Eq a => a -> a -> Bool
== Ordering
EQ

instance Ord LevelPrime where
   compare :: LevelPrime -> LevelPrime -> Ordering
compare (LevelPrime (T
a,T
_,T
_)) (LevelPrime (T
b,T
_,T
_))  =
      forall a. Ord a => a -> a -> Ordering
compare (T -> Int
Permutation.size T
a) (T -> Int
Permutation.size T
b)

{-
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.
-}
levelPrime :: Integer -> LevelPrime
levelPrime :: Integer -> LevelPrime
levelPrime Integer
n =
   let perm :: T
perm = Int -> T
Permutation.multiplicative forall a b. (a -> b) -> a -> b
$ forall a b. (C a, C b) => a -> b
fromIntegral Integer
n
   in  (T, T, T) -> LevelPrime
LevelPrime
          (T
perm, T -> T
Permutation.reverse T
perm, T -> T
Permutation.inverse T
perm)


data LevelCachePrime sig y =
   LevelCachePrime (Permutation.T, Permutation.T) (sig y)
      deriving (Int -> LevelCachePrime sig y -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCachePrime sig y -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCachePrime sig y] -> ShowS
forall (sig :: * -> *) y.
Show (sig y) =>
LevelCachePrime sig y -> [Char]
showList :: [LevelCachePrime sig y] -> ShowS
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCachePrime sig y] -> ShowS
show :: LevelCachePrime sig y -> [Char]
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCachePrime sig y -> [Char]
showsPrec :: Int -> LevelCachePrime sig y -> ShowS
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCachePrime sig y -> ShowS
Show)

levelCachePrime ::
   (Element y, SigG.Transform sig y) =>
   LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y
levelCachePrime :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y
levelCachePrime (LevelPrime (T
perm, T
rev, T
inv)) Cache sig y
subCache y
z sig y
sig =
   forall (sig :: * -> *) y. (T, T) -> sig y -> LevelCachePrime sig y
LevelCachePrime (T
rev, T
inv)
      ((\sig y
zs -> forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
recipInteger sig y
zs) sig y
zs) forall a b. (a -> b) -> a -> b
$
       forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache forall a b. (a -> b) -> a -> b
$
       forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm forall a b. (a -> b) -> a -> b
$
       forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig forall a b. (a -> b) -> a -> b
$
       forall a. (a -> a) -> a -> T a
SigS.iterate (y
zforall a. C a => a -> a -> a
*) y
z)

{- |
Rader's algorithm for prime length signals.
-}
transformPrime ::
   (Element y, SigG.Transform sig y) =>
   LevelCachePrime sig y -> (Cache sig y, Cache sig y) -> sig y -> sig y
transformPrime :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelCachePrime sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
transformPrime (LevelCachePrime (T
rev, T
inv) sig y
zs) (Cache sig y, Cache sig y)
subCaches =
   forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL (forall a. HasCallStack => [Char] -> a
error [Char]
"transformPrime: empty signal") forall a b. (a -> b) -> a -> b
$
   \y
x0 sig y
rest ->
      forall (sig :: * -> *) y.
(Transform0 sig, Storage (sig y)) =>
y -> sig y -> sig y
SigG.cons (forall (sig :: * -> *) y s.
(Read0 sig, Storage (sig y)) =>
(s -> y -> s) -> s -> sig y -> s
SigG.foldL forall a. C a => a -> a -> a
(+) y
x0 sig y
rest) forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y0 y1.
(Transform0 sig, Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
SigG.map (y
x0forall a. C a => a -> a -> a
+) forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
inv forall a b. (a -> b) -> a -> b
$
      forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveSpectrumCyclicCache (Cache sig y, Cache sig y)
subCaches sig y
zs forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
rev sig y
rest

{-
Cyclic.reverse xs = shiftR 1 (reverse xs)
Cyclic.reverse (xs <*> ys) = Cyclic.reverse xs <*> Cyclic.reverse ys
Cyclic.reverse (Cyclic.reverse xs) = xs

We could move the 'Cyclic.reverse' over to the z-vector,
but then we would have to reverse again after convolution.

zs <*> Cyclic.reverse rest
 = Cyclic.reverse (Cyclic.reverse zs <*> rest)
-}

{-
This uses Cyclic.filter instead of Cyclic.convolve.
This is simpler, but Fourier.convolveCyclic is a bit simpler than Fourier.filterCyclic,
since it does not need to reverse an operand.
-}
_transformPrimeAlt ::
   (Ring.C y, SigG.Transform sig y) =>
   LevelPrime -> y -> sig y -> sig y
_transformPrimeAlt :: forall y (sig :: * -> *).
(C y, Transform sig y) =>
LevelPrime -> y -> sig y -> sig y
_transformPrimeAlt (LevelPrime (T
perm, T
_, T
inv)) y
z =
   forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL (forall a. HasCallStack => [Char] -> a
error [Char]
"transformPrime: empty signal") forall a b. (a -> b) -> a -> b
$
   \y
x0 sig y
rest ->
      forall (sig :: * -> *) y.
(Transform0 sig, Storage (sig y)) =>
y -> sig y -> sig y
SigG.cons (forall (sig :: * -> *) y s.
(Read0 sig, Storage (sig y)) =>
(s -> y -> s) -> s -> sig y -> s
SigG.foldL forall a. C a => a -> a -> a
(+) y
x0 sig y
rest) forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y0 y1.
(Transform0 sig, Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
SigG.map (y
x0forall a. C a => a -> a -> a
+) forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
inv forall a b. (a -> b) -> a -> b
$
      forall (sig :: * -> *) y.
(Transform sig y, C y) =>
sig y -> sig y -> sig y
Cyclic.filterNaive
         (forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm sig y
rest)
         (forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm (forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
rest (forall a. (a -> a) -> a -> T a
SigS.iterate (y
zforall a. C a => a -> a -> a
*) y
z)))



{- |
Filter window stored as spectrum
such that it can be applied efficiently to long signals.
-}
data Window sig y =
   Window Int (Cache sig y, Cache sig y) (sig y)
   deriving (Int -> Window sig y -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Window sig y -> ShowS
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Window sig y] -> ShowS
forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Window sig y -> [Char]
showList :: [Window sig y] -> ShowS
$cshowList :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Window sig y] -> ShowS
show :: Window sig y -> [Char]
$cshow :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Window sig y -> [Char]
showsPrec :: Int -> Window sig y -> ShowS
$cshowsPrec :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Window sig y -> ShowS
Show)


window ::
   (Element y, SigG.Transform sig y) =>
   sig y -> Window sig y
window :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> Window sig y
window sig y
x =
   if forall sig. Read sig => sig -> Bool
CutG.null sig y
x
     then forall (sig :: * -> *) y.
Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
Window Int
0 (forall (sig :: * -> *) y. Cache sig y
CacheIdentity, forall (sig :: * -> *) y. Cache sig y
CacheIdentity) forall sig. Monoid sig => sig
CutG.empty
     else
       let size :: Int
size  = forall sig. Read sig => sig -> Int
CutG.length sig y
x
           size2 :: Int
size2 = Int
2 forall a. C a => a -> a -> a
* forall a. (C a, Bits a) => a -> a
NumberTheory.ceilingPowerOfTwo Int
size
           padded :: sig y
padded =
              forall sig. Transform sig => Int -> sig -> sig
SigG.take Int
size2 forall a b. (a -> b) -> a -> b
$
              forall sig. Monoid sig => sig -> sig -> sig
CutG.append sig y
x forall a b. (a -> b) -> a -> b
$
                 let pad :: sig y
pad = forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
x forall a b. (a -> b) -> a -> b
$ forall a. a -> T a
SigS.repeat forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
addId sig y
x
                 in  forall sig. Monoid sig => sig -> sig -> sig
CutG.append sig y
pad (forall sig. Monoid sig => sig -> sig -> sig
SigG.append sig y
pad sig y
pad)
           caches :: (Cache sig y, Cache sig y)
caches@(Cache sig y
cache, Cache sig y
_cacheInv) =
              forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> (Cache sig y, Cache sig y)
cacheDuplex sig y
padded
       in  forall (sig :: * -> *) y.
Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
Window
              (Int
size2forall a. C a => a -> a -> a
-Int
sizeforall a. C a => a -> a -> a
+Int
1)
              (Cache sig y, Cache sig y)
caches
              (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache forall a b. (a -> b) -> a -> b
$
               forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
recipInteger sig y
padded) sig y
padded)

{- |
Efficient convolution of a large filter window
with a probably infinite signal.
-}
convolveWithWindow ::
   (Element y, SigG.Transform sig y) =>
   Window sig y -> sig y -> sig y
convolveWithWindow :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Window sig y -> sig y -> sig y
convolveWithWindow (Window Int
blockSize (Cache sig y, Cache sig y)
caches sig y
spectrum) sig y
b =
   if Int
blockSizeforall a. Eq a => a -> a -> Bool
==forall a. C a => a
zero
     then forall sig. Monoid sig => sig
CutG.empty
     else
       let windowSize :: Int
windowSize = forall sig. Read sig => sig -> Int
SigG.length sig y
spectrum forall a. C a => a -> a -> a
- Int
blockSize
       in  forall x acc. (x -> acc -> acc) -> acc -> T x -> acc
SigS.foldR (forall a (sig :: * -> *).
(C a, Transform sig a) =>
Int -> sig a -> sig a -> sig a
FiltNRG.addShiftedSimple Int
blockSize) forall sig. Monoid sig => sig
CutG.empty forall a b. (a -> b) -> a -> b
$
           forall a b. (a -> b) -> T a -> T b
SigS.map
              (\sig y
block ->
                 forall sig. Transform sig => Int -> sig -> sig
SigG.take (Int
windowSize forall a. C a => a -> a -> a
+ forall sig. Read sig => sig -> Int
SigG.length sig y
block) forall a b. (a -> b) -> a -> b
$
                 forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveSpectrumCyclicCache (Cache sig y, Cache sig y)
caches sig y
spectrum forall a b. (a -> b) -> a -> b
$
                 forall a b c. (a -> b -> c) -> b -> a -> c
flip forall sig. Monoid sig => sig -> sig -> sig
CutG.append
                    {-
                    The last block may be shorter than blockSize
                    and thus needs more padding.
                    -}
                    (forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
spectrum forall a b. (a -> b) -> a -> b
$ forall a. a -> T a
SigS.repeat forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
addId sig y
b) forall a b. (a -> b) -> a -> b
$
                 sig y
block) forall a b. (a -> b) -> a -> b
$
           forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical Int
blockSize sig y
b


{- |
Signal must have equal size and must not be empty.
-}
convolveCyclic ::
   (Element y, SigG.Transform sig y) =>
   sig y -> sig y -> sig y
convolveCyclic :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> sig y -> sig y
convolveCyclic sig y
x =
   let len :: Integer
len = forall a b. (C a, C b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall sig. Read sig => sig -> Int
SigG.length sig y
x
       ((Direction, y)
z,(Direction, y)
zInv) =
          forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
x
   in  forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveCyclicCache
          (forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan (Integer -> Plan
plan Integer
len) (Direction, y)
z sig y
x,
           forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan (Integer -> Plan
plan Integer
len) (Direction, y)
zInv sig y
x)
          sig y
x

convolveCyclicCache ::
   (Element y, SigG.Transform sig y) =>
   (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveCyclicCache :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveCyclicCache (Cache sig y, Cache sig y)
caches sig y
x =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveSpectrumCyclicCache (Cache sig y, Cache sig y)
caches forall a b. (a -> b) -> a -> b
$
   forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
recipInteger sig y
x) forall a b. (a -> b) -> a -> b
$ forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (forall a b. (a, b) -> a
fst (Cache sig y, Cache sig y)
caches) sig y
x

{- |
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 ::
   (Element y, SigG.Transform sig y) =>
   (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveSpectrumCyclicCache :: forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveSpectrumCyclicCache (Cache sig y
cache,Cache sig y
cacheInv) sig y
x sig y
y =
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cacheInv forall a b. (a -> b) -> a -> b
$
   forall (sig :: * -> *) a b c.
(Read sig a, Transform sig b, Transform sig c) =>
(a -> b -> c) -> sig a -> sig b -> sig c
SigG.zipWith forall a. C a => a -> a -> a
(*) sig y
x forall a b. (a -> b) -> a -> b
$
   forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache sig y
y

{-
Test:

let xs = [0,1,0,0,0,0 :: Complex.T Double]; z = fst $ conjugatePrimitiveRootsOfUnity xs in print (transformNaive z xs) >> print (transformCoprime z (2,3) xs)
-}