{- |
module: Factor.Bz
description: The Berlekamp-Zassenhaus algorithm for factoring in Z[x]
license: MIT

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

import qualified Data.List as List

import Factor.Gfpx (Gfpx)
import qualified Factor.Gfpx as Gfpx
import Factor.Prime (Prime)
import qualified Factor.Prime as Prime
import Factor.Util
import Factor.Zx (Zx)
import qualified Factor.Zx as Zx

-------------------------------------------------------------------------------
-- The Landau-Mignotte bound on the absolute value of factor coefficients
--
-- https://www.math.fsu.edu/~hoeij/papers/issac11/A.pdf
-------------------------------------------------------------------------------

factorCoeffBound :: Zx -> Integer
factorCoeffBound :: Zx -> Integer
factorCoeffBound Zx
f | Zx -> Bool
Zx.isZero Zx
f =
    [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"Z[x] factor coefficient bound is undefined for zero"
factorCoeffBound Zx
f =
    Integer -> Integer -> Integer
nthRoot Integer
2 (Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (Integer
2 Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
d) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
m
  where
    d :: Integer
d = Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Zx -> Int
Zx.degree Zx
f)
    m :: Integer
m = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (((Int, Integer) -> Integer) -> [(Int, Integer)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer)
-> ((Int, Integer) -> Integer) -> (Int, Integer) -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Integer) -> Integer
forall a b. (a, b) -> b
snd) (Zx -> [(Int, Integer)]
Zx.monomials Zx
f))

-------------------------------------------------------------------------------
-- Finding factors modulo a suitable prime p
-------------------------------------------------------------------------------

monicGfpx :: Prime -> Zx -> Gfpx
monicGfpx :: Integer -> Zx -> Gfpx
monicGfpx Integer
p = (Integer, Gfpx) -> Gfpx
forall a b. (a, b) -> b
snd ((Integer, Gfpx) -> Gfpx) -> (Zx -> (Integer, Gfpx)) -> Zx -> Gfpx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Gfpx -> (Integer, Gfpx)
Gfpx.constMonic Integer
p (Gfpx -> (Integer, Gfpx)) -> (Zx -> Gfpx) -> Zx -> (Integer, Gfpx)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Zx -> Gfpx
Gfpx.fromZx Integer
p

suitablePrime :: Zx -> Prime
suitablePrime :: Zx -> Integer
suitablePrime Zx
f = [Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter Integer -> Bool
suitable [Integer]
Prime.primes
  where
    suitable :: Integer -> Bool
suitable Integer
p = Bool -> Bool
not (Integer -> Integer -> Bool
divides Integer
p Integer
c) Bool -> Bool -> Bool
&& Integer -> Gfpx -> Bool
Gfpx.squareFree Integer
p (Integer -> Zx -> Gfpx
monicGfpx Integer
p Zx
f)
    c :: Integer
c = Zx -> Integer
Zx.leadingCoeff Zx
f

-------------------------------------------------------------------------------
-- Hensel lifting factors modulo p^k
--
-- https://www.csd.uwo.ca/~mmorenom/CS874/Lectures/Newton2Hensel.html/node17.html
-------------------------------------------------------------------------------

henselLiftQuadratic :: Zx -> (Integer,(Gfpx,Gfpx,Gfpx,Gfpx)) ->
                       (Integer,(Gfpx,Gfpx,Gfpx,Gfpx))
henselLiftQuadratic :: Zx
-> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
-> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
henselLiftQuadratic Zx
f (Integer
m,(Gfpx
g,Gfpx
h,Gfpx
s,Gfpx
t)) = (Integer
m',(Gfpx
g',Gfpx
h',Gfpx
s',Gfpx
t'))
  where
    m' :: Integer
m' = Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
m
    e :: Gfpx
e = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.subtract Integer
m' (Integer -> Zx -> Gfpx
Gfpx.fromZx Integer
m' Zx
f) (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
g Gfpx
h)
    (Gfpx
q,Gfpx
r) = Integer -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
Gfpx.division Integer
m' (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
s Gfpx
e) Gfpx
h
    g' :: Gfpx
g' = Integer -> [Gfpx] -> Gfpx
Gfpx.sum Integer
m' [Gfpx
g, Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
t Gfpx
e, Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
q Gfpx
g]
    h' :: Gfpx
h' = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.add Integer
m' Gfpx
h Gfpx
r
    b :: Gfpx
b = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.subtract Integer
m'
          (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.add Integer
m' (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
s Gfpx
g') (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
t Gfpx
h'))
          Gfpx
Gfpx.one
    (Gfpx
c,Gfpx
d) = Integer -> Gfpx -> Gfpx -> (Gfpx, Gfpx)
Gfpx.division Integer
m' (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
s Gfpx
b) Gfpx
h'
    s' :: Gfpx
s' = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.subtract Integer
m' Gfpx
s Gfpx
d
    t' :: Gfpx
t' = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.subtract Integer
m'
           (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
t (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.subtract Integer
m' Gfpx
Gfpx.one Gfpx
b))
           (Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
m' Gfpx
c Gfpx
g')

henselLiftModulus :: Zx -> Prime -> (Int,Integer)
henselLiftModulus :: Zx -> Integer -> (Int, Integer)
henselLiftModulus Zx
f Integer
p =
    [(Int, Integer)] -> (Int, Integer)
forall a. [a] -> a
head ([(Int, Integer)] -> (Int, Integer))
-> [(Int, Integer)] -> (Int, Integer)
forall a b. (a -> b) -> a -> b
$ ((Int, Integer) -> Bool) -> [(Int, Integer)] -> [(Int, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
List.dropWhile (\(Int
_,Integer
pk) -> Integer
pk Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
b) ([(Int, Integer)] -> [(Int, Integer)])
-> [(Int, Integer)] -> [(Int, Integer)]
forall a b. (a -> b) -> a -> b
$
    ((Int, Integer) -> (Int, Integer))
-> (Int, Integer) -> [(Int, Integer)]
forall a. (a -> a) -> a -> [a]
iterate (\(Int
k,Integer
pk) -> (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Integer
pk Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
pk)) (Int
0,Integer
p)
  where
    b :: Integer
b = Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Zx -> Integer
factorCoeffBound Zx
f Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1

henselLiftFactors :: Zx -> Prime -> Int -> [Gfpx] -> [Gfpx]
henselLiftFactors :: Zx -> Integer -> Int -> [Gfpx] -> [Gfpx]
henselLiftFactors Zx
f0 Integer
p Int
k = [Gfpx] -> Zx -> [Gfpx] -> [Gfpx]
go [] Zx
f0
  where
    go :: [Gfpx] -> Zx -> [Gfpx] -> [Gfpx]
go [Gfpx]
_ Zx
_ [] = [Char] -> [Gfpx]
forall a. HasCallStack => [Char] -> a
error [Char]
"no factors to Hensel lift"
    go [Gfpx]
_ Zx
_ [Gfpx
_] = [Char] -> [Gfpx]
forall a. HasCallStack => [Char] -> a
error [Char]
"single factor to Hensel lift"
    go [Gfpx]
acc Zx
f (Gfpx
g : [Gfpx]
gs) =
        if [Gfpx] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Gfpx]
gs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 then [Gfpx] -> [Gfpx]
forall a. [a] -> [a]
reverse [Gfpx]
acc [Gfpx] -> [Gfpx] -> [Gfpx]
forall a. [a] -> [a] -> [a]
++ [Gfpx
gk,Gfpx
hk]
        else [Gfpx] -> Zx -> [Gfpx] -> [Gfpx]
go (Gfpx
gk Gfpx -> [Gfpx] -> [Gfpx]
forall a. a -> [a] -> [a]
: [Gfpx]
acc) (Gfpx -> Zx
Gfpx.toZx Gfpx
hk) [Gfpx]
gs
      where
        h :: Gfpx
h = Integer -> [Gfpx] -> Gfpx
Gfpx.product Integer
p [Gfpx]
gs
        (Gfpx
_,(Gfpx
s,Gfpx
t)) = Integer -> Gfpx -> Gfpx -> (Gfpx, (Gfpx, Gfpx))
Gfpx.egcd Integer
p Gfpx
g Gfpx
h
        (Integer
_,(Gfpx
gk,Gfpx
hk,Gfpx
_,Gfpx
_)) = ((Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
 -> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx)))
-> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
-> [(Integer, (Gfpx, Gfpx, Gfpx, Gfpx))]
forall a. (a -> a) -> a -> [a]
iterate (Zx
-> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
-> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
henselLiftQuadratic Zx
f) (Integer
p,(Gfpx
g,Gfpx
h,Gfpx
s,Gfpx
t)) [(Integer, (Gfpx, Gfpx, Gfpx, Gfpx))]
-> Int -> (Integer, (Gfpx, Gfpx, Gfpx, Gfpx))
forall a. [a] -> Int -> a
!! Int
k

-------------------------------------------------------------------------------
-- Recombining factors mod p^k to find true factors in Z[x]
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials
-- https://www.math.fsu.edu/~hoeij/papers/issac11/A.pdf
-------------------------------------------------------------------------------

recombineFactors :: Zx -> Integer -> [Gfpx] -> [Zx]
recombineFactors :: Zx -> Integer -> [Gfpx] -> [Zx]
recombineFactors Zx
f Integer
pk =
    \[Gfpx]
hs -> case [Gfpx] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Gfpx]
hs of
             Int
