{- |
module: Factor.Gfpx
description: The polynomial ring GF(p)[x]
license: MIT

maintainer: Joe Leslie-Hurd <joe@gilith.com>
stability: provisional
portability: portable
-}
module Factor.Gfpx
where

import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap
import qualified Data.List as List
import System.Random (RandomGen)

import Factor.Prime (Gfp,Prime)
import qualified Factor.Prime as Prime
import Factor.Term (Term(..),Var)
import qualified Factor.Term as Term
import Factor.Util
import Factor.Zx (Zx)
import qualified Factor.Zx as Zx

-------------------------------------------------------------------------------
-- The polynomial ring GF(p)[x]
-------------------------------------------------------------------------------

data Gfpx =
    Gfpx
      {Gfpx -> Int
degree :: Int,
       Gfpx -> IntMap Gfp
coeffMap :: IntMap Gfp}
  deriving (Gfpx -> Gfpx -> Bool
(Gfpx -> Gfpx -> Bool) -> (Gfpx -> Gfpx -> Bool) -> Eq Gfpx
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Gfpx -> Gfpx -> Bool
$c/= :: Gfpx -> Gfpx -> Bool
== :: Gfpx -> Gfpx -> Bool
$c== :: Gfpx -> Gfpx -> Bool
Eq,Eq Gfpx
Eq Gfpx
-> (Gfpx -> Gfpx -> Ordering)
-> (Gfpx -> Gfpx -> Bool)
-> (Gfpx -> Gfpx -> Bool)
-> (Gfpx -> Gfpx -> Bool)
-> (Gfpx -> Gfpx -> Bool)
-> (Gfpx -> Gfpx -> Gfpx)
-> (Gfpx -> Gfpx -> Gfpx)
-> Ord Gfpx
Gfpx -> Gfpx -> Bool
Gfpx -> Gfpx -> Ordering
Gfpx -> Gfpx -> Gfpx
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Gfpx -> Gfpx -> Gfpx
$cmin :: Gfpx -> Gfpx -> Gfpx
max :: Gfpx -> Gfpx -> Gfpx
$cmax :: Gfpx -> Gfpx -> Gfpx
>= :: Gfpx -> Gfpx -> Bool
$c>= :: Gfpx -> Gfpx -> Bool
> :: Gfpx -> Gfpx -> Bool
$c> :: Gfpx -> Gfpx -> Bool
<= :: Gfpx -> Gfpx -> Bool
$c<= :: Gfpx -> Gfpx -> Bool
< :: Gfpx -> Gfpx -> Bool
$c< :: Gfpx -> Gfpx -> Bool
compare :: Gfpx -> Gfpx -> Ordering
$ccompare :: Gfpx -> Gfpx -> Ordering
$cp1Ord :: Eq Gfpx
Ord)

instance Show Gfpx where
  show :: Gfpx -> String
show = Zx -> String
forall a. Show a => a -> String
show (Zx -> String) -> (Gfpx -> Zx) -> Gfpx -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> Zx
toZx

valid :: Prime -> Gfpx -> Bool
valid :: Gfp -> Gfpx -> Bool
valid Gfp
p Gfpx
f =
    Zx -> Bool
Zx.valid (Gfpx -> Zx
toZx Gfpx
f) Bool -> Bool -> Bool
&&
    (Gfp -> Bool) -> [Gfp] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Gfp -> Gfp -> Bool
Prime.valid Gfp
p) (IntMap Gfp -> [Gfp]
forall a. IntMap a -> [a]
IntMap.elems (Gfpx -> IntMap Gfp
coeffMap Gfpx
f))

fromNormCoeffMap :: IntMap Gfp -> Gfpx
fromNormCoeffMap :: IntMap Gfp -> Gfpx
fromNormCoeffMap IntMap Gfp
c | IntMap Gfp -> Bool
forall a. IntMap a -> Bool
IntMap.null IntMap Gfp
c = Gfpx
zero
fromNormCoeffMap IntMap Gfp
c | Bool
otherwise = Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = Int
d, coeffMap :: IntMap Gfp
coeffMap = IntMap Gfp
c}
  where (Int
d,Gfp
_) = IntMap Gfp -> (Int, Gfp)
forall a. IntMap a -> (Int, a)
IntMap.findMax IntMap Gfp
c

fromCoeffMap :: IntMap Gfp -> Gfpx
fromCoeffMap :: IntMap Gfp -> Gfpx
fromCoeffMap = IntMap Gfp -> Gfpx
fromNormCoeffMap (IntMap Gfp -> Gfpx)
-> (IntMap Gfp -> IntMap Gfp) -> IntMap Gfp -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Gfp -> Bool) -> IntMap Gfp -> IntMap Gfp
forall a. (a -> Bool) -> IntMap a -> IntMap a
IntMap.filter (Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
(/=) Gfp
0)

-------------------------------------------------------------------------------
-- Polynomial operations
-------------------------------------------------------------------------------

isZero :: Gfpx -> Bool
isZero :: Gfpx -> Bool
isZero Gfpx
f = Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0

isOne :: Gfpx -> Bool
isOne :: Gfpx -> Bool
isOne Gfpx
f = Gfpx
f Gfpx -> Gfpx -> Bool
forall a. Eq a => a -> a -> Bool
== Gfpx
one

isConstant :: Gfpx -> Bool
isConstant :: Gfpx -> Bool
isConstant Gfpx
f = Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0

isLinear :: Gfpx -> Bool
isLinear :: Gfpx -> Bool
isLinear Gfpx
f = Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1

isMonic :: Gfpx -> Bool
isMonic :: Gfpx -> Bool
isMonic Gfpx
f = Gfpx -> Gfp
leadingCoeff Gfpx
f Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
1

constMonic :: Prime -> Gfpx -> (Gfp,Gfpx)
constMonic :: Gfp -> Gfpx -> (Gfp, Gfpx)
constMonic Gfp
p Gfpx
f =
    case Gfpx -> Gfp
leadingCoeff Gfpx
f of
      Gfp
0 -> String -> (Gfp, Gfpx)
forall a. HasCallStack => String -> a
error String
"constMonic: zero polynomial"
      Gfp
1 -> (Gfp
1,Gfpx
f)
      Gfp
c -> (Gfp
c, Gfp -> Gfp -> Gfpx -> Gfpx
multiplyConstant Gfp
p (Gfp -> Gfp -> Gfp
Prime.invert Gfp
p Gfp
c) Gfpx
f)

powerCoeff :: Gfpx -> Int -> Gfp
powerCoeff :: Gfpx -> Int -> Gfp
powerCoeff Gfpx
f Int
i = Gfp -> Int -> IntMap Gfp -> Gfp
forall a. a -> Int -> IntMap a -> a
IntMap.findWithDefault Gfp
0 Int
i (Gfpx -> IntMap Gfp
coeffMap Gfpx
f)

constantCoeff :: Gfpx -> Gfp
constantCoeff :: Gfpx -> Gfp
constantCoeff Gfpx
f = Gfpx -> Int -> Gfp
powerCoeff Gfpx
f Int
0

linearCoeff :: Gfpx -> Gfp
linearCoeff :: Gfpx -> Gfp
linearCoeff Gfpx
f = Gfpx -> Int -> Gfp
powerCoeff Gfpx
f Int
1

leadingCoeff :: Gfpx -> Gfp
leadingCoeff :: Gfpx -> Gfp
leadingCoeff Gfpx
f = Gfpx -> Int -> Gfp
powerCoeff Gfpx
f (Gfpx -> Int
degree Gfpx
f)

monomials :: Gfpx -> [(Int,Gfp)]
monomials :: Gfpx -> [(Int, Gfp)]
monomials = IntMap Gfp -> [(Int, Gfp)]
forall a. IntMap a -> [(Int, a)]
IntMap.toAscList (IntMap Gfp -> [(Int, Gfp)])
-> (Gfpx -> IntMap Gfp) -> Gfpx -> [(Int, Gfp)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> IntMap Gfp
coeffMap

lengthMonomials :: Gfpx -> Int
lengthMonomials :: Gfpx -> Int
lengthMonomials = IntMap Gfp -> Int
forall a. IntMap a -> Int
IntMap.size (IntMap Gfp -> Int) -> (Gfpx -> IntMap Gfp) -> Gfpx -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> IntMap Gfp
coeffMap

filterMonomials :: (Int -> Gfp -> Bool) -> Gfpx -> Gfpx
filterMonomials :: (Int -> Gfp -> Bool) -> Gfpx -> Gfpx
filterMonomials Int -> Gfp -> Bool
p = IntMap Gfp -> Gfpx
fromNormCoeffMap (IntMap Gfp -> Gfpx) -> (Gfpx -> IntMap Gfp) -> Gfpx -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Gfp -> Bool) -> IntMap Gfp -> IntMap Gfp
forall a. (Int -> a -> Bool) -> IntMap a -> IntMap a
IntMap.filterWithKey Int -> Gfp -> Bool
p (IntMap Gfp -> IntMap Gfp)
-> (Gfpx -> IntMap Gfp) -> Gfpx -> IntMap Gfp
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> IntMap Gfp
coeffMap

constant :: Gfp -> Gfpx
constant :: Gfp -> Gfpx
constant = Int -> Gfp -> Gfpx
monomial Int
0

variable :: Gfpx
variable :: Gfpx
variable = Int -> Gfp -> Gfpx
monomial Int
1 Gfp
1

monomial :: Int -> Gfp -> Gfpx
monomial :: Int -> Gfp -> Gfpx
monomial Int
_ Gfp
0 = Gfpx
zero
monomial Int
d Gfp
x = Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = Int
d, coeffMap :: IntMap Gfp
coeffMap = Int -> Gfp -> IntMap Gfp
forall a. Int -> a -> IntMap a
IntMap.singleton Int
d Gfp
x}

simpleRoot :: Prime -> Gfp -> Gfpx
simpleRoot :: Gfp -> Gfp -> Gfpx
simpleRoot Gfp
p Gfp
r = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p Gfpx
variable (Gfp -> Gfpx
constant Gfp
r)

