{- |
Permutations of signals as needed for Fast Fourier transforms.
Most functions are independent of the Signal framework.
We could move them as well to Synthesizer.Basic.
-}
module Synthesizer.Generic.Permutation (
   T,
   apply,
   size,
   transposition,
   skewGrid,
   skewGridInv,
   skewGridCRT,
   skewGridCRTInv,
   multiplicative,
   inverse,
   reverse,
   ) where

import qualified Synthesizer.Basic.NumberTheory as NumberTheory

import qualified Synthesizer.Generic.Signal as SigG
import qualified Synthesizer.State.Signal as SigS

import qualified Data.StorableVector.ST.Strict as SVST
import qualified Data.StorableVector as SV

import qualified Algebra.PrincipalIdealDomain as PID

import Prelude hiding (reverse, )



type T = SV.Vector Int

apply ::
   (SigG.Transform sig y) =>
   T -> sig y -> sig y
apply :: forall (sig :: * -> *) y. Transform sig y => T -> sig y -> sig y
apply T
p sig y
xs =
   forall (sig :: * -> *) a b.
(Transform sig a, Transform sig b) =>
sig a -> T b -> sig b
SigG.takeStateMatch sig y
xs forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> T a -> T b
SigS.map (forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> Int -> y
SigG.index sig y
xs) forall a b. (a -> b) -> a -> b
$
   forall a. Storable a => Vector a -> T a
SigS.fromStrictStorableSignal T
p


size :: T -> Int
size :: T -> Int
size = forall a. Vector a -> Int
SV.length


{- |
> inverse (transposition n m) = transposition m n
-}
transposition ::
   Int -> Int -> T
transposition :: Int -> Int -> T
transposition Int
n Int
m =
   forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN (Int
nforall a. Num a => a -> a -> a
*Int
m)
      (\(Int
i,Int
j,Int
k0) -> forall a. a -> Maybe a
Just (Int
i,
         case forall a. Enum a => a -> a
pred Int
k0 of
            Int
0  -> let j1 :: Int
j1 = Int
jforall a. Num a => a -> a -> a
+Int
1 in (Int
j1, Int
j1, Int
m)
            Int
k1 -> (Int
iforall a. Num a => a -> a -> a
+Int
n, Int
j, Int
k1)))
      (Int
0,Int
0,Int
m)


{-
In general the inverse of a skewGrid
does not look like even a generalized skewGrid.
E.g. @inverse $ skewGrid 3 4@.
-}
skewGrid ::
   Int -> Int -> T
skewGrid :: Int -> Int -> T
skewGrid Int
n Int
m =
   let len :: Int
len = Int
nforall a. Num a => a -> a -> a
*Int
m
   in  forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN Int
len
          (\(Int
i0,Int
k0) -> forall a. a -> Maybe a
Just (Int
i0,
             let k1 :: Int
k1 = forall a. Enum a => a -> a
pred Int
k0
                 i1 :: Int
i1 = Int
i0forall a. Num a => a -> a -> a
+Int
n
             in  if Int
k1forall a. Eq a => a -> a -> Bool
==Int
0
                   then (forall a. Integral a => a -> a -> a
mod (Int
i1forall a. Num a => a -> a -> a
+Int
m) Int
len, Int
m)
                   else (forall a. Integral a => a -> a -> a
mod Int
i1 Int
len, Int
k1)))
          (Int
0,Int
m)

{- |
> inverse (skewGrid n m) == skewGridInv n m

In general the inverse of a skewGrid
cannot be expressed like skewGrid or skewGridCRT.
E.g. @inverse $ skewGrid 3 4@.
-}
skewGridInv ::
   Int -> Int -> T
skewGridInv :: Int -> Int -> T
skewGridInv Int
n Int
m =
   forall a. Storable a => [a] -> Vector a
SV.pack forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> [a] -> [b]
map (\Int
k ->
      let Just (Int
i,Int
j) = forall a. C a => a -> a -> a -> Maybe (a, a)
PID.diophantine Int
k Int
n Int
m
      in  forall a. Integral a => a -> a -> a
mod Int
i Int
m forall a. Num a => a -> a -> a
+ forall a. Integral a => a -> a -> a
mod Int
j Int
n forall a. Num a => a -> a -> a
* Int
m) forall a b. (a -> b) -> a -> b
$
   forall a. Int -> [a] -> [a]
take (Int
nforall a. Num a => a -> a -> a
*Int
m) forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Int
1forall a. Num a => a -> a -> a
+) Int
0

skewGridCRT ::
   Int -> Int -> T
skewGridCRT :: Int -> Int -> T
skewGridCRT Int
n Int
m =
   let len :: Int
