{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -Wall -fno-warn-orphans #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.Polynomial.Factorization.FiniteField
-- Copyright   :  (c) Masahiro Sakai 2012-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Factoriation of polynomial over a finite field.
--
-- References:
--
-- * <http://en.wikipedia.org/wiki/Factorization_of_polynomials_over_a_finite_field_and_irreducibility_tests>
--
-- * <http://en.wikipedia.org/wiki/Berlekamp%27s_algorithm>
--
-- * Martin Kreuzer and Lorenzo Robbiano. Computational Commutative Algebra 1. Springer Verlag, 2000.
--
-----------------------------------------------------------------------------
module ToySolver.Data.Polynomial.Factorization.FiniteField
  ( factor
  , sqfree
  , berlekamp
  , basisOfBerlekampSubalgebra
  ) where

import Control.Exception (assert)
import Data.FiniteField
import Data.List
import Data.Ord
import Data.Set (Set)
import qualified Data.Set as Set
import GHC.TypeLits

import ToySolver.Data.Polynomial.Base (Polynomial, UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial.Base as P
import qualified ToySolver.Data.Polynomial.GroebnerBasis as GB

instance KnownNat p => P.Factor (UPolynomial (PrimeField p)) where
  factor :: UPolynomial (PrimeField p)
-> [(UPolynomial (PrimeField p), Integer)]
factor = forall k.
(Ord k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
factor

instance KnownNat p => P.SQFree (UPolynomial (PrimeField p)) where
  sqfree :: UPolynomial (PrimeField p)
-> [(UPolynomial (PrimeField p), Integer)]
sqfree = forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree

factor :: forall k. (Ord k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)]
factor :: forall k.
(Ord k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
factor UPolynomial k
f = do
  (UPolynomial k
g,Integer
n) <- forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree UPolynomial k
f
  if forall t. Degree t => t -> Integer
P.deg UPolynomial k
g forall a. Ord a => a -> a -> Bool
> Integer
0
    then do
      UPolynomial k
h <- forall k.
(Eq k, Ord k, FiniteField k) =>
UPolynomial k -> [UPolynomial k]
berlekamp UPolynomial k
g
      forall (m :: * -> *) a. Monad m => a -> m a
return (UPolynomial k
h,Integer
n)
    else do
      forall (m :: * -> *) a. Monad m => a -> m a
return (UPolynomial k
g,Integer
n)

-- | Square-free decomposition of univariate polynomials over a finite field.
sqfree :: forall k. (Eq k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)]
sqfree :: forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree UPolynomial k
f
  | k
c forall a. Eq a => a -> a -> Bool
== k
1    = forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' UPolynomial k
f
  | Bool
otherwise = (forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
c, Integer
1) forall a. a -> [a] -> [a]
: forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' (forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff (forall a. Fractional a => a -> a -> a
/k
c) UPolynomial k
f)
  where
    c :: k
c = forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat UPolynomial k
f

sqfree' :: forall k. (Eq k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' :: forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' UPolynomial k
0 = []
sqfree' UPolynomial k
f
  | UPolynomial k
g forall a. Eq a => a -> a -> Bool
== UPolynomial k
0    = [(UPolynomial k
h, Integer
nforall a. Num a => a -> a -> a
*Integer
p) | (UPolynomial k
h,Integer
n) <- forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' (forall k. (Eq k, FiniteField k) => UPolynomial k -> UPolynomial k
polyPthRoot UPolynomial k
f)]
  | Bool
otherwise = forall {k}.
(FiniteField k, Eq k) =>
Integer
-> UPolynomial k
-> UPolynomial k
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go Integer
1 UPolynomial k
c0 UPolynomial k
w0 []
  where
    p :: Integer
p = forall k. FiniteField k => k -> Integer
char (forall a. HasCallStack => a
undefined :: k)
    g :: UPolynomial k
g = forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> v -> Polynomial k v
P.deriv UPolynomial k
f X
X
    c0 :: UPolynomial k
c0 = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.gcd UPolynomial k
f UPolynomial k
g
    w0 :: UPolynomial k
w0 = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.div UPolynomial k
f UPolynomial k
c0
    go :: Integer
-> UPolynomial k
-> UPolynomial k
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go !Integer
i UPolynomial k
c UPolynomial k
w ![(UPolynomial k, Integer)]
result
      | UPolynomial k
w forall a. Eq a => a -> a -> Bool
== UPolynomial k
1    =
          if UPolynomial k
c forall a. Eq a => a -> a -> Bool
== UPolynomial k
1
          then [(UPolynomial k, Integer)]
result
          else [(UPolynomial k, Integer)]
result forall a. [a] -> [a] -> [a]
++ [(UPolynomial k
h, Integer
nforall a. Num a => a -> a -> a
*Integer
p) | (UPolynomial k
h,Integer
n) <- forall k.
(Eq k, FiniteField k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfree' (forall k. (Eq k, FiniteField k) => UPolynomial k -> UPolynomial k
polyPthRoot UPolynomial k
c)]
      | Bool
