{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE BangPatterns #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.Polynomial.Factorization.Kronecker
-- Copyright   :  (c) Masahiro Sakai 2012-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Factoriation of integer-coefficient polynomial using Kronecker's method.
--
-- References:
--
-- * <http://en.wikipedia.org/wiki/Polynomial_factorization>
--
-----------------------------------------------------------------------------
module ToySolver.Data.Polynomial.Factorization.Kronecker
  ( factor
  ) where

import Data.List
import Data.MultiSet (MultiSet)
import qualified Data.MultiSet as MultiSet
import Data.Numbers.Primes (primes)
import Data.Ratio

import ToySolver.Data.Polynomial.Base (Polynomial, UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial.Base as P
import qualified ToySolver.Data.Polynomial.Interpolation.Lagrange as Interpolation
import ToySolver.Internal.Util (isInteger)

factor :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]
factor :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]
factor UPolynomial Integer
0 = [(UPolynomial Integer
0,Integer
1)]
factor UPolynomial Integer
1 = []
factor UPolynomial Integer
p | UPolynomial Integer -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = [(UPolynomial Integer
p,Integer
1)]
factor UPolynomial Integer
p = [(Integer -> UPolynomial Integer
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Integer
c, Integer
1) | Integer
c Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
1] [(UPolynomial Integer, Integer)]
-> [(UPolynomial Integer, Integer)]
-> [(UPolynomial Integer, Integer)]
forall a. [a] -> [a] -> [a]
++ [(UPolynomial Integer
q, Occur -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Occur
m) | (UPolynomial Integer
q,Occur
m) <- MultiSet (UPolynomial Integer) -> [(UPolynomial Integer, Occur)]
forall a. MultiSet a -> [(a, Occur)]
MultiSet.toOccurList MultiSet (UPolynomial Integer)
qs]
  where
    (Integer
c,MultiSet (UPolynomial Integer)
qs) = (Integer, MultiSet (UPolynomial Integer))
-> (Integer, MultiSet (UPolynomial Integer))
normalize (UPolynomial Integer -> Integer
forall k v. (ContPP k, Ord v) => Polynomial k v -> k
P.cont UPolynomial Integer
p, UPolynomial Integer -> MultiSet (UPolynomial Integer)
factor' (UPolynomial Integer -> Polynomial (PPCoeff Integer) X
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Integer
p))

normalize :: (Integer, MultiSet (UPolynomial Integer)) -> (Integer, MultiSet (UPolynomial Integer))
normalize :: (Integer, MultiSet (UPolynomial Integer))
-> (Integer, MultiSet (UPolynomial Integer))
normalize (Integer
c,MultiSet (UPolynomial Integer)
ps) = [(UPolynomial Integer, Occur)]
-> Integer
-> MultiSet (UPolynomial Integer)
-> (Integer, MultiSet (UPolynomial Integer))
forall a.
(Num a, Ord a) =>
[(Polynomial a X, Occur)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go (MultiSet (UPolynomial Integer) -> [(UPolynomial Integer, Occur)]
forall a. MultiSet a -> [(a, Occur)]
MultiSet.toOccurList MultiSet (UPolynomial Integer)
ps) Integer
c MultiSet (UPolynomial Integer)
forall a. MultiSet a
MultiSet.empty
  where
    go :: [(Polynomial a X, Occur)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [] !a
c !MultiSet (Polynomial a X)
qs = (a
c, MultiSet (Polynomial a X)
qs)
    go ((Polynomial a X
p,Occur
m) : [(Polynomial a X, Occur)]
ps) !a
c !MultiSet (Polynomial a X)
qs
      | Polynomial a X -> Integer
forall t. Degree t => t -> Integer
P.deg Polynomial a X
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = [(Polynomial a X, Occur)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Occur)]
ps (a
c a -> a -> a
forall a. Num a => a -> a -> a
* (Monomial X -> Polynomial a X -> a
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (X -> Monomial X
forall a v. Var a v => v -> a
P.var X
X) Polynomial a X
p) a -> Occur -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ Occur
m) MultiSet (Polynomial a X)
qs
      | MonomialOrder X -> Polynomial a X -> a
forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat Polynomial a X
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = [(Polynomial a X, Occur)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Occur)]
ps (a
c a -> a -> a
forall a. Num a => a -> a -> a
* (-a
1)a -> Occur -> a
forall a b. (Num a, Integral b) => a -> b -> a
^Occur
m) (Polynomial a X
-> Occur -> MultiSet (Polynomial a X) -> MultiSet (Polynomial a X)
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany (-Polynomial a X
p) Occur
m MultiSet (Polynomial a X)
qs)
      | Bool
