{- |
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 = a -> a
forall a. C a => a -> a
recip (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (sig (T a) -> Int
forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig)) a -> a -> T a
forall a. a -> a -> T a
+: a
forall a. C a => a
zero
   addId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
addId sig (T a)
_sig = T a
forall a. C a => a
zero
   multId :: forall (sig :: * -> *). Read sig (T a) => sig (T a) -> T a
multId sig (T a)
_sig = T a
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, T a -> T a
forall a. C a => T a -> T a
Complex.conjugate T a
x)) (T a -> (T a, T a)) -> T a -> (T a, T a)
forall a b. (a -> b) -> a -> b
$
      case sig (T a) -> Int
forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig of
         Int
1 -> T a
forall a. C a => a
one
         Int
2 -> T a -> T a
forall a. C a => a -> a
negate T a
forall a. C a => a
one
         Int
3 -> (a -> a
forall a. C a => a -> a
negate a
forall a. C a => a
one a -> a -> T a
forall a. a -> a -> T a
+: a -> a
forall a. C a => a -> a
sqrt a
3) T a -> T a -> T a
forall a. C a => a -> a -> a
/ T a
2
         Int
4 -> a
forall a. C a => a
zero a -> a -> T a
forall a. a -> a -> T a
+: a
forall a. C a => a
one
         Int
5 ->
            let sqrt5 :: a
sqrt5 = a -> a
forall a. C a => a -> a
sqrt a
5
            in  ((a
sqrt5 a -> a -> a
forall a. C a => a -> a -> a
- a
1) a -> a -> T a
forall a. a -> a -> T a
+: a -> a
forall a. C a => a -> a
sqrt a
2 a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
sqrt(a
5 a -> a -> a
forall a. C a => a -> a -> a
+ a
sqrt5)) T a -> T a -> T a
forall a. C a => a -> a -> a
/ T a
4
         Int
6 -> (a
forall a. C a => a
one a -> a -> T a
forall a. a -> a -> T a
+: a -> a
forall a. C a => a -> a
sqrt a
3) T a -> T a -> T a
forall a. C a => a -> a -> a
/ T a
2
         Int
8 -> a -> T a -> T a
forall a. C a => a -> T a -> T a
Complex.scale (a -> a
forall a. C a => a -> a
sqrt a
2 a -> a -> a
forall a. C a => a -> a -> a
/ a
2) (a
forall a. C a => a
one a -> a -> T a
forall a. a -> a -> T a
+: a
forall a. C a => a
one)
         Int
12 -> (a -> a
forall a. C a => a -> a
sqrt a
3 a -> a -> T a
forall a. a -> a -> T a
+: a
forall a. C a => a
one) T a -> T a -> T a
forall a. C a => a -> a -> a
/ T a
2
         Int
n -> a -> T a
forall a. C a => a -> T a
Complex.cis (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
forall a. C a => a
pi a -> a -> a
forall a. C a => a -> a -> a
/ Int -> 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 =
      T a -> T a
forall a. C a => a -> a
recip (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (sig (T a) -> Int
forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig) a -> a -> T a
forall a. C a => a -> a -> T a
/: T a -> a
forall a. T a -> a
RC.modulus (sig (T a) -> T a
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 = a
forall a. C a => a
zero a -> a -> T a
forall a. C a => a -> a -> T a
/: T a -> a
forall a. T a -> a
RC.modulus (sig (T a) -> T a
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 = a
forall a. C a => a
one a -> a -> T a
forall a. C a => a -> a -> T a
/: T a -> a
forall a. T a -> a
RC.modulus (sig (T a) -> T a
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 = T a -> a
forall a. T a -> a
RC.modulus (sig (T a) -> T a
forall (sig :: * -> *) y. Read sig y => sig y -> y
head sig (T a)
sig)
          order :: Order
order@(NumberTheory.Order Integer
expo) =
             a -> Order
forall a. PrimitiveRoot a => a -> Order
NumberTheory.maximumOrderOfPrimitiveRootsOfUnity a
modu
          a
r:[a]
_ = a -> Order -> [a]
forall a. (PrimitiveRoot a, Eq a) => a -> Order -> [a]
NumberTheory.primitiveRootsOfUnity a
modu Order
order
          n :: Integer
n = Integer -> Integer -> Integer
forall a. (C a, C a) => a -> a -> a
Integral.divChecked Integer
expo (Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (sig (T a) -> Int
forall sig. Read sig => sig -> Int
SigG.length sig (T a)
sig))
          z :: T a
z = (a
r a -> a -> T a
forall a. C a => a -> a -> T a
/: a
modu) T a -> Integer -> T a
forall a. C a => a -> Integer -> a
^ Integer
n
      in  (T a
z, T a -> T a
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 =
   y -> (y -> T y -> y) -> T y -> y
forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL ([Char] -> y
forall a. HasCallStack => [Char] -> a
error [Char]
"Generic.Signal.head: empty signal") y -> T y -> y
forall a b. a -> b -> a
const (T y -> y) -> (sig y -> T y) -> sig y -> y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   sig y -> T y
forall y. Storage (sig y) => sig y -> T y
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) = sig y -> (y, y)
forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> (y, y)
forall (sig :: * -> *). 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 =
   Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (sig y -> Cache sig y
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 =
   Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (sig y -> Cache sig y
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 =
   Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache (Plan -> (Direction, y) -> sig y -> Cache sig y
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 -> y -> sig y -> sig y
forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
transform2 y
zs sig y
xs
            LevelCache3 (y, y)
zs -> (y, y) -> sig y -> sig y
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 -> (y, y, y) -> sig y -> sig y
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 -> (y, y, y, y) -> sig y -> sig y
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 ->
         LevelCacheNaive y -> sig y -> sig y
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 ->
         LevelCacheRadix2 sig y -> Cache sig y -> sig y -> sig y
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 ->
         LevelCachePrime sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
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 ->
         LevelCoprime -> (Cache sig y, Cache sig y) -> sig y -> sig y
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 ->
         LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> sig y -> sig y
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]
(Int -> Plan -> ShowS)
-> (Plan -> [Char]) -> ([Plan] -> ShowS) -> Show Plan
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Plan -> ShowS
showsPrec :: Int -> Plan -> ShowS
$cshow :: Plan -> [Char]
show :: Plan -> [Char]
$cshowList :: [Plan] -> ShowS
showList :: [Plan] -> ShowS
Show)

{-
efficient swallow comparison
only correct for Plans generated by 'plan'.
-}
instance Eq Plan where
   Plan
p0 == :: Plan -> Plan -> Bool
== Plan
p1  =  Plan -> Plan -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Plan
p0 Plan
p1 Ordering -> Ordering -> Bool
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) -> LevelSmall -> LevelSmall -> Ordering
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
_) -> LevelRadix2 -> LevelRadix2 -> Ordering
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
_) -> LevelPrime -> LevelPrime -> Ordering
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)
_) -> LevelCoprime -> LevelCoprime -> Ordering
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)
_) -> LevelComposite -> LevelComposite -> Ordering
forall a. Ord a => a -> a -> Ordering
compare LevelComposite
l0 LevelComposite
l1