0 -> []
             Int
1 -> [Zx
f]
             Int
_ -> [(Gfpx, [Gfpx], [Gfpx])] -> [Zx]
go [(Integer -> Zx -> Gfpx
Gfpx.fromZx Integer
pk (Integer -> Zx
Zx.constant Integer
c), [], [Gfpx]
hs)]
  where
    c :: Integer
c = Zx -> Integer
Zx.leadingCoeff Zx
f
    d :: Int
d = Zx -> Int
Zx.degree Zx
f
    db :: Int
db = Int
d Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2

    -- Check every subset of hs with degree sum at most db
    -- If s is a subset of t, then check s before t
    go :: [(Gfpx, [Gfpx], [Gfpx])] -> [Zx]
go [] = [Zx
f]
    go ((Gfpx
_,[Gfpx]
_,[]) : [(Gfpx, [Gfpx], [Gfpx])]
wl) = [(Gfpx, [Gfpx], [Gfpx])] -> [Zx]
go [(Gfpx, [Gfpx], [Gfpx])]
wl
    go ((Gfpx
g, [Gfpx]
gs, Gfpx
h : [Gfpx]
hs) : [(Gfpx, [Gfpx], [Gfpx])]
wl) = Gfpx
-> Gfpx
-> [Gfpx]
-> [Gfpx]
-> (Gfpx, [Gfpx], [Gfpx])
-> [(Gfpx, [Gfpx], [Gfpx])]
-> [Zx]
check Gfpx
g Gfpx
h [Gfpx]
gs [Gfpx]
hs (Gfpx
g, Gfpx
h Gfpx -> [Gfpx] -> [Gfpx]
forall a. a -> [a] -> [a]
: [Gfpx]
gs, [Gfpx]
hs) [(Gfpx, [Gfpx], [Gfpx])]
wl

    check :: Gfpx
-> Gfpx
-> [Gfpx]
-> [Gfpx]
-> (Gfpx, [Gfpx], [Gfpx])
-> [(Gfpx, [Gfpx], [Gfpx])]
-> [Zx]
check Gfpx
g Gfpx
h [Gfpx]
gs [Gfpx]
hs (Gfpx, [Gfpx], [Gfpx])
w [(Gfpx, [Gfpx], [Gfpx])]
wl =
        if Gfpx -> Int
Gfpx.degree Gfpx
g Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Gfpx -> Int
Gfpx.degree Gfpx
h Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
db then [(Gfpx, [Gfpx], [Gfpx])] -> [Zx]
go ((Gfpx, [Gfpx], [Gfpx])
w (Gfpx, [Gfpx], [Gfpx])
-> [(Gfpx, [Gfpx], [Gfpx])] -> [(Gfpx, [Gfpx], [Gfpx])]
forall a. a -> [a] -> [a]
: [(Gfpx, [Gfpx], [Gfpx])]
wl)
        else if Zx -> Zx -> Bool
Zx.divides Zx
gz Zx
f then Zx
gz Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: Zx -> Integer -> [Gfpx] -> [Zx]
recombineFactors Zx
f' Integer
pk ([Gfpx]
gs [Gfpx] -> [Gfpx] -> [Gfpx]
forall a. [a] -> [a] -> [a]
++ [Gfpx]
hs)
        else [(Gfpx, [Gfpx], [Gfpx])] -> [Zx]
go ((Gfpx, [Gfpx], [Gfpx])
w (Gfpx, [Gfpx], [Gfpx])
-> [(Gfpx, [Gfpx], [Gfpx])] -> [(Gfpx, [Gfpx], [Gfpx])]
forall a. a -> [a] -> [a]
: (Gfpx
g',[Gfpx]
gs,[Gfpx]
hs) (Gfpx, [Gfpx], [Gfpx])
-> [(Gfpx, [Gfpx], [Gfpx])] -> [(Gfpx, [Gfpx], [Gfpx])]
forall a. a -> [a] -> [a]
: [(Gfpx, [Gfpx], [Gfpx])]
wl)
      where
        g' :: Gfpx
g' = Integer -> Gfpx -> Gfpx -> Gfpx
Gfpx.multiply Integer
pk Gfpx
g Gfpx
h
        gz :: Zx
