module ToySolver.Data.Polynomial.Factorization.Hensel.Internal
( hensel
, cabook_proposition_5_10
, cabook_proposition_5_11
) where
import Control.Exception (assert)
import Data.FiniteField
import Data.Proxy
import GHC.TypeLits
import ToySolver.Data.Polynomial.Base (UPolynomial)
import qualified ToySolver.Data.Polynomial.Base as P
hensel :: forall p. KnownNat p => UPolynomial Integer -> [UPolynomial (PrimeField p)] -> Integer -> [UPolynomial Integer]
hensel f fs1 k
| k <= 0 = error "ToySolver.Data.Polynomial.Factorization.Hensel.hensel: k <= 0"
| otherwise = assert precondition $ go 1 (map (P.mapCoeff Data.FiniteField.toInteger) fs1)
where
precondition =
P.mapCoeff fromInteger f == product fs1 &&
P.deg f == P.deg (product fs1)
p :: Integer
p = natVal (Proxy :: Proxy p)
go :: Integer -> [UPolynomial Integer] -> [UPolynomial Integer]
go !i fs
| i==k = assert (check i fs) $ fs
| otherwise = assert (check i fs) $ go (i+1) (lift i fs)
check :: Integer -> [UPolynomial Integer] -> Bool
check k' fs =
and
[ P.mapCoeff (`mod` pk) f == P.mapCoeff (`mod` pk) (product fs)
, fs1 == map (P.mapCoeff fromInteger) fs
, and [P.deg fi1 == P.deg fik | (fi1, fik) <- zip fs1 fs]
]
where
pk = p ^ k'
lift :: Integer -> [UPolynomial Integer] -> [UPolynomial Integer]
lift k' fs = fs'
where
pk = p^k'
pk1 = p^(k'+1)
h :: UPolynomial Integer
h = P.mapCoeff (\c -> (c `mod` pk1) `div` pk) (f product fs)
hs :: [UPolynomial (PrimeField p)]
hs = cabook_proposition_5_11 (map (P.mapCoeff fromInteger) fs) (P.mapCoeff fromInteger h)
fs' :: [UPolynomial Integer]
fs' = [ P.mapCoeff (`mod` pk1) (fi + P.constant pk * P.mapCoeff Data.FiniteField.toInteger hi)
| (fi, hi) <- zip fs hs ]
cabook_proposition_5_10
:: forall k. (Num k, Fractional k, Eq k)
=> [UPolynomial k] -> [UPolynomial k]
cabook_proposition_5_10 fs = normalize (go fs)
where
check :: [UPolynomial k] -> [UPolynomial k] -> Bool
check es fs' = sum [ei * (product fs' `P.div` fi) | (ei,fi) <- zip es fs'] == 1
go :: [UPolynomial k] -> [UPolynomial k]
go [] = error "cabook_proposition_5_10: empty list"
go [fi] = assert (check [1] [fi]) [1]
go (fi : fs') =
case P.exgcd (product fs') fi of
(g,ei,v) ->
assert (g == 1) $
let es' = go fs'
es = ei : map (v*) es'
in assert (check es (fi : fs')) es
normalize :: [UPolynomial k] -> [UPolynomial k]
normalize es = assert (check es2 fs) es2
where
es2 = zipWith P.mod es fs
cabook_proposition_5_11
:: forall k. (Fractional k, Ord k)
=> [UPolynomial k] -> UPolynomial k -> [UPolynomial k]
cabook_proposition_5_11 fs g =
assert (P.deg g <= P.deg (product fs)) $
assert (P.deg c <= 0) $
assert (check es2) $
es2
where
es = map (g*) $ cabook_proposition_5_10 fs
c = sum [ei `P.div` fi | (ei,fi) <- zip es fs]
es2 = case zipWith P.mod es fs of
e2' : es2' -> e2' + c * head fs : es2'
check :: [UPolynomial k] -> Bool
check es' =
sum [ei * (product fs `P.div` fi) | (ei,fi) <- zip es' fs] == g &&
and [P.deg ei <= P.deg fi | (ei,fi) <- zip es' fs]