plan :: Integer -> Plan
plan :: Integer -> Plan
plan Integer
n =
   State PlanMap Plan -> PlanMap -> Plan
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 =
   [(Integer, Plan)] -> PlanMap
forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList ([(Integer, Plan)] -> PlanMap) -> [(Integer, Plan)] -> PlanMap
forall a b. (a -> b) -> a -> b
$ [Integer] -> [Plan] -> [(Integer, Plan)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] ([Plan] -> [(Integer, Plan)]) -> [Plan] -> [(Integer, Plan)]
forall a b. (a -> b) -> a -> b
$
   Plan
PlanIdentity Plan -> [Plan] -> [Plan]
forall a. a -> [a] -> [a]
:
   Plan
PlanIdentity Plan -> [Plan] -> [Plan]
forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level2 Plan -> [Plan] -> [Plan]
forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level3 Plan -> [Plan] -> [Plan]
forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level4 Plan -> [Plan] -> [Plan]
forall a. a -> [a] -> [a]
:
   LevelSmall -> Plan
PlanSmall LevelSmall
Level5 Plan -> [Plan] -> [Plan]
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 Integer -> Integer -> (Integer, Integer)
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) (Plan -> Plan) -> State PlanMap Plan -> State PlanMap Plan
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 ((Integer, Integer) -> Bool)
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
a,Integer
b) -> Integer
aInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
1 Bool -> Bool -> Bool
&& Integer -> Integer -> Integer
forall a. C a => a -> a -> a
gcd Integer
a Integer
b Integer -> Integer -> Bool
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) ((Plan, Plan) -> Plan)
-> StateT PlanMap Identity (Plan, Plan) -> State PlanMap Plan
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                   (Integer, Integer) -> StateT PlanMap Identity (Plan, Plan)
planWithMapUpdate2 (Integer, Integer)
q2
                [(Integer, Integer)]
_ ->
                   let ((Integer, Integer)
q2 : [(Integer, Integer)]
_) = [(Integer, Integer)]
facs
                   in  if (Integer, Integer) -> Integer
forall a b. (a, b) -> a
fst (Integer, Integer)
q2 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1
                         then LevelPrime -> Plan -> Plan
PlanPrime (Integer -> LevelPrime
levelPrime (Integer -> LevelPrime) -> Integer -> LevelPrime
forall a b. (a -> b) -> a -> b
$ (Integer, Integer) -> Integer
forall a b. (a, b) -> b
snd (Integer, Integer)
q2) (Plan -> Plan) -> State PlanMap Plan -> State PlanMap Plan
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                              Integer -> State PlanMap Plan
planWithMapUpdate (Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1)
                         else LevelComposite -> (Plan, Plan) -> Plan
PlanComposite ((Integer, Integer) -> LevelComposite
levelComposite (Integer, Integer)
q2) ((Plan, Plan) -> Plan)
-> StateT PlanMap Identity (Plan, Plan) -> State PlanMap Plan
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
                              (Integer, Integer) -> StateT PlanMap Identity (Plan, Plan)
planWithMapUpdate2 (Integer, Integer)
q2

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

planWithMapUpdate2 :: (Integer, Integer) -> State.State PlanMap (Plan, Plan)
planWithMapUpdate2 :: (Integer, Integer) -> StateT PlanMap Identity (Plan, Plan)
planWithMapUpdate2 =
   (State PlanMap Plan
 -> State PlanMap Plan -> StateT PlanMap Identity (Plan, Plan))
-> (State PlanMap Plan, State PlanMap Plan)
-> StateT PlanMap Identity (Plan, Plan)
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Plan -> Plan -> (Plan, Plan))
-> State PlanMap Plan
-> State PlanMap Plan
-> StateT PlanMap Identity (Plan, Plan)
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (,)) ((State PlanMap Plan, State PlanMap Plan)
 -> StateT PlanMap Identity (Plan, Plan))