evaluate :: Prime -> Gfpx -> Gfp -> Gfp
evaluate :: Gfp -> Gfpx -> Gfp -> Gfp
evaluate Gfp
p Gfpx
f Gfp
x = Int -> (Int, Gfp) -> Gfp
forall a. Integral a => a -> (a, Gfp) -> Gfp
align Int
0 ((Int, Gfp) -> Gfp) -> (Int, Gfp) -> Gfp
forall a b. (a -> b) -> a -> b
$ (Int -> Gfp -> (Int, Gfp) -> (Int, Gfp))
-> (Int, Gfp) -> IntMap Gfp -> (Int, Gfp)
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Gfp -> (Int, Gfp) -> (Int, Gfp)
forall a. Integral a => a -> Gfp -> (a, Gfp) -> (a, Gfp)
fma (Int
0,Gfp
0) (IntMap Gfp -> (Int, Gfp)) -> IntMap Gfp -> (Int, Gfp)
forall a b. (a -> b) -> a -> b
$ Gfpx -> IntMap Gfp
coeffMap Gfpx
f
  where
    fma :: a -> Gfp -> (a, Gfp) -> (a, Gfp)
fma a
i Gfp
c (a, Gfp)
z = (a
i, Gfp -> Gfp -> Gfp -> Gfp
Prime.add Gfp
p (a -> (a, Gfp) -> Gfp
forall a. Integral a => a -> (a, Gfp) -> Gfp
align a
i (a, Gfp)
z) Gfp
c)

    align :: a -> (a, Gfp) -> Gfp
align a
i (a
d,Gfp
z) =
        if a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then Gfp
z else Gfp -> Gfp -> Gfp -> Gfp -> Gfp
Prime.multiplyExp Gfp
p Gfp
z Gfp
x (a -> Gfp
forall a. Integral a => a -> Gfp
toInteger a
k)
      where
        k :: a
k = a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
i

derivative :: Prime -> Gfpx -> Gfpx
derivative :: Gfp -> Gfpx -> Gfpx
derivative Gfp
p (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
cm}) =
    if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 then Gfpx
zero
    else Int -> Gfpx -> Gfpx
multiplyPower (-Int
1) (Gfpx -> Gfpx) -> Gfpx -> Gfpx
forall a b. (a -> b) -> a -> b
$ IntMap Gfp -> Gfpx
fromNormCoeffMap (IntMap Gfp -> Gfpx) -> IntMap Gfp -> Gfpx
forall a b. (a -> b) -> a -> b
$ IntMap Gfp -> IntMap Gfp
deriv IntMap Gfp
cm
  where
    deriv :: IntMap Gfp -> IntMap Gfp
deriv = (Int -> Gfp -> Maybe Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (Int -> a -> Maybe b) -> IntMap a -> IntMap b
IntMap.mapMaybeWithKey Int -> Gfp -> Maybe Gfp
multDeg
    multDeg :: Int -> Gfp -> Maybe Gfp
multDeg Int
i Gfp
c = if Gfp
c' Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0 then Maybe Gfp
forall a. Maybe a
Nothing else Gfp -> Maybe Gfp
forall a. a -> Maybe a
Just Gfp
c'
      where c' :: Gfp
c' = Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p (Gfp -> Int -> Gfp
Prime.fromInt Gfp
p Int
i) Gfp
c

fromCoeff :: [Gfp] -> Gfpx
fromCoeff :: [Gfp] -> Gfpx
fromCoeff = IntMap Gfp -> Gfpx
fromCoeffMap (IntMap Gfp -> Gfpx) -> ([Gfp] -> IntMap Gfp) -> [Gfp] -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Int, Gfp)] -> IntMap Gfp
forall a. [(Int, a)] -> IntMap a
IntMap.fromList ([(Int, Gfp)] -> IntMap Gfp)
-> ([Gfp] -> [(Int, Gfp)]) -> [Gfp] -> IntMap Gfp
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Int] -> [Gfp] -> [(Int, Gfp)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..]

fromZx :: Prime -> Zx -> Gfpx
fromZx :: Gfp -> Zx -> Gfpx
fromZx Gfp
p = IntMap Gfp -> Gfpx
fromNormCoeffMap (IntMap Gfp -> Gfpx) -> (Zx -> IntMap Gfp) -> Zx -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Gfp -> Maybe Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (a -> Maybe b) -> IntMap a -> IntMap b
IntMap.mapMaybe Gfp -> Maybe Gfp
f (IntMap Gfp -> IntMap Gfp)
-> (Zx -> IntMap Gfp) -> Zx -> IntMap Gfp
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> IntMap Gfp
Zx.coeffMap
  where
    f :: Gfp -> Maybe Gfp
f Gfp
c = if Gfp
x Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0 then Maybe Gfp
forall a. Maybe a
Nothing else Gfp -> Maybe Gfp
forall a. a -> Maybe a
Just Gfp
x
      where x :: Gfp
x = Gfp -> Gfp -> Gfp
Prime.fromInteger Gfp
p Gfp
c

toZx :: Gfpx -> Zx
toZx :: Gfpx -> Zx
toZx (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
c}) = Zx :: Int -> IntMap Gfp -> Zx
Zx.Zx {degree :: Int
Zx.degree = Int
d, coeffMap :: IntMap Gfp
Zx.coeffMap = IntMap Gfp
c}

toSmallestZx :: Prime -> Gfpx -> Zx
toSmallestZx :: Gfp -> Gfpx -> Zx
toSmallestZx Gfp
p (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
c}) =
    Zx :: Int -> IntMap Gfp -> Zx
Zx.Zx {degree :: Int
Zx.degree = Int
d, coeffMap :: IntMap Gfp
Zx.coeffMap = (Gfp -> Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (a -> b) -> IntMap a -> IntMap b
IntMap.map (Gfp -> Gfp -> Gfp
Prime.toSmallestInteger Gfp
p) IntMap Gfp
c}

uniform :: RandomGen r => Prime -> Int -> r -> (Gfpx,r)
uniform :: Gfp -> Int -> r -> (Gfpx, r)
uniform Gfp
_ Int
d r
r | Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = (Gfpx
zero,r
r)
uniform Gfp
p Int
d r
r = ([Gfp] -> Gfpx
fromCoeff [Gfp]
l, r
r')
  where ([Gfp]
l,r
r') = (r -> (Gfp, r)) -> Int -> r -> ([Gfp], r)
forall b a. (b -> (a, b)) -> Int -> b -> ([a], b)
unfoldrN (Gfp -> r -> (Gfp, r)
forall r. RandomGen r => Gfp -> r -> (Gfp, r)
Prime.uniform Gfp
p) (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) r
r

polyToTerm :: Prime -> Var -> Gfpx -> Term
polyToTerm :: Gfp -> String -> Gfpx -> Term
polyToTerm Gfp
p String
v Gfpx
f = String -> Zx -> Term
Zx.toTerm String
v (Gfp -> Gfpx -> Zx
toSmallestZx Gfp
p Gfpx
f)

toTerm :: Prime -> Var -> Gfpx -> Term
toTerm :: Gfp -> String -> Gfpx -> Term
toTerm Gfp
p String
v Gfpx
f = Term -> Gfp -> Term
Term.modulo (Gfp -> String -> Gfpx -> Term
polyToTerm Gfp
p String
v Gfpx
f) Gfp
p

-------------------------------------------------------------------------------
-- Ring operations
-------------------------------------------------------------------------------

zero :: Gfpx
zero :: Gfpx
zero = Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = -Int
1, coeffMap :: IntMap Gfp
coeffMap = IntMap Gfp
forall a. IntMap a
IntMap.empty}

one :: Gfpx
one :: Gfpx
one = Gfp -> Gfpx
constant Gfp
1

negate :: Prime -> Gfpx -> Gfpx
negate :: Gfp -> Gfpx -> Gfpx
negate Gfp
_ Gfpx
f | Gfpx -> Bool
isZero Gfpx
f = Gfpx
zero
negate Gfp
p Gfpx
f = Gfpx
f {coeffMap :: IntMap Gfp
coeffMap = (Gfp -> Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (a -> b) -> IntMap a -> IntMap b
IntMap.map (Gfp -> Gfp -> Gfp
Prime.negate Gfp
p) (Gfpx -> IntMap Gfp
coeffMap Gfpx
f)}

add :: Prime -> Gfpx -> Gfpx -> Gfpx
add :: Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
_ Gfpx
f Gfpx
g | Gfpx -> Bool
isZero Gfpx
f = Gfpx
g
add Gfp
_ Gfpx
f Gfpx
g | Gfpx -> Bool
isZero Gfpx
g = Gfpx
f
add Gfp
p Gfpx
f Gfpx
g | Bool
otherwise = IntMap Gfp -> Gfpx
fromNormCoeffMap IntMap Gfp
c
  where
    c :: IntMap Gfp
c = (Int -> Gfp -> Gfp -> Maybe Gfp)
-> (IntMap Gfp -> IntMap Gfp)
-> (IntMap Gfp -> IntMap Gfp)
-> IntMap Gfp
-> IntMap Gfp
-> IntMap Gfp
forall a b c.
(Int -> a -> b -> Maybe c)
-> (IntMap a -> IntMap c)
-> (IntMap b -> IntMap c)
-> IntMap a
-> IntMap b
-> IntMap c
IntMap.mergeWithKey Int -> Gfp -> Gfp -> Maybe Gfp
forall p. p -> Gfp -> Gfp -> Maybe Gfp
addCoeff IntMap Gfp -> IntMap Gfp
forall a. a -> a
id IntMap Gfp -> IntMap Gfp
forall a. a -> a
id (Gfpx -> IntMap Gfp
coeffMap Gfpx
f) (Gfpx -> IntMap Gfp
coeffMap Gfpx
g)

    addCoeff :: p -> Gfp -> Gfp -> Maybe Gfp
addCoeff p
_ Gfp
fx Gfp
gx = if Gfp
x Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0 then Maybe Gfp
forall a. Maybe a
Nothing else Gfp -> Maybe Gfp
forall a. a -> Maybe a
Just Gfp
x
      where x :: Gfp
x = Gfp -> Gfp -> Gfp -> Gfp
Prime.add Gfp
p Gfp
fx Gfp
gx

sum :: Prime -> [Gfpx] -> Gfpx
sum :: Gfp -> [Gfpx] -> Gfpx
sum Gfp
p = (Gfpx -> Gfpx -> Gfpx) -> Gfpx -> [Gfpx] -> Gfpx
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p) Gfpx
zero

subtract :: Prime -> Gfpx -> Gfpx -> Gfpx
subtract :: Gfp -> Gfpx -> Gfpx -> Gfpx
subtract Gfp
p Gfpx
f Gfpx
g = Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx
Factor.Gfpx.negate Gfp
p Gfpx
g)

multiply :: Prime -> Gfpx -> Gfpx -> Gfpx
multiply :: Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
f Gfpx
g | Gfpx -> Int
lengthMonomials Gfpx
g Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Gfpx -> Int
lengthMonomials Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
g Gfpx
f
multiply Gfp
_ Gfpx
f Gfpx
_ | Gfpx -> Bool
isZero Gfpx
f = Gfpx
zero
multiply Gfp
p Gfpx
f Gfpx
g | Bool
otherwise = (Int -> Gfp -> Gfpx -> Gfpx) -> Gfpx -> IntMap Gfp -> Gfpx
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Gfp -> Gfpx -> Gfpx
fma Gfpx
zero (Gfpx -> IntMap Gfp
coeffMap Gfpx
f)
  where fma :: Int -> Gfp -> Gfpx -> Gfpx
fma Int
i Gfp
c Gfpx
z = Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p (Int -> Gfpx -> Gfpx
multiplyPower Int
i (Gfp -> Gfp -> Gfpx -> Gfpx
multiplyConstant Gfp
p Gfp
c Gfpx
g)) Gfpx
z

square :: Prime -> Gfpx -> Gfpx
square :: Gfp -> Gfpx -> Gfpx
square Gfp
p Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
f Gfpx
f

cube :: Prime -> Gfpx -> Gfpx
cube :: Gfp -> Gfpx -> Gfpx
cube Gfp
p Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx
square Gfp
p Gfpx
f)

