{-# 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 | forall t. Degree t => t -> Integer
P.deg UPolynomial Integer
p forall a. Eq a => a -> a -> Bool
== Integer
0 = [(UPolynomial Integer
p,Integer
1)]
factor UPolynomial Integer
p = [(forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Integer
c, Integer
1) | Integer
c forall a. Eq a => a -> a -> Bool
/= Integer
1] forall a. [a] -> [a] -> [a]
++ [(UPolynomial Integer
q, forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m) | (UPolynomial Integer
q,Int
m) <- forall a. MultiSet a -> [(a, Int)]
MultiSet.toOccurList MultiSet (UPolynomial Integer)
qs]
  where
    (Integer
c,MultiSet (UPolynomial Integer)
qs) = (Integer, MultiSet (UPolynomial Integer))
-> (Integer, MultiSet (UPolynomial Integer))
normalize (forall k v. (ContPP k, Ord v) => Polynomial k v -> k
P.cont UPolynomial Integer
p, UPolynomial Integer -> MultiSet (UPolynomial Integer)
factor' (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) = forall {a}.
(Num a, Ord a) =>
[(Polynomial a X, Int)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go (forall a. MultiSet a -> [(a, Int)]
MultiSet.toOccurList MultiSet (UPolynomial Integer)
ps) Integer
c forall a. MultiSet a
MultiSet.empty
  where
    go :: [(Polynomial a X, Int)]
-> 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,Int
m) : [(Polynomial a X, Int)]
ps) !a
c !MultiSet (Polynomial a X)
qs
      | forall t. Degree t => t -> Integer
P.deg Polynomial a X
p forall a. Eq a => a -> a -> Bool
== Integer
0 = [(Polynomial a X, Int)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Int)]
ps (a
c forall a. Num a => a -> a -> a
* (forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (forall a v. Var a v => v -> a
P.var X
X) Polynomial a X
p) forall a b. (Num a, Integral b) => a -> b -> a
^ Int
m) MultiSet (Polynomial a X)
qs
      | forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat Polynomial a X
p forall a. Ord a => a -> a -> Bool
< a
0 = [(Polynomial a X, Int)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Int)]
ps (a
c forall a. Num a => a -> a -> a
* (-a
1)forall a b. (Num a, Integral b) => a -> b -> a
^Int
m) (forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany (-Polynomial a X
p) Int
m MultiSet (Polynomial a X)
qs)
      | Bool
otherwise = [(Polynomial a X, Int)]
-> a -> MultiSet (Polynomial a X) -> (a, MultiSet (Polynomial a X))
go [(Polynomial a X, Int)]
ps a
c (forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany Polynomial a X
p Int
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 (forall a. a -> MultiSet a
MultiSet.singleton UPolynomial Integer
p) 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
      | 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' (forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
p Int
m MultiSet (UPolynomial Integer)
ret)
            Just (UPolynomial Integer
q1,UPolynomial Integer
q2) ->
              MultiSet (UPolynomial Integer)