-> ((Integer, Integer) -> (State PlanMap Plan, State PlanMap Plan))
-> (Integer, Integer)
-> StateT PlanMap Identity (Plan, Plan)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Integer -> State PlanMap Plan, Integer -> State PlanMap Plan)
-> (Integer, Integer) -> (State PlanMap Plan, State PlanMap Plan)
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
[Cache sig y] -> ShowS
Cache sig y -> [Char]
(Int -> Cache sig y -> ShowS)
-> (Cache sig y -> [Char])
-> ([Cache sig y] -> ShowS)
-> Show (Cache sig y)
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]
$cshowsPrec :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Cache sig y -> ShowS
showsPrec :: Int -> Cache sig y -> ShowS
$cshow :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Cache sig y -> [Char]
show :: Cache sig y -> [Char]
$cshowList :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Cache sig y] -> ShowS
showList :: [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 =
   Plan -> (Direction, y) -> sig y -> Cache sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan
      (Integer -> Plan
plan (Integer -> Plan) -> Integer -> Plan
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
xs)
      (((Direction, y), (Direction, y)) -> (Direction, y)
forall a b. (a, b) -> a
fst (((Direction, y), (Direction, y)) -> (Direction, y))
-> ((Direction, y), (Direction, y)) -> (Direction, y)
forall a b. (a -> b) -> a -> b
$ sig y -> ((Direction, y), (Direction, y))
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 =
   Plan -> (Direction, y) -> sig y -> Cache sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Plan -> (Direction, y) -> sig y -> Cache sig y
cacheFromPlan
      (Integer -> Plan
plan (Integer -> Plan) -> Integer -> Plan
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
xs)
      (((Direction, y), (Direction, y)) -> (Direction, y)
forall a b. (a, b) -> b
snd (((Direction, y), (Direction, y)) -> (Direction, y))
-> ((Direction, y), (Direction, y)) -> (Direction, y)
forall a b. (a -> b) -> a -> b
$ sig y -> ((Direction, y), (Direction, y))
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 (Integer -> Plan) -> Integer -> Plan
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
xs
       ((Direction, y)
z,(Direction, y)
zInv) = sig y -> ((Direction, y), (Direction, y))
forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
xs
   in  State (CacheMap sig y) (Cache sig y, Cache sig y)
-> CacheMap sig y -> (Cache sig y, Cache sig y)
forall s a. State s a -> s -> a
State.evalState
          ((Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> State (CacheMap sig y) (Cache sig y, Cache sig y)
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)) (CacheMap sig y -> (Cache sig y, Cache sig y))
-> CacheMap sig y -> (Cache sig y, Cache sig y)
forall a b. (a -> b) -> a -> b
$
       CacheMap sig y
forall k a. Map k a
Map.empty


data Direction = Forward | Backward
   deriving (Int -> Direction -> ShowS
[Direction] -> ShowS
Direction -> [Char]
(Int -> Direction -> ShowS)
-> (Direction -> [Char])
-> ([Direction] -> ShowS)
-> Show Direction
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Direction -> ShowS
showsPrec :: Int -> Direction -> ShowS
$cshow :: Direction -> [Char]
show :: Direction -> [Char]
$cshowList :: [Direction] -> ShowS
showList :: [Direction] -> ShowS
Show, Direction -> Direction -> Bool
(Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool) -> Eq Direction
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Direction -> Direction -> Bool
== :: Direction -> Direction -> Bool
$c/= :: Direction -> Direction -> Bool
/= :: Direction -> Direction -> Bool
Eq, Eq Direction
Eq Direction =>
(Direction -> Direction -> Ordering)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Bool)
-> (Direction -> Direction -> Direction)
-> (Direction -> Direction -> Direction)
-> Ord 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
$ccompare :: Direction -> Direction -> Ordering
compare :: Direction -> Direction -> Ordering
$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
>= :: Direction -> Direction -> Bool
$cmax :: Direction -> Direction -> Direction
max :: Direction -> Direction -> Direction
$cmin :: Direction -> Direction -> Direction
min :: Direction -> Direction -> Direction
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 =
   State (CacheMap sig y) (Cache sig y)
-> CacheMap sig y -> Cache sig y
forall s a. State s a -> s -> a
State.evalState (Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
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) (CacheMap sig y -> Cache sig y) -> CacheMap sig y -> Cache sig y
forall a b. (a -> b) -> a -> b
$
   CacheMap sig y
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 -> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a. a -> StateT (CacheMap sig y) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Cache sig y -> State (CacheMap sig y) (Cache sig y))
-> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a b. (a -> b) -> a -> b
$ Cache sig y
forall (sig :: * -> *) y. Cache sig y
CacheIdentity
      PlanSmall LevelSmall
size -> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a. a -> StateT (CacheMap sig y) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Cache sig y -> State (CacheMap sig y) (Cache sig y))
-> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a b. (a -> b) -> a -> b
$ LevelCacheSmall y -> Cache sig y
forall (sig :: * -> *) y. LevelCacheSmall y -> Cache sig y
CacheSmall (LevelCacheSmall y -> Cache sig y)
-> LevelCacheSmall y -> Cache sig y
forall a b. (a -> b) -> a -> b
$
         case LevelSmall
size of
            LevelSmall
Level2 -> y -> LevelCacheSmall y
forall y. y -> LevelCacheSmall y
LevelCache2 (y -> LevelCacheSmall y) -> y -> LevelCacheSmall y
forall a b. (a -> b) -> a -> b
$ y -> y
forall y. C y => y -> y
cache2 y
z
            LevelSmall
Level3 -> (y, y) -> LevelCacheSmall y
forall y. (y, y) -> LevelCacheSmall y
LevelCache3 ((y, y) -> LevelCacheSmall y) -> (y, y) -> LevelCacheSmall y
forall a b. (a -> b) -> a -> b
$ y -> (y, y)
forall y. C y => y -> (y, y)
cache3 y
z
            LevelSmall
Level4 -> (y, y, y) -> LevelCacheSmall y
forall y. (y, y, y) -> LevelCacheSmall y
LevelCache4 ((y, y, y) -> LevelCacheSmall y) -> (y, y, y) -> LevelCacheSmall y
forall a b. (a -> b) -> a -> b
$ y -> (y, y, y)
forall y. C y => y -> (y, y, y)
cache4 y
z
            LevelSmall
Level5 -> (y, y, y, y) -> LevelCacheSmall y
forall y. (y, y, y, y) -> LevelCacheSmall y
LevelCache5 ((y, y, y, y) -> LevelCacheSmall y)
-> (y, y, y, y) -> LevelCacheSmall y
forall a b. (a -> b) -> a -> b
$ y -> (y, y, y, y)
forall y. C y => y -> (y, y, y, y)
cache5 y
z
      Plan
PlanNaive ->
         Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a. a -> StateT (CacheMap sig y) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Cache sig y -> State (CacheMap sig y) (Cache sig y))
-> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a b. (a -> b) -> a -> b
$ LevelCacheNaive y -> Cache sig y
forall (sig :: * -> *) y. LevelCacheNaive y -> Cache sig y
CacheNaive (LevelCacheNaive y -> Cache sig y)
-> LevelCacheNaive y -> Cache sig y
forall a b. (a -> b) -> a -> b
$ y -> LevelCacheNaive y
forall y. y -> LevelCacheNaive y
LevelCacheNaive y
z
      PlanRadix2 level :: LevelRadix2
level@(LevelRadix2 Int
n2) Plan
subPlan ->
         let subxs :: sig y
subxs = Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take Int
n2 sig y
xs
         in  LevelCacheRadix2 sig y -> Cache sig y -> Cache sig y