otherwise = [(Polynomial a X, Occur)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Occur)]
ps a
c (Polynomial a X
-> Occur -> MultiSet (Polynomial a X) -> MultiSet (Polynomial a X)
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany Polynomial a X
p Occur
m MultiSet (Polynomial a X)
qs)

factor' :: UPolynomial Integer -> MultiSet (UPolynomial Integer)
factor' :: UPolynomial Integer -> MultiSet (UPolynomial Integer)
factor' UPolynomial Integer
p = MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
go (UPolynomial Integer -> MultiSet (UPolynomial Integer)
forall a. a -> MultiSet a
MultiSet.singleton UPolynomial Integer
p) MultiSet (UPolynomial Integer)
forall a. MultiSet a
MultiSet.empty
  where
    go :: MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
go MultiSet (UPolynomial Integer)
ps MultiSet (UPolynomial Integer)
ret
      | MultiSet (UPolynomial Integer) -> Bool
forall a. MultiSet a -> Bool
MultiSet.null MultiSet (UPolynomial Integer)
ps = MultiSet (UPolynomial Integer)
ret
      | Bool
otherwise =
          case UPolynomial Integer
-> Maybe (UPolynomial Integer, UPolynomial Integer)
factor2 UPolynomial Integer
p of
            Maybe (UPolynomial Integer, UPolynomial Integer)
Nothing ->
              MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
go MultiSet (UPolynomial Integer)
ps' (UPolynomial Integer
-> Occur
-> MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer)
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
p Occur
m MultiSet (UPolynomial Integer)
ret)
            Just (UPolynomial Integer
q1,UPolynomial Integer
q2) ->
              MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
go (UPolynomial Integer
-> Occur
-> MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer)
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
q1 Occur
m (MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer))
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
forall a b. (a -> b) -> a -> b
$ UPolynomial Integer
-> Occur
-> MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer)
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
q2 Occur
m MultiSet (UPolynomial Integer)
ps') MultiSet (UPolynomial Integer)
ret
          where
            p :: UPolynomial Integer
p   = MultiSet (UPolynomial Integer) -> UPolynomial Integer
forall a. MultiSet a -> a
MultiSet.findMin MultiSet (UPolynomial Integer)
ps
            m :: Occur
m   = UPolynomial Integer -> MultiSet (UPolynomial Integer) -> Occur
forall a. Ord a => a -> MultiSet a -> Occur
MultiSet.occur UPolynomial Integer
p MultiSet (UPolynomial Integer)
ps
            ps' :: MultiSet (UPolynomial Integer)
ps' = UPolynomial Integer
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
forall a. Ord a => a -> MultiSet a -> MultiSet a
MultiSet.deleteAll UPolynomial Integer
p MultiSet (UPolynomial Integer)
ps

factor2 :: UPolynomial Integer -> Maybe (UPolynomial Integer, UPolynomial Integer)
factor2 :: UPolynomial Integer
-> Maybe (UPolynomial Integer, UPolynomial Integer)
factor2 UPolynomial Integer
p | UPolynomial Integer
p UPolynomial Integer -> UPolynomial Integer -> Bool
forall a. Eq a => a -> a -> Bool
== X -> UPolynomial Integer
forall a v. Var a v => v -> a
P.var X
X = Maybe (UPolynomial Integer, UPolynomial Integer)
forall a. Maybe a
Nothing
factor2 UPolynomial Integer
p =
  case ((Integer, Integer) -> Bool)
