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

f :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f :: UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)
f UPolynomial k
p1 UPolynomial k
p2 = Bool -> (UPolynomial k, Integer) -> (UPolynomial k, Integer)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (UPolynomial k
p1 UPolynomial k -> UPolynomial k -> Bool
forall a. Eq a => a -> a -> Bool
== UPolynomial k
p2 UPolynomial k -> Integer -> UPolynomial k
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
m UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* UPolynomial k
q) ((UPolynomial k, Integer) -> (UPolynomial k, Integer))
-> (UPolynomial k, Integer) -> (UPolynomial k, Integer)
forall a b. (a -> b) -> a -> b
$ (UPolynomial k, Integer)
result
  where
    result :: (UPolynomial k, Integer)
result@(UPolynomial k
q, Integer
m) = Integer -> UPolynomial k -> (UPolynomial k, Integer)
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 UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
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
mb -> b -> b
forall 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 = UPolynomial Rational -> [(UPolynomial Rational, Integer)]
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 = Rational
-> [(UPolynomial Integer, Integer)]
-> [(UPolynomial Rational, Integer)]
-> [(UPolynomial Integer, Integer)]
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 [] (UPolynomial Rational -> [(UPolynomial Rational, Integer)]
forall a. SQFree a => a -> [(a, Integer)]
P.sqfree ((Integer -> Rational)
-> UPolynomial Integer -> UPolynomial Rational
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff Integer -> Rational
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 [] =
        Bool -> [(Polynomial k v, b)] -> [(Polynomial k v, b)]
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Ratio k -> k
forall a. Ratio a -> a
denominator Ratio k
u k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
1) ([(Polynomial k v, b)] -> [(Polynomial k v, b)])
-> [(Polynomial k v, b)] -> [(Polynomial k v, b)]
forall a b. (a -> b) -> a -> b
$
          [(k -> Polynomial k v
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (Ratio k -> k
forall a. Ratio a -> a
numerator Ratio k
u), b
1) | Ratio k
u Ratio k -> Ratio k -> Bool
forall a. Eq a => a -> a -> Bool
/= Ratio k
1] [(Polynomial k v, b)]
-> [(Polynomial k v, b)] -> [(Polynomial k v, b)]
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)
        | Polynomial (Ratio k) v -> Integer
forall t. Degree t => t -> Integer
P.deg Polynomial (Ratio k) v
g Integer -> Integer -> Bool
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 Ratio k -> Ratio k -> Ratio k
forall a. Num a => a -> a -> a
* Monomial v -> Polynomial (Ratio k) v -> Ratio k
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff Monomial v
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 Ratio k -> Ratio k -> Ratio k
forall a. Num a => a -> a -> a
* (Polynomial (Ratio k) v -> Ratio k
forall k v. (ContPP k, Ord v) => Polynomial k v -> k
P.cont Polynomial (Ratio k) v
g)Ratio k -> b -> Ratio k
forall a b. (Num a, Integral b) => a -> b -> a
^b
n) ((Polynomial (Ratio k) v -> Polynomial (PPCoeff (Ratio k)) v
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp Polynomial (Ratio k) v
g, b
n) (Polynomial k v, b)
-> [(Polynomial k v, b)] -> [(Polynomial k v, b)]
forall a. a -> [a] -> [a]
: [(Polynomial k v, b)]
ys) [(Polynomial (Ratio k) v, b)]
xs