-- | Low level stuff

{-# LANGUAGE BangPatterns, DataKinds, KindSignatures, GADTs #-}
{-# LANGUAGE ScopedTypeVariables, TypeApplications #-}

module Math.FiniteField.GaloisField.Small.Internal where

--------------------------------------------------------------------------------

import Data.Bits
import Data.Int
import Data.Word
import Data.List

import GHC.TypeNats (Nat)

import Control.Monad
import Control.Monad.ST

import Foreign.C
import Foreign.Ptr
import Foreign.Storable
import Foreign.Marshal

import Data.Vector.Unboxed (Vector, MVector)
import qualified Data.Vector.Unboxed         as Vec
import qualified Data.Vector.Generic.Mutable as MVec

import System.IO.Unsafe as Unsafe
import Control.Monad.ST.Unsafe (unsafeIOToST)      -- debugging only

import Math.FiniteField.Primes
import Math.FiniteField.TypeLevel
import Math.FiniteField.Conway
import Math.FiniteField.Conway.Internal
import Math.FiniteField.Misc
 
import qualified Math.FiniteField.PrimeField.Small.Raw as Raw

--------------------------------------------------------------------------------  

type P = Word64
type M = Int
type C = Ptr Word32
type F = Vector Word32

--------------------------------------------------------------------------------  

neg :: P -> F -> F
neg :: Word64 -> F -> F
neg !Word64
p !F
xs = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
Vec.map Word32 -> Word32
neg1 F
xs where
  neg1 :: Word32 -> Word32
neg1 !Word32
x = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64
Raw.neg Word64
p (Word32 -> Word64
w32toW64 Word32
x))

add :: P -> F -> F -> F
add :: Word64 -> F -> F -> F
add !Word64
p !F
xs !F
ys = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
Vec.zipWith Word32 -> Word32 -> Word32
add1 F
xs F
ys where 
  add1 :: Word32 -> Word32 -> Word32
add1 !Word32
x !Word32
y = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.add Word64
p (Word32 -> Word64
w32toW64 Word32
x) (Word32 -> Word64
w32toW64 Word32
y))

sub :: P -> F -> F -> F
sub :: Word64 -> F -> F -> F
sub !Word64
p !F
xs !F
ys = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
Vec.zipWith Word32 -> Word32 -> Word32
sub1 F
xs F
ys where
  sub1 :: Word32 -> Word32 -> Word32
sub1 !Word32
x !Word32
y = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.sub Word64
p (Word32 -> Word64
w32toW64 Word32
x) (Word32 -> Word64
w32toW64 Word32
y))

mul :: C -> F -> F -> F
mul :: C -> F -> F -> F
mul !C
ptr !F
xs !F
ys = forall a. IO a -> a
Unsafe.unsafePerformIO (C -> F -> F -> IO F
mulIO C
ptr F
xs F
ys)

inv :: P -> C -> F -> F
inv :: Word64 -> C -> F -> F
inv !Word64
p !C
ptr !F
xs = Word64 -> C -> F -> F
invEuclid Word64
p C
ptr F
xs

div :: P -> C -> F -> F -> F
div :: Word64 -> C -> F -> F -> F
div !Word64
p !C
ptr !F
xs !F
ys = C -> F -> F -> F
mul C
ptr F
xs (Word64 -> C -> F -> F
inv Word64
p C
ptr F
ys)

--------------------------------------------------------------------------------  
-- * multiplication

mulIO :: C -> F -> F -> IO F
mulIO :: C -> F -> F -> IO F
mulIO !C
ptr !F
xs !F
ys = 
  do
    !Word32
p0 <- forall a. Storable a => Ptr a -> IO a
peek C
ptr             :: IO Word32
    !Word32
m0 <- forall a. Storable a => Ptr a -> IO a
peek (forall a b. Ptr a -> Int -> Ptr b
plusPtr C
ptr Int
4) :: IO Word32
    let !p :: Word64
p  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Word32
p0   :: Word64
        !m :: Int
m  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Word32
m0   :: Int
        !m1 :: Int
m1 = Int
m forall a. Num a => a -> a -> a
- Int
1             :: Int
    F