-> [(Integer, Integer)] -> Maybe (Integer, Integer)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (\(Integer
_,Integer
yi) -> Integer
yiInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0) [(Integer, Integer)]
vs of
    Just (Integer
xi,Integer
_) ->
      let q1 :: UPolynomial Integer
q1 = UPolynomial Integer
x UPolynomial Integer -> UPolynomial Integer -> UPolynomial Integer
forall a. Num a => a -> a -> a
- Integer -> UPolynomial Integer
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Integer
xi
          q2 :: UPolynomial Rational
q2 = UPolynomial Rational
p' UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` (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. Num a => Integer -> a
fromInteger UPolynomial Integer
q1
      in (UPolynomial Integer, UPolynomial Integer)
-> Maybe (UPolynomial Integer, UPolynomial Integer)
forall a. a -> Maybe a
Just (UPolynomial Integer
q1, UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Rational
q2)
    Maybe (Integer, Integer)
Nothing ->
      let qs :: [UPolynomial Rational]
qs = ([(Rational, Rational)] -> UPolynomial Rational)
-> [[(Rational, Rational)]] -> [UPolynomial Rational]
forall a b. (a -> b) -> [a] -> [b]
map [(Rational, Rational)] -> UPolynomial Rational
forall k. (Eq k, Fractional k) => [(k, k)] -> UPolynomial k
Interpolation.interpolate ([[(Rational, Rational)]] -> [UPolynomial Rational])
-> [[(Rational, Rational)]] -> [UPolynomial Rational]
forall a b. (a -> b) -> a -> b
$
                  [[(Rational, Rational)]] -> [[(Rational, Rational)]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[(Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
xi, Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
z) | Integer
z <- Integer -> [Integer]
factors Integer
yi] | (Integer
xi,Integer
yi) <- [(Integer, Integer)]
vs]
          zs :: [(UPolynomial Rational, UPolynomial Rational)]
zs = [ (UPolynomial Rational
q1,UPolynomial Rational
q2)
               | UPolynomial Rational
q1 <- [UPolynomial Rational]
qs, UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q1 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0, UPolynomial Rational -> Bool
isUPolyZ UPolynomial Rational
q1
               , let (UPolynomial Rational
q2,UPolynomial Rational
r) = UPolynomial Rational
p' UPolynomial Rational
-> UPolynomial Rational
-> (UPolynomial Rational, UPolynomial Rational)
forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
`P.divMod` UPolynomial Rational
q1
               , UPolynomial Rational
r UPolynomial Rational -> UPolynomial Rational -> Bool
forall a. Eq a => a -> a -> Bool
== UPolynomial Rational
0, UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q2 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0, UPolynomial Rational -> Bool
isUPolyZ UPolynomial Rational
q2
               ]
      in case [(UPolynomial Rational, UPolynomial Rational)]
zs of
           [] -> Maybe (UPolynomial Integer, UPolynomial Integer)
forall a. Maybe a
Nothing
           (UPolynomial Rational
q1,UPolynomial Rational
q2):[(UPolynomial Rational, UPolynomial Rational)]
_ -> (UPolynomial Integer, UPolynomial Integer)
-> Maybe (UPolynomial Integer, UPolynomial Integer)
forall a. a -> Maybe a
Just (UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Rational
q1, UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Rational
q2)
  where
    n :: Integer
n = UPolynomial Integer -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Integer
p Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2
    xs :: [Integer]
xs = Occur -> [Integer] -> [Integer]
forall a. Occur -> [a] -> [a]
take (Integer -> Occur
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n Occur -> Occur -> Occur
forall a. Num a => a -> a -> a
+ Occur
1) [Integer]
xvalues
    vs :: [(Integer, Integer)]
vs = [(Integer
x, (X -> Integer) -> UPolynomial Integer -> Integer
forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> Integer
x) UPolynomial Integer
p) | Integer
x <- [Integer]
xs]
    x :: UPolynomial Integer
