{- |
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 =
   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
xs (T y -> sig y) -> T y -> sig y
forall a b. (a -> b) -> a -> b
$
   (Int -> y) -> T Int -> T y
forall a b. (a -> b) -> T a -> T b
SigS.map (sig y -> Int -> y
forall y. Storage (sig y) => sig y -> Int -> y
forall (sig :: * -> *) y.
(Read0 sig, Storage (sig y)) =>
sig y -> Int -> y
SigG.index sig y
xs) (T Int -> T y) -> T Int -> T y
forall a b. (a -> b) -> a -> b
$
   T -> T Int
forall a. Storable a => Vector a -> T a
SigS.fromStrictStorableSignal T
p


size :: T -> Int
size :: T -> Int
size = T -> Int
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 =
   (T, Maybe (Int, Int, Int)) -> T
forall a b. (a, b) -> a
fst ((T, Maybe (Int, Int, Int)) -> T)
-> (T, Maybe (Int, Int, Int)) -> T
forall a b. (a -> b) -> a -> b
$ Int
-> ((Int, Int, Int) -> Maybe (Int, (Int, Int, Int)))
-> (Int, Int, Int)
-> (T, Maybe (Int, Int, Int))
forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m)
      (\(Int
i,Int
j,Int
k0) -> (Int, (Int, Int, Int)) -> Maybe (Int, (Int, Int, Int))
forall a. a -> Maybe a
Just (Int
i,
         case Int -> Int
forall a. Enum a => a -> a
pred Int
k0 of
            Int
0  -> let j1 :: Int
j1 = Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 in (Int
j1, Int
j1, Int
m)
            Int
k1 -> (Int
iInt -> Int -> Int
forall 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
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m
   in  (T, Maybe (Int, Int)) -> T
forall a b. (a, b) -> a
fst ((T, Maybe (Int, Int)) -> T) -> (T, Maybe (Int, Int)) -> T
forall a b. (a -> b) -> a -> b
$ Int
-> ((Int, Int) -> Maybe (Int, (Int, Int)))
-> (Int, Int)
-> (T, Maybe (Int, Int))
forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN Int
len
          (\(Int
i0,Int
k0) -> (Int, (Int, Int)) -> Maybe (Int, (Int, Int))
forall a. a -> Maybe a
Just (Int
i0,
             let k1 :: Int
k1 = Int -> Int
forall a. Enum a => a -> a
pred Int
k0
                 i1 :: Int
i1 = Int
i0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n
             in  if Int
k1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0
                   then (Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod (Int
i1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
m) Int
len, Int
m)
                   else (Int -> Int -> Int
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 =
   [Int] -> T
forall a. Storable a => [a] -> Vector a
SV.pack ([Int] -> T) -> [Int] -> T
forall a b. (a -> b) -> a -> b
$
   (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
k ->
      let Just (Int
i,Int
j) = Int -> Int -> Int -> Maybe (Int, Int)
forall a. C a => a -> a -> a -> Maybe (a, a)
PID.diophantine Int
k Int
n Int
m
      in  Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
i Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
j Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$
   Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int
1Int -> Int -> Int
forall 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
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m
       (Int
ni,Int
mi) = (Int, (Int, Int)) -> (Int, Int)
forall a b. (a, b) -> b
snd ((Int, (Int, Int)) -> (Int, Int))
-> (Int, (Int, Int)) -> (Int, Int)
forall a b. (a -> b) -> a -> b
$ Int -> Int -> (Int, (Int, Int))
forall a. C a => a -> a -> (a, (a, a))
PID.extendedGCD Int
n Int
m
   in  (T, Maybe (Int, Int)) -> T
forall a b. (a, b) -> a
fst ((T, Maybe (Int, Int)) -> T) -> (T, Maybe (Int, Int)) -> T
forall a b. (a -> b) -> a -> b
$ Int
-> ((Int, Int) -> Maybe (Int, (Int, Int)))
-> (Int, Int)
-> (T, Maybe (Int, Int))
forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN Int
len
          (\(Int
i0,Int
k0) -> (Int, (Int, Int)) -> Maybe (Int, (Int, Int))
forall a. a -> Maybe a
Just (Int
i0,
             let k1 :: Int
k1 = Int -> Int
forall a. Enum a => a -> a
pred Int
k0
                 i1 :: Int
i1 = Int
i0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
niInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n
             in  if Int
k1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0
                   then (Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod (Int
i1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
miInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m) Int
len, Int
m)
                   else (Int -> Int -> Int
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 =
   (T, [Int]) -> T
forall a b. (a, b) -> a
fst ((T, [Int]) -> T) -> (T, [Int]) -> T
forall a b. (a -> b) -> a -> b
$ Int -> [Int] -> (T, [Int])
forall a. Storable a => Int -> [a] -> (Vector a, [a])
SV.packN (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m) ([Int] -> (T, [Int])) -> [Int] -> (T, [Int])
forall a b. (a -> b) -> a -> b
$
   (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
k Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
k Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$
   (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int
1Int -> Int -> Int
forall 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 = Int -> Integer
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 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
forall a. Bounded a => a
maxBound :: Int)
         then [Char] -> T
forall a. HasCallStack => [Char] -> a
error [Char]
"signal too long for Int indexing"
         else (T, Maybe Int) -> T
forall a b. (a, b) -> a
fst ((T, Maybe Int) -> T) -> (T, Maybe Int) -> T
forall a b. (a -> b) -> a -> b
$ Int -> (Int -> Maybe (Int, Int)) -> Int -> (T, Maybe Int)
forall b a.
Storable b =>
Int -> (a -> Maybe (b, a)) -> a -> (Vector b, Maybe a)
SV.unfoldrN (Int
niInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
                 (\Int
x -> (Int, Int) -> Maybe (Int, Int)
forall a. a -> Maybe a
Just (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
gen Int -> Int -> Int
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 s. ST s (Vector s Int)) -> T
forall e. Storable e => (forall s. ST s (Vector s e)) -> Vector e
SVST.runSTVector
      (do Vector s Int
inv <- Int -> ST s (Vector s Int)
forall e s. Storable e => Int -> ST s (Vector s e)
SVST.new_ (T -> Int
forall a. Vector a -> Int
SV.length T
perm)
          T (ST s ()) -> ST s ()
forall (m :: * -> *) a. Monad m => T (m a) -> m ()
SigS.sequence_ (T (ST s ()) -> ST s ()) -> T (ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$
             (Int -> Int -> ST s ()) -> T Int -> T Int -> T (ST s ())
forall a b c. (a -> b -> c) -> T a -> T b -> T c
SigS.zipWith (Vector s Int -> Int -> Int -> ST s ()
forall e s. Storable e => Vector s e -> Int -> e -> ST s ()
SVST.write Vector s Int
inv)
                (T -> T Int
forall a. Storable a => Vector a -> T a
SigS.fromStrictStorableSignal T
perm)
                ((Int -> Int) -> Int -> T Int
forall a. (a -> a) -> a -> T a
SigS.iterate (Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+) Int
0)
          Vector s Int -> ST s (Vector s Int)
forall a. a -> ST s a
forall (m :: * -> *) a. Monad m => a -> m a
return Vector s Int
inv)

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