{-# 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 :: forall k.
(Eq k, Fractional k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfreeChar0 UPolynomial k
0 = []
sqfreeChar0 UPolynomial k
p = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [UPolynomial k
qforall a b. (Num a, Integral b) => a -> b -> a
^Integer
m | (UPolynomial k
q,Integer
m) <- [(UPolynomial k, Integer)]
result] forall a. Eq a => a -> a -> Bool
== UPolynomial k
p) forall a b. (a -> b) -> a -> b
$ [(UPolynomial k, Integer)]
result
  where
    result :: [(UPolynomial k, Integer)]
result = forall {k}.
(Eq k, Fractional k) =>
UPolynomial k
-> UPolynomial k
-> Integer
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go UPolynomial k
p (UPolynomial k
p forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.gcd UPolynomial k
p (forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> v -> Polynomial k v
P.deriv UPolynomial k
p X
X)) Integer
0 []
    go :: UPolynomial k
-> UPolynomial k
-> Integer
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go UPolynomial k
p UPolynomial k
flat !Integer
m [(UPolynomial k, Integer)]
result
      | forall t. Degree t => t -> Integer
P.deg UPolynomial k
flat forall a. Ord a => a -> a -> Bool
<= Integer
0 = [(UPolynomial k
p,Integer
1) | UPolynomial k
p forall a. Eq a => a -> a -> Bool
/= UPolynomial k
1] forall a. [a] -> [a] -> [a]
++ forall a. [a] -> [a]
reverse [(UPolynomial k, Integer)]
result
      | Bool
otherwise     = UPolynomial k
-> UPolynomial k
-> Integer
-> [(UPolynomial k, Integer)]
-> [(UPolynomial k, Integer)]
go UPolynomial k
p' UPolynomial k
flat' Integer
m' ((UPolynomial k
flat forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` UPolynomial k
flat', Integer
m') forall a. a -> [a] -> [a]
: [(UPolynomial k, Integer)]
result)
          where
            (UPolynomial k
p',Integer
n) = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f UPolynomial k
p UPolynomial k
flat
            flat' :: UPolynomial k
flat'  = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
P.gcd UPolynomial k
p' UPolynomial k
flat
            m' :: Integer
m' = Integer
m forall a. Num a => a -> a -> a
+ Integer
n

f :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f :: forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f UPolynomial k
p1 UPolynomial k
p2 = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (UPolynomial k
p1 forall a. Eq a => a -> a -> Bool
== UPolynomial k
p2 forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
m forall a. Num a => a -> a -> a
* UPolynomial k
q) forall a b. (a -> b) -> a -> b
$ (UPolynomial k, Integer)
result
  where
    result :: (UPolynomial k, Integer)
result@(UPolynomial k
q, Integer
m) = forall {b}. Num b => b -> UPolynomial k -> (UPolynomial k, b)
go Integer
0 UPolynomial k
p1
    go :: b -> UPolynomial k -> (UPolynomial k, b)
go !b
m UPolynomial k
p =
      case UPolynomial k
p forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
`P.divMod` UPolynomial k
p2 of
        (UPolynomial k
q, UPolynomial k
0) -> b -> UPolynomial k -> (UPolynomial k, b)
go (b
mforall a. Num a => a -> a -> a
+b
1) UPolynomial k
q
        (UPolynomial k, UPolynomial k)
_ -> (UPolynomial k
p, b
m)


instance P.SQFree (UPolynomial Rational) where
  sqfree :: UPolynomial Rational -> [(UPolynomial Rational, Integer)]
sqfree = forall k.
(Eq k, Fractional k) =>
UPolynomial k -> [(UPolynomial k, Integer)]
sqfreeChar0

instance P.SQFree (UPolynomial Integer) where
  sqfree :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]
sqfree UPolynomial Integer
0 = [(UPolynomial Integer
0,Integer
1)]
  sqfree UPolynomial Integer
f = forall {k} {b} {v}.
(Integral k, Integral b, Ord v) =>
Ratio k
-> [(Polynomial k v, b)]
-> [(Polynomial (Ratio k) v, b)]
-> [(Polynomial k v, b)]
go Rational
1 [] (forall a. SQFree a => a -> [(a, Integer)]
P.sqfree (forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff forall a b. (Integral a, Num b) => a -> b
fromIntegral UPolynomial Integer
f))
    where
      go :: Ratio k
-> [(Polynomial k v, b)]
-> [(Polynomial (Ratio k) v, b)]
-> [(Polynomial k v, b)]
go !Ratio k
u [(Polynomial k v, b)]
ys [] =
        forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall a. Ratio a -> a
denominator Ratio k
u forall a. Eq a => a -> a -> Bool
== k
1) forall a b. (a -> b) -> a -> b
$
          [(forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (forall a. Ratio a -> a
numerator Ratio k
u), b
1) | Ratio k
u forall a. Eq a => a -> a -> Bool
/= Ratio k
1] forall a. [a] -> [a] -> [a]
++ [(Polynomial k v, b)]
ys
      go !Ratio k
u [(Polynomial k v, b)]
ys ((Polynomial (Ratio k) v
g,b
n):[(Polynomial (Ratio k) v, b)]
xs)
        | forall t. Degree t => t -> Integer
P.deg Polynomial (Ratio k) v
g forall a. Ord a => a -> a -> Bool
<= Integer
0 = Ratio k
-> [(Polynomial k v, b)]
-> [(Polynomial (Ratio k) v, b)]
-> [(Polynomial k v, b)]
go (Ratio k
u forall a. Num a => a -> a -> a
* forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff forall v. Monomial v
P.mone Polynomial (Ratio k) v
g) [(Polynomial k v, b)]
ys [(Polynomial (Ratio k) v, b)]
xs
        | Bool
otherwise    = Ratio k
-> [(Polynomial k v, b)]
-> [(Polynomial (Ratio k) v, b)]
-> [(Polynomial k v, b)]
go (Ratio k
u forall a. Num a => a -> a -> a
* (forall k v. (ContPP k, Ord v) => Polynomial k v -> k
P.cont Polynomial (Ratio k) v
g)forall a b. (Num a, Integral b) => a -> b -> a
^b
n) ((forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp Polynomial (Ratio k) v
g, b
n) forall a. a -> [a] -> [a]
: [(Polynomial k v, b)]
ys) [(Polynomial (Ratio k) v, b)]
xs