x = X -> UPolynomial Integer
forall a v. Var a v => v -> a
P.var X
X
    p' :: UPolynomial Rational
    p' :: UPolynomial Rational
p' = (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. Num a => Integer -> a
fromInteger UPolynomial Integer
p

isUPolyZ :: UPolynomial Rational -> Bool
isUPolyZ :: UPolynomial Rational -> Bool
isUPolyZ UPolynomial Rational
p = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Rational -> Bool
forall a. RealFrac a => a -> Bool
isInteger Rational
c | (Rational
c,Monomial X
_) <- UPolynomial Rational -> [(Rational, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial Rational
p]

-- [0, 1, -1, 2, -2, 3, -3 ..]
xvalues :: [Integer]
xvalues :: [Integer]
xvalues = Integer
0 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
interleave [Integer
1,Integer
2..] [-Integer
1,-Integer
2..]

interleave :: [a] -> [a] -> [a]
interleave :: [a] -> [a] -> [a]
interleave [a]
xs [] = [a]
xs
interleave [] [a]
ys     = [a]
ys
interleave (a
x:[a]
xs) [a]
ys = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
interleave [a]
ys [a]
xs

factors :: Integer -> [Integer]
factors :: Integer -> [Integer]
factors Integer
0 = []
factors Integer
x = [Integer]
xs [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Integer
forall a. Num a => a -> a
negate [Integer]
xs
  where
    ps :: MultiSet Integer
ps = Integer -> MultiSet Integer
primeFactors (Integer -> Integer
forall a. Num a => a -> a
abs Integer
x)
    xs :: [Integer]
xs = ([Integer] -> Integer) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$ [[Integer]] -> [[Integer]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [Occur -> [Integer] -> [Integer]
forall a. Occur -> [a] -> [a]
take (Occur
nOccur -> Occur -> Occur
forall a. Num a => a -> a -> a
+Occur
1) ((Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
pInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*) Integer
1) | (Integer
p,Occur
n) <- MultiSet Integer -> [(Integer, Occur)]
forall a. MultiSet a -> [(a, Occur)]
MultiSet.toOccurList MultiSet Integer
ps]

primeFactors :: Integer -> MultiSet Integer
primeFactors :: Integer -> MultiSet Integer
primeFactors Integer
0 = MultiSet Integer
forall a. MultiSet a
MultiSet.empty
primeFactors Integer
n = Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
go Integer
n [Integer]
forall int. Integral int => [int]
primes MultiSet Integer
forall a. MultiSet a
MultiSet.empty
  where
    go :: Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
    go :: Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
go Integer
1 ![Integer]
_ !MultiSet Integer
result = MultiSet Integer
result
    go Integer
n (Integer
p:[Integer]
ps) !MultiSet Integer
result
      | Integer
pInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
n   = Integer -> MultiSet Integer -> MultiSet Integer
forall a. Ord a => a -> MultiSet a -> MultiSet a
MultiSet.insert Integer
n MultiSet Integer
result
      | Bool
otherwise =
          case Integer -> Integer -> (Occur, Integer)
f Integer
p Integer
n of
            (Occur
m,Integer
n') -> Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
go Integer
n' [Integer]
ps (Integer -> Occur -> MultiSet Integer -> MultiSet Integer
forall a. Ord a => a -> Occur -> MultiSet a -> MultiSet a
MultiSet.insertMany Integer
p Occur
m MultiSet Integer
result)

    f :: Integer -> Integer -> (Int, Integer)
    f :: Integer -> Integer -> (Occur, Integer)
f Integer
p = Occur -> Integer -> (Occur, Integer)
forall a. Num a => a -> Integer -> (a, Integer)
go2 Occur
0
      where
        go2 :: a -> Integer -> (a, Integer)
go2 !a
m !Integer
n
          | Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = a -> Integer -> (a, Integer)
go2 (a
ma -> a -> a
forall a. Num a => a -> a -> a
+a
1) (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
p)
          | Bool
otherwise = (a
m, Integer
n)