module Data.Polynomial.Factorization.SquareFree
( sqfree
) where
import Control.Exception
import Data.Polynomial
sqfree :: (Eq k, Fractional k) => UPolynomial k -> [(UPolynomial k, Integer)]
sqfree 0 = []
sqfree p = assert (product [q^m | (q,m) <- result] == p) $ result
where
result = go p (p `polyDiv` polyGCD p (deriv p X)) 0 []
go p flat !m result
| deg flat <= 0 = [(p,1) | p /= 1] ++ reverse result
| otherwise = go p' flat' m' ((flat `polyDiv` flat', m') : result)
where
(p',n) = f p flat
flat' = polyGCD p' flat
m' = m + n
f :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f p1 p2 = assert (p1 == p2 ^ m * q) $ result
where
result@(q, m) = go 0 p1
go !m p =
case p `polyDivMod` p2 of
(q, 0) -> go (m+1) q
_ -> (p, m)