product :: Prime -> [Gfpx] -> Gfpx
product :: Gfp -> [Gfpx] -> Gfpx
product Gfp
p = (Gfpx -> Gfpx -> Gfpx) -> Gfpx -> [Gfpx] -> Gfpx
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p) Gfpx
one

multiplyConstant :: Prime -> Gfp -> Gfpx -> Gfpx
multiplyConstant :: Gfp -> Gfp -> Gfpx -> Gfpx
multiplyConstant Gfp
_ Gfp
0 Gfpx
_ = Gfpx
zero
multiplyConstant Gfp
_ Gfp
1 Gfpx
f = Gfpx
f
multiplyConstant Gfp
p Gfp
c Gfpx
f = IntMap Gfp -> Gfpx
fromNormCoeffMap (IntMap Gfp -> Gfpx) -> IntMap Gfp -> Gfpx
forall a b. (a -> b) -> a -> b
$ (Gfp -> Maybe Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (a -> Maybe b) -> IntMap a -> IntMap b
IntMap.mapMaybe Gfp -> Maybe Gfp
m (IntMap Gfp -> IntMap Gfp) -> IntMap Gfp -> IntMap Gfp
forall a b. (a -> b) -> a -> b
$ Gfpx -> IntMap Gfp
coeffMap Gfpx
f
  where m :: Gfp -> Maybe Gfp
m Gfp
x = if Gfp
y Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0 then Maybe Gfp
forall a. Maybe a
Nothing else Gfp -> Maybe Gfp
forall a. a -> Maybe a
Just Gfp
y where y :: Gfp
y = Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p Gfp
c Gfp
x

multiplyPower :: Int -> Gfpx -> Gfpx
multiplyPower :: Int -> Gfpx -> Gfpx
multiplyPower Int
0 Gfpx
f = Gfpx
f
multiplyPower Int
i (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
c}) =
    if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then Gfpx
zero
    else Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i, coeffMap :: IntMap Gfp
coeffMap = (Int -> Int) -> IntMap Gfp -> IntMap Gfp
forall a. (Int -> Int) -> IntMap a -> IntMap a
IntMap.mapKeysMonotonic (Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
i) IntMap Gfp
c}

multiplyExp :: Prime -> Gfpx -> Gfpx -> Integer -> Gfpx
multiplyExp :: Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
multiplyExp Gfp
_ Gfpx
z Gfpx
_ Gfp
_ | Gfpx -> Bool
isZero Gfpx
z = Gfpx
zero
multiplyExp Gfp
_ Gfpx
z Gfpx
_ Gfp
0 = Gfpx
z
multiplyExp Gfp
_ Gfpx
_ Gfpx
f Gfp
_ | Gfpx -> Bool
isZero Gfpx
f = Gfpx
zero
multiplyExp Gfp
_ Gfpx
z Gfpx
f Gfp
_ | Gfpx -> Bool
isOne Gfpx
f = Gfpx
z
multiplyExp Gfp
p Gfpx
z0 Gfpx
f0 Gfp
k0 = Gfpx -> Gfpx -> Gfp -> Gfpx
forall t. Integral t => Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z0 Gfpx
f0 Gfp
k0
  where
    go :: Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z Gfpx
_ t
0 = Gfpx
z
    go Gfpx
z Gfpx
f t
k = Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z' Gfpx
f' t
k'
      where
        z' :: Gfpx
z' = if t -> Bool
forall a. Integral a => a -> Bool
even t
k then Gfpx
z else Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
z Gfpx
f
        f' :: Gfpx
f' = Gfp -> Gfpx -> Gfpx
square Gfp
p Gfpx
f
        k' :: t
k' = t
k t -> t -> t
forall a. Integral a => a -> a -> a
`div` t
2

exp :: Prime -> Gfpx -> Integer -> Gfpx
exp :: Gfp -> Gfpx -> Gfp -> Gfpx
exp Gfp
p = Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
multiplyExp Gfp
p Gfpx
one

-------------------------------------------------------------------------------
-- Polynomial composition
-------------------------------------------------------------------------------

compose :: Prime -> Gfpx -> Gfpx -> Gfpx
compose :: Gfp -> Gfpx -> Gfpx -> Gfpx
compose Gfp
p Gfpx
f Gfpx
g = Int -> (Int, Gfpx) -> Gfpx
forall a. Integral a => a -> (a, Gfpx) -> Gfpx
align Int
0 ((Int, Gfpx) -> Gfpx) -> (Int, Gfpx) -> Gfpx
forall a b. (a -> b) -> a -> b
$ (Int -> Gfp -> (Int, Gfpx) -> (Int, Gfpx))
-> (Int, Gfpx) -> IntMap Gfp -> (Int, Gfpx)
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Gfp -> (Int, Gfpx) -> (Int, Gfpx)
forall a. Integral a => a -> Gfp -> (a, Gfpx) -> (a, Gfpx)
fma (Int
0,Gfpx
zero) (IntMap Gfp -> (Int, Gfpx)) -> IntMap Gfp -> (Int, Gfpx)
forall a b. (a -> b) -> a -> b
$ Gfpx -> IntMap Gfp
coeffMap Gfpx
f
  where
    fma :: a -> Gfp -> (a, Gfpx) -> (a, Gfpx)
fma a
i Gfp
c (a, Gfpx)
z = (a
i, Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p (a -> (a, Gfpx) -> Gfpx
forall a. Integral a => a -> (a, Gfpx) -> Gfpx
align a
i (a, Gfpx)
z) (Gfp -> Gfpx
constant Gfp
c))

    align :: a -> (a, Gfpx) -> Gfpx
align a
i (a
d,Gfpx
z) = if a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then Gfpx
z else Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
multiplyExp Gfp
p Gfpx
z Gfpx
g (a -> Gfp
forall a. Integral a => a -> Gfp
toInteger a
k)
      where k :: a
k = a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
i

-------------------------------------------------------------------------------
-- Division
-------------------------------------------------------------------------------

division :: Prime -> Gfpx -> Gfpx -> (Gfpx,Gfpx)
division :: Gfp -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
division Gfp
p Gfpx
f0 Gfpx
g = if Int
gd Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then String -> (Gfpx, Gfpx)
forall a. HasCallStack => String -> a
error String
"Gfpx.division by zero" else Gfpx -> Gfpx -> (Gfpx, Gfpx)
go Gfpx
zero Gfpx
f0
  where
    go :: Gfpx -> Gfpx -> (Gfpx, Gfpx)
go Gfpx
q Gfpx
f = if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then (Gfpx
q,Gfpx
f) else Gfpx -> Gfpx -> (Gfpx, Gfpx)
go Gfpx
q' Gfpx
f'
      where
        fd :: Int
fd = Gfpx -> Int
degree Gfpx
f
        d :: Int
d = Int
fd Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
gd
        xd :: Gfpx
xd = Int -> Gfp -> Gfpx
monomial Int
d (Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p (Gfpx -> Int -> Gfp
powerCoeff Gfpx
f Int
fd) Gfp
gx)
        -- f - x^d*g = q*g + r ==> f = (q + x^d)*g + r
        q' :: Gfpx
q' = Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p Gfpx
q Gfpx
xd
        f' :: Gfpx
f' = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
xd Gfpx
g)
    gd :: Int
gd = Gfpx -> Int
degree Gfpx
g
    gx :: Gfp
gx = Gfp -> Gfp -> Gfp
Prime.invert Gfp
p (Gfpx -> Int -> Gfp
powerCoeff Gfpx
g Int
gd)

quotient :: Prime -> Gfpx -> Gfpx -> Gfpx
quotient :: Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
f Gfpx
g = (Gfpx, Gfpx) -> Gfpx
forall a b. (a, b) -> a
fst ((Gfpx, Gfpx) -> Gfpx) -> (Gfpx, Gfpx) -> Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
Factor.Gfpx.division Gfp
p Gfpx
f Gfpx
g

remainder :: Prime -> Gfpx -> Gfpx -> Gfpx
remainder :: Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p Gfpx
f Gfpx
g = (Gfpx, Gfpx) -> Gfpx
forall a b. (a, b) -> b
snd ((Gfpx, Gfpx) -> Gfpx) -> (Gfpx, Gfpx) -> Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
Factor.Gfpx.division Gfp
p Gfpx
f Gfpx
g

divides :: Prime -> Gfpx -> Gfpx -> Bool
divides :: Gfp -> Gfpx -> Gfpx -> Bool
divides Gfp
p Gfpx
f Gfpx
g = Gfpx -> Bool
isZero Gfpx
g Bool -> Bool -> Bool
|| (Bool -> Bool
not (Gfpx -> Bool
isZero Gfpx
f) Bool -> Bool -> Bool
&& Gfpx -> Bool
isZero (Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p Gfpx
g Gfpx
f))

properDivisor :: Prime -> Gfpx -> Gfpx -> Bool
properDivisor :: Gfp -> Gfpx -> Gfpx -> Bool
properDivisor Gfp
p Gfpx
f Gfpx
g =
    Gfp -> Gfpx -> Gfpx -> Bool
Factor.Gfpx.divides Gfp
p Gfpx
f Gfpx
g Bool -> Bool -> Bool
&&
    Bool -> Bool