forall (sig :: * -> *) y.
LevelCacheRadix2 sig y -> Cache sig y -> Cache sig y
CacheRadix2 (LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y
levelCacheRadix2 LevelRadix2
level y
z sig y
subxs) (Cache sig y -> Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
             Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
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
zy -> y -> y
forall 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 ->
            LevelCachePrime sig y -> (Cache sig y, Cache sig y) -> Cache sig y
forall (sig :: * -> *) y.
LevelCachePrime sig y -> (Cache sig y, Cache sig y) -> Cache sig y
CachePrime
               (LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y
levelCachePrime LevelPrime
level ((Cache sig y, Cache sig y) -> Cache sig y
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)
         ((Cache sig y, Cache sig y) -> Cache sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         let subxs :: sig y
subxs = Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take (T -> Int
Permutation.size T
perm) sig y
xs
         in  (Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
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)
                (sig y -> ((Direction, y), (Direction, y))
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 ->
         LevelCoprime -> (Cache sig y, Cache sig y) -> Cache sig y
forall (sig :: * -> *) y.
LevelCoprime -> (Cache sig y, Cache sig y) -> Cache sig y
CacheCoprime LevelCoprime
level ((Cache sig y, Cache sig y) -> Cache sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         (Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
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
zy -> Integer -> y
forall a. C a => a -> Integer -> a
^Integer
m), (Direction
d,y
zy -> Integer -> y
forall a. C a => a -> Integer -> a
^Integer
n))
            (Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take (Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n) sig y
xs, Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take (Integer -> Int
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 ->
         LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> Cache sig y
forall (sig :: * -> *) y.
LevelCacheComposite sig y
-> (Cache sig y, Cache sig y) -> Cache sig y
CacheComposite (LevelComposite -> y -> sig y -> LevelCacheComposite sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
LevelComposite -> y -> sig y -> LevelCacheComposite sig y
levelCacheComposite LevelComposite
level y
z sig y
xs) ((Cache sig y, Cache sig y) -> Cache sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
         (Plan, Plan)
-> ((Direction, y), (Direction, y))
-> (sig y, sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
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
zy -> Integer -> y
forall a. C a => a -> Integer -> a
^Integer
m), (Direction
d,y
zy -> Integer -> y
forall a. C a => a -> Integer -> a
^Integer
n))
            (Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take (Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n) sig y
xs, Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
CutG.take (Integer -> Int
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, (Direction, y) -> Direction
forall a b. (a, b) -> a
fst (Direction, y)
z)
   Maybe (Cache sig y)
item <- (CacheMap sig y -> Maybe (Cache sig y))
-> StateT (CacheMap sig y) Identity (Maybe (Cache sig y))
forall (m :: * -> *) s a. Monad m => (s -> a) -> StateT s m a
State.gets ((Plan, Direction) -> CacheMap sig y -> Maybe (Cache sig y)
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 -> Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a. a -> StateT (CacheMap sig y) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Cache sig y
c
      Maybe (Cache sig y)
Nothing ->
         Plan
-> (Direction, y) -> sig y -> State (CacheMap sig y) (Cache sig y)
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 State (CacheMap sig y) (Cache sig y)
-> (Cache sig y -> State (CacheMap sig y) (Cache sig y))
-> State (CacheMap sig y) (Cache sig y)
forall a b.
StateT (CacheMap sig y) Identity a
-> (a -> StateT (CacheMap sig y) Identity b)
-> StateT (CacheMap sig y) Identity b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Cache sig y
m ->
         (CacheMap sig y -> CacheMap sig y)
-> StateT (CacheMap sig y) Identity ()
forall (m :: * -> *) s. Monad m => (s -> s) -> StateT s m ()
State.modify ((Plan, Direction)
-> Cache sig y -> CacheMap sig y -> CacheMap sig y
forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert (Plan, Direction)
key Cache sig y
m) StateT (CacheMap sig y) Identity ()
-> State (CacheMap sig y) (Cache sig y)
-> State (CacheMap sig y) (Cache sig y)
forall a b.
StateT (CacheMap sig y) Identity a
-> StateT (CacheMap sig y) Identity b
-> StateT (CacheMap sig y) Identity b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>>
         Cache sig y -> State (CacheMap sig y) (Cache sig y)
forall a. a -> StateT (CacheMap sig y) Identity a
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) =
   (Cache sig y -> Cache sig y -> (Cache sig y, Cache sig y))
-> StateT (CacheMap sig y) Identity (Cache sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y)
-> StateT (CacheMap sig y) Identity (Cache sig y, Cache sig y)
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (,)
      (Plan
-> (Direction, y)
-> sig y
-> StateT (CacheMap sig y) Identity (Cache sig y)
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)
      (Plan
-> (Direction, y)
-> sig y
-> StateT (CacheMap sig y) Identity (Cache sig y)
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
[LevelCacheNaive y] -> ShowS
LevelCacheNaive y -> [Char]
(Int -> LevelCacheNaive y -> ShowS)
-> (LevelCacheNaive y -> [Char])
-> ([LevelCacheNaive y] -> ShowS)
-> Show (LevelCacheNaive y)
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
$cshowsPrec :: forall y. Show y => Int -> LevelCacheNaive y -> ShowS
showsPrec :: Int -> LevelCacheNaive y -> ShowS
$cshow :: forall y. Show y => LevelCacheNaive y -> [Char]
show :: LevelCacheNaive y -> [Char]
$cshowList :: forall y. Show y => [LevelCacheNaive y] -> ShowS
showList :: [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 =
   sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
   (y -> y) -> T y -> T y
forall a b. (a -> b) -> T a -> T b
SigS.map
      (T y -> T y -> y
forall a. C a => T a -> T a -> a
scalarProduct1 (sig y -> T y
forall y. Storage (sig y) => sig y -> T y
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
sig) (T y -> y) -> (y -> T y) -> y -> y
forall b c a. (b -> c) -> (a -> b) -> a -> c
. sig y -> y -> T y
forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> y -> T y
powers sig y
sig)
      (sig y -> y -> T y
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 =
   (a -> a -> a) -> T a -> a
forall x. (x -> x -> x) -> T x -> x
SigS.foldL1 a -> a -> a
forall a. C a => a -> a -> a
(+) (T a -> a) -> T a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> T a -> T a -> T a
forall a b c. (a -> b -> c) -> T a -> T b -> T c
SigS.zipWith a -> a -> a
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 =
   sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
   y -> T y -> T y
forall y. C y => y -> T y -> T y
Ana.chirpTransform y
z (T y -> T y) -> T y -> T y
forall a b. (a -> b) -> a -> b
$ sig y -> T y
forall y. Storage (sig y) => sig y -> T y
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 = (y -> y) -> y -> T y
forall a. (a -> a) -> a -> T a
SigS.iterate (y
cy -> y -> y
forall a. C a => a -> a -> a
*) (y -> T y) -> y -> T y
forall a b. (a -> b) -> a -> b
$ sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
multId sig y
sig


data LevelSmall = Level2 | Level3 | Level4 | Level5
   deriving (Int -> LevelSmall -> ShowS
[LevelSmall] -> ShowS
LevelSmall -> [Char]
(Int -> LevelSmall -> ShowS)
-> (LevelSmall -> [Char])
-> ([LevelSmall] -> ShowS)
-> Show LevelSmall
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LevelSmall -> ShowS
showsPrec :: Int -> LevelSmall -> ShowS
$cshow :: LevelSmall -> [Char]
show :: LevelSmall -> [Char]
$cshowList :: [LevelSmall] -> ShowS
showList :: [LevelSmall] -> ShowS
Show, LevelSmall -> LevelSmall -> Bool
(LevelSmall -> LevelSmall -> Bool)
-> (LevelSmall -> LevelSmall -> Bool) -> Eq LevelSmall
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LevelSmall -> LevelSmall -> Bool
== :: LevelSmall -> LevelSmall -> Bool
$c/= :: LevelSmall -> LevelSmall -> Bool
/= :: LevelSmall -> LevelSmall -> Bool
Eq, Eq LevelSmall
Eq LevelSmall =>
(LevelSmall -> LevelSmall -> Ordering)
-> (LevelSmall -> LevelSmall -> Bool)
-> (LevelSmall -> LevelSmall -> Bool)
-> (LevelSmall -> LevelSmall -> Bool)
-> (LevelSmall -> LevelSmall -> Bool)
-> (LevelSmall -> LevelSmall -> LevelSmall)
-> (LevelSmall -> LevelSmall -> LevelSmall)
-> Ord 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
$ccompare :: LevelSmall -> LevelSmall -> Ordering
compare :: LevelSmall -> LevelSmall -> Ordering
$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
>= :: LevelSmall -> LevelSmall -> Bool
$cmax :: LevelSmall -> LevelSmall -> LevelSmall
max :: LevelSmall -> LevelSmall -> LevelSmall
$cmin :: LevelSmall -> LevelSmall -> LevelSmall
min :: LevelSmall -> LevelSmall -> LevelSmall
Ord, Int -> LevelSmall
LevelSmall -> Int
LevelSmall -> [LevelSmall]
LevelSmall -> LevelSmall
LevelSmall -> LevelSmall -> [LevelSmall]
LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall]
(LevelSmall -> LevelSmall)
-> (LevelSmall -> LevelSmall)
-> (Int -> LevelSmall)
-> (LevelSmall -> Int)
-> (LevelSmall -> [LevelSmall])
-> (LevelSmall -> LevelSmall -> [LevelSmall])
-> (LevelSmall -> LevelSmall -> [LevelSmall])
-> (LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall])
-> Enum 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
$csucc :: LevelSmall -> LevelSmall
succ :: LevelSmall -> LevelSmall
$cpred :: LevelSmall -> LevelSmall
pred :: LevelSmall -> LevelSmall
$ctoEnum :: Int -> LevelSmall
toEnum :: Int -> LevelSmall
$cfromEnum :: LevelSmall -> Int
fromEnum :: LevelSmall -> Int
$cenumFrom :: LevelSmall -> [LevelSmall]
enumFrom :: LevelSmall -> [LevelSmall]
$cenumFromThen :: LevelSmall -> LevelSmall -> [LevelSmall]
enumFromThen :: LevelSmall -> LevelSmall -> [LevelSmall]
$cenumFromTo :: LevelSmall -> LevelSmall -> [LevelSmall]
enumFromTo :: LevelSmall -> LevelSmall -> [LevelSmall]
$cenumFromThenTo :: LevelSmall -> LevelSmall -> LevelSmall -> [LevelSmall]
enumFromThenTo :: LevelSmall -> LevelSmall -> 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
[LevelCacheSmall y] -> ShowS
LevelCacheSmall y -> [Char]
(Int -> LevelCacheSmall y -> ShowS)
-> (LevelCacheSmall y -> [Char])
-> ([LevelCacheSmall y] -> ShowS)
-> Show (LevelCacheSmall y)
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
$cshowsPrec :: forall y. Show y => Int -> LevelCacheSmall y -> ShowS
showsPrec :: Int -> LevelCacheSmall y -> ShowS
$cshow :: forall y. Show y => LevelCacheSmall y -> [Char]
show :: LevelCacheSmall y -> [Char]
$cshowList :: forall y. Show y => [LevelCacheSmall y] -> ShowS
showList :: [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
zy -> y -> y
forall 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
zy -> y -> y
forall a. C a => a -> a -> a
*y
z in (y
z,y
z2,y
zy -> y -> y
forall 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
zy -> y -> y
forall a. C a => a -> a -> a
*y
z in (y
z,y
z2,y
zy -> y -> y
forall a. C a => a -> a -> a
*y
z2,y
z2y -> y -> y
forall 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]
_ = sig y -> [y]
forall y. Storage (sig y) => sig y -> [y]
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
   in  sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       [y] -> T y
forall y. [y] -> T y
SigS.fromList [y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
x1, y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
zy -> y -> y
forall 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]
_ = sig y -> [y]
forall y. Storage (sig y) => sig y -> [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)) = (y, y) -> (y, y) -> ((y, y), (y, y))
forall y. C y => Pair y -> Pair y -> (Pair y, Pair y)
Cyclic.sumAndConvolvePair (y
x1,y
x2) (y
z,y
z2)
   in  sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       [y] -> T y
forall y. [y] -> T y
SigS.fromList [y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
s, y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
zx1, y
x0y -> y -> y
forall 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]
_ = sig y -> [y]
forall y. Storage (sig y) => sig y -> [y]
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> [y]
SigG.toList sig y
sig
       x02a :: y
x02a = y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
x2; x02b :: y
x02b = y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
z2y -> y -> y
forall a. C a => a -> a -> a
*y
x2
       x13a :: y
x13a = y
x1y -> y -> y
forall a. C a => a -> a -> a
+y
x3; x13b :: y
x13b = y
x1y -> y -> y
forall a. C a => a -> a -> a
+y
z2y -> y -> y
forall a. C a => a -> a -> a
*y
x3
   in  sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       [y] -> T y
forall y. [y] -> T y
SigS.fromList [y
x02ay -> y -> y
forall a. C a => a -> a -> a
+   y
x13a, y
x02by -> y -> y
forall a. C a => a -> a -> a
+y
z y -> y -> y
forall a. C a => a -> a -> a
*y
x13b,
                      y
x02ay -> y -> y
forall a. C a => a -> a -> a
+y
z2y -> y -> y
forall a. C a => a -> a -> a
*y
x13a, y
x02by -> y -> y
forall a. C a => a -> a -> a
+y
z3y -> y -> y
forall 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]
_ = sig y -> [y]
forall y. Storage (sig y) => sig y -> [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)) =
          (y, y, y, y) -> (y, y, y, y) -> ((y, y), (y, y, y, y))
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  sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       [y] -> T y
forall y. [y] -> T y
SigS.fromList [y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
s, y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
d1, y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
d2, y
x0y -> y -> y
forall a. C a => a -> a -> a
+y
d3, y
x0y -> y -> y
forall 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]
(Int -> LevelRadix2 -> ShowS)
-> (LevelRadix2 -> [Char])
-> ([LevelRadix2] -> ShowS)
-> Show LevelRadix2
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LevelRadix2 -> ShowS
showsPrec :: Int -> LevelRadix2 -> ShowS
$cshow :: LevelRadix2 -> [Char]
show :: LevelRadix2 -> [Char]
$cshowList :: [LevelRadix2] -> ShowS
showList :: [LevelRadix2] -> ShowS
Show, LevelRadix2 -> LevelRadix2 -> Bool
(LevelRadix2 -> LevelRadix2 -> Bool)
-> (LevelRadix2 -> LevelRadix2 -> Bool) -> Eq LevelRadix2
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LevelRadix2 -> LevelRadix2 -> Bool
== :: LevelRadix2 -> LevelRadix2 -> Bool
$c/= :: LevelRadix2 -> LevelRadix2 -> Bool
/= :: LevelRadix2 -> LevelRadix2 -> Bool
Eq, Eq LevelRadix2
Eq LevelRadix2 =>
(LevelRadix2 -> LevelRadix2 -> Ordering)
-> (LevelRadix2 -> LevelRadix2 -> Bool)
-> (LevelRadix2 -> LevelRadix2 -> Bool)
-> (LevelRadix2 -> LevelRadix2 -> Bool)
-> (LevelRadix2 -> LevelRadix2 -> Bool)
-> (LevelRadix2 -> LevelRadix2 -> LevelRadix2)
-> (LevelRadix2 -> LevelRadix2 -> LevelRadix2)
-> Ord 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
$ccompare :: LevelRadix2 -> LevelRadix2 -> Ordering
compare :: LevelRadix2 -> LevelRadix2 -> Ordering
$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
>= :: LevelRadix2 -> LevelRadix2 -> Bool
$cmax :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
max :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
$cmin :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
min :: LevelRadix2 -> LevelRadix2 -> LevelRadix2
Ord)

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


data LevelCacheRadix2 sig y =
   LevelCacheRadix2 Int (sig y)
   deriving (Int -> LevelCacheRadix2 sig y -> ShowS
[LevelCacheRadix2 sig y] -> ShowS
LevelCacheRadix2 sig y -> [Char]
(Int -> LevelCacheRadix2 sig y -> ShowS)
-> (LevelCacheRadix2 sig y -> [Char])
-> ([LevelCacheRadix2 sig y] -> ShowS)
-> Show (LevelCacheRadix2 sig y)
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]
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheRadix2 sig y -> ShowS
showsPrec :: Int -> LevelCacheRadix2 sig y -> ShowS
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheRadix2 sig y -> [Char]
show :: LevelCacheRadix2 sig y -> [Char]
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheRadix2 sig y] -> ShowS
showList :: [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 =
   Int -> sig y -> LevelCacheRadix2 sig y
forall (sig :: * -> *) y. Int -> sig y -> LevelCacheRadix2 sig y
LevelCacheRadix2 Int
n2
      (sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$ sig y -> y -> T y
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) = Int -> sig y -> (sig y, sig y)
forall sig. Transform sig => Int -> sig -> (sig, sig)
SigG.splitAt Int
n2 sig y
sig
       fs0 :: sig y
fs0 = Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$ (y -> y -> y) -> sig y -> sig y -> sig y
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 y -> y -> y
forall a. C a => a -> a -> a
(+) sig y
xs0 sig y
xs1
       fs1 :: sig y
fs1 = Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
                (y -> y -> y -> y) -> sig y -> sig y -> sig y -> sig y
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
wy -> y -> y
forall a. C a => a -> a -> a
*(y
x0y -> y -> y
forall a. C a => a -> a -> a
-y
x1))
                   sig y
twiddle sig y
xs0 sig y
xs1
   in  sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       T y -> T y -> T y
forall y. T y -> T y -> T y
SigS.interleave (sig y -> T y
forall y. Storage (sig y) => sig y -> T y
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState sig y
fs0) (sig y -> T y
forall y. Storage (sig y) => sig y -> T y
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]
(Int -> LevelComposite -> ShowS)
-> (LevelComposite -> [Char])
-> ([LevelComposite] -> ShowS)
-> Show LevelComposite
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LevelComposite -> ShowS
showsPrec :: Int -> LevelComposite -> ShowS
$cshow :: LevelComposite -> [Char]
show :: LevelComposite -> [Char]
$cshowList :: [LevelComposite] -> ShowS
showList :: [LevelComposite] -> ShowS
Show)

instance Eq LevelComposite where
   LevelComposite
a == :: LevelComposite -> LevelComposite -> Bool
== LevelComposite
b  =  LevelComposite -> LevelComposite -> Ordering
forall a. Ord a => a -> a -> Ordering
compare LevelComposite
a LevelComposite
b Ordering -> Ordering -> Bool
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)
_)  =
      (Integer, Integer) -> (Integer, Integer) -> Ordering
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 = Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n
       mi :: Int
mi = Integer -> Int
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
[LevelCacheComposite sig y] -> ShowS
LevelCacheComposite sig y -> [Char]
(Int -> LevelCacheComposite sig y -> ShowS)
-> (LevelCacheComposite sig y -> [Char])
-> ([LevelCacheComposite sig y] -> ShowS)
-> Show (LevelCacheComposite sig y)
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]
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCacheComposite sig y -> ShowS
showsPrec :: Int -> LevelCacheComposite sig y -> ShowS
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCacheComposite sig y -> [Char]
show :: LevelCacheComposite sig y -> [Char]
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCacheComposite sig y] -> ShowS
showList :: [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 =
   (Integer, Integer) -> (T, T) -> sig y -> LevelCacheComposite sig y
forall (sig :: * -> *) y.
(Integer, Integer) -> (T, T) -> sig y -> LevelCacheComposite sig y
LevelCacheComposite (Integer
n,Integer
m) (T, T)
transpose (sig y -> LevelCacheComposite sig y)
-> sig y -> LevelCacheComposite sig y
forall a b. (a -> b) -> a -> b
$
   sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
   (((Integer, y, y) -> (y, (Integer, y, y)))
 -> (Integer, y, y) -> T y)
-> (Integer, y, y)
-> ((Integer, y, y) -> (y, (Integer, y, y)))
-> T y
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Integer, y, y) -> (y, (Integer, y, y))) -> (Integer, y, y) -> T y
forall acc y. (acc -> (y, acc)) -> acc -> T y
SigS.generateInfinite (Integer
n, sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
multId sig y
sig, sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
multId sig y
sig) (((Integer, y, y) -> (y, (Integer, y, y))) -> T y)
-> ((Integer, y, y) -> (y, (Integer, y, y))) -> T y
forall a b. (a -> b) -> a -> b
$ \(Integer
i,y
zi,y
zij) ->
   (y
zij,
    case Integer -> Integer
forall a. Enum a => a -> a
pred Integer
i of
      Integer
0 -> (Integer
n, y
ziy -> y -> y
forall a. C a => a -> a -> a
*y
z, sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
multId sig y
sig)
      Integer
i1 -> (Integer
i1, y
zi, y
zijy -> y -> y
forall 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 =
   T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeMN (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       sig y -> T (sig y) -> sig y
forall (sig :: * -> *) y.
Transform sig y =>
sig y -> T (sig y) -> sig y
concatRechunk sig y
sig (T (sig y) -> sig y) -> (sig y -> T (sig y)) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       (sig y -> sig y) -> T (sig y) -> T (sig y)
forall a b. (a -> b) -> T a -> T b
SigS.map (Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCacheM) (T (sig y) -> T (sig y))
-> (sig y -> T (sig y)) -> sig y -> T (sig y)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Int -> sig y -> T (sig y)
forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
m) (sig y -> T (sig y)) -> (sig y -> sig y) -> sig y -> T (sig y)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeNM (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
--       concatRechunk sig .
       (y -> y -> y) -> sig y -> sig y -> sig y
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 y -> y -> y
forall a. C a => a -> a -> a
(*) sig y
twiddle (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       T (sig y) -> sig y
forall m. Monoid m => T m -> m
SigS.fold (T (sig y) -> sig y) -> (sig y -> T (sig y)) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       (sig y -> sig y) -> T (sig y) -> T (sig y)
forall a b. (a -> b) -> T a -> T b
SigS.map (Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCacheN) (T (sig y) -> T (sig y))
-> (sig y -> T (sig y)) -> sig y -> T (sig y)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Int -> sig y -> T (sig y)
forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n) (sig y -> T (sig y)) -> (sig y -> sig y) -> sig y -> T (sig y)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transposeMN (sig y -> sig y) -> sig y -> sig y
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]
(Int -> LevelCoprime -> ShowS)
-> (LevelCoprime -> [Char])
-> ([LevelCoprime] -> ShowS)
-> Show LevelCoprime
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LevelCoprime -> ShowS
showsPrec :: Int -> LevelCoprime -> ShowS
$cshow :: LevelCoprime -> [Char]
show :: LevelCoprime -> [Char]
$cshowList :: [LevelCoprime] -> ShowS
showList :: [LevelCoprime] -> ShowS
Show)

instance Eq LevelCoprime where
   LevelCoprime
a == :: LevelCoprime -> LevelCoprime -> Bool
== LevelCoprime
b  =  LevelCoprime -> LevelCoprime -> Ordering
forall a. Ord a => a -> a -> Ordering
compare LevelCoprime
a LevelCoprime
b Ordering -> Ordering -> Bool
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)
_)  =
      (Integer, Integer) -> (Integer, Integer) -> Ordering
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 = Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n
       mi :: Int