len = Int
nforall a. Num a => a -> a -> a
*Int
m
       (Int
ni,Int
mi) = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall a. C a => a -> a -> (a, (a, a))
PID.extendedGCD Int
n Int
m
   in  forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN Int
len
          (\(Int
i0,Int
k0) -> forall a. a -> Maybe a
Just (Int
i0,
             let k1 :: Int
k1 = forall a. Enum a => a -> a
pred Int
k0
                 i1 :: Int
i1 = Int
i0forall a. Num a => a -> a -> a
+Int
niforall a. Num a => a -> a -> a
*Int
n
             in  if Int
k1forall a. Eq a => a -> a -> Bool
==Int
0
                   then (forall a. Integral a => a -> a -> a
mod (Int
i1forall a. Num a => a -> a -> a
+Int
miforall a. Num a => a -> a -> a
*Int
m) Int
len, Int
m)
                   else (forall a. Integral a => a -> a -> a
mod Int
i1 Int
len, Int
k1)))
          (Int
0,Int
m)

skewGridCRTInv ::
   Int -> Int -> T
skewGridCRTInv :: Int -> Int -> T
skewGridCRTInv Int
n Int
m =
   forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Int -> [a] -> (Vector a, [a])
SV.packN (Int
nforall a. Num a => a -> a -> a
*Int
m) forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> forall a. Integral a => a -> a -> a
mod Int
k Int
m forall a. Num a => a -> a -> a
+ forall a. Integral a => a -> a -> a
mod Int
k Int
n forall a. Num a => a -> a -> a
* Int
m) forall a b. (a -> b) -> a -> b
$
   forall a. (a -> a) -> a -> [a]
iterate (Int
1forall a. Num a => a -> a -> a
+) Int
0


{- |
Beware of 0-based indices stored in the result vector.
-}
multiplicative :: Int -> T
multiplicative :: Int -> T
multiplicative Int
ni =
   let n :: Integer
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ni
       gen :: Integer
gen = Integer -> Integer
NumberTheory.multiplicativeGenerator Integer
n
   in  {-
       Since 'gen' is usually 2 or 3,
       the error should occur really only for huge signals.
       -}
       if Integer
gen forall a. Num a => a -> a -> a
* Integer
n forall a. Ord a => a -> a -> Bool
> forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a. Bounded a => a
maxBound :: Int)
         then forall a. HasCallStack => [Char] -> a
error [Char]
"signal too long for Int indexing"
         else forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN (Int
niforall a. Num a => a -> a -> a
-Int
1)
                 (\Int
x -> forall a. a -> Maybe a
Just (Int
xforall a. Num a => a -> a -> a
-Int
1, forall a. Integral a => a -> a -> a
mod (forall a. Num a => Integer -> a
fromInteger Integer
gen forall a. Num a => a -> a -> a
* Int
x) Int
ni)) Int
1

{- |
We only need to compute the inverse permutation explicitly,
because not all signal structures support write to arbitrary indices,
thus Generic.Write does not support it.
For strict StorableVector it would be more efficient
to build the vector directly.

It holds:

> inverse . inverse == id
-}
inverse :: T -> T
inverse :: T -> T
inverse T
perm =
   forall e. Storable e => (forall s. ST s (Vector s e)) -> Vector e
SVST.runSTVector
      (do Vector s Int
inv <- forall e s. Storable e => Int -> ST s (Vector s e)
SVST.new_ (forall a. Vector a -> Int
SV.length T
perm)
          forall (m :: * -> *) a. Monad m => T (m a) -> m ()
SigS.sequence_ forall a b. (a -> b) -> a -> b
$
             forall a b c. (a -> b -> c) -> T a -> T b -> T c
SigS.zipWith (forall e s. Storable e => Vector s e -> Int -> e -> ST s ()
SVST.write Vector s Int
inv)
                (forall a. Storable a => Vector a -> T a
SigS.fromStrictStorableSignal T
perm)
                (forall a. (a -> a) -> a -> T a
SigS.iterate (Int
1forall a. Num a => a -> a -> a
+) Int
0)
          forall (m :: * -> *) a. Monad m => a -> m a
return Vector s Int
inv)

reverse :: T -> T
reverse :: T -> T
reverse T
perm =
   forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN (forall a. Vector a -> Int
SV.length T
perm)
      (\Maybe Int
mn -> forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$
         case Maybe Int
mn of
            Maybe Int
Nothing -> (forall a. Storable a => Vector a -> a
SV.head T
perm, forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall a. Vector a -> Int
SV.length T
perm)
            Just Int
n ->
               let n1 :: Int
n1 = Int
nforall a. Num a => a -> a -> a
-Int
1
               in  (forall a. Storable a => Vector a -> Int -> a
SV.index T
perm Int
n1, forall a. a -> Maybe a
Just Int
n1))
      forall a. Maybe a
Nothing