otherwise = Integer
-> UPolynomial k
-> UPolynomial k
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go (Integer
iforall a. Num a => a -> a -> a
+Integer
1) UPolynomial k
c' UPolynomial k
w' [(UPolynomial k, Integer)]
result'
          where
            y :: UPolynomial k
y  = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.gcd UPolynomial k
w UPolynomial k
c
            z :: UPolynomial k
z  = UPolynomial k
w forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` UPolynomial k
y
            c' :: UPolynomial k
c' = UPolynomial k
c forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` UPolynomial k
y
            w' :: UPolynomial k
w' = UPolynomial k
y
            result' :: [(UPolynomial k, Integer)]
result' = [(UPolynomial k
z,Integer
i) | UPolynomial k
z forall a. Eq a => a -> a -> Bool
/= UPolynomial k
1] forall a. [a] -> [a] -> [a]
++ [(UPolynomial k, Integer)]
result

polyPthRoot :: forall k. (Eq k, FiniteField k) => UPolynomial k -> UPolynomial k
polyPthRoot :: forall k. (Eq k, FiniteField k) => UPolynomial k -> UPolynomial k
polyPthRoot Polynomial k X
f = forall a. HasCallStack => Bool -> a -> a
assert (forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> v -> Polynomial k v
P.deriv Polynomial k X
f X
X forall a. Eq a => a -> a -> Bool
== Polynomial k X
0) forall a b. (a -> b) -> a -> b
$
  forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [(forall k. FiniteField k => k -> k
pthRoot k
c, forall {t}. Degree t => t -> Monomial X
g Monomial X
mm) | (k
c,Monomial X
mm) <- forall k v. Polynomial k v -> [Term k v]
P.terms Polynomial k X
f]
  where
    p :: Integer
p = forall k. FiniteField k => k -> Integer
char (forall a. HasCallStack => a
undefined :: k)
    g :: t -> Monomial X
g t
mm = forall a v. Var a v => v -> a
P.var X
X forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (forall t. Degree t => t -> Integer
P.deg t
mm forall a. Integral a => a -> a -> a
`div` Integer
p)

-- | Berlekamp algorithm for polynomial factorization.
--
-- Input polynomial is assumed to be monic and square-free.
berlekamp :: forall k. (Eq k, Ord k, FiniteField k) => UPolynomial k -> [UPolynomial k]
berlekamp :: forall k.
(Eq k, Ord k, FiniteField k) =>
UPolynomial k -> [UPolynomial k]
berlekamp UPolynomial k
f = Set (UPolynomial k) -> [UPolynomial k] -> [UPolynomial k]
go (forall a. a -> Set a
Set.singleton UPolynomial k
f) [UPolynomial k]
basis
  where
    go :: Set (UPolynomial k) -> [UPolynomial k] -> [UPolynomial k]
    go :: Set (UPolynomial k) -> [UPolynomial k] -> [UPolynomial k]
go Set (UPolynomial k)
_ [] = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$ [Char]
"ToySolver.Data.Polynomial.Factorization.FiniteField.berlekamp: should not happen"
    go Set (UPolynomial k)
fs (UPolynomial k
b:[UPolynomial k]
bs)
      | forall a. Set a -> Int
Set.size Set (UPolynomial k)
fs forall a. Eq a => a -> a -> Bool
== Int
r = forall a. Set a -> [a]
Set.toList Set (UPolynomial k)
fs
      | Bool
otherwise = Set (UPolynomial k) -> [UPolynomial k] -> [UPolynomial k]
go (forall (f :: * -> *) a. (Foldable f, Ord a) => f (Set a) -> Set a
Set.unions [UPolynomial k -> Set (UPolynomial k)
func UPolynomial k
fi | UPolynomial k
fi <- forall a. Set a -> [a]
Set.toList Set (UPolynomial k)
fs]) [UPolynomial k]
bs
        where
          func :: UPolynomial k -> Set (UPolynomial k)
func UPolynomial k
fi = forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$ [UPolynomial k]
hs2 forall a. [a] -> [a] -> [a]
++ [UPolynomial k]
hs1
            where
              hs1 :: [UPolynomial k]
hs1 = [UPolynomial k
h | k
k <- forall k. FiniteField k => [k]
allValues, let h :: UPolynomial k
h = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.gcd UPolynomial k
fi (UPolynomial k
b forall a. Num a => a -> a -> a
- forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
k), forall t. Degree t => t -> Integer
P.deg UPolynomial k
h forall a. Ord a => a -> a -> Bool
> Integer
0]
              hs2 :: [UPolynomial k]
hs2 = if forall t. Degree t => t -> Integer
P.deg UPolynomial k
g forall a. Ord a => a -> a -> Bool
> Integer
0 then [UPolynomial k
g] else []
                where
                  g :: UPolynomial k