-> MultiSet (UPolynomial Integer) -> MultiSet (UPolynomial Integer)
go (forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
q1 Int
m forall a b. (a -> b) -> a -> b
$ forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany UPolynomial Integer
q2 Int
m MultiSet (UPolynomial Integer)
ps') MultiSet (UPolynomial Integer)
ret
          where
            p :: UPolynomial Integer
p   = forall a. MultiSet a -> a
MultiSet.findMin MultiSet (UPolynomial Integer)
ps
            m :: Int
m   = forall a. Ord a => a -> MultiSet a -> Int
MultiSet.occur UPolynomial Integer
p MultiSet (UPolynomial Integer)
ps
            ps' :: MultiSet (UPolynomial Integer)
ps' = 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 forall a. Eq a => a -> a -> Bool
== forall a v. Var a v => v -> a
P.var X
X = forall a. Maybe a
Nothing
factor2 UPolynomial Integer
p =
  case forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find (\(Integer
_,Integer
yi) -> Integer
yiforall a. Eq a => a -> a -> Bool
==Integer
0) [(Integer, Integer)]
vs of
    Just (Integer
xi,Integer
_) ->
      let q1 :: UPolynomial Integer
q1 = UPolynomial Integer
x forall a. Num a => a -> a -> a
- forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Integer
xi
          q2 :: UPolynomial Rational
q2 = UPolynomial Rational
p' forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.div` forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff forall a. Num a => Integer -> a
fromInteger UPolynomial Integer
q1
      in forall a. a -> Maybe a
Just (UPolynomial Integer
q1, 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 = forall a b. (a -> b) -> [a] -> [b]
map forall k. (Eq k, Fractional k) => [(k, k)] -> UPolynomial k
Interpolation.interpolate forall a b. (a -> b) -> a -> b
$
                  forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[(forall a. Num a => Integer -> a
fromInteger Integer
xi, 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, forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q1 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' forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)
`P.divMod` UPolynomial Rational
q1
               , UPolynomial Rational
r forall a. Eq a => a -> a -> Bool
== UPolynomial Rational
0, forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q2 forall a. Ord a => a -> a -> Bool
> Integer
0, UPolynomial Rational -> Bool
isUPolyZ UPolynomial Rational
q2
               ]
      in case [(UPolynomial Rational, UPolynomial Rational)]
zs of
           [] -> forall a. Maybe a
Nothing
           (UPolynomial Rational
q1,UPolynomial Rational
q2):[(UPolynomial Rational, UPolynomial Rational)]
_ -> forall a. a -> Maybe a
Just (forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Rational
q1, forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp UPolynomial Rational
q2)
  where
    n :: Integer
n = forall t. Degree t => t -> Integer
P.deg UPolynomial Integer
p forall a. Integral a => a -> a -> a
`div` Integer
2
    xs :: [Integer]
xs = forall a. Int -> [a] -> [a]
take (forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n forall a. Num a => a -> a -> a
+ Int
1) [Integer]
xvalues
    vs :: [(Integer, Integer)]
vs = [(Integer
x, 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 = forall a v. Var a v => v -> a
P.var X
X
    p' :: UPolynomial Rational
    p' :: UPolynomial Rational
p' = forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff forall a. Num a => Integer -> a
fromInteger UPolynomial Integer
p

isUPolyZ :: UPolynomial Rational -> Bool
isUPolyZ :: UPolynomial Rational -> Bool
isUPolyZ UPolynomial Rational
p = forall (t :: * -> *). Foldable t => t Bool -> Bool
and [forall a. RealFrac a => a -> Bool
isInteger Rational
c | (Rational
c,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 forall a. a -> [a] -> [a]
: forall a. [a] -> [a] -> [a]
interleave [Integer
1,Integer
2..] [-Integer
1,-Integer
2..]

interleave :: [a] -> [a] -> [a]
interleave :: forall a. [a] -> [a] -> [a]
interleave [a]
xs [] = [a]
xs
interleave [] [a]
ys     = [a]
ys
interleave (a
x:[a]
xs) [a]
ys = a
x forall 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 forall a. [a] -> [a] -> [a]
++ forall a b. (a -> b) -> [a] -> [b]
map forall a. Num a => a -> a
negate [Integer]
xs
  where
    ps :: MultiSet Integer
ps = Integer -> MultiSet Integer
primeFactors (forall a. Num a => a -> a
abs Integer
x)
    xs :: [Integer]
xs = forall a b. (a -> b) -> [a] -> [b]
map forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [forall a. Int -> [a] -> [a]
take (Int
nforall a. Num a => a -> a -> a
+Int
1) (forall a. (a -> a) -> a -> [a]
iterate (Integer
pforall a. Num a => a -> a -> a
*) Integer
1) | (Integer
p,Int
n) <- forall a. MultiSet a -> [(a, Int)]
MultiSet.toOccurList MultiSet Integer
ps]

primeFactors :: Integer -> MultiSet Integer
primeFactors :: Integer -> MultiSet Integer
primeFactors Integer
0 = forall a. MultiSet a
MultiSet.empty
primeFactors Integer
n = Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
go Integer
n forall int. Integral int => [int]
primes 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
pforall a. Num a => a -> a -> a
*Integer
p forall a. Ord a => a -> a -> Bool
> Integer
n   = forall a. Ord a => a -> MultiSet a -> MultiSet a
MultiSet.insert Integer
n MultiSet Integer
result
      | Bool
otherwise =
          case Integer -> Integer -> (Int, Integer)
f Integer
p Integer
n of
            (Int
m,Integer
n') -> Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer
go Integer
n' [Integer]
ps (forall a. Ord a => a -> Int -> MultiSet a -> MultiSet a
MultiSet.insertMany Integer
p Int
m MultiSet Integer
result)

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