module ToySolver.Data.Polynomial.Factorization.SquareFree
( sqfreeChar0
) where
import Control.Exception
import Data.Ratio
import ToySolver.Data.Polynomial.Base (UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial.Base as P
sqfreeChar0 :: (Eq k, Fractional k) => UPolynomial k -> [(UPolynomial k, Integer)]
sqfreeChar0 0 = []
sqfreeChar0 p = assert (product [q^m | (q,m) <- result] == p) $ result
where
result = go p (p `P.div` P.gcd p (P.deriv p X)) 0 []
go p flat !m result
| P.deg flat <= 0 = [(p,1) | p /= 1] ++ reverse result
| otherwise = go p' flat' m' ((flat `P.div` flat', m') : result)
where
(p',n) = f p flat
flat' = P.gcd 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 `P.divMod` p2 of
(q, 0) -> go (m+1) q
_ -> (p, m)
instance P.SQFree (UPolynomial Rational) where
sqfree = sqfreeChar0
instance P.SQFree (UPolynomial Integer) where
sqfree 0 = [(0,1)]
sqfree f = go 1 [] (P.sqfree (P.mapCoeff fromIntegral f))
where
go !u ys [] =
assert (denominator u == 1) $
[(P.constant (numerator u), 1) | u /= 1] ++ ys
go !u ys ((g,n):xs)
| P.deg g <= 0 = go (u * P.coeff P.mone g) ys xs
| otherwise = go (u * (P.cont g)^n) ((P.mapCoeff numerator (P.pp g), n) : ys) xs