g = UPolynomial k
fi forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [UPolynomial k]
hs1
    basis :: [UPolynomial k]
basis = forall k.
(Ord k, FiniteField k) =>
UPolynomial k -> [UPolynomial k]
basisOfBerlekampSubalgebra UPolynomial k
f
    r :: Int
r     = forall (t :: * -> *) a. Foldable t => t a -> Int
length [UPolynomial k]
basis

basisOfBerlekampSubalgebra :: forall k. (Ord k, FiniteField k) => UPolynomial k -> [UPolynomial k]
basisOfBerlekampSubalgebra :: forall k.
(Ord k, FiniteField k) =>
UPolynomial k -> [UPolynomial k]
basisOfBerlekampSubalgebra UPolynomial k
f =
  forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (forall a. a -> Down a
Down forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Degree t => t -> Integer
P.deg) forall a b. (a -> b) -> a -> b
$
    forall a b. (a -> b) -> [a] -> [b]
map (forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
P.toMonic MonomialOrder X
P.nat) forall a b. (a -> b) -> a -> b
$
      [UPolynomial k]
basis
  where
    q :: Integer
q    = forall k. FiniteField k => k -> Integer
order (forall a. HasCallStack => a
undefined :: k)
    d :: Integer
d    = forall t. Degree t => t -> Integer
P.deg UPolynomial k
f
    x :: UPolynomial k
x    = forall a v. Var a v => v -> a
P.var X
X

    qs :: [UPolynomial k]
    qs :: [UPolynomial k]
qs = [(UPolynomial k
xforall a b. (Num a, Integral b) => a -> b -> a
^(Integer
qforall a. Num a => a -> a -> a
*Integer
i)) forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.mod` UPolynomial k
f | Integer
i <- [Integer
0 .. Integer
d forall a. Num a => a -> a -> a
- Integer
1]]

    gb :: [Polynomial k Int]
    gb :: [Polynomial k Int]
gb = forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
GB.basis forall v. Ord v => MonomialOrder v
P.grlex [Polynomial k Int
p3 | (Polynomial k Int
p3,Monomial X
_) <- forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial (Polynomial k Int)
p2]

    p1 :: Polynomial k Int
    p1 :: Polynomial k Int
p1 = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall a v. Var a v => v -> a
P.var Int
i forall a. Num a => a -> a -> a
* (forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
qi (\X
X -> forall a v. Var a v => v -> a
P.var (-Int
1)) forall a. Num a => a -> a -> a
- (forall a v. Var a v => v -> a
P.var (-Int
1) forall a b. (Num a, Integral b) => a -> b -> a
^ Int
i)) | (Int
i, UPolynomial k
qi) <- forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] [UPolynomial k]
qs]
    p2 :: UPolynomial (Polynomial k Int)
    p2 :: UPolynomial (Polynomial k Int)
p2 = forall k v.
(Ord k, Num k, Ord v) =>
Polynomial k v -> v -> UPolynomial (Polynomial k v)
P.toUPolynomialOf Polynomial k Int
p1 (-Int
1)

    es :: [(Int, Polynomial k Int)]
es  = [(Int
i, forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce forall v. Ord v => MonomialOrder v
P.grlex (forall a v. Var a v => v -> a
P.var Int
i) [Polynomial k Int]
gb) | Int
i <- [Int
0 .. forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
d forall a. Num a => a -> a -> a
- Int
1]]
    vs1 :: [Int]
vs1 = [Int
i           | (Int
i, Polynomial k Int
gi_def) <- [(Int, Polynomial k Int)]
es, Polynomial k Int
gi_def forall a. Eq a => a -> a -> Bool
== forall a v. Var a v => v -> a
P.var Int
i]
    vs2 :: [(Int, Polynomial k Int)]
vs2 = [(Int
i, Polynomial k Int
gi_def) | (Int
i, Polynomial k Int
gi_def) <- [(Int, Polynomial k Int)]
es, Polynomial k Int
gi_def forall a. Eq a => a -> a -> Bool
/= forall a v. Var a v => v -> a
P.var Int
i]

    basis :: [UPolynomial k]
    basis :: [UPolynomial k]
basis = [ UPolynomial k
xforall a b. (Num a, Integral b) => a -> b -> a
^Int
i forall a. Num a => a -> a -> a
+ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (forall {a} {a}. (Eq a, Num a) => a -> a -> a
delta Int
i) Polynomial k Int
gj_def) forall a. Num a => a -> a -> a
* UPolynomial k
xforall a b. (Num a, Integral b) => a -> b -> a
^Int
j | (Int
j, Polynomial k Int
gj_def) <- [(Int, Polynomial k Int)]
vs2] | Int
i <- [Int]
vs1 ]
      where
        delta :: a -> a -> a
delta a
i a
k
          | a
kforall a. Eq a => a -> a -> Bool
==a
i      = a
1
          | Bool
otherwise = a
0