gz = Zx -> Zx
Zx.primitive (Zx -> Zx) -> Zx -> Zx
forall a b. (a -> b) -> a -> b
$ Integer -> Gfpx -> Zx
Gfpx.toSmallestZx Integer
pk Gfpx
g'
        f' :: Zx
f' = Zx -> Zx -> Zx
Zx.quotient Zx
f Zx
gz

-------------------------------------------------------------------------------
-- Factorization using the Berlekamp-Zassenhaus algorithm
--
-- https://en.wikipedia.org/wiki/Factorization_of_polynomials
-- https://www.math.fsu.edu/~hoeij/papers/issac11/A.pdf
-------------------------------------------------------------------------------

factorSquareFree :: Zx -> [Zx]
factorSquareFree :: Zx -> [Zx]
factorSquareFree Zx
f | Zx -> Bool
Zx.isOne Zx
f = []
factorSquareFree Zx
f | Zx -> Bool
Zx.isLinear Zx
f = [Zx
f]
factorSquareFree Zx
f =
    if [Gfpx] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Gfpx]
gs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 then [Zx
f] else Zx -> Integer -> [Gfpx] -> [Zx]
recombineFactors Zx
f Integer
pk [Gfpx]
hs
  where
    p :: Integer
p = Zx -> Integer
suitablePrime Zx
f
    g :: Gfpx
g = Integer -> Zx -> Gfpx
monicGfpx Integer
p Zx
f
    gs :: [Gfpx]
gs = Integer -> Gfpx -> [Gfpx]
Gfpx.factorSquareFreeBerlekamp Integer
p Gfpx
g
    (Int
k,Integer
pk) = Zx -> Integer -> (Int, Integer)
henselLiftModulus Zx
f Integer
p
    hs :: [Gfpx]
hs = Zx -> Integer -> Int -> [Gfpx] -> [Gfpx]
henselLiftFactors (Gfpx -> Zx
Gfpx.toZx (Integer -> Zx -> Gfpx
monicGfpx Integer
pk Zx
f)) Integer
p Int
k [Gfpx]
gs

factorPrimitive :: Zx -> [(Zx,Integer)]
factorPrimitive :: Zx -> [(Zx, Integer)]
factorPrimitive = [[(Zx, Integer)]] -> [(Zx, Integer)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[(Zx, Integer)]] -> [(Zx, Integer)])
-> (Zx -> [[(Zx, Integer)]]) -> Zx -> [(Zx, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Zx -> [(Zx, Integer)])
-> [Integer] -> [Zx] -> [[(Zx, Integer)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Zx -> [(Zx, Integer)]
forall b. b -> Zx -> [(Zx, b)]
fsf [Integer
1..] ([Zx] -> [[(Zx, Integer)]])
-> (Zx -> [Zx]) -> Zx -> [[(Zx, Integer)]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> [Zx]
Zx.squareFreeDecomposition
  where fsf :: b -> Zx -> [(Zx, b)]
fsf b
k = (Zx -> (Zx, b)) -> [Zx] -> [(Zx, b)]
forall a b. (a -> b) -> [a] -> [b]
map (\Zx
g -> (Zx
g,b
k)) ([Zx] -> [(Zx, b)]) -> (Zx -> [Zx]) -> Zx -> [(Zx, b)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> [Zx]
factorSquareFree

factor :: Zx -> (Integer,[(Zx,Integer)])
factor :: Zx -> (Integer, [(Zx, Integer)])
factor Zx
f = (Integer
c, Zx -> [(Zx, Integer)]
factorPrimitive Zx
g)
  where (Integer
c,Zx
g) = Zx -> (Integer, Zx)
Zx.contentPrimitive Zx
f

-------------------------------------------------------------------------------
-- Irreducibility test
-------------------------------------------------------------------------------

irreducible :: Zx -> Bool
irreducible :: Zx -> Bool
irreducible Zx
f =
    Zx -> Bool
Zx.isPrimitive Zx
f Bool -> Bool -> Bool
&&
    Zx -> Bool
Zx.squareFree Zx
f Bool -> Bool -> Bool
&&
    [Zx] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Zx -> [Zx]
factorSquareFree Zx
f) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1

{-
ghci Factor/Bz.hs
:break factorSquareFree
factor $ Zx.fromCoeff [4,47,-2,-23,18,10]
factor $ Zx.swinnertonDyer !! 4
factor $ Zx.multiply (Zx.swinnertonDyer !! 3) (Zx.swinnertonDyer !! 4)
factor $ Zx.multiply (Zx.fromCoeff [2,-2]) (Zx.fromCoeff [-3,3])
:show bindings
-}