mi = Integer -> Int
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 =
          sig y -> T (sig y) -> sig y
forall (sig :: * -> *) y.
Transform sig y =>
sig y -> T (sig y) -> sig y
concatRechunk sig y
sig (T (sig y) -> sig y) -> (sig y -> T (sig y)) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          (sig y -> sig y) -> T (sig y) -> T (sig y)
forall a b. (a -> b) -> T a -> T b
SigS.map (Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache) (T (sig y) -> T (sig y))
-> (sig y -> T (sig y)) -> sig y -> T (sig y)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          Int -> sig y -> T (sig y)
forall sig. Transform sig => Int -> sig -> T sig
SigG.sliceVertical (a -> Int
forall a b. (C a, C b) => a -> b
fromIntegral a
j) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$ sig y
sig
   in  T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
gridInv (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Cache sig y -> Integer -> sig y -> sig y
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 (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
transpose (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       Cache sig y -> Integer -> sig y -> sig y
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 (sig y -> sig y) -> (sig y -> sig y) -> sig y -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
       T -> sig y -> sig y
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 =
   sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
pattern (T y -> sig y) -> (T (sig y) -> T y) -> T (sig y) -> sig y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   sig y -> T y
forall y. Storage (sig y) => sig y -> T y
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> T y
SigG.toState (sig y -> T y) -> (T (sig y) -> sig y) -> T (sig y) -> T y
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   T (sig y) -> sig y
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]
(Int -> LevelPrime -> ShowS)
-> (LevelPrime -> [Char])
-> ([LevelPrime] -> ShowS)
-> Show LevelPrime
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LevelPrime -> ShowS
showsPrec :: Int -> LevelPrime -> ShowS
$cshow :: LevelPrime -> [Char]
show :: LevelPrime -> [Char]
$cshowList :: [LevelPrime] -> ShowS
showList :: [LevelPrime] -> ShowS
Show)

instance Eq LevelPrime where
   LevelPrime
a == :: LevelPrime -> LevelPrime -> Bool
== LevelPrime
b  =  LevelPrime -> LevelPrime -> Ordering
forall a. Ord a => a -> a -> Ordering
compare LevelPrime
a LevelPrime
b Ordering -> Ordering -> Bool
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
_))  =
      Int -> Int -> Ordering
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 (Int -> T) -> Int -> T
forall a b. (a -> b) -> a -> b
$ Integer -> Int
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
[LevelCachePrime sig y] -> ShowS
LevelCachePrime sig y -> [Char]
(Int -> LevelCachePrime sig y -> ShowS)
-> (LevelCachePrime sig y -> [Char])
-> ([LevelCachePrime sig y] -> ShowS)
-> Show (LevelCachePrime sig y)
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]
$cshowsPrec :: forall (sig :: * -> *) y.
Show (sig y) =>
Int -> LevelCachePrime sig y -> ShowS
showsPrec :: Int -> LevelCachePrime sig y -> ShowS
$cshow :: forall (sig :: * -> *) y.
Show (sig y) =>
LevelCachePrime sig y -> [Char]
show :: LevelCachePrime sig y -> [Char]
$cshowList :: forall (sig :: * -> *) y.
Show (sig y) =>
[LevelCachePrime sig y] -> ShowS
showList :: [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 =
   (T, T) -> sig y -> LevelCachePrime sig y
forall (sig :: * -> *) y. (T, T) -> sig y -> LevelCachePrime sig y
LevelCachePrime (T
rev, T
inv)
      ((\sig y
zs -> y -> sig y -> sig y
forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
recipInteger sig y
zs) sig y
zs) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
       Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
subCache (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
       T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
       sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
sig (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
       (y -> y) -> y -> T y
forall a. (a -> a) -> a -> T a
SigS.iterate (y
zy -> y -> y
forall 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 =
   sig y -> (y -> sig y -> sig y) -> sig y -> sig y
forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL ([Char] -> sig y
forall a. HasCallStack => [Char] -> a
error [Char]
"transformPrime: empty signal") ((y -> sig y -> sig y) -> sig y -> sig y)
-> (y -> sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
   \y
x0 sig y
rest ->
      y -> sig y -> sig y
forall y. Storage (sig y) => y -> sig y -> sig y
forall (sig :: * -> *) y.
(Transform0 sig, Storage (sig y)) =>
y -> sig y -> sig y
SigG.cons ((y -> y -> y) -> y -> sig y -> y
forall y s. Storage (sig y) => (s -> y -> s) -> s -> sig y -> s
forall (sig :: * -> *) y s.
(Read0 sig, Storage (sig y)) =>
(s -> y -> s) -> s -> sig y -> s
SigG.foldL y -> y -> y
forall a. C a => a -> a -> a
(+) y
x0 sig y
rest) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      (y -> y) -> sig y -> sig y
forall y0 y1.
(Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
forall (sig :: * -> *) y0 y1.
(Transform0 sig, Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
SigG.map (y
x0y -> y -> y
forall a. C a => a -> a -> a
+) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
inv (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
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 (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      T -> sig y -> sig y
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 =
   sig y -> (y -> sig y -> sig y) -> sig y -> sig y
forall (sig :: * -> *) y a.
Transform sig y =>
a -> (y -> sig y -> a) -> sig y -> a
SigG.switchL ([Char] -> sig y
forall a. HasCallStack => [Char] -> a
error [Char]
"transformPrime: empty signal") ((y -> sig y -> sig y) -> sig y -> sig y)
-> (y -> sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
   \y
x0 sig y
rest ->
      y -> sig y -> sig y
forall y. Storage (sig y) => y -> sig y -> sig y
forall (sig :: * -> *) y.
(Transform0 sig, Storage (sig y)) =>
y -> sig y -> sig y
SigG.cons ((y -> y -> y) -> y -> sig y -> y
forall y s. Storage (sig y) => (s -> y -> s) -> s -> sig y -> s
forall (sig :: * -> *) y s.
(Read0 sig, Storage (sig y)) =>
(s -> y -> s) -> s -> sig y -> s
SigG.foldL y -> y -> y
forall a. C a => a -> a -> a
(+) y
x0 sig y
rest) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      (y -> y) -> sig y -> sig y
forall y0 y1.
(Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
forall (sig :: * -> *) y0 y1.
(Transform0 sig, Storage (sig y0), Storage (sig y1)) =>
(y0 -> y1) -> sig y0 -> sig y1
SigG.map (y
x0y -> y -> y
forall a. C a => a -> a -> a
+) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
inv (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
      sig y -> sig y -> sig y
forall (sig :: * -> *) y.
(Transform sig y, C y) =>
sig y -> sig y -> sig y
Cyclic.filterNaive
         (T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm sig y
rest)
         (T -> sig y -> sig y
forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
Permutation.apply T
perm (sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
rest ((y -> y) -> y -> T y
forall a. (a -> a) -> a -> T a
SigS.iterate (y
zy -> y -> y
forall 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
[Window sig y] -> ShowS
Window sig y -> [Char]
(Int -> Window sig y -> ShowS)
-> (Window sig y -> [Char])
-> ([Window sig y] -> ShowS)
-> Show (Window sig y)
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]
$cshowsPrec :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Int -> Window sig y -> ShowS
showsPrec :: Int -> Window sig y -> ShowS
$cshow :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
Window sig y -> [Char]
show :: Window sig y -> [Char]
$cshowList :: forall (sig :: * -> *) y.
(Show y, Show (sig y)) =>
[Window sig y] -> ShowS
showList :: [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 sig y -> Bool
forall sig. Read sig => sig -> Bool
CutG.null sig y
x
     then Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
forall (sig :: * -> *) y.
Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
Window Int
0 (Cache sig y
forall (sig :: * -> *) y. Cache sig y
CacheIdentity, Cache sig y
forall (sig :: * -> *) y. Cache sig y
CacheIdentity) sig y
forall sig. Monoid sig => sig
CutG.empty
     else
       let size :: Int
size  = sig y -> Int
forall sig. Read sig => sig -> Int
CutG.length sig y
x
           size2 :: Int
size2 = Int
2 Int -> Int -> Int
forall a. C a => a -> a -> a
* Int -> Int
forall a. (C a, Bits a) => a -> a
NumberTheory.ceilingPowerOfTwo Int
size
           padded :: sig y
padded =
              Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
SigG.take Int
size2 (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
              sig y -> sig y -> sig y
forall sig. Monoid sig => sig -> sig -> sig
CutG.append sig y
x (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
                 let pad :: sig y
pad = sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
x (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$ y -> T y
forall a. a -> T a
SigS.repeat (y -> T y) -> y -> T y
forall a b. (a -> b) -> a -> b
$ sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
addId sig y
x
                 in  sig y -> sig y -> sig y
forall sig. Monoid sig => sig -> sig -> sig
CutG.append sig y
pad (sig y -> sig y -> sig y
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) =
              sig y -> (Cache sig y, Cache sig y)
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
sig y -> (Cache sig y, Cache sig y)
cacheDuplex sig y
padded
       in  Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
forall (sig :: * -> *) y.
Int -> (Cache sig y, Cache sig y) -> sig y -> Window sig y
Window
              (Int
size2Int -> Int -> Int
forall a. C a => a -> a -> a
-Int
sizeInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1)
              (Cache sig y, Cache sig y)
caches
              (Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cache (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
               y -> sig y -> sig y
forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). 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
blockSizeInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
forall a. C a => a
zero
     then sig y
forall sig. Monoid sig => sig
CutG.empty
     else
       let windowSize :: Int
windowSize = sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
spectrum Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
blockSize
       in  (sig y -> sig y -> sig y) -> sig y -> T (sig y) -> sig y
forall x acc. (x -> acc -> acc) -> acc -> T x -> acc
SigS.foldR (Int -> sig y -> sig y -> sig y
forall a (sig :: * -> *).
(C a, Transform sig a) =>
Int -> sig a -> sig a -> sig a
FiltNRG.addShiftedSimple Int
blockSize) sig y
forall sig. Monoid sig => sig
CutG.empty (T (sig y) -> sig y) -> T (sig y) -> sig y
forall a b. (a -> b) -> a -> b
$
           (sig y -> sig y) -> T (sig y) -> T (sig y)
forall a b. (a -> b) -> T a -> T b
SigS.map
              (\sig y
block ->
                 Int -> sig y -> sig y
forall sig. Transform sig => Int -> sig -> sig
SigG.take (Int
windowSize Int -> Int -> Int
forall a. C a => a -> a -> a
+ sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
block) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
                 (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
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 (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
                 (sig y -> sig y -> sig y) -> sig y -> sig y -> sig y
forall a b c. (a -> b -> c) -> b -> a -> c
flip sig y -> sig y -> sig y
forall sig. Monoid sig => sig -> sig -> sig
CutG.append
                    {-
                    The last block may be shorter than blockSize
                    and thus needs more padding.
                    -}
                    (sig y -> T y -> sig y
forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
spectrum (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$ y -> T y
forall a. a -> T a
SigS.repeat (y -> T y) -> y -> T y
forall a b. (a -> b) -> a -> b
$ sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
addId sig y
b) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
                 sig y
block) (T (sig y) -> T (sig y)) -> T (sig y) -> T (sig y)
forall a b. (a -> b) -> a -> b
$
           Int -> sig y -> T (sig y)
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 = Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ sig y -> Int
forall sig. Read sig => sig -> Int
SigG.length sig y
x
       ((Direction, y)
z,(Direction, y)
zInv) =
          sig y -> ((Direction, y), (Direction, y))
forall y (sig :: * -> *).
(Element y, Read sig y) =>
sig y -> ((Direction, y), (Direction, y))
directionPrimitiveRootsOfUnity sig y
x
   in  (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
(Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
convolveCyclicCache
          (Plan -> (Direction, y) -> sig y -> Cache sig y
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,
           Plan -> (Direction, y) -> sig y -> Cache sig y
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 =
   (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y
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 -> sig y -> sig y) -> sig y -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
   y -> sig y -> sig y
forall y (sig :: * -> *).
(C y, Transform sig y) =>
y -> sig y -> sig y
FiltNRG.amplify (sig y -> y
forall y (sig :: * -> *). (Element y, Read sig y) => sig y -> y
forall (sig :: * -> *). Read sig y => sig y -> y
recipInteger sig y
x) (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$ Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache ((Cache sig y, Cache sig y) -> Cache sig y
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 =
   Cache sig y -> sig y -> sig y
forall y (sig :: * -> *).
(Element y, Transform sig y) =>
Cache sig y -> sig y -> sig y
transformWithCache Cache sig y
cacheInv (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
   (y -> y -> y) -> sig y -> sig y -> sig y
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 y -> y -> y
forall a. C a => a -> a -> a
(*) sig y
x (sig y -> sig y) -> sig y -> sig y
forall a b. (a -> b) -> a -> b
$
   Cache sig y -> sig y -> sig y
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)
-}