cvec <- forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN (Int
mforall a. Num a => a -> a -> a
+Int
1) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray (Int
mforall a. Num a => a -> a -> a
+Int
1) (forall a b. Ptr a -> Int -> Ptr b
plusPtr C
ptr Int
8)
    MVector RealWorld Word32
mvec <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
Vec.unsafeThaw (Word64 -> F -> F -> F
mulPoly Word64
p F
xs F
ys)
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Int -> [Int]
downIndices Int
m1) forall a b. (a -> b) -> a -> b
$ \ !Int
k -> do
      Word32
c0 <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MVec.unsafeRead MVector RealWorld Word32
mvec (Int
mforall a. Num a => a -> a -> a
+Int
kforall a. Num a => a -> a -> a
-Int
1)
      -- print (k,c0)
      let !c :: Word64
c = forall a b. (Integral a, Num b) => a -> b
fromIntegral Word32
c0 :: Word64
      forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Word64
c forall a. Eq a => a -> a -> Bool
/= Word64
0) forall a b. (a -> b) -> a -> b
$ do
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
m] forall a b. (a -> b) -> a -> b
$ \ !Int
j -> do     -- NOTE: this could be [0..m1], a bit faster, but while debugging i leave it as is
          let !y :: Word64
y = Word32 -> Word64
w32toW64 (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
cvec Int
j)
          forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> (a -> a) -> Int -> m ()
MVec.unsafeModify MVector RealWorld Word32
mvec (\ !Word32
z -> Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.sub Word64
p (Word32 -> Word64
w32toW64 Word32
z) (Word64 -> Word64 -> Word64 -> Word64
Raw.mul Word64
p Word64
c Word64
y))) (Int
jforall a. Num a => a -> a -> a
+Int
kforall a. Num a => a -> a -> a
-Int
1)
    -- print =<< Vec.freeze mvec
    [Word32]
xs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
0..Int
mforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \ !Int
i -> forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MVec.unsafeRead MVector RealWorld Word32
mvec Int
i
    forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN Int
m [Word32]
xs

--------------------------------------------------------------------------------  e
-- * polynomials

-- -- | Multiplication of polynomials. Note: we assume both have length @m@
-- mulVecPolyFix  :: P -> M -> Vector Word32 -> Vector Word32 -> Vector Word32
-- mulVecPolyFIx !p !m !xs !ys = 
--   Vec.create $ do
--     mvec <- MVec.replicate mbig 0
--     forM_ [0..m-1] $ \ !i -> do
--       let !x = w32toW64 (Vec.unsafeIndex xs i)
--       forM_ [0..m-1] $ \ !j -> do
--         let !y = w32toW64 (Vec.unsafeIndex ys j)
--         MVec.unsafeModify mvec (\ !z -> w64toW32 (Raw.add p (w32toW64 z) (Raw.mul p x y))) (i+j)
--     return mvec
--   where
--     !mbig = m+m-1

resizePoly :: Int -> Vector Word32 -> Vector Word32
resizePoly :: Int -> F -> F
resizePoly Int
m F
vec = forall a. Unbox a => (forall s. ST s (MVector s a)) -> Vector a
Vec.create forall a b. (a -> b) -> a -> b
$ do
  let !l :: Int
l = forall a. Ord a => a -> a -> a
min Int
m (forall a. Unbox a => Vector a -> Int
Vec.length F
vec) forall a. Num a => a -> a -> a
- Int
1
  MVector s Word32
mvec <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
Int -> a -> m (v (PrimState m) a)
MVec.replicate Int
m Word32
0
  forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
l] forall a b. (a -> b) -> a -> b
$ \ !Int
i -> forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
MVec.unsafeWrite MVector s Word32
mvec Int
i (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
vec Int
i)
  forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Word32
mvec

-- | The vectors can have different lengths
addPoly :: P -> Vector Word32 -> Vector Word32 -> Vector Word32
addPoly :: Word64 -> F -> F -> F
addPoly !Word64
p !F
xs !F
ys 
  | forall a. Unbox a => Vector a -> Int
Vec.length F
xs forall a. Eq a => a -> a -> Bool
== forall a. Unbox a => Vector a -> Int
Vec.length F
ys = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
Vec.zipWith Word32 -> Word32 -> Word32
add1 F
xs F
ys 
  | Bool
otherwise = forall a. Unbox a => [a] -> Vector a
Vec.fromList (forall a b c. a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith Word32
0 Word32
0 Word32 -> Word32 -> Word32
add1 (forall a. Unbox a => Vector a -> [a]
Vec.toList F
xs) (forall a. Unbox a => Vector a -> [a]
Vec.toList F
ys))
  where
    add1 :: Word32 -> Word32 -> Word32
add1 !Word32
x !Word32
y = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.add Word64
p (Word32 -> Word64
w32toW64 Word32
x) (Word32 -> Word64
w32toW64 Word32
y))

-- | The vectors can have different lengths
subPoly :: P -> Vector Word32 -> Vector Word32 -> Vector Word32
subPoly :: Word64 -> F -> F -> F
subPoly !Word64
p !F
xs !F
ys 
  | forall a. Unbox a => Vector a -> Int
Vec.length F
xs forall a. Eq a => a -> a -> Bool
== forall a. Unbox a => Vector a -> Int
Vec.length F
ys = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
Vec.zipWith Word32 -> Word32 -> Word32
sub1 F
xs F
ys 
  | Bool
otherwise = forall a. Unbox a => [a] -> Vector a
Vec.fromList (forall a b c. a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith Word32
0 Word32
0 Word32 -> Word32 -> Word32
sub1 (forall a. Unbox a => Vector a -> [a]
Vec.toList F
xs) (forall a. Unbox a => Vector a -> [a]
Vec.toList F
ys))
  where
    sub1 :: Word32 -> Word32 -> Word32
sub1 !Word32
x !Word32
y = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.sub Word64
p (Word32 -> Word64
w32toW64 Word32
x) (Word32 -> Word64
w32toW64 Word32
y))

scalePoly :: P -> Word32 -> Vector Word32 -> Vector Word32
scalePoly :: Word64 -> Word32 -> F -> F
scalePoly !Word64
p !Word32
s !F
ys = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
Vec.map Word32 -> Word32
mul1 F
ys where
  mul1 :: Word32 -> Word32
mul1 !Word32
y = Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.mul Word64
p Word64
s64 (Word32 -> Word64
w32toW64 Word32
y))
  s64 :: Word64
s64 = (Word32 -> Word64
w32toW64 Word32
s)

mulPoly  :: P -> Vector Word32 -> Vector Word32 -> Vector Word32
mulPoly :: Word64 -> F -> F -> F
mulPoly !Word64
p !F
xs !F
ys =
  forall a. Unbox a => (forall s. ST s (MVector s a)) -> Vector a
Vec.create forall a b. (a -> b) -> a -> b
$ do
    MVector s Word32
mvec <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
Int -> a -> m (v (PrimState m) a)
MVec.replicate Int
mbig Word32
0
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
m1forall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \ !Int
i -> do
      let !x :: Word64
x = Word32 -> Word64
w32toW64 (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
xs Int
i)
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
m2forall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \ !Int
j -> do
        let !y :: Word64
y = Word32 -> Word64
w32toW64 (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
ys Int
j)
        forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> (a -> a) -> Int -> m ()
MVec.unsafeModify MVector s Word32
mvec (\ !Word32
z -> Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.add Word64
p (Word32 -> Word64
w32toW64 Word32
z) (Word64 -> Word64 -> Word64 -> Word64
Raw.mul Word64
p Word64
x Word64
y))) (Int
iforall a. Num a => a -> a -> a
+Int
j)
    forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Word32
mvec
  where
    !m1 :: Int
m1 = forall a. Unbox a => Vector a -> Int
Vec.length F
xs
    !m2 :: Int
m2 = forall a. Unbox a => Vector a -> Int
Vec.length F
ys
    !mbig :: Int
mbig = Int
m1forall a. Num a => a -> a -> a
+Int
m2forall a. Num a => a -> a -> a
-Int
1

zeroPoly :: Int -> Vector Word32 
zeroPoly :: Int -> F
zeroPoly Int
n = forall a. Unbox a => Int -> a -> Vector a
Vec.replicate Int
n Word32
0

onePoly :: Int -> Vector Word32
onePoly :: Int -> F
onePoly Int
n = forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN Int
n (Word32
1 forall a. a -> [a] -> [a]
: forall a. Int -> a -> [a]
replicate (Int
nforall a. Num a => a -> a -> a
-Int
1) Word32
0)

isZeroPoly :: Vector Word32 -> Bool
isZeroPoly :: F -> Bool
isZeroPoly F
v = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (forall a. Eq a => a -> a -> Bool
==Word32
0) (forall a. Unbox a => Vector a -> [a]
Vec.toList F
v)

polyDegree :: Vector Word32 -> Int
polyDegree :: F -> Int
polyDegree F
v = Int -> Int
go (Int
nforall a. Num a => a -> a -> a
-Int
1) where
  !n :: Int
n = forall a. Unbox a => Vector a -> Int
Vec.length F
v
  go :: Int -> Int
go !Int
i = if Int
i forall a. Ord a => a -> a -> Bool
< Int
0 
    then (-Int
1) 
    else if forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
v Int
i forall a. Eq a => a -> a -> Bool
/= Word32
0 
      then Int
i
      else Int -> Int
go (Int
iforall a. Num a => a -> a -> a
-Int
1)

polyLongDivision :: P -> Vector Word32 -> Vector Word32 -> (Vector Word32, Vector Word32)
polyLongDivision :: Word64 -> F -> F -> (F, F)
polyLongDivision !Word64
p !F
as !F
bs = forall a. (forall s. ST s a) -> a
runST forall s. ST s (F, F)
action where
  action :: forall s. ST s (Vector Word32, Vector Word32) 
  action :: forall s. ST s (F, F)
action = do
    let !m :: Int
m = forall a. Unbox a => Vector a -> Int
Vec.length F
as
    let !d :: Int
d = F -> Int
polyDegree F
bs
        !lcf0 :: Word32
lcf0 = forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
bs Int
d   -- leading coefficient
        !lcf :: Word64
lcf  = Word32 -> Word64
w32toW64 Word32
lcf0
    MVector s Word32
rem  <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
Vec.thaw F
as
    MVector s Word32
quot <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
Int -> a -> m (v (PrimState m) a)
MVec.replicate @(ST s) @MVector Int
m (Word32
0 :: Word32) 
   
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Int -> [Int]
downIndices (Int
mforall a. Num a => a -> a -> a
-Int
d)) forall a b. (a -> b) -> a -> b
$ \ !Int
k -> do
      Word32
c0 <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MVec.unsafeRead MVector s Word32
rem (Int
dforall a. Num a => a -> a -> a
+Int
kforall a. Num a => a -> a -> a
-Int
1)
      -- print (k,c0)
      forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Word32
c0 forall a. Eq a => a -> a -> Bool
/= Word32
0) forall a b. (a -> b) -> a -> b
$ do
        let !s :: Word64
s = Word64 -> Word64 -> Word64 -> Word64
Raw.div Word64
p (Word32 -> Word64
w32toW64 Word32
c0) Word64
lcf
        -- Vec.freeze rem >>= \r -> unsafeIOToST $ print (k,c0,s,r)
        forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
MVec.unsafeWrite MVector s Word32
quot (Int
kforall a. Num a => a -> a -> a
-Int
1) (Word64 -> Word32
w64toW32 Word64
s)
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
d] forall a b. (a -> b) -> a -> b
$ \ !Int
j -> do     -- NOTE: this could be [0..d-1], a bit faster, but while debugging i leave it as is
          let !y :: Word64
y = Word32 -> Word64
w32toW64 (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
bs Int
j)
          forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> (a -> a) -> Int -> m ()
MVec.unsafeModify MVector s Word32
rem (\ !Word32
z -> Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64 -> Word64
Raw.sub Word64
p (Word32 -> Word64
w32toW64 Word32
z) (Word64 -> Word64 -> Word64 -> Word64
Raw.mul Word64
p Word64
s Word64
y))) (Int
jforall a. Num a => a -> a -> a
+Int
kforall a. Num a => a -> a -> a
-Int
1)
  
    [Word32]
fquot <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
0..Int
mforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \ !Int
i -> forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MVec.unsafeRead MVector s Word32
quot Int
i
    [Word32]
frem  <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
0..Int
dforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \ !Int
i -> forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MVec.unsafeRead MVector s Word32
rem  Int
i
    forall (m :: * -> *) a. Monad m => a -> m a
return (forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN Int
m [Word32]
fquot, forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN Int
d [Word32]
frem)

downIndices :: Int -> [Int]
downIndices :: Int -> [Int]
downIndices !Int
k = if Int
k forall a. Ord a => a -> a -> Bool
<= Int
0 then [] else Int
k forall a. a -> [a] -> [a]
: Int -> [Int]
downIndices (Int
kforall a. Num a => a -> a -> a
-Int
1)

polyQuotient :: P -> Vector Word32 -> Vector Word32 -> Vector Word32
polyQuotient :: Word64 -> F -> F -> F
polyQuotient Word64
p F
xs F
ys = forall a b. (a, b) -> a
fst (Word64 -> F -> F -> (F, F)
polyLongDivision Word64
p F
xs F
ys)

--------------------------------------------------------------------------------  

getCPoly :: Ptr Word32 -> Vector Word32
getCPoly :: C -> F
getCPoly C
ptr = forall a. IO a -> a
Unsafe.unsafePerformIO (C -> IO F
getCPolyIO C
ptr)

getCPolyIO :: Ptr Word32 -> IO (Vector Word32)
getCPolyIO :: C -> IO F
getCPolyIO C
ptr = do
  !Word32
m0 <- forall a. Storable a => Ptr a -> IO a
peek (forall a b. Ptr a -> Int -> Ptr b
plusPtr C
ptr Int
4) :: IO Word32
  let !m :: Int
m  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Word32
m0   :: Int
  forall a. Unbox a => Int -> [a] -> Vector a
Vec.fromListN (Int
mforall a. Num a => a -> a -> a
+Int
1) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray (Int
mforall a. Num a => a -> a -> a
+Int
1) (forall a b. Ptr a -> Int -> Ptr b
plusPtr C
ptr Int
8)

-- | Based on <https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Modular_integers>
invEuclid :: P -> C -> Vector Word32 -> Vector Word32
invEuclid :: Word64 -> C -> F -> F
invEuclid !Word64
p !C
ptr !F
a = F -> F -> F -> F -> F
worker (Int -> F
zeroPoly Int
n) (Int -> F
onePoly Int
n) (C -> F
getCPoly C
ptr) F
a where

  n :: Int
n = forall a. Unbox a => Vector a -> Int
Vec.length F
a

  worker :: Vector Word32 -> Vector Word32 -> Vector Word32 -> Vector Word32 -> Vector Word32
  worker :: F -> F -> F -> F -> F
worker !F
t !F
newt !F
r !F
newr = case F -> Bool
isZeroPoly F
newr of

    Bool
False -> let quot :: F
quot = Word64 -> F -> F -> F
polyQuotient Word64
p F
r F
newr 
                 r' :: F
r'   = Word64 -> F -> F -> F
subPoly Word64
p F
r (Word64 -> F -> F -> F
mulPoly Word64
p F
quot F
newr)
                 t' :: F
t'   = Word64 -> F -> F -> F
subPoly Word64
p F
t (Word64 -> F -> F -> F
mulPoly Word64
p F
quot F
newt)
             in  F -> F -> F -> F -> F
worker F
newt F
t' F
newr F
r'

    Bool
True  -> if (F -> Int
polyDegree F
r forall a. Ord a => a -> a -> Bool
> Int
0)
      then Int -> F
zeroPoly (forall a. Unbox a => Vector a -> Int
Vec.length F
a)
      else let r0 :: Word64
r0 = Word32 -> Word64
w32toW64 (forall a. Unbox a => Vector a -> Int -> a
Vec.unsafeIndex F
r Int
0)
           in  Int -> F -> F
resizePoly Int
n (Word64 -> Word32 -> F -> F
scalePoly Word64
p (Word64 -> Word32
w64toW32 (Word64 -> Word64 -> Word64
Raw.inv Word64
p Word64
r0)) F
t)

--------------------------------------------------------------------------------