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
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))
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
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
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
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
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
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