content: {-# OPTIONS_HADDOCK show-extensions #-}
  {-# LANGUAGE BangPatterns #-}
  {-# LANGUAGE FlexibleInstances #-}
  {-# LANGUAGE TypeSynonymInstances #-}
  -----------------------------------------------------------------------------
  -- |
  -- Module      :  ToySolver.Data.Polynomial.Factorization.SquareFree
  -- Copyright   :  (c) Masahiro Sakai 2013
  -- License     :  BSD-style
  --
  -- Maintainer  :  masahiro.sakai@gmail.com
  -- Stability   :  provisional
  -- Portability :  non-portable
  --
  -- References:
  --
  -- * <http://www.math.kobe-u.ac.jp/Asir/ca.pdf>
  --
  -----------------------------------------------------------------------------
  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
  
  -- | Square-free decomposition of univariate polynomials over a field of characteristic 0.
  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.pp g, n) : ys) xs
  