not (Gfpx -> Bool
isConstant Gfpx
f) Bool -> Bool -> Bool
&&
    Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Gfpx -> Int
degree Gfpx
g

-------------------------------------------------------------------------------
-- Every Euclidean domain a admits the definition of a greatest common
-- divisor function
--
--   egcd :: a -> a -> (a,(a,a))
--
-- such that if (g,(s,t)) = egcd x y then:
--
--   1. divides g x
--   2. divides g y
--   3. s*x + t*y == g
-------------------------------------------------------------------------------

egcd :: Prime -> Gfpx -> Gfpx -> (Gfpx,(Gfpx,Gfpx))
egcd :: Gfp -> Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
egcd Gfp
p = Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
go
  where
    go :: Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
go Gfpx
f Gfpx
g | Gfpx -> Bool
isZero Gfpx
g =
        if Gfpx -> Bool
isZero Gfpx
f then (Gfpx
zero,(Gfpx
zero,Gfpx
zero)) else (Gfpx
h, (Gfp -> Gfpx
constant Gfp
x, Gfpx
zero))
      where
        x :: Gfp
x = Gfp -> Gfp -> Gfp
Prime.invert Gfp
p (Gfpx -> Gfp
leadingCoeff Gfpx
f)
        h :: Gfpx
h = Gfp -> Gfp -> Gfpx -> Gfpx
multiplyConstant Gfp
p Gfp
x Gfpx
f
    go Gfpx
f Gfpx
g | Bool
otherwise =
        (Gfpx
h, (Gfpx
t, Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p Gfpx
s (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
q Gfpx
t)))
      where
        (Gfpx
q,Gfpx
r) = Gfp -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
Factor.Gfpx.division Gfp
p Gfpx
f Gfpx
g
        (Gfpx
h,(Gfpx
s,Gfpx
t)) = Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
go Gfpx
g Gfpx
r

gcd :: Prime -> Gfpx -> Gfpx -> Gfpx
gcd :: Gfp -> Gfpx -> Gfpx -> Gfpx
gcd Gfp
p Gfpx
f Gfpx
g = (Gfpx, (Gfpx, Gfpx)) -> Gfpx
forall a b. (a, b) -> a
fst ((Gfpx, (Gfpx, Gfpx)) -> Gfpx) -> (Gfpx, (Gfpx, Gfpx)) -> Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
Factor.Gfpx.egcd Gfp
p Gfpx
f Gfpx
g

chineseRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx -> Gfpx -> Gfpx
chineseRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx -> Gfpx
chineseRemainder Gfp
p Gfpx
f Gfpx
g =
    \Gfpx
x Gfpx
y -> Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
x Gfpx
tg) (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
y Gfpx
sf)) Gfpx
fg
  where
    (Gfpx
_,(Gfpx
s,Gfpx
t)) = Gfp -> Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
Factor.Gfpx.egcd Gfp
p Gfpx
f Gfpx
g
    fg :: Gfpx
fg = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
f Gfpx
g
    sf :: Gfpx
sf = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
s Gfpx
f
    tg :: Gfpx
tg = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
t Gfpx
g

-------------------------------------------------------------------------------
-- Ring operations modulo a nonzero polynomial f
-------------------------------------------------------------------------------

multiplyRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder Gfp
p Gfpx
f Gfpx
g Gfpx
h = Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
g Gfpx
h) Gfpx
f

squareRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx
squareRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx
squareRemainder Gfp
p Gfpx
f Gfpx
g = Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder Gfp
p Gfpx
f Gfpx
g Gfpx
g

multiplyExpRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx -> Integer -> Gfpx
multiplyExpRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfp -> Gfpx
multiplyExpRemainder Gfp
_ Gfpx
_ Gfpx
z Gfpx
_ Gfp
_ | Gfpx -> Bool
isZero Gfpx
z = Gfpx
zero
multiplyExpRemainder Gfp
p Gfpx
f Gfpx
z Gfpx
_ Gfp
0 = Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p Gfpx
z Gfpx
f
multiplyExpRemainder Gfp
_ Gfpx
_ Gfpx
_ Gfpx
g Gfp
_ | Gfpx -> Bool
isZero Gfpx
g = Gfpx
zero
multiplyExpRemainder Gfp
p Gfpx
f Gfpx
z Gfpx
g Gfp
_ | Gfpx -> Bool
isOne Gfpx
g = Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p Gfpx
z Gfpx
f
multiplyExpRemainder Gfp
p Gfpx
f Gfpx
z0 Gfpx
g0 Gfp
k0 = Gfpx -> Gfpx -> Gfp -> Gfpx
forall t. Integral t => Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z0 Gfpx
g0 Gfp
k0
  where
    go :: Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z Gfpx
_ t
0 = Gfpx
z
    go Gfpx
z Gfpx
g t
k = Gfpx -> Gfpx -> t -> Gfpx
go Gfpx
z' Gfpx
g' t
k'
      where
        z' :: Gfpx
z' = if t -> Bool
forall a. Integral a => a -> Bool
even t
k then Gfpx
z else Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder Gfp
p Gfpx
f Gfpx
z Gfpx
g
        g' :: Gfpx
g' = Gfp -> Gfpx -> Gfpx -> Gfpx
squareRemainder Gfp
p Gfpx
f Gfpx
g
        k' :: t
k' = t
k t -> t -> t
forall a. Integral a => a -> a -> a
`div` t
2

expRemainder :: Prime -> Gfpx -> Gfpx -> Integer -> Gfpx
expRemainder :: Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfp -> Gfpx
multiplyExpRemainder Gfp
p Gfpx
f Gfpx
one

composeRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx -> Gfpx
composeRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
composeRemainder Gfp
p Gfpx
f Gfpx
g Gfpx
h = Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfpx
compose Gfp
p Gfpx
g Gfpx
h) Gfpx
f

invertRemainderF :: Prime -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
invertRemainderF :: Gfp -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
invertRemainderF Gfp
_ Gfpx
_ Gfpx
h | Gfpx -> Bool
isZero Gfpx
h = String -> Factor Gfpx Gfpx
forall a. HasCallStack => String -> a
error String
"cannot invert zero polynomial"
invertRemainderF Gfp
p Gfpx
f Gfpx
h =
    if Gfpx -> Bool
isOne Gfpx
g then Gfpx -> Factor Gfpx Gfpx
forall a b. b -> Either a b
Right Gfpx
t else Gfpx -> Factor Gfpx Gfpx
forall a b. a -> Either a b
Left Gfpx
g
  where
    (Gfpx
g,(Gfpx
_,Gfpx
t)) = Gfp -> Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
Factor.Gfpx.egcd Gfp
p Gfpx
f Gfpx
h

invertRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx
invertRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx
invertRemainder Gfp
p Gfpx
f Gfpx
h = Factor Gfpx Gfpx -> Gfpx
forall f a. Show f => Factor f a -> a
runFactor (Factor Gfpx Gfpx -> Gfpx) -> Factor Gfpx Gfpx -> Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
invertRemainderF Gfp
p Gfpx
f Gfpx
h

divideRemainderF :: Prime -> Gfpx -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
divideRemainderF :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
divideRemainderF Gfp
p Gfpx
f Gfpx
g Gfpx
h = do
    Gfpx
i <- Gfp -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
invertRemainderF Gfp
p Gfpx
f Gfpx
h
    Gfpx -> Factor Gfpx Gfpx
forall (m :: * -> *) a. Monad m => a -> m a
return (Gfpx -> Factor Gfpx Gfpx) -> Gfpx -> Factor Gfpx Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder Gfp
p Gfpx
f Gfpx
g Gfpx
i

divideRemainder :: Prime -> Gfpx -> Gfpx -> Gfpx -> Gfpx
divideRemainder :: Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
divideRemainder Gfp
p Gfpx
f Gfpx
g Gfpx
h = Factor Gfpx Gfpx -> Gfpx
forall f a. Show f => Factor f a -> a
runFactor (Factor Gfpx Gfpx -> Gfpx) -> Factor Gfpx Gfpx -> Gfpx
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Gfpx -> Gfpx -> Factor Gfpx Gfpx
divideRemainderF Gfp
p Gfpx
f Gfpx
g Gfpx
h

-------------------------------------------------------------------------------
-- Finding all roots of a polynomial [Briggs1998, sec 4.2]
-------------------------------------------------------------------------------

roots :: Prime -> Gfpx -> [Gfp]
roots :: Gfp -> Gfpx -> [Gfp]
roots Gfp
p | Gfp
p Gfp -> Gfp -> Bool
forall a. Ord a => a -> a -> Bool
< Gfp
3 = \Gfpx
f -> (Gfp -> Bool) -> [Gfp] -> [Gfp]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Gfp
x -> Gfp -> Gfpx -> Gfp -> Gfp
evaluate Gfp
p Gfpx
f Gfp
x Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0) [Gfp
0..(Gfp
pGfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
-Gfp
1)]
roots Gfp
p | Bool
otherwise =
    \Gfpx
f -> if Gfpx -> Bool
isLinear Gfpx
f then Gfpx -> [Gfp]
lin Gfpx
f
          else [Gfp] -> [Gfp]
forall a. Ord a => [a] -> [a]
List.sort (Gfpx -> [Gfp]
go (Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfpx -> Gfpx
xp Gfpx
f)))
  where
    go :: Gfpx -> [Gfp]
go Gfpx
f | Gfpx -> Bool
isLinear Gfpx
f = Gfpx -> [Gfp]
lin Gfpx
f
    go Gfpx
f | Gfpx -> Gfp
constantCoeff Gfpx
f Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
== Gfp
0 = Gfp
0 Gfp -> [Gfp] -> [Gfp]
forall a. a -> [a] -> [a]
: Gfpx -> [Gfp]
go (Int -> Gfpx -> Gfpx
multiplyPower (-Int
1) Gfpx
f)
    go Gfpx
f | Bool
otherwise =
        if Int
0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
d1 Bool -> Bool -> Bool
&& Int
d1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Gfpx -> Int
degree Gfpx
f then Gfpx -> [Gfp]
go Gfpx
f1 [Gfp] -> [Gfp] -> [Gfp]
forall a. [a] -> [a] -> [a]
++ Gfpx -> [Gfp]
go Gfpx
f2
        else (Gfp -> Gfp) -> [Gfp] -> [Gfp]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp -> Gfp -> Gfp -> Gfp
Prime.add Gfp
p Gfp
1) (Gfpx -> [Gfp]
go (Gfp -> Gfpx -> Gfpx -> Gfpx
compose Gfp
p Gfpx
f Gfpx
x1))
      where
        d1 :: Int
d1 = Gfpx -> Int
degree Gfpx
f1
        f1 :: Gfpx
f1 = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfpx -> Gfpx
xp1 Gfpx
f)
        f2 :: Gfpx
f2 = Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
f Gfpx
f1

    lin :: Gfpx -> [Gfp]
lin Gfpx
f | Gfpx -> Bool
isZero Gfpx
f = [Gfp
0..(Gfp
pGfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
-Gfp
1)]
    lin Gfpx
f | Gfpx -> Bool
isConstant Gfpx
f = []
    lin Gfpx
f | Bool
otherwise = [Gfp -> Gfp -> Gfp -> Gfp
Prime.divide Gfp
p (Gfp -> Gfp -> Gfp
Prime.negate Gfp
p Gfp
b) Gfp
a]  -- ax + b = 0
      where
        a :: Gfp
a = Gfpx -> Gfp
linearCoeff Gfpx
f
        b :: Gfp
b = Gfpx -> Gfp
constantCoeff Gfpx
f

    -- x^p - x == product [ (x - i) | 0 <= i < p ]
    xp :: Gfpx -> Gfpx
xp Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
variable Gfp
p) Gfpx
variable

    -- x^p - x == x * (x^((p-1)/2) + 1) * (x^((p-1)/2) - 1)
    xp1 :: Gfpx -> Gfpx
xp1 Gfpx
f = Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
variable ((Gfp
p Gfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
- Gfp
1) Gfp -> Gfp -> Gfp
forall a. Integral a => a -> a -> a
`div` Gfp
2)) Gfpx
one

    -- x + 1
    x1 :: Gfpx
x1 = Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p (Int -> Gfp -> Gfpx
monomial Int
1 Gfp
1) Gfpx
one

totallySplits :: Zx -> Prime -> Maybe [Gfp]
totallySplits :: Zx -> Gfp -> Maybe [Gfp]
totallySplits Zx
f Gfp
p = if [Gfp] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Gfp]
rs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Zx -> Int
Zx.degree Zx
f then [Gfp] -> Maybe [Gfp]
forall a. a -> Maybe a
Just [Gfp]
rs else Maybe [Gfp]
forall a. Maybe a
Nothing
  where rs :: [Gfp]
rs = Gfp -> Gfpx -> [Gfp]
roots Gfp
p (Gfp -> Zx -> Gfpx
fromZx Gfp
p Zx
f)

-------------------------------------------------------------------------------
-- A monic polynomial f of degree d is irreducible in GF(p)[x] if:
--
-- 1. f divides x^(p^d) - x
-- 2. For all prime divisors e of d, gcd (f, x^(p^(d/e)) - x) == 1
-------------------------------------------------------------------------------

irreducible :: Prime -> Gfpx -> Bool
irreducible :: Gfp -> Gfpx -> Bool
irreducible Gfp
p Gfpx
f =
    Gfpx -> Bool
isLinear Gfpx
f Bool -> Bool -> Bool
||
    ((Gfp -> Bool) -> [Gfp] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Gfpx -> Bool
isOne (Gfpx -> Bool) -> (Gfp -> Gfpx) -> Gfp -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfpx -> Gfpx) -> (Gfp -> Gfpx) -> Gfp -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfp -> Gfpx
xpde) [Gfp]
el Bool -> Bool -> Bool
&&
     Gfpx -> Bool
isZero (Gfp -> Gfpx
xpde Gfp
1))
  where
    d :: Gfp
d = Int -> Gfp
forall a. Integral a => a -> Gfp
toInteger (Gfpx -> Int
degree Gfpx
f)
    el :: [Gfp]
el = [Gfp] -> [Gfp]
forall a. [a] -> [a]
reverse (((Gfp, Gfp) -> Gfp) -> [(Gfp, Gfp)] -> [Gfp]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp, Gfp) -> Gfp
forall a b. (a, b) -> a
fst (([(Gfp, Gfp)], Gfp) -> [(Gfp, Gfp)]
forall a b. (a, b) -> a
fst (Gfp -> ([(Gfp, Gfp)], Gfp)
Prime.trialDivision Gfp
d)))

    -- x^(p^(d/e)) - x
    xpde :: Gfp -> Gfpx
xpde Gfp
e = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
variable Gfp
pde) Gfpx
variable
      where pde :: Gfp
pde = Gfp
p Gfp -> Gfp -> Gfp
forall a b. (Num a, Integral b) => a -> b -> a
^ (Gfp
d Gfp -> Gfp -> Gfp
forall a. Integral a => a -> a -> a
`div` Gfp
e)

-------------------------------------------------------------------------------
-- Given a polynomial f in Z[x], use Hensel's Lemma to lift a root r
-- modulo p to a sequence of roots r_k modulo p^k.
-------------------------------------------------------------------------------

liftRoot :: Zx -> Prime -> Gfp -> [(Integer,Integer)]
liftRoot :: Zx -> Gfp -> Gfp -> [(Gfp, Gfp)]
liftRoot Zx
f Gfp
p Gfp
r = Gfp -> Gfp -> [(Gfp, Gfp)]
go Gfp
p Gfp
r
  where
     a :: Gfp
a = Gfp -> Gfp -> Gfp
Prime.invert Gfp
p (Gfp -> Gfpx -> Gfp -> Gfp
evaluate Gfp
p (Gfp -> Gfpx -> Gfpx
derivative Gfp
p (Gfp -> Zx -> Gfpx
fromZx Gfp
p Zx
f)) Gfp
r)
     go :: Gfp -> Gfp -> [(Gfp, Gfp)]
go Gfp
pk Gfp
rk = (Gfp
pk,Gfp
rk) (Gfp, Gfp) -> [(Gfp, Gfp)] -> [(Gfp, Gfp)]
forall a. a -> [a] -> [a]
: Gfp -> Gfp -> [(Gfp, Gfp)]
go Gfp
pk' Gfp
rk'
       where
          pk' :: Gfp
pk' = Gfp
pk Gfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
* Gfp
p
          rk' :: Gfp
rk' = (Gfp
rk Gfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
- Gfp -> Gfpx -> Gfp -> Gfp
evaluate Gfp
pk' (Gfp -> Zx -> Gfpx
fromZx Gfp
pk' Zx
f) Gfp
rk Gfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
* Gfp
a) Gfp -> Gfp -> Gfp
forall a. Integral a => a -> a -> a
`mod` Gfp
pk'

-------------------------------------------------------------------------------
-- The Frobenius endomorphism
-------------------------------------------------------------------------------

frobenius :: Prime -> Gfpx -> Gfpx
frobenius :: Gfp -> Gfpx -> Gfpx
frobenius Gfp
_ Gfpx
f | Gfpx -> Bool
isConstant Gfpx
f = Gfpx
f
frobenius Gfp
p (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
m}) =
    Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = Int
d', coeffMap :: IntMap Gfp
coeffMap = IntMap Gfp
m'}
  where
    d' :: Int
d' = Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p'
    m' :: IntMap Gfp
m' = (Int -> Int) -> IntMap Gfp -> IntMap Gfp
forall a. (Int -> Int) -> IntMap a -> IntMap a
IntMap.mapKeysMonotonic (Int -> Int -> Int
forall a. Num a => a -> a -> a
(*) Int
p') IntMap Gfp
m
    p' :: Int
p' = Gfp -> Int
forall a. Num a => Gfp -> a
fromInteger Gfp
p

frobeniusRange :: Prime -> Gfpx -> Bool
frobeniusRange :: Gfp -> Gfpx -> Bool
frobeniusRange Gfp
p = (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Int -> Bool
pDiv ([Int] -> Bool) -> (Gfpx -> [Int]) -> Gfpx -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IntMap Gfp -> [Int]
forall a. IntMap a -> [Int]
IntMap.keys (IntMap Gfp -> [Int]) -> (Gfpx -> IntMap Gfp) -> Gfpx -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> IntMap Gfp
coeffMap
  where
    p' :: Int
p' = Gfp -> Int
forall a. Num a => Gfp -> a
fromInteger Gfp
p
    pDiv :: Int -> Bool
pDiv Int
i = Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
p' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0

-- Only defined for polynomials satisfying frobeniusRange
frobeniusInverse :: Prime -> Gfpx -> Gfpx
frobeniusInverse :: Gfp -> Gfpx -> Gfpx
frobeniusInverse Gfp
_ Gfpx
f | Gfpx -> Bool
isConstant Gfpx
f = Gfpx
f
frobeniusInverse Gfp
p (Gfpx {degree :: Gfpx -> Int
degree = Int
d, coeffMap :: Gfpx -> IntMap Gfp
coeffMap = IntMap Gfp
m}) =
    Gfpx :: Int -> IntMap Gfp -> Gfpx
Gfpx {degree :: Int
degree = Int
d', coeffMap :: IntMap Gfp
coeffMap = IntMap Gfp
m'}
  where
    d' :: Int
d' = Int
d Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
p'
    m' :: IntMap Gfp
m' = (Int -> Int) -> IntMap Gfp -> IntMap Gfp
forall a. (Int -> Int) -> IntMap a -> IntMap a
IntMap.mapKeysMonotonic Int -> Int
pDiv IntMap Gfp
m
    p' :: Int
p' = Gfp -> Int
forall a. Num a => Gfp -> a
fromInteger Gfp
p
    pDiv :: Int -> Int
pDiv Int
i = if Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
p' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
p'
             else String -> Int
forall a. HasCallStack => String -> a
error (String -> Int) -> String -> Int
forall a b. (a -> b) -> a -> b
$ String
"power does not divide " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Gfp -> String
forall a. Show a => a -> String
show Gfp
p

-------------------------------------------------------------------------------
-- Square-free decomposition
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields
-------------------------------------------------------------------------------

squareFree :: Prime -> Gfpx -> Bool
squareFree :: Gfp -> Gfpx -> Bool
squareFree Gfp
p Gfpx
f =
    if Gfpx -> Bool
isZero Gfpx
f then String -> Bool
forall a. HasCallStack => String -> a
error String
"GF(p)[x] square-free property not defined for zero"
    else Gfpx -> Bool
isLinear Gfpx
f Bool -> Bool -> Bool
|| Gfpx -> Bool
isConstant (Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx
derivative Gfp
p Gfpx
f))

squareFreeDecomposition :: Prime -> Gfpx -> [Gfpx]
squareFreeDecomposition :: Gfp -> Gfpx -> [Gfpx]
squareFreeDecomposition Gfp
_ Gfpx
f | Gfpx -> Bool
isZero Gfpx
f =
    String -> [Gfpx]
forall a. HasCallStack => String -> a
error String
"GF(p)[x] square-free decomposition not defined for zero"
squareFreeDecomposition Gfp
_ Gfpx
f | Gfpx -> Bool
isLinear Gfpx
f =
    [Gfpx
f]
squareFreeDecomposition Gfp
p Gfpx
f | Bool -> Bool
not (Gfpx -> Bool
isMonic Gfpx
f) =
    case Gfp -> Gfpx -> [Gfpx]
squareFreeDecomposition Gfp
p Gfpx
g of
      [] -> String -> [Gfpx]
forall a. HasCallStack => String -> a
error String
"GF(p)[x] square-free decomposition returned an empty list"
      Gfpx
a1 : [Gfpx]
al -> Gfp -> Gfp -> Gfpx -> Gfpx
multiplyConstant Gfp
p Gfp
c Gfpx
a1 Gfpx -> [Gfpx] -> [Gfpx]
forall a. a -> [a] -> [a]
: [Gfpx]
al
  where
    (Gfp
c,Gfpx
g) = Gfp -> Gfpx -> (Gfp, Gfpx)
constMonic Gfp
p Gfpx
f
squareFreeDecomposition Gfp
p Gfpx
f0 =
    IntMap Gfpx -> [Gfpx]
output (IntMap Gfpx -> [Gfpx]) -> IntMap Gfpx -> [Gfpx]
forall a b. (a -> b) -> a -> b
$ Gfpx -> IntMap Gfpx
decompose Gfpx
f0
  where
    decompose :: Gfpx -> IntMap Gfpx
decompose Gfpx
f =
        if Gfpx -> Bool
isOne Gfpx
cp then IntMap Gfpx
m
        else if Gfpx -> Bool
isLinear Gfpx
cp then Int -> Gfpx -> IntMap Gfpx -> IntMap Gfpx
insert Int
p' Gfpx
cp IntMap Gfpx
m
        else IntMap Gfpx -> IntMap Gfpx -> IntMap Gfpx
forall a. IntMap a -> IntMap a -> IntMap a
IntMap.union IntMap Gfpx
m (IntMap Gfpx -> IntMap Gfpx
forall a. IntMap a -> IntMap a
frob (Gfpx -> IntMap Gfpx
decompose Gfpx
cp))
      where
        c :: Gfpx
c = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx
derivative Gfp
p Gfpx
f)
        w :: Gfpx
w = Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
f Gfpx
c
        (Gfpx
cp,IntMap Gfpx
m) = Gfpx -> Gfpx -> Int -> IntMap Gfpx -> (Gfpx, IntMap Gfpx)
loop Gfpx
c Gfpx
w Int
1 IntMap Gfpx
forall a. IntMap a
IntMap.empty

    loop :: Gfpx -> Gfpx -> Int -> IntMap Gfpx -> (Gfpx, IntMap Gfpx)
loop Gfpx
c Gfpx
w Int
_ IntMap Gfpx
m | Gfpx -> Bool
isOne Gfpx
w = (Gfp -> Gfpx -> Gfpx
frobeniusInverse Gfp
p Gfpx
c, IntMap Gfpx
m)
    loop Gfpx
c Gfpx
w Int
k IntMap Gfpx
m = Gfpx -> Gfpx -> Int -> IntMap Gfpx -> (Gfpx, IntMap Gfpx)
loop (Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
c Gfpx
y) Gfpx
y (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Int -> Gfpx -> IntMap Gfpx -> IntMap Gfpx
insert Int
k Gfpx
a IntMap Gfpx
m)
      where
        y :: Gfpx
y = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
w Gfpx
c
        a :: Gfpx
a = Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
w Gfpx
y

    p' :: Int
p' = Gfp -> Int
forall a. Num a => Gfp -> a
fromInteger Gfp
p

    frob :: IntMap a -> IntMap a
frob = (Int -> Int) -> IntMap a -> IntMap a
forall a. (Int -> Int) -> IntMap a -> IntMap a
IntMap.mapKeysMonotonic (Int -> Int -> Int
forall a. Num a => a -> a -> a
(*) Int
p')

    insert :: Int -> Gfpx -> IntMap Gfpx -> IntMap Gfpx
insert Int
_ Gfpx
a IntMap Gfpx
m | Gfpx -> Bool
isOne Gfpx
a = IntMap Gfpx
m
    insert Int
k Gfpx
a IntMap Gfpx
m = Int -> Gfpx -> IntMap Gfpx -> IntMap Gfpx
forall a. Int -> a -> IntMap a -> IntMap a
IntMap.insert Int
k Gfpx
a IntMap Gfpx
m

    output :: IntMap Gfpx -> [Gfpx]
output IntMap Gfpx
m = (Int -> Gfpx) -> [Int] -> [Gfpx]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> Gfpx -> Int -> IntMap Gfpx -> Gfpx
forall a. a -> Int -> IntMap a -> a
IntMap.findWithDefault Gfpx
one Int
k IntMap Gfpx
m) [Int
1..Int
d]
      where (Int
d,Gfpx
_) = IntMap Gfpx -> (Int, Gfpx)
forall a. IntMap a -> (Int, a)
IntMap.findMax IntMap Gfpx
m

squareFreeRecomposition :: Prime -> [Gfpx] -> Gfpx
squareFreeRecomposition :: Gfp -> [Gfpx] -> Gfpx
squareFreeRecomposition Gfp
_ [] = String -> Gfpx
forall a. HasCallStack => String -> a
error String
"GF(p)[x] square-free recomposition of empty list"
squareFreeRecomposition Gfp
p [Gfpx]
al = (Gfpx, Gfpx) -> Gfpx
forall a b. (a, b) -> a
fst ((Gfpx, Gfpx) -> Gfpx) -> (Gfpx, Gfpx) -> Gfpx
forall a b. (a -> b) -> a -> b
$ (Gfpx -> (Gfpx, Gfpx) -> (Gfpx, Gfpx))
-> (Gfpx, Gfpx) -> [Gfpx] -> (Gfpx, Gfpx)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Gfpx -> (Gfpx, Gfpx) -> (Gfpx, Gfpx)
mult (Gfpx
one,Gfpx
one) [Gfpx]
al
  where mult :: Gfpx -> (Gfpx, Gfpx) -> (Gfpx, Gfpx)
mult Gfpx
a (Gfpx
f,Gfpx
g) = (Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
g' Gfpx
f, Gfpx
g') where g' :: Gfpx
g' = Gfp -> Gfpx -> Gfpx -> Gfpx
multiply Gfp
p Gfpx
a Gfpx
g

-------------------------------------------------------------------------------
-- Factorization using Berlekamp's algorithm
--
-- Input:  Polynomial f assumed to be monic and square-free
-- Output: List of irreducible factors of f
--
-- https://en.wikipedia.org/wiki/Berlekamp%27s_algorithm
-------------------------------------------------------------------------------

matrixBerlekamp :: Prime -> Gfpx -> [[Gfp]]
matrixBerlekamp :: Gfp -> Gfpx -> [[Gfp]]
matrixBerlekamp Gfp
p Gfpx
f = [[Gfp]] -> [[Gfp]]
forall a. [[a]] -> [[a]]
List.transpose [[Gfp]]
rows
  where
    rows :: [[Gfp]]
rows = ((Int, Gfpx) -> [Gfp]) -> [(Int, Gfpx)] -> [[Gfp]]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Gfpx) -> [Gfp]
row ([Int] -> [Gfpx] -> [(Int, Gfpx)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int] -> [Int]
forall a. [a] -> [a]
tail [Int]
il) ((Gfpx -> Gfpx) -> Gfpx -> [Gfpx]
forall a. (a -> a) -> a -> [a]
iterate (Gfp -> Gfpx -> Gfpx -> Gfpx -> Gfpx
multiplyRemainder Gfp
p Gfpx
f Gfpx
xp) Gfpx
xp))
    row :: (Int, Gfpx) -> [Gfp]
row (Int
i,Gfpx
xpi) = (Int -> Gfp) -> [Int] -> [Gfp]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Gfpx -> Int -> Gfp
entry Int
i Gfpx
xpi) [Int]
il
    entry :: Int -> Gfpx -> Int -> Gfp
entry Int
i Gfpx
xpi Int
j = if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Gfp -> Gfp -> Gfp -> Gfp
Prime.subtract Gfp
p Gfp
x Gfp
1 else Gfp
x
      where x :: Gfp
x = Gfpx -> Int -> Gfp
powerCoeff Gfpx
xpi Int
j
    il :: [Int]
il = [Int
0..(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)]
    n :: Int
n = Gfpx -> Int
degree Gfpx
f
    xp :: Gfpx
xp = Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
variable Gfp
p

nullBerlekamp :: Prime -> [[Gfp]] -> Maybe [Gfp]
nullBerlekamp :: Gfp -> [[Gfp]] -> Maybe [Gfp]
nullBerlekamp Gfp
p = [[Gfp]] -> Maybe [Gfp]
go
  where
    go :: [[Gfp]] -> Maybe [Gfp]
go ([] : [[Gfp]]
_) = Maybe [Gfp]
forall a. Maybe a
Nothing
    go [[Gfp]]
rows =
        case [(Gfp, [Gfp])]
nr of
          [] -> [Gfp] -> Maybe [Gfp]
forall a. a -> Maybe a
Just (Gfp
1 Gfp -> [Gfp] -> [Gfp]
forall a. a -> [a] -> [a]
: (Gfp -> Gfp) -> [Gfp] -> [Gfp]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp -> Gfp -> Gfp
forall a b. a -> b -> a
const Gfp
0) ((Gfp, [Gfp]) -> [Gfp]
forall a b. (a, b) -> b
snd ([(Gfp, [Gfp])] -> (Gfp, [Gfp])
forall a. [a] -> a
head [(Gfp, [Gfp])]
zr)))
          (Gfp
xh,[Gfp]
xt) : [(Gfp, [Gfp])]
yl ->
              case [[Gfp]] -> Maybe [Gfp]
go ([[Gfp]]
yl' [[Gfp]] -> [[Gfp]] -> [[Gfp]]
forall a. [a] -> [a] -> [a]
++ ((Gfp, [Gfp]) -> [Gfp]) -> [(Gfp, [Gfp])] -> [[Gfp]]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp, [Gfp]) -> [Gfp]
forall a b. (a, b) -> b
snd [(Gfp, [Gfp])]
zr) of
                Maybe [Gfp]
Nothing -> Maybe [Gfp]
forall a. Maybe a
Nothing
                Just [Gfp]
v -> [Gfp] -> Maybe [Gfp]
forall a. a -> Maybe a
Just (Gfp -> Gfp -> Gfp
Prime.negate Gfp
p ([Gfp] -> [Gfp] -> Gfp
dotprod [Gfp]
xt' [Gfp]
v) Gfp -> [Gfp] -> [Gfp]
forall a. a -> [a] -> [a]
: [Gfp]
v)
            where
              xh' :: Gfp
xh' = Gfp -> Gfp -> Gfp
Prime.invert Gfp
p Gfp
xh
              xt' :: [Gfp]
xt' = (Gfp -> Gfp) -> [Gfp] -> [Gfp]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p Gfp
xh') [Gfp]
xt
              yl' :: [[Gfp]]
yl' = ((Gfp, [Gfp]) -> [Gfp]) -> [(Gfp, [Gfp])] -> [[Gfp]]
forall a b. (a -> b) -> [a] -> [b]
map (Gfp, [Gfp]) -> [Gfp]
subx [(Gfp, [Gfp])]
yl
              subx :: (Gfp, [Gfp]) -> [Gfp]
subx (Gfp
yh,[Gfp]
yt) = (Gfp -> Gfp -> Gfp) -> [Gfp] -> [Gfp] -> [Gfp]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Gfp -> Gfp -> Gfp -> Gfp
cancel Gfp
yh) [Gfp]
yt [Gfp]
xt'
      where
        ([(Gfp, [Gfp])]
zr,[(Gfp, [Gfp])]
nr) = ((Gfp, [Gfp]) -> Bool)
-> [(Gfp, [Gfp])] -> ([(Gfp, [Gfp])], [(Gfp, [Gfp])])
forall a. (a -> Bool) -> [a] -> ([a], [a])
List.partition (Gfp -> Gfp -> Bool
forall a. Eq a => a -> a -> Bool
(==) Gfp
0 (Gfp -> Bool) -> ((Gfp, [Gfp]) -> Gfp) -> (Gfp, [Gfp]) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Gfp, [Gfp]) -> Gfp
forall a b. (a, b) -> a
fst) (([Gfp] -> (Gfp, [Gfp])) -> [[Gfp]] -> [(Gfp, [Gfp])]
forall a b. (a -> b) -> [a] -> [b]
map [Gfp] -> (Gfp, [Gfp])
forall a. [a] -> (a, [a])
uncons [[Gfp]]
rows)

    cancel :: Gfp -> Gfp -> Gfp -> Gfp
cancel Gfp
a Gfp
b Gfp
c = Gfp -> Gfp -> Gfp -> Gfp
Prime.subtract Gfp
p Gfp
b (Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p Gfp
a Gfp
c)
    dotprod :: [Gfp] -> [Gfp] -> Gfp
dotprod [Gfp]
u [Gfp]
v = Gfp -> [Gfp] -> Gfp
Prime.sum Gfp
p ((Gfp -> Gfp -> Gfp) -> [Gfp] -> [Gfp] -> [Gfp]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Gfp -> Gfp -> Gfp -> Gfp
Prime.multiply Gfp
p) [Gfp]
u [Gfp]
v)
    uncons :: [a] -> (a, [a])
uncons [a]
l = ([a] -> a
forall a. [a] -> a
head [a]
l, [a] -> [a]
forall a. [a] -> [a]
tail [a]
l)

splitBerlekamp :: Prime -> Gfpx -> Maybe (Gfpx,Gfpx)
splitBerlekamp :: Gfp -> Gfpx -> Maybe (Gfpx, Gfpx)
splitBerlekamp Gfp
_ Gfpx
f | Bool -> Bool
not (Gfpx -> Bool
isMonic Gfpx
f) =
    String -> Maybe (Gfpx, Gfpx)
forall a. HasCallStack => String -> a
error String
"GF(p)[x] Berlekamp factorization requires monic input"
splitBerlekamp Gfp
_ Gfpx
f | Gfpx -> Bool
isLinear Gfpx
f = Maybe (Gfpx, Gfpx)
forall a. Maybe a
Nothing
splitBerlekamp Gfp
p Gfpx
f =
    case Gfp -> [[Gfp]] -> Maybe [Gfp]
nullBerlekamp Gfp
p [[Gfp]]
m of
      Maybe [Gfp]
Nothing -> Maybe (Gfpx, Gfpx)
forall a. Maybe a
Nothing
      Just [Gfp]
v -> (Gfpx, Gfpx) -> Maybe (Gfpx, Gfpx)
forall a. a -> Maybe a
Just (Gfpx
h, Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
f Gfpx
h)
        where
          h :: Gfpx
h = [Gfpx] -> Gfpx
forall a. [a] -> a
head ([Gfpx] -> Gfpx) -> [Gfpx] -> Gfpx
forall a b. (a -> b) -> a -> b
$ (Gfpx -> Bool) -> [Gfpx] -> [Gfpx]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> (Gfpx -> Bool) -> Gfpx -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfpx -> Bool
isOne) ((Gfp -> Gfpx) -> [Gfp] -> [Gfpx]
forall a b. (a -> b) -> [a] -> [b]
map Gfp -> Gfpx
gi [Gfp
0..(Gfp
pGfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
-Gfp
1)])
          -- Property: Product { gi i | 0 <= i < p } = f
          gi :: Gfp -> Gfpx
gi = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfpx -> Gfpx) -> (Gfp -> Gfpx) -> Gfp -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfp -> Gfpx -> Gfpx -> Gfpx
add Gfp
p Gfpx
g (Gfpx -> Gfpx) -> (Gfp -> Gfpx) -> Gfp -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfp -> Gfpx
constant
          -- Property: expRemainder p f g p == g
          g :: Gfpx
g = [Gfp] -> Gfpx
fromCoeff (Gfp
0 Gfp -> [Gfp] -> [Gfp]
forall a. a -> [a] -> [a]
: [Gfp]
v)
  where
    m :: [[Gfp]]
m = Gfp -> Gfpx -> [[Gfp]]
matrixBerlekamp Gfp
p Gfpx
f

-------------------------------------------------------------------------------
-- Distinct degree factorization
--
-- Input:  Polynomial f assumed to be monic and square-free
-- Output: List of (g,d) where g is product of all monic irreducible degree d
--         factors of f
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields
-------------------------------------------------------------------------------

factorDistinctDegree :: Prime -> Gfpx -> [(Gfpx,Int)]
factorDistinctDegree :: Gfp -> Gfpx -> [(Gfpx, Int)]
factorDistinctDegree Gfp
_ Gfpx
f | Gfpx -> Bool
isZero Gfpx
f =
    String -> [(Gfpx, Int)]
forall a. HasCallStack => String -> a
error String
"GF(p)[x] distinct degree factorization not defined for zero"
factorDistinctDegree Gfp
_ Gfpx
f | Bool -> Bool
not (Gfpx -> Bool
isMonic Gfpx
f) =
    String -> [(Gfpx, Int)]
forall a. HasCallStack => String -> a
error String
"GF(p)[x] distinct degree factorization requires monic input"
factorDistinctDegree Gfp
p Gfpx
f0 =
    Int -> Gfpx -> Gfpx -> [(Gfpx, Int)]
go Int
1 Gfpx
f0 Gfpx
variable
  where
    -- Invariant: degree f >= 2*d ==> x == remainder p (x^(p^(d-1))) f
    go :: Int -> Gfpx -> Gfpx -> [(Gfpx, Int)]
go Int
d Gfpx
f Gfpx
_ | Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d = if Gfpx -> Bool
isOne Gfpx
f then [] else [(Gfpx
f, Gfpx -> Int
degree Gfpx
f)]
    go Int
d Gfpx
f Gfpx
x | Bool
otherwise = if Gfpx -> Bool
isOne Gfpx
g then Int -> Gfpx -> Gfpx -> [(Gfpx, Int)]
go Int
d' Gfpx
f Gfpx
x' else (Gfpx
g,Int
d) (Gfpx, Int) -> [(Gfpx, Int)] -> [(Gfpx, Int)]
forall a. a -> [a] -> [a]
: Int -> Gfpx -> Gfpx -> [(Gfpx, Int)]
go Int
d' Gfpx
f' Gfpx
x''
      where
        x' :: Gfpx
x' = Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
x Gfp
p
        g :: Gfpx
g = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
f (Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p Gfpx
x' Gfpx
variable)
        d' :: Int
d' = Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        f' :: Gfpx
f' = Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
f Gfpx
g
        x'' :: Gfpx
x'' = Gfp -> Gfpx -> Gfpx -> Gfpx
remainder Gfp
p Gfpx
x' Gfpx
f'

-------------------------------------------------------------------------------
-- Equal-degree factorization
--
-- Input: Polynomial f assumed to be product of distinct monic irreducible
--        degree d polynomials
-- Output: List of degree d factors of f
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields
-------------------------------------------------------------------------------

factorEqualDegreeBerlekamp :: Prime -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp :: Gfp -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp Gfp
_ Gfpx
f Int
d | Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d = [Gfpx
f]
factorEqualDegreeBerlekamp Gfp
p Gfpx
f Int
d =
    case Gfp -> Gfpx -> Maybe (Gfpx, Gfpx)
splitBerlekamp Gfp
p Gfpx
f of
      Just (Gfpx
g,Gfpx
h) -> Gfp -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp Gfp
p Gfpx
g Int
d [Gfpx] -> [Gfpx] -> [Gfpx]
forall a. [a] -> [a] -> [a]
++
                    Gfp -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp Gfp
p Gfpx
h Int
d
      Maybe (Gfpx, Gfpx)
Nothing -> String -> [Gfpx]
forall a. HasCallStack => String -> a
error String
"Berlekamp could not split polynomial"

factorEqualDegree :: RandomGen r => Prime -> Gfpx -> Int -> r -> ([Gfpx],r)
factorEqualDegree :: Gfp -> Gfpx -> Int -> r -> ([Gfpx], r)
factorEqualDegree Gfp
p Gfpx
f Int
d r
r | Gfp
p Gfp -> Gfp -> Bool
forall a. Ord a => a -> a -> Bool
< Gfp
3 = (Gfp -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp Gfp
p Gfpx
f Int
d, r
r)
factorEqualDegree Gfp
_ Gfpx
f Int
d r
r | Gfpx -> Int
degree Gfpx
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d = ([Gfpx
f],r
r)
factorEqualDegree Gfp
p Gfpx
f Int
d r
r0 = [Gfpx] -> [Gfpx] -> r -> ([Gfpx], r)
forall t. RandomGen t => [Gfpx] -> [Gfpx] -> t -> ([Gfpx], t)
go [] [Gfpx
f] r
r0
  where
    go :: [Gfpx] -> [Gfpx] -> t -> ([Gfpx], t)
go [Gfpx]
dl [] t
r = ([Gfpx]
dl,t
r)
    go [Gfpx]
dl [Gfpx]
hl t
r = [Gfpx] -> [Gfpx] -> t -> ([Gfpx], t)
go ([Gfpx]
dl' [Gfpx] -> [Gfpx] -> [Gfpx]
forall a. [a] -> [a] -> [a]
++ [Gfpx]
dl) [Gfpx]
hl' t
r'
      where
        (Gfpx
u,t
r') = Gfp -> Int -> t -> (Gfpx, t)
forall r. RandomGen r => Gfp -> Int -> r -> (Gfpx, r)
uniform Gfp
p Int
n t
r
        v :: Gfpx
v = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.subtract Gfp
p (Gfp -> Gfpx -> Gfpx -> Gfp -> Gfpx
expRemainder Gfp
p Gfpx
f Gfpx
u Gfp
k) Gfpx
one
        ([Gfpx]
dl',[Gfpx]
hl') = (Gfpx -> Bool) -> [Gfpx] -> ([Gfpx], [Gfpx])
forall a. (a -> Bool) -> [a] -> ([a], [a])
List.partition Gfpx -> Bool
degreeD ([Gfpx] -> ([Gfpx], [Gfpx])) -> [Gfpx] -> ([Gfpx], [Gfpx])
forall a b. (a -> b) -> a -> b
$ (Gfpx -> [Gfpx]) -> [Gfpx] -> [Gfpx]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Gfpx -> Gfpx -> [Gfpx]
split Gfpx
v) [Gfpx]
hl

    split :: Gfpx -> Gfpx -> [Gfpx]
split Gfpx
v Gfpx
h =
        if Gfpx -> Bool
isOne Gfpx
g Bool -> Bool -> Bool
|| Gfpx -> Int
degree Gfpx
g Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Gfpx -> Int
degree Gfpx
h then [Gfpx
h] else [Gfpx
g, Gfp -> Gfpx -> Gfpx -> Gfpx
quotient Gfp
p Gfpx
h Gfpx
g]
      where
        g :: Gfpx
g = Gfp -> Gfpx -> Gfpx -> Gfpx
Factor.Gfpx.gcd Gfp
p Gfpx
v Gfpx
h

    degreeD :: Gfpx -> Bool
degreeD Gfpx
h = Gfpx -> Int
degree Gfpx
h Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d

    n :: Int
n = Gfpx -> Int
degree Gfpx
f
    k :: Gfp
k = (Gfp
p Gfp -> Int -> Gfp
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
d Gfp -> Gfp -> Gfp
forall a. Num a => a -> a -> a
- Gfp
1) Gfp -> Gfp -> Gfp
forall a. Integral a => a -> a -> a
`div` Gfp
2

-------------------------------------------------------------------------------
-- Square-free factorization
--
-- Input:  Polynomial f assumed to be monic and square-free
-- Output: List of monic irreducible factors of f
-------------------------------------------------------------------------------

factorSquareFreeBerlekamp :: Prime -> Gfpx -> [Gfpx]
factorSquareFreeBerlekamp :: Gfp -> Gfpx -> [Gfpx]
factorSquareFreeBerlekamp Gfp
_ Gfpx
f | Gfpx -> Bool
isOne Gfpx
f = []
factorSquareFreeBerlekamp Gfp
_ Gfpx
f | Gfpx -> Bool
isLinear Gfpx
f = [Gfpx
f]
factorSquareFreeBerlekamp Gfp
p Gfpx
f = ((Gfpx, Int) -> [Gfpx]) -> [(Gfpx, Int)] -> [Gfpx]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Gfpx, Int) -> [Gfpx]
fed (Gfp -> Gfpx -> [(Gfpx, Int)]
factorDistinctDegree Gfp
p Gfpx
f)
  where fed :: (Gfpx, Int) -> [Gfpx]
fed = (Gfpx -> Int -> [Gfpx]) -> (Gfpx, Int) -> [Gfpx]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Gfpx -> Int -> [Gfpx]) -> (Gfpx, Int) -> [Gfpx])
-> (Gfpx -> Int -> [Gfpx]) -> (Gfpx, Int) -> [Gfpx]
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> Int -> [Gfpx]
factorEqualDegreeBerlekamp Gfp
p

factorSquareFree :: RandomGen r => Prime -> Gfpx -> r -> ([Gfpx],r)
factorSquareFree :: Gfp -> Gfpx -> r -> ([Gfpx], r)
factorSquareFree Gfp
_ Gfpx
f r
r | Gfpx -> Bool
isOne Gfpx
f = ([],r
r)
factorSquareFree Gfp
_ Gfpx
f r
r | Gfpx -> Bool
isLinear Gfpx
f = ([Gfpx
f],r
r)
factorSquareFree Gfp
p Gfpx
f r
r0 = ((Gfpx, Int) -> ([Gfpx], r) -> ([Gfpx], r))
-> ([Gfpx], r) -> [(Gfpx, Int)] -> ([Gfpx], r)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Gfpx, Int) -> ([Gfpx], r) -> ([Gfpx], r)
forall b. RandomGen b => (Gfpx, Int) -> ([Gfpx], b) -> ([Gfpx], b)
fed ([],r
r0) (Gfp -> Gfpx -> [(Gfpx, Int)]
factorDistinctDegree Gfp
p Gfpx
f)
  where
    fed :: (Gfpx, Int) -> ([Gfpx], b) -> ([Gfpx], b)
fed (Gfpx
h,Int
d) ([Gfpx]
gl,b
r) = ([Gfpx]
gl' [Gfpx] -> [Gfpx] -> [Gfpx]
forall a. [a] -> [a] -> [a]
++ [Gfpx]
gl, b
r')
      where ([Gfpx]
gl',b
r') = Gfp -> Gfpx -> Int -> b -> ([Gfpx], b)
forall r. RandomGen r => Gfp -> Gfpx -> Int -> r -> ([Gfpx], r)
factorEqualDegree Gfp
p Gfpx
h Int
d b
r

-------------------------------------------------------------------------------
-- Factorization of monic polynomials
--
-- Input:  Monic polynomial f
-- Output: List of monic irreducible polynomials g_i and positive integers k_i
--         such that f = product_i g_i^k_i
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields
-------------------------------------------------------------------------------

factorMonicBerlekamp :: Prime -> Gfpx -> [(Gfpx,Integer)]
factorMonicBerlekamp :: Gfp -> Gfpx -> [(Gfpx, Gfp)]
factorMonicBerlekamp Gfp
p = [[(Gfpx, Gfp)]] -> [(Gfpx, Gfp)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[(Gfpx, Gfp)]] -> [(Gfpx, Gfp)])
-> (Gfpx -> [[(Gfpx, Gfp)]]) -> Gfpx -> [(Gfpx, Gfp)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Gfp -> Gfpx -> [(Gfpx, Gfp)])
-> [Gfp] -> [Gfpx] -> [[(Gfpx, Gfp)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Gfp -> Gfpx -> [(Gfpx, Gfp)]
forall b. b -> Gfpx -> [(Gfpx, b)]
fsf [Gfp
1..] ([Gfpx] -> [[(Gfpx, Gfp)]])
-> (Gfpx -> [Gfpx]) -> Gfpx -> [[(Gfpx, Gfp)]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Gfp -> Gfpx -> [Gfpx]
squareFreeDecomposition Gfp
p
  where fsf :: b -> Gfpx -> [(Gfpx, b)]
fsf b
k Gfpx
h = (Gfpx -> (Gfpx, b)) -> [Gfpx] -> [(Gfpx, b)]
forall a b. (a -> b) -> [a] -> [b]
map (\Gfpx
g -> (Gfpx
g,b
k)) ([Gfpx] -> [(Gfpx, b)]) -> [Gfpx] -> [(Gfpx, b)]
forall a b. (a -> b) -> a -> b
$ Gfp -> Gfpx -> [Gfpx]
factorSquareFreeBerlekamp Gfp
p Gfpx
h

factorMonic :: RandomGen r => Prime -> Gfpx -> r -> ([(Gfpx,Integer)],r)
factorMonic :: Gfp -> Gfpx -> r -> ([(Gfpx, Gfp)], r)
factorMonic Gfp
p Gfpx
f r
r0 = ((Gfp, Gfpx) -> ([(Gfpx, Gfp)], r) -> ([(Gfpx, Gfp)], r))
-> ([(Gfpx, Gfp)], r) -> [(Gfp, Gfpx)] -> ([(Gfpx, Gfp)], r)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Gfp, Gfpx) -> ([(Gfpx, Gfp)], r) -> ([(Gfpx, Gfp)], r)
forall b b.
RandomGen b =>
(b, Gfpx) -> ([(Gfpx, b)], b) -> ([(Gfpx, b)], b)
fsf ([],r
r0) ([Gfp] -> [Gfpx] -> [(Gfp, Gfpx)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Gfp
1..] (Gfp -> Gfpx -> [Gfpx]
squareFreeDecomposition Gfp
p Gfpx
f))
  where
    fsf :: (b, Gfpx) -> ([(Gfpx, b)], b) -> ([(Gfpx, b)], b)
fsf (b
k,Gfpx
h) ([(Gfpx, b)]
gkl,b
r) = ((Gfpx -> (Gfpx, b)) -> [Gfpx] -> [(Gfpx, b)]
forall a b. (a -> b) -> [a] -> [b]
map (\Gfpx
g -> (Gfpx
g,b
k)) [Gfpx]
gl [(Gfpx, b)] -> [(Gfpx, b)] -> [(Gfpx, b)]
forall a. [a] -> [a] -> [a]
++ [(Gfpx, b)]
gkl, b
r')
      where ([Gfpx]
gl,b
r') = Gfp -> Gfpx -> b -> ([Gfpx], b)
forall r. RandomGen r => Gfp -> Gfpx -> r -> ([Gfpx], r)
factorSquareFree Gfp
p Gfpx
h b
r