{-# LANGUAGE BangPatterns, DataKinds, ExtendedDefaultRules          #-}
{-# LANGUAGE FlexibleContexts, GADTs, MultiParamTypeClasses         #-}
{-# LANGUAGE NoImplicitPrelude, OverloadedLabels, OverloadedStrings #-}
{-# LANGUAGE PatternSynonyms, PolyKinds, RecordWildCards            #-}
{-# LANGUAGE ScopedTypeVariables, TupleSections, TypeApplications   #-}
{-# LANGUAGE ViewPatterns                                           #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
module Algebra.Ring.Polynomial.Factorise
  ( -- * Factorisation
    factorise, factorQBigPrime, factorHensel,
    isReducible,
    -- * Internal helper functions
    distinctDegFactor,
    equalDegreeSplitM, equalDegreeFactorM,
    henselStep, multiHensel, clearDenom,
    squareFreePart, pthRoot,
    yun, squareFreeDecompFiniteField
  ) where
import           Algebra.Arithmetic                 hiding (modPow)
import           Algebra.Field.Prime
import           Algebra.Prelude.Core
import           Algebra.Ring.Euclidean.Quotient
import           Algebra.Ring.Polynomial.Univariate
import           Control.Applicative                ((<|>))
import           Control.Lens                       (both, ifoldMap, (%~), (&))
import           Control.Monad                      (guard, replicateM)
import           Control.Monad.Loops                (untilJust)
import           Control.Monad.Random               (MonadRandom, Random,
                                                     getRandom, getRandomR)
import qualified Data.DList                         as DL
import qualified Data.FMList                        as FML
import           Data.IntMap                        (IntMap)
import qualified Data.IntMap.Strict                 as IM
import           Data.Maybe                         (fromJust)
import           Data.Monoid                        ((<>))
import           Data.Proxy                         (Proxy (..))
import qualified Data.Set                           as S
import qualified Data.Traversable                   as F
import qualified Data.Vector                        as V
import           Math.NumberTheory.Primes           (Prime, precPrime)
import           Math.NumberTheory.Primes           (Prime (unPrime))
import qualified Math.NumberTheory.Primes           as PRIMES
import           Numeric.Decidable.Zero             (isZero)
import           Numeric.Domain.GCD                 (gcd, lcm)
import qualified Numeric.Field.Fraction             as F
import qualified Prelude                            as P


-- | @distinctDegFactor f@ computes the distinct-degree decomposition of the given
--   square-free polynomial over finite field @f@.
distinctDegFactor :: forall k. (Eq k, FiniteField k)
                  => Unipol k     -- ^ Square-free polynomial over finite field.
                  -> [(Natural, Unipol k)]   -- ^ Distinct-degree decomposition.
distinctDegFactor :: Unipol k -> [(Natural, Unipol k)]
distinctDegFactor Unipol k
f0 = [Natural] -> [Unipol k] -> [(Natural, Unipol k)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Natural
1..] ([Unipol k] -> [(Natural, Unipol k)])
-> [Unipol k] -> [(Natural, Unipol k)]
forall a b. (a -> b) -> a -> b
$ ([Unipol k] -> [Unipol k])
-> Unipol k -> Unipol k -> [Unipol k] -> [Unipol k]
go [Unipol k] -> [Unipol k]
forall a. a -> a
id (Ordinal (Arity (Unipol k)) -> Unipol k
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol k))
forall (n :: Nat). (0 < n) => Ordinal n
OZ :: Unipol k) Unipol k
f0 []
  where
    go :: ([Unipol k] -> [Unipol k])
-> Unipol k -> Unipol k -> [Unipol k] -> [Unipol k]
go [Unipol k] -> [Unipol k]
gs Unipol k
h Unipol k
f =
      let h' :: Unipol k
h' = Unipol k -> Natural -> Unipol k -> Unipol k
forall poly. Euclidean poly => poly -> Natural -> poly -> poly
modPow Unipol k
h (Proxy k -> Natural
forall k (proxy :: * -> *). FiniteField k => proxy k -> Natural
order (Proxy k
forall k (t :: k). Proxy t
Proxy :: Proxy k)) Unipol k
f
          g' :: Unipol k
g' = Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd (Unipol k
h' Unipol k -> Unipol k -> Unipol k
forall r. Group r => r -> r -> r
- Ordinal (Arity (Unipol k)) -> Unipol k
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol k))
0) Unipol k
f
          f' :: Unipol k
f' = Unipol k
f Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
g'
          gs' :: [Unipol k] -> [Unipol k]
gs' = [Unipol k] -> [Unipol k]
gs ([Unipol k] -> [Unipol k])
-> ([Unipol k] -> [Unipol k]) -> [Unipol k] -> [Unipol k]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Unipol k
g' Unipol k -> [Unipol k] -> [Unipol k]
forall a. a -> [a] -> [a]
:)
      in if Unipol k
f' Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
== Unipol k
forall r. Unital r => r
one
         then [Unipol k] -> [Unipol k]
gs'
         else ([Unipol k] -> [Unipol k])
-> Unipol k -> Unipol k -> [Unipol k] -> [Unipol k]
go [Unipol k] -> [Unipol k]
gs' Unipol k
h' Unipol k
f'

-- | @'modPow' f n g@ computes \(f^n \mod{g}\).
modPow :: (Euclidean poly)
       => poly -> Natural -> poly -> poly
modPow :: poly -> Natural -> poly -> poly
modPow poly
a Natural
p poly
f = poly -> (forall r. Reifies r poly => Quotient r poly) -> poly
forall a.
Euclidean a =>
a -> (forall r. Reifies r a => Quotient r a) -> a
withQuotient poly
f ((forall r. Reifies r poly => Quotient r poly) -> poly)
-> (forall r. Reifies r poly => Quotient r poly) -> poly
forall a b. (a -> b) -> a -> b
$
  Quotient r poly -> Natural -> Quotient r poly
forall r. Multiplicative r => r -> Natural -> r
repeatedSquare (poly -> Quotient r poly
forall k (r :: k) a.
(Reifies r a, Euclidean a) =>
a -> Quotient r a
quotient poly
a) Natural
p

traceCharTwo :: (Unital m, Monoidal m) => Natural -> m -> m
traceCharTwo :: Natural -> m -> m
traceCharTwo Natural
m m
a = (m -> m -> m) -> m -> [m] -> m
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' m -> m -> m
forall r. Additive r => r -> r -> r
(+) m
forall m. Monoidal m => m
zero ([m] -> m) -> [m] -> m
forall a b. (a -> b) -> a -> b
$ Int -> [m] -> [m]
forall a. Int -> [a] -> [a]
take (Natural -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
m) ([m] -> [m]) -> [m] -> [m]
forall a b. (a -> b) -> a -> b
$ (m -> m) -> m -> [m]
forall a. (a -> a) -> a -> [a]
iterate (\(!m
x) ->m
xm -> m -> m
forall r. Multiplicative r => r -> r -> r
*m
x) m
a

equalDegreeSplitM
  :: forall k m. (MonadRandom m, Random k, CoeffRing k,  FiniteField k)
  => Unipol k
  -> Natural
  -> m (Maybe (Unipol k))
equalDegreeSplitM :: Unipol k -> Natural -> m (Maybe (Unipol k))
equalDegreeSplitM Unipol k
f Natural
d
  | Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Natural -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0 = Maybe (Unipol k) -> m (Maybe (Unipol k))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Unipol k)
forall a. Maybe a
Nothing
  | Bool
otherwise = do
    let q :: Natural
q = Natural -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Natural -> Natural) -> Natural -> Natural
forall a b. (a -> b) -> a -> b
$ Proxy k -> Natural
forall k (proxy :: * -> *). FiniteField k => proxy k -> Natural
order (Proxy k
forall k (t :: k). Proxy t
Proxy :: Proxy k)
    Int
e <- (Int, Int) -> m Int
forall (m :: * -> *) a. (MonadRandom m, Random a) => (a, a) -> m a
getRandomR (Int
1, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
P.- Int
1)
    [k]
cs <- Int -> m k -> m [k]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM (Int -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
e) m k
forall (m :: * -> *) a. (MonadRandom m, Random a) => m a
getRandom
    let a :: Unipol k
a = Ordinal (Arity (Unipol k)) -> Unipol k
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol k))
0 Unipol k -> Natural -> Unipol k
forall r. Unital r => r -> Natural -> r
^ Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
e Unipol k -> Unipol k -> Unipol k
forall r. Additive r => r -> r -> r
+
            [Unipol k] -> Unipol k
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum ((Unipol k -> Unipol k -> Unipol k)
-> [Unipol k] -> [Unipol k] -> [Unipol k]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Unipol k -> Unipol k -> Unipol k
forall r. Multiplicative r => r -> r -> r
(*) ((k -> Unipol k) -> [k] -> [Unipol k]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map k -> Unipol k
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff [k]
cs) [Ordinal (Arity (Unipol k)) -> Unipol k
forall poly. IsPolynomial poly => Ordinal (Arity poly) -> poly
var Ordinal (Arity (Unipol k))
0 Unipol k -> Natural -> Unipol k
forall r. Unital r => r -> Natural -> r
^ Natural
l | Natural
l <-[Natural
0..]])
        g1 :: Unipol k
g1 = Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd Unipol k
a Unipol k
f
    Maybe (Unipol k) -> m (Maybe (Unipol k))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Unipol k) -> m (Maybe (Unipol k)))
-> Maybe (Unipol k) -> m (Maybe (Unipol k))
forall a b. (a -> b) -> a -> b
$ (Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Unipol k
g1 Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= Unipol k
forall r. Unital r => r
one) Maybe () -> Maybe (Unipol k) -> Maybe (Unipol k)
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Unipol k -> Maybe (Unipol k)
forall (m :: * -> *) a. Monad m => a -> m a
return Unipol k
g1)
         Maybe (Unipol k) -> Maybe (Unipol k) -> Maybe (Unipol k)
forall (f :: * -> *) a. Alternative f => f a -> f a -> f a
<|> do let b :: Unipol k
b | Unipol k -> Natural
forall r. Characteristic r => Unipol r -> Natural
charUnipol Unipol k
f Natural -> Natural -> Bool
forall a. Eq a => a -> a -> Bool
== Natural
2  =
                          Unipol k
-> (forall r. Reifies r (Unipol k) => Quotient r (Unipol k))
-> Unipol k
forall a.
Euclidean a =>
a -> (forall r. Reifies r a => Quotient r a) -> a
withQuotient Unipol k
f ((forall r. Reifies r (Unipol k) => Quotient r (Unipol k))
 -> Unipol k)
-> (forall r. Reifies r (Unipol k) => Quotient r (Unipol k))
-> Unipol k
forall a b. (a -> b) -> a -> b
$
                          Natural -> Quotient r (Unipol k) -> Quotient r (Unipol k)
forall m. (Unital m, Monoidal m) => Natural -> m -> m
traceCharTwo (Unipol k -> Natural
forall r. FiniteField r => Unipol r -> Natural
powerUnipol Unipol k
fNatural -> Natural -> Natural
forall r. Multiplicative r => r -> r -> r
*Natural
d) (Unipol k -> Quotient r (Unipol k)
forall k (r :: k) a.
(Reifies r a, Euclidean a) =>
a -> Quotient r a
quotient Unipol k
a)
                      | Bool
otherwise = Unipol k -> Natural -> Unipol k -> Unipol k
forall poly. Euclidean poly => poly -> Natural -> poly -> poly
modPow Unipol k
a (Natural -> Natural
forall a. Enum a => a -> a
pred (Natural
qNatural -> Natural -> Natural
forall r. Unital r => r -> Natural -> r
^Natural
d) Natural -> Natural -> Natural
forall a. Integral a => a -> a -> a
`div`Natural
2) Unipol k
f
                    g2 :: Unipol k
g2 = Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd (Unipol k
b Unipol k -> Unipol k -> Unipol k
forall r. Group r => r -> r -> r
- Unipol k
forall r. Unital r => r
one) Unipol k
f
                Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Unipol k
g2 Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= Unipol k
forall r. Unital r => r
one Bool -> Bool -> Bool
&& Unipol k
g2 Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= Unipol k
f)
                Unipol k -> Maybe (Unipol k)
forall (m :: * -> *) a. Monad m => a -> m a
return Unipol k
g2
  where n :: Int
n = Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol k
f

equalDegreeFactorM :: (Eq k, FiniteField k, MonadRandom m, Random k)
                   => Unipol k -> Natural -> m [Unipol k]
equalDegreeFactorM :: Unipol k -> Natural -> m [Unipol k]
equalDegreeFactorM Unipol k
f Natural
d = Unipol k -> m ([Unipol k] -> [Unipol k])
go Unipol k
f m ([Unipol k] -> [Unipol k])
-> (([Unipol k] -> [Unipol k]) -> m [Unipol k]) -> m [Unipol k]
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \[Unipol k] -> [Unipol k]
a -> [Unipol k] -> m [Unipol k]
forall (m :: * -> *) a. Monad m => a -> m a
return ([Unipol k] -> [Unipol k]
a [])
  where
    go :: Unipol k -> m ([Unipol k] -> [Unipol k])
go Unipol k
h
      | Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol k
h Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = ([Unipol k] -> [Unipol k]) -> m ([Unipol k] -> [Unipol k])
forall (m :: * -> *) a. Monad m => a -> m a
return [Unipol k] -> [Unipol k]
forall a. a -> a
id
      | Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol k
h) Natural -> Natural -> Bool
forall a. Eq a => a -> a -> Bool
== Natural
d = ([Unipol k] -> [Unipol k]) -> m ([Unipol k] -> [Unipol k])
forall (m :: * -> *) a. Monad m => a -> m a
return (Unipol k
hUnipol k -> [Unipol k] -> [Unipol k]
forall a. a -> [a] -> [a]
:)
      | Bool
otherwise = do
          Unipol k
g <- m (Maybe (Unipol k)) -> m (Unipol k)
forall (m :: * -> *) a. Monad m => m (Maybe a) -> m a
untilJust (Unipol k -> Natural -> m (Maybe (Unipol k))
forall k (m :: * -> *).
(MonadRandom m, Random k, CoeffRing k, FiniteField k) =>
Unipol k -> Natural -> m (Maybe (Unipol k))
equalDegreeSplitM Unipol k
h Natural
d)
          [Unipol k] -> [Unipol k]
l <- Unipol k -> m ([Unipol k] -> [Unipol k])
go Unipol k
g
          [Unipol k] -> [Unipol k]
r <- Unipol k -> m ([Unipol k] -> [Unipol k])
go (Unipol k
h Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
g)
          ([Unipol k] -> [Unipol k]) -> m ([Unipol k] -> [Unipol k])
forall (m :: * -> *) a. Monad m => a -> m a
return (([Unipol k] -> [Unipol k]) -> m ([Unipol k] -> [Unipol k]))
-> ([Unipol k] -> [Unipol k]) -> m ([Unipol k] -> [Unipol k])
forall a b. (a -> b) -> a -> b
$ [Unipol k] -> [Unipol k]
l ([Unipol k] -> [Unipol k])
-> ([Unipol k] -> [Unipol k]) -> [Unipol k] -> [Unipol k]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Unipol k] -> [Unipol k]
r

factorSquareFree :: (Eq k, Random k, FiniteField k, MonadRandom m)
                 => Unipol k -> m [Unipol k]
factorSquareFree :: Unipol k -> m [Unipol k]
factorSquareFree Unipol k
f =
   [[Unipol k]] -> [Unipol k]
forall w. Monoid w => [w] -> w
concat ([[Unipol k]] -> [Unipol k]) -> m [[Unipol k]] -> m [Unipol k]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((Natural, Unipol k) -> m [Unipol k])
-> [(Natural, Unipol k)] -> m [[Unipol k]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM ((Natural -> Unipol k -> m [Unipol k])
-> (Natural, Unipol k) -> m [Unipol k]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Natural -> Unipol k -> m [Unipol k])
 -> (Natural, Unipol k) -> m [Unipol k])
-> (Natural -> Unipol k -> m [Unipol k])
-> (Natural, Unipol k)
-> m [Unipol k]
forall a b. (a -> b) -> a -> b
$ (Unipol k -> Natural -> m [Unipol k])
-> Natural -> Unipol k -> m [Unipol k]
forall a b c. (a -> b -> c) -> b -> a -> c
flip Unipol k -> Natural -> m [Unipol k]
forall k (m :: * -> *).
(Eq k, FiniteField k, MonadRandom m, Random k) =>
Unipol k -> Natural -> m [Unipol k]
equalDegreeFactorM) (((Natural, Unipol k) -> Bool)
-> [(Natural, Unipol k)] -> [(Natural, Unipol k)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= Unipol k
forall r. Unital r => r
one) (Unipol k -> Bool)
-> ((Natural, Unipol k) -> Unipol k) -> (Natural, Unipol k) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Natural, Unipol k) -> Unipol k
forall a b. (a, b) -> b
snd) ([(Natural, Unipol k)] -> [(Natural, Unipol k)])
-> [(Natural, Unipol k)] -> [(Natural, Unipol k)]
forall a b. (a -> b) -> a -> b
$ Unipol k -> [(Natural, Unipol k)]
forall k.
(Eq k, FiniteField k) =>
Unipol k -> [(Natural, Unipol k)]
distinctDegFactor Unipol k
f)

squareFreePart :: (Eq k, FiniteField k)
               => Unipol k -> Unipol k
squareFreePart :: Unipol k -> Unipol k
squareFreePart Unipol k
f =
  let !n :: Natural
n = Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol k
f
      u :: Unipol k
u  = Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd Unipol k
f (Ordinal (Arity (Unipol k)) -> Unipol k -> Unipol k
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity (Unipol k))
0 Unipol k
f)
      v :: Unipol k
v  = Unipol k
f Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
u
      f' :: Unipol k
f' = Unipol k
u Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd Unipol k
u (Unipol k
v Unipol k -> Natural -> Unipol k
forall r. Unital r => r -> Natural -> r
^ Natural
n)
  in if Unipol k
f' Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
== Unipol k
forall r. Unital r => r
one
     then Unipol k
v
     else Unipol k
v Unipol k -> Unipol k -> Unipol k
forall r. Multiplicative r => r -> r -> r
* Unipol k -> Unipol k
forall k. (Eq k, FiniteField k) => Unipol k -> Unipol k
squareFreePart (Unipol k -> Unipol k
forall r. (CoeffRing r, FiniteField r) => Unipol r -> Unipol r
pthRoot Unipol k
f')

-- | Yun's method to compute square-free decomposition of
--   a univariate polynomial over field of characteristic \(0\).
--
--   Use 'squareFreeDecompFiniteField' for finite fields.
yun :: (CoeffRing r, Field r)
    => Unipol r -> IntMap (Unipol r)
yun :: Unipol r -> IntMap (Unipol r)
yun Unipol r
f = let f' :: Unipol r
f' = Ordinal (Arity (Unipol r)) -> Unipol r -> Unipol r
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity (Unipol r))
forall (n :: Nat). (0 < n) => Ordinal n
OZ Unipol r
f
            u :: Unipol r
u  = Unipol r -> Unipol r -> Unipol r
forall d. GCDDomain d => d -> d -> d
gcd Unipol r
f Unipol r
f'
        in Int
-> IntMap (Unipol r) -> Unipol r -> Unipol r -> IntMap (Unipol r)
forall t.
(IsOrderedPolynomial t, Euclidean t, (1 <=? Arity t) ~ 'True) =>
Int -> IntMap t -> t -> t -> IntMap t
go Int
1 IntMap (Unipol r)
forall a. IntMap a
IM.empty (Unipol r
f Unipol r -> Unipol r -> Unipol r
forall d. Euclidean d => d -> d -> d
`quot` Unipol r
u) (Unipol r
f' Unipol r -> Unipol r -> Unipol r
forall d. Euclidean d => d -> d -> d
`quot` Unipol r
u)
  where
    go :: Int -> IntMap t -> t -> t -> IntMap t
go !Int
i IntMap t
dic t
v t
w =
      let t :: t
t  = t
w t -> t -> t
forall r. Group r => r -> r -> r
- Ordinal (Arity t) -> t -> t
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity t)
forall (n :: Nat). (0 < n) => Ordinal n
OZ t
v
          h :: t
h  = t -> t -> t
forall d. GCDDomain d => d -> d -> d
gcd t
v t
t
          v' :: t
v' = t
v t -> t -> t
forall d. Euclidean d => d -> d -> d
`quot` t
h
          w' :: t
w' = t
t t -> t -> t
forall d. Euclidean d => d -> d -> d
`quot` t
h
          dic' :: IntMap t
dic' = Int -> t -> IntMap t -> IntMap t
forall a. Int -> a -> IntMap a -> IntMap a
IM.insert Int
i t
h IntMap t
dic
      in if t
v' t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
forall r. Unital r => r
one
         then IntMap t
dic'
         else Int -> IntMap t -> t -> t -> IntMap t
go (Int
iInt -> Int -> Int
forall r. Additive r => r -> r -> r
+Int
1) IntMap t
dic' t
v' t
w'

charUnipol :: forall r. Characteristic r => Unipol r -> Natural
charUnipol :: Unipol r -> Natural
charUnipol Unipol r
_ = Proxy r -> Natural
forall r (proxy :: * -> *). Characteristic r => proxy r -> Natural
char (Proxy r
forall k (t :: k). Proxy t
Proxy :: Proxy r)

powerUnipol :: forall r. FiniteField r => Unipol r -> Natural
powerUnipol :: Unipol r -> Natural
powerUnipol Unipol r
_ = Proxy r -> Natural
forall k (proxy :: * -> *). FiniteField k => proxy k -> Natural
power (Proxy r
forall k (t :: k). Proxy t
Proxy :: Proxy r)

pthRoot :: forall r. (CoeffRing r, FiniteField r) => Unipol r -> Unipol r
pthRoot :: Unipol r -> Unipol r
pthRoot Unipol r
f =
  let !p :: Natural
p = Unipol r -> Natural
forall r. Characteristic r => Unipol r -> Natural
charUnipol Unipol r
f
      !q :: Natural
q = Proxy r -> Natural
forall k (proxy :: * -> *). FiniteField k => proxy k -> Natural
order (Proxy r -> Natural) -> Proxy r -> Natural
forall a b. (a -> b) -> a -> b
$ Proxy r
forall k (t :: k). Proxy t
Proxy @r
      !s :: Natural
s = Natural
q Natural -> Natural -> Natural
forall a. Integral a => a -> a -> a
`div` Natural
p
  in [Unipol r] -> Unipol r
forall (f :: * -> *) m. (Foldable f, Monoidal m) => f m -> m
sum  [ (r
ar -> Natural -> r
forall r. Unital r => r -> Natural -> r
^Natural
s) r -> Unipol r -> Unipol r
forall r m. Module (Scalar r) m => r -> m -> m
.*. IsLabel "x" (Unipol r)
Unipol r
#xUnipol r -> Natural -> Unipol r
forall r. Unital r => r -> Natural -> r
^Natural
i
          | (Int
ip, r
a) <- IntMap r -> [(Int, r)]
forall a. IntMap a -> [(Int, a)]
IM.toList (IntMap r -> [(Int, r)]) -> IntMap r -> [(Int, r)]
forall a b. (a -> b) -> a -> b
$ Unipol r -> IntMap r
forall k. Unipol k -> IntMap k
termsUnipol Unipol r
f
          , (Natural
i, Natural
0) <- [Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ip Natural -> Natural -> (Natural, Natural)
forall a. Integral a => a -> a -> (a, a)
`P.divMod` Natural
p]
          ]

-- | Square-free decomposition of polynomials over a finite field.
--
--   Use 'yun' for a polynomials over a field of characteristic zero.
squareFreeDecompFiniteField
  :: (Eq k, FiniteField k)
  => Unipol k -> IntMap (Unipol k)
squareFreeDecompFiniteField :: Unipol k -> IntMap (Unipol k)
squareFreeDecompFiniteField Unipol k
f =
  let dcmp :: IntMap (Unipol k)
dcmp = Unipol k -> IntMap (Unipol k)
forall r. (CoeffRing r, Field r) => Unipol r -> IntMap (Unipol r)
yun Unipol k
f
      h' :: Unipol k
h'   = Mult (Unipol k) -> Unipol k
forall a. Mult a -> a
runMult (Mult (Unipol k) -> Unipol k) -> Mult (Unipol k) -> Unipol k
forall a b. (a -> b) -> a -> b
$ (Int -> Unipol k -> Mult (Unipol k))
-> IntMap (Unipol k) -> Mult (Unipol k)
forall i (f :: * -> *) m a.
(FoldableWithIndex i f, Monoid m) =>
(i -> a -> m) -> f a -> m
ifoldMap (\Int
i Unipol k
g -> Unipol k -> Mult (Unipol k)
forall a. a -> Mult a
Mult (Unipol k -> Mult (Unipol k)) -> Unipol k -> Mult (Unipol k)
forall a b. (a -> b) -> a -> b
$ Unipol k
g Unipol k -> Natural -> Unipol k
forall r. Unital r => r -> Natural -> r
^ Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) IntMap (Unipol k)
dcmp
      fhInv :: Unipol k
fhInv = Unipol k
f Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
h'
      p :: Int
p    = Natural -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Natural -> Int) -> Natural -> Int
forall a b. (a -> b) -> a -> b
$ Unipol k -> Natural
forall r. Characteristic r => Unipol r -> Natural
charUnipol Unipol k
f
      powers :: IntMap (Unipol k)
powers =
        (Int -> Int) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. (Int -> Int) -> IntMap a -> IntMap a
IM.mapKeysMonotonic (Int
pInt -> Int -> Int
forall r. Multiplicative r => r -> r -> r
*) (IntMap (Unipol k) -> IntMap (Unipol k))
-> IntMap (Unipol k) -> IntMap (Unipol k)
forall a b. (a -> b) -> a -> b
$
        Unipol k -> IntMap (Unipol k)
forall k. (Eq k, FiniteField k) => Unipol k -> IntMap (Unipol k)
squareFreeDecompFiniteField (Unipol k -> IntMap (Unipol k)) -> Unipol k -> IntMap (Unipol k)
forall a b. (a -> b) -> a -> b
$ Unipol k -> Unipol k
forall r. (CoeffRing r, FiniteField r) => Unipol r -> Unipol r
pthRoot Unipol k
fhInv
      ((IntMap (Unipol k)
pns, IntMap (Unipol k)
pnjs), IntMap (Unipol k)
js) =
        ((IntMap (Unipol k), IntMap (Unipol k))
 -> Int
 -> Unipol k
 -> ((IntMap (Unipol k), IntMap (Unipol k)), Unipol k))
-> (IntMap (Unipol k), IntMap (Unipol k))
-> IntMap (Unipol k)
-> ((IntMap (Unipol k), IntMap (Unipol k)), IntMap (Unipol k))
forall a b c.
(a -> Int -> b -> (a, c)) -> a -> IntMap b -> (a, IntMap c)
IM.mapAccumWithKey
            (\(IntMap (Unipol k)
ps', IntMap (Unipol k)
pnjs') Int
j Unipol k
g0 ->
                let (Unipol k
g', IntMap (Unipol k, Unipol k)
pnjsps'') =
                      (Unipol k -> Unipol k -> (Unipol k, (Unipol k, Unipol k)))
-> Unipol k
-> IntMap (Unipol k)
-> (Unipol k, IntMap (Unipol k, Unipol k))
forall a b c. (a -> b -> (a, c)) -> a -> IntMap b -> (a, IntMap c)
IM.mapAccum
                        (\Unipol k
g Unipol k
h ->
                          let u :: Unipol k
u = Unipol k
g Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
`gcd` Unipol k
h
                              h' :: Unipol k
h' = Unipol k
h Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
u
                          in (Unipol k
g Unipol k -> Unipol k -> Unipol k
forall d. Euclidean d => d -> d -> d
`quot` Unipol k
u, (Unipol k
u, Unipol k
h'))
                        )
                      Unipol k
g0 IntMap (Unipol k)
ps'
                    pnjs'' :: IntMap (Unipol k)
pnjs'' = (Unipol k -> Bool) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter ((Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0) (Int -> Bool) -> (Unipol k -> Int) -> Unipol k -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree')
                      (IntMap (Unipol k) -> IntMap (Unipol k))
-> IntMap (Unipol k) -> IntMap (Unipol k)
forall a b. (a -> b) -> a -> b
$ (Int -> Int) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. (Int -> Int) -> IntMap a -> IntMap a
IM.mapKeysMonotonic (Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
j)
                      (IntMap (Unipol k) -> IntMap (Unipol k))
-> IntMap (Unipol k) -> IntMap (Unipol k)
forall a b. (a -> b) -> a -> b
$ (Unipol k, Unipol k) -> Unipol k
forall a b. (a, b) -> a
fst ((Unipol k, Unipol k) -> Unipol k)
-> IntMap (Unipol k, Unipol k) -> IntMap (Unipol k)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> IntMap (Unipol k, Unipol k)
pnjsps''
                    ps'' :: IntMap (Unipol k)
ps'' = (Unipol k -> Bool) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter ((Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0) (Int -> Bool) -> (Unipol k -> Int) -> Unipol k -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree')
                        (IntMap (Unipol k) -> IntMap (Unipol k))
-> IntMap (Unipol k) -> IntMap (Unipol k)
forall a b. (a -> b) -> a -> b
$ (Unipol k, Unipol k) -> Unipol k
forall a b. (a, b) -> b
snd ((Unipol k, Unipol k) -> Unipol k)
-> IntMap (Unipol k, Unipol k) -> IntMap (Unipol k)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> IntMap (Unipol k, Unipol k)
pnjsps''
                in (( IntMap (Unipol k)
ps'', IntMap (Unipol k)
pnjs' IntMap (Unipol k) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. Semigroup a => a -> a -> a
<> IntMap (Unipol k)
pnjs''), Unipol k
g')
            )
            (IntMap (Unipol k)
powers, IntMap (Unipol k)
forall a. IntMap a
IM.empty)
            (IntMap (Unipol k)
 -> ((IntMap (Unipol k), IntMap (Unipol k)), IntMap (Unipol k)))
-> IntMap (Unipol k)
-> ((IntMap (Unipol k), IntMap (Unipol k)), IntMap (Unipol k))
forall a b. (a -> b) -> a -> b
$ (IntMap (Unipol k), IntMap (Unipol k)) -> IntMap (Unipol k)
forall a b. (a, b) -> a
fst (Int -> IntMap (Unipol k) -> (IntMap (Unipol k), IntMap (Unipol k))
forall a. Int -> IntMap a -> (IntMap a, IntMap a)
IM.split Int
p IntMap (Unipol k)
dcmp)
  in if (Int, Unipol k) -> Int
forall a b. (a, b) -> a
fst (IntMap (Unipol k) -> (Int, Unipol k)
forall a. IntMap a -> (Int, a)
IM.findMax IntMap (Unipol k)
dcmp) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
p Bool -> Bool -> Bool
&& Unipol k
fhInv Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
== Unipol k
forall r. Unital r => r
one
    then IntMap (Unipol k)
dcmp
    else [IntMap (Unipol k)] -> IntMap (Unipol k)
forall (f :: * -> *) a. Foldable f => f (IntMap a) -> IntMap a
IM.unions [(Unipol k -> Bool) -> IntMap (Unipol k) -> IntMap (Unipol k)
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter ((Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
0) (Int -> Bool) -> (Unipol k -> Int) -> Unipol k -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree') IntMap (Unipol k)
js, IntMap (Unipol k)
pnjs, IntMap (Unipol k)
pns]


newtype GCD a = GCD { GCD a -> a
runGCD :: a }

instance Euclidean a => P.Semigroup (GCD a) where
  <> :: GCD a -> GCD a -> GCD a
(<>) = (a -> GCD a
forall a. a -> GCD a
GCD (a -> GCD a) -> (GCD a -> a) -> GCD a -> GCD a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.) ((GCD a -> a) -> GCD a -> GCD a)
-> (GCD a -> GCD a -> a) -> GCD a -> GCD a -> GCD a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall d. GCDDomain d => d -> d -> d
gcd (a -> a -> a) -> (GCD a -> a) -> GCD a -> GCD a -> a
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` GCD a -> a
forall a. GCD a -> a
runGCD)
instance Euclidean a => Monoid (GCD a) where
  mappend :: GCD a -> GCD a -> GCD a
mappend = GCD a -> GCD a -> GCD a
forall a. Semigroup a => a -> a -> a
(<>)
  mempty :: GCD a
mempty = a -> GCD a
forall a. a -> GCD a
GCD a
forall m. Monoidal m => m
zero

newtype LCM a = LCM { LCM a -> a
runLCM :: a }

instance Euclidean a => P.Semigroup (LCM a) where
  <> :: LCM a -> LCM a -> LCM a
(<>) = (a -> LCM a
forall a. a -> LCM a
LCM (a -> LCM a) -> (LCM a -> a) -> LCM a -> LCM a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.) ((LCM a -> a) -> LCM a -> LCM a)
-> (LCM a -> LCM a -> a) -> LCM a -> LCM a -> LCM a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall d. GCDDomain d => d -> d -> d
lcm (a -> a -> a) -> (LCM a -> a) -> LCM a -> LCM a -> a
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` LCM a -> a
forall a. LCM a -> a
runLCM)
instance Euclidean a => Monoid (LCM a) where
  mappend :: LCM a -> LCM a -> LCM a
mappend = LCM a -> LCM a -> LCM a
forall a. Semigroup a => a -> a -> a
(<>)
  mempty :: LCM a
mempty = a -> LCM a
forall a. a -> LCM a
LCM a
forall r. Unital r => r
one

-- | Factorise a polynomial over finite field using Cantor-Zassenhaus algorithm
factorise :: (MonadRandom m, CoeffRing k, FiniteField k, Random k)
          => Unipol k -> m [(Unipol k, Natural)]
factorise :: Unipol k -> m [(Unipol k, Natural)]
factorise Unipol k
f0
  | Unipol k -> Bool
forall r. DecidableZero r => r -> Bool
isZero Unipol k
f0 = [(Unipol k, Natural)] -> m [(Unipol k, Natural)]
forall (f :: * -> *) a. Applicative f => a -> f a
pure [(Unipol k
forall m. Monoidal m => m
zero, Natural
1)]
  | Bool
otherwise = do
  let f :: Unipol k
f = Unipol k -> Unipol k
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize Unipol k
f0
      c :: Coefficient (Unipol k)
c = Unipol k -> Coefficient (Unipol k)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol k
f0
  [(Unipol k, Natural)]
monoFacts <- [[(Unipol k, Natural)]] -> [(Unipol k, Natural)]
forall w. Monoid w => [w] -> w
concat
    ([[(Unipol k, Natural)]] -> [(Unipol k, Natural)])
-> m [[(Unipol k, Natural)]] -> m [(Unipol k, Natural)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((Int, Unipol k) -> m [(Unipol k, Natural)])
-> [(Int, Unipol k)] -> m [[(Unipol k, Natural)]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\(Int
r, Unipol k
h) -> (Unipol k -> (Unipol k, Natural))
-> [Unipol k] -> [(Unipol k, Natural)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (,Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r) ([Unipol k] -> [(Unipol k, Natural)])
-> m [Unipol k] -> m [(Unipol k, Natural)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Unipol k -> m [Unipol k]
forall k (m :: * -> *).
(Eq k, Random k, FiniteField k, MonadRandom m) =>
Unipol k -> m [Unipol k]
factorSquareFree Unipol k
h) (IntMap (Unipol k) -> [(Int, Unipol k)]
forall a. IntMap a -> [(Int, a)]
IM.toList (IntMap (Unipol k) -> [(Int, Unipol k)])
-> IntMap (Unipol k) -> [(Int, Unipol k)]
forall a b. (a -> b) -> a -> b
$  Unipol k -> IntMap (Unipol k)
forall k. (Eq k, FiniteField k) => Unipol k -> IntMap (Unipol k)
squareFreeDecompFiniteField Unipol k
f)
  [(Unipol k, Natural)] -> m [(Unipol k, Natural)]
forall (f :: * -> *) a. Applicative f => a -> f a
pure ([(Unipol k, Natural)] -> m [(Unipol k, Natural)])
-> [(Unipol k, Natural)] -> m [(Unipol k, Natural)]
forall a b. (a -> b) -> a -> b
$ if k
Coefficient (Unipol k)
c k -> k -> Bool
forall a. Eq a => a -> a -> Bool
== k
forall r. Unital r => r
one
    then [(Unipol k, Natural)]
monoFacts
    else (Coefficient (Unipol k) -> Unipol k
forall poly. IsPolynomial poly => Coefficient poly -> poly
injectCoeff Coefficient (Unipol k)
c, Natural
1) (Unipol k, Natural)
-> [(Unipol k, Natural)] -> [(Unipol k, Natural)]
forall a. a -> [a] -> [a]
: [(Unipol k, Natural)]
monoFacts

clearDenom :: (CoeffRing a, Euclidean a)
           => Unipol (Fraction a) -> (a, Unipol a)
clearDenom :: Unipol (Fraction a) -> (a, Unipol a)
clearDenom Unipol (Fraction a)
f =
  let g :: a
g = (Fraction a -> a -> a) -> a -> Map (Monomial 1) (Fraction a) -> a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (a -> a -> a
forall d. GCDDomain d => d -> d -> d
lcm (a -> a -> a) -> (Fraction a -> a) -> Fraction a -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Fraction a -> a
forall t. Fraction t -> t
denominator) a
forall r. Unital r => r
one (Map (Monomial 1) (Fraction a) -> a)
-> Map (Monomial 1) (Fraction a) -> a
forall a b. (a -> b) -> a -> b
$ Unipol (Fraction a)
-> Map
     (Monomial (Arity (Unipol (Fraction a))))
     (Coefficient (Unipol (Fraction a)))
forall poly.
IsPolynomial poly =>
poly -> Map (Monomial (Arity poly)) (Coefficient poly)
terms' Unipol (Fraction a)
f
  in (a
g, (Fraction a -> a) -> Unipol (Fraction a) -> Unipol a
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Fraction a -> a
forall t. Fraction t -> t
numerator (Fraction a -> a) -> (Fraction a -> Fraction a) -> Fraction a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((a
g a -> a -> Fraction a
forall d. GCDDomain d => d -> d -> Fraction d
F.% a
forall r. Unital r => r
one)Fraction a -> Fraction a -> Fraction a
forall r. Multiplicative r => r -> r -> r
*)) Unipol (Fraction a)
f)

-- | Factorise the given primitive and square-freeinteger-coefficient polynomial,
--   choosing a large enough prime.
factorQBigPrime :: (MonadRandom m)
               => Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
factorQBigPrime :: Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
factorQBigPrime = (Unipol Integer -> m [Unipol Integer])
-> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
forall (m :: * -> *).
MonadRandom m =>
(Unipol Integer -> m [Unipol Integer])
-> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
wrapSQFFactor Unipol Integer -> m [Unipol Integer]
forall (m :: * -> *).
MonadRandom m =>
Unipol Integer -> m [Unipol Integer]
factorSqFreeQBP

-- | Factorise the given
--   interger-coefficient polynomial by Hensel lifting.
factorHensel :: (MonadRandom m)
             => Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
factorHensel :: Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
factorHensel = (Unipol Integer -> m [Unipol Integer])
-> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
forall (m :: * -> *).
MonadRandom m =>
(Unipol Integer -> m [Unipol Integer])
-> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
wrapSQFFactor Unipol Integer -> m [Unipol Integer]
forall (m :: * -> *).
MonadRandom m =>
Unipol Integer -> m [Unipol Integer]
factorHenselSqFree

wrapSQFFactor :: (MonadRandom m)
              => (Unipol Integer -> m [Unipol Integer])
              -> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
wrapSQFFactor :: (Unipol Integer -> m [Unipol Integer])
-> Unipol Integer -> m (Integer, IntMap (Set (Unipol Integer)))
wrapSQFFactor Unipol Integer -> m [Unipol Integer]
fac Unipol Integer
f0
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = (Integer, IntMap (Set (Unipol Integer)))
-> m (Integer, IntMap (Set (Unipol Integer)))
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f0, IntMap (Set (Unipol Integer))
forall a. IntMap a
IM.empty)
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Integer, IntMap (Set (Unipol Integer)))
-> m (Integer, IntMap (Set (Unipol Integer)))
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Unipol Integer -> Coefficient (Unipol Integer)
forall poly.
(IsPolynomial poly, Euclidean (Coefficient poly)) =>
poly -> Coefficient poly
content Unipol Integer
f0, Int -> Set (Unipol Integer) -> IntMap (Set (Unipol Integer))
forall a. Int -> a -> IntMap a
IM.singleton Int
1 (Set (Unipol Integer) -> IntMap (Set (Unipol Integer)))
-> Set (Unipol Integer) -> IntMap (Set (Unipol Integer))
forall a b. (a -> b) -> a -> b
$ Unipol Integer -> Set (Unipol Integer)
forall a. a -> Set a
S.singleton (Unipol Integer -> Set (Unipol Integer))
-> Unipol Integer -> Set (Unipol Integer)
forall a b. (a -> b) -> a -> b
$ Unipol Integer -> Unipol Integer
forall poly.
(Euclidean (Coefficient poly), IsPolynomial poly) =>
poly -> poly
pp Unipol Integer
f0)
  | Bool
otherwise =  do
    let (Unipol Integer
g, Integer
c) = (Unipol Integer -> Unipol Integer
forall poly.
(Euclidean (Coefficient poly), IsPolynomial poly) =>
poly -> poly
pp Unipol Integer
f0, Unipol Integer -> Coefficient (Unipol Integer)
forall poly.
(IsPolynomial poly, Euclidean (Coefficient poly)) =>
poly -> Coefficient poly
content Unipol Integer
f0)
    IntMap (Integer, [Unipol Integer])
ts0 <- (Unipol (Fraction Integer) -> m (Integer, [Unipol Integer]))
-> IntMap (Unipol (Fraction Integer))
-> m (IntMap (Integer, [Unipol Integer]))
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
F.mapM ((Unipol Integer -> m [Unipol Integer])
-> (Integer, Unipol Integer) -> m (Integer, [Unipol Integer])
forall (f :: * -> *) t a t1.
Functor f =>
(t -> f a) -> (t1, t) -> f (t1, a)
secondM Unipol Integer -> m [Unipol Integer]
fac ((Integer, Unipol Integer) -> m (Integer, [Unipol Integer]))
-> (Unipol (Fraction Integer) -> (Integer, Unipol Integer))
-> Unipol (Fraction Integer)
-> m (Integer, [Unipol Integer])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol (Fraction Integer) -> (Integer, Unipol Integer)
forall a.
(CoeffRing a, Euclidean a) =>
Unipol (Fraction a) -> (a, Unipol a)
clearDenom) (Unipol (Fraction Integer) -> IntMap (Unipol (Fraction Integer))
forall r. (CoeffRing r, Field r) => Unipol r -> IntMap (Unipol r)
yun (Unipol (Fraction Integer) -> IntMap (Unipol (Fraction Integer)))
-> Unipol (Fraction Integer) -> IntMap (Unipol (Fraction Integer))
forall a b. (a -> b) -> a -> b
$ Unipol (Fraction Integer) -> Unipol (Fraction Integer)
forall poly.
(Field (Coefficient poly), IsOrderedPolynomial poly) =>
poly -> poly
monoize (Unipol (Fraction Integer) -> Unipol (Fraction Integer))
-> Unipol (Fraction Integer) -> Unipol (Fraction Integer)
forall a b. (a -> b) -> a -> b
$ (Integer -> Fraction Integer)
-> Unipol Integer -> Unipol (Fraction Integer)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Integer -> Integer -> Fraction Integer
forall d. GCDDomain d => d -> d -> Fraction d
F.% Integer
1) Unipol Integer
g)
    let anss :: [(Int, (Integer, [Unipol Integer]))]
anss = IntMap (Integer, [Unipol Integer])
-> [(Int, (Integer, [Unipol Integer]))]
forall a. IntMap a -> [(Int, a)]
IM.toList IntMap (Integer, [Unipol Integer])
ts0
        k :: Integer
k = Integer
c Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
g Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` [Integer] -> Integer
forall (f :: * -> *) r. (Foldable f, Unital r) => f r -> r
product (((Int, (Integer, [Unipol Integer])) -> Integer)
-> [(Int, (Integer, [Unipol Integer]))] -> [Integer]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (\(Int
p,(Integer
i,[Unipol Integer]
_)) -> Integer
i Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
P.^ Int
p) [(Int, (Integer, [Unipol Integer]))]
anss)
    (Integer, IntMap (Set (Unipol Integer)))
-> m (Integer, IntMap (Set (Unipol Integer)))
forall (m :: * -> *) a. Monad m => a -> m a
return (Integer
k, [(Int, Set (Unipol Integer))] -> IntMap (Set (Unipol Integer))
forall a. [(Int, a)] -> IntMap a
IM.fromList ([(Int, Set (Unipol Integer))] -> IntMap (Set (Unipol Integer)))
-> [(Int, Set (Unipol Integer))] -> IntMap (Set (Unipol Integer))
forall a b. (a -> b) -> a -> b
$ ((Int, (Integer, [Unipol Integer])) -> (Int, Set (Unipol Integer)))
-> [(Int, (Integer, [Unipol Integer]))]
-> [(Int, Set (Unipol Integer))]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (((Integer, [Unipol Integer]) -> Set (Unipol Integer))
-> (Int, (Integer, [Unipol Integer]))
-> (Int, Set (Unipol Integer))
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second (((Integer, [Unipol Integer]) -> Set (Unipol Integer))
 -> (Int, (Integer, [Unipol Integer]))
 -> (Int, Set (Unipol Integer)))
-> ((Integer, [Unipol Integer]) -> Set (Unipol Integer))
-> (Int, (Integer, [Unipol Integer]))
-> (Int, Set (Unipol Integer))
forall a b. (a -> b) -> a -> b
$ [Unipol Integer] -> Set (Unipol Integer)
forall a. Ord a => [a] -> Set a
S.fromList ([Unipol Integer] -> Set (Unipol Integer))
-> ((Integer, [Unipol Integer]) -> [Unipol Integer])
-> (Integer, [Unipol Integer])
-> Set (Unipol Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, [Unipol Integer]) -> [Unipol Integer]
forall a b. (a, b) -> b
snd) [(Int, (Integer, [Unipol Integer]))]
anss)
  where
    n :: Int
n = Unipol Integer -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Integer
f0

secondM :: Functor f => (t -> f a) -> (t1, t) -> f (t1, a)
secondM :: (t -> f a) -> (t1, t) -> f (t1, a)
secondM t -> f a
f (t1
a, t
b)= (t1
a,) (a -> (t1, a)) -> f a -> f (t1, a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> f a
f t
b

-- | Factors square-free primitive integer-coefficient polynomial,
--   using big enough prime.
factorSqFreeQBP :: (MonadRandom m)
                => Unipol Integer -> m [Unipol Integer]
{-# INLINE factorSqFreeQBP #-}
factorSqFreeQBP :: Unipol Integer -> m [Unipol Integer]
factorSqFreeQBP =
  (Prime Integer -> Unipol Integer -> BasicConsts)
-> (BasicConsts -> Unipol Integer -> Maybe (Integer, Integer))
-> (BasicConsts -> Prime Integer -> Unipol Integer -> Integer)
-> (BasicConsts
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
forall (m :: * -> *) a.
MonadRandom m =>
(Prime Integer -> Unipol Integer -> a)
-> (a -> Unipol Integer -> Maybe (Integer, Integer))
-> (a -> Prime Integer -> Unipol Integer -> Integer)
-> (a
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
factorPrimitiveIntegerPolWith
    ((Unipol Integer -> BasicConsts)
-> Prime Integer -> Unipol Integer -> BasicConsts
forall a b. a -> b -> a
const Unipol Integer -> BasicConsts
calcBasicConsts)
    (\BasicConsts{Double
Integer
lc :: BasicConsts -> Integer
kB :: BasicConsts -> Double
kA :: BasicConsts -> Integer
lc :: Integer
kB :: Double
kA :: Integer
..} -> Maybe (Integer, Integer)
-> Unipol Integer -> Maybe (Integer, Integer)
forall a b. a -> b -> a
const (Maybe (Integer, Integer)
 -> Unipol Integer -> Maybe (Integer, Integer))
-> Maybe (Integer, Integer)
-> Unipol Integer
-> Maybe (Integer, Integer)
forall a b. (a -> b) -> a -> b
$
        (Integer, Integer) -> Maybe (Integer, Integer)
forall a. a -> Maybe a
Just (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ Double
2Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
*Double
kB, Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ Double
4Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
*Double
kB)
    )
    ((Prime Integer -> Unipol Integer -> Integer)
-> BasicConsts -> Prime Integer -> Unipol Integer -> Integer
forall a b. a -> b -> a
const ((Prime Integer -> Unipol Integer -> Integer)
 -> BasicConsts -> Prime Integer -> Unipol Integer -> Integer)
-> (Prime Integer -> Unipol Integer -> Integer)
-> BasicConsts
-> Prime Integer
-> Unipol Integer
-> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Unipol Integer -> Integer
forall a b. a -> b -> a
const (Integer -> Unipol Integer -> Integer)
-> (Prime Integer -> Integer)
-> Prime Integer
-> Unipol Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Prime Integer -> Integer
forall a. Prime a -> a
unPrime)
    ((BasicConsts
  -> Prime Integer
  -> Unipol Integer
  -> [Unipol Integer]
  -> [Unipol Integer])
 -> Unipol Integer -> m [Unipol Integer])
-> (BasicConsts
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ (Prime Integer
 -> Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
-> BasicConsts
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
forall a b. a -> b -> a
const ((Prime Integer
  -> Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
 -> BasicConsts
 -> Prime Integer
 -> Unipol Integer
 -> [Unipol Integer]
 -> [Unipol Integer])
-> (Prime Integer
    -> Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
-> BasicConsts
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ (Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
forall a b. a -> b -> a
const ((Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
 -> Prime Integer
 -> Unipol Integer
 -> [Unipol Integer]
 -> [Unipol Integer])
-> (Unipol Integer -> [Unipol Integer] -> [Unipol Integer])
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ ([Unipol Integer] -> [Unipol Integer])
-> Unipol Integer -> [Unipol Integer] -> [Unipol Integer]
forall a b. a -> b -> a
const [Unipol Integer] -> [Unipol Integer]
forall a. a -> a
id

data BasicConsts =
  BasicConsts
  { BasicConsts -> Integer
kA :: !Integer
  , BasicConsts -> Double
kB :: !Double
  , BasicConsts -> Integer
lc :: !Integer
  }

data HenselConsts = HenselConsts
  { HenselConsts -> BasicConsts
basicConsts :: !BasicConsts
  , HenselConsts -> Natural
primePower  :: Natural
  }

factorHenselSqFree
  :: MonadRandom m
  => Unipol Integer -> m [Unipol Integer]
{-# INLINE factorHenselSqFree #-}
factorHenselSqFree :: Unipol Integer -> m [Unipol Integer]
factorHenselSqFree =
  (Prime Integer -> Unipol Integer -> HenselConsts)
-> (HenselConsts -> Unipol Integer -> Maybe (Integer, Integer))
-> (HenselConsts -> Prime Integer -> Unipol Integer -> Integer)
-> (HenselConsts
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
forall (m :: * -> *) a.
MonadRandom m =>
(Prime Integer -> Unipol Integer -> a)
-> (a -> Unipol Integer -> Maybe (Integer, Integer))
-> (a -> Prime Integer -> Unipol Integer -> Integer)
-> (a
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
factorPrimitiveIntegerPolWith
    Prime Integer -> Unipol Integer -> HenselConsts
calcHenselConsts
    ((Unipol Integer -> Maybe (Integer, Integer))
-> HenselConsts -> Unipol Integer -> Maybe (Integer, Integer)
forall a b. a -> b -> a
const ((Unipol Integer -> Maybe (Integer, Integer))
 -> HenselConsts -> Unipol Integer -> Maybe (Integer, Integer))
-> (Unipol Integer -> Maybe (Integer, Integer))
-> HenselConsts
-> Unipol Integer
-> Maybe (Integer, Integer)
forall a b. (a -> b) -> a -> b
$ Maybe (Integer, Integer)
-> Unipol Integer -> Maybe (Integer, Integer)
forall a b. a -> b -> a
const Maybe (Integer, Integer)
forall a. Maybe a
Nothing)
    (\HenselConsts{Natural
BasicConsts
primePower :: Natural
basicConsts :: BasicConsts
primePower :: HenselConsts -> Natural
basicConsts :: HenselConsts -> BasicConsts
..} (Prime Integer -> Integer
forall a. Prime a -> a
unPrime -> Integer
p) Unipol Integer
_ ->
        Integer
p Integer -> Natural -> Integer
forall r. Unital r => r -> Natural -> r
^ Natural
primePower
    )
    ((HenselConsts
  -> Prime Integer
  -> Unipol Integer
  -> [Unipol Integer]
  -> [Unipol Integer])
 -> Unipol Integer -> m [Unipol Integer])
-> (HenselConsts
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ \HenselConsts{basicConsts :: HenselConsts -> BasicConsts
basicConsts = BasicConsts{Double
Integer
lc :: Integer
kB :: Double
kA :: Integer
lc :: BasicConsts -> Integer
kB :: BasicConsts -> Double
kA :: BasicConsts -> Integer
..}, Natural
primePower :: Natural
primePower :: HenselConsts -> Natural
..} (Prime Integer -> Integer
forall a. Prime a -> a
unPrime -> Integer
p) Unipol Integer
f [Unipol Integer]
gs ->
        Natural
-> Int -> Unipol Integer -> [Unipol Integer] -> [Unipol Integer]
multiHensel (Integer -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
p) (Natural -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Natural
primePower) Unipol Integer
f [Unipol Integer]
gs

{-# INLINE calcHenselConsts #-}
calcHenselConsts :: Prime Integer -> Unipol Integer -> HenselConsts
calcHenselConsts :: Prime Integer -> Unipol Integer -> HenselConsts
calcHenselConsts Prime Integer
p Unipol Integer
f =
  let primePower :: Natural
primePower = Double -> Natural
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Natural) -> Double -> Natural
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. Floating a => a -> a -> a
logBase (Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Double) -> Integer -> Double
forall a b. (a -> b) -> a -> b
$ Prime Integer -> Integer
forall a. Prime a -> a
unPrime Prime Integer
p) (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
2 Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
* Double
kB Double -> Double -> Double
forall r. Additive r => r -> r -> r
+ Double
1
      basicConsts :: BasicConsts
basicConsts@BasicConsts{Double
Integer
lc :: Integer
kA :: Integer
kB :: Double
lc :: BasicConsts -> Integer
kB :: BasicConsts -> Double
kA :: BasicConsts -> Integer
..} = Unipol Integer -> BasicConsts
calcBasicConsts Unipol Integer
f
  in HenselConsts :: BasicConsts -> Natural -> HenselConsts
HenselConsts{Natural
BasicConsts
basicConsts :: BasicConsts
primePower :: Natural
primePower :: Natural
basicConsts :: BasicConsts
..}

calcBasicConsts :: Unipol Integer -> BasicConsts
calcBasicConsts :: Unipol Integer -> BasicConsts
calcBasicConsts Unipol Integer
f =
  let lc :: Coefficient (Unipol Integer)
lc = Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f
      kA :: Norm (Coefficient (Unipol Integer))
kA = Unipol Integer -> Norm (Coefficient (Unipol Integer))
forall poly.
(IsPolynomial poly, Normed (Coefficient poly)) =>
poly -> Norm (Coefficient poly)
maxNorm Unipol Integer
f
      n :: Int
n  = Unipol Integer -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Integer
f
      kB :: Double
kB = Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall r. Additive r => r -> r -> r
+ Double
1) Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
* Double
2 Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
P.^ Int
n Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
* Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer
Norm (Coefficient (Unipol Integer))
kA Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Integer
Coefficient (Unipol Integer)
lc)
  in BasicConsts :: Integer -> Double -> Integer -> BasicConsts
BasicConsts{Double
Integer
Coefficient (Unipol Integer)
Norm (Coefficient (Unipol Integer))
kB :: Double
kA :: Norm (Coefficient (Unipol Integer))
lc :: Coefficient (Unipol Integer)
lc :: Integer
kB :: Double
kA :: Integer
..}

factorPrimitiveIntegerPolWith
  :: MonadRandom m
  => (Prime Integer -> Unipol Integer -> a)
      -- ^ Preprocess
  -> (a -> Unipol Integer -> Maybe (Integer, Integer))
      -- ^ Possible searching range for primes
  -> (a -> Prime Integer -> Unipol Integer -> Integer)
      -- ^ Modulus for final coefficients
  -> (a -> Prime Integer -> Unipol Integer
      -> [Unipol Integer] -> [Unipol Integer])
      -- ^ Post processor for factor; identity in big-prime, multiHensel for hensel.
  -> Unipol Integer -> m [Unipol Integer]
factorPrimitiveIntegerPolWith :: (Prime Integer -> Unipol Integer -> a)
-> (a -> Unipol Integer -> Maybe (Integer, Integer))
-> (a -> Prime Integer -> Unipol Integer -> Integer)
-> (a
    -> Prime Integer
    -> Unipol Integer
    -> [Unipol Integer]
    -> [Unipol Integer])
-> Unipol Integer
-> m [Unipol Integer]
factorPrimitiveIntegerPolWith Prime Integer -> Unipol Integer -> a
pre a -> Unipol Integer -> Maybe (Integer, Integer)
mRange a -> Prime Integer -> Unipol Integer -> Integer
modulus a
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
postProcess Unipol Integer
f =
  let lc :: Coefficient (Unipol Integer)
lc = Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f
      info :: a
info = Prime Integer -> Unipol Integer -> a
pre Prime Integer
p Unipol Integer
f
      range :: [Prime Integer]
range = [Prime Integer]
-> ((Integer, Integer) -> [Prime Integer])
-> Maybe (Integer, Integer)
-> [Prime Integer]
forall b a. b -> (a -> b) -> Maybe a -> b
maybe [Prime Integer]
forall a. Integral a => [Prime a]
PRIMES.primes (\(Integer
l,Integer
r) -> [Integer -> Prime Integer
forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
precPrime Integer
l .. Integer -> Prime Integer
forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
precPrime Integer
r])
          (Maybe (Integer, Integer) -> [Prime Integer])
-> Maybe (Integer, Integer) -> [Prime Integer]
forall a b. (a -> b) -> a -> b
$ a -> Unipol Integer -> Maybe (Integer, Integer)
mRange a
info Unipol Integer
f
      Just Prime Integer
p = (Prime Integer -> Bool) -> [Prime Integer] -> Maybe (Prime Integer)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find Prime Integer -> Bool
isGoodPrime [Prime Integer]
range
  in Integer
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> m [Unipol Integer])
-> m [Unipol Integer]
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField (Prime Integer -> Integer
forall a. Prime a -> a
unPrime Prime Integer
p) ((forall (p :: Nat).
  KnownNat p =>
  Proxy (F p) -> m [Unipol Integer])
 -> m [Unipol Integer])
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> m [Unipol Integer])
-> m [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
fp -> do
    let lc' :: F p
lc' = Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
fp Integer
Coefficient (Unipol Integer)
lc
        f0 :: Unipol (F p)
f0 = (Integer -> F p) -> Unipol Integer -> Unipol (F p)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol ((F p -> F p -> F p
forall r. Division r => r -> r -> r
/F p
lc') (F p -> F p) -> (Integer -> F p) -> Integer -> F p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
fp) Unipol Integer
f
    [Unipol Integer]
fps <- ((Unipol (F p), Natural) -> Unipol Integer)
-> [(Unipol (F p), Natural)] -> [Unipol Integer]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Integer -> Unipol Integer -> Unipol Integer
normalizeMod (Prime Integer -> Integer
forall a. Prime a -> a
unPrime Prime Integer
p)
              (Unipol Integer -> Unipol Integer)
-> ((Unipol (F p), Natural) -> Unipol Integer)
-> (Unipol (F p), Natural)
-> Unipol Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (F p -> Integer) -> Unipol (F p) -> Unipol Integer
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol F p -> Integer
forall k (p :: k). F p -> Integer
naturalRepr (Unipol (F p) -> Unipol Integer)
-> ((Unipol (F p), Natural) -> Unipol (F p))
-> (Unipol (F p), Natural)
-> Unipol Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Unipol (F p), Natural) -> Unipol (F p)
forall a b. (a, b) -> a
fst)
        ([(Unipol (F p), Natural)] -> [Unipol Integer])
-> m [(Unipol (F p), Natural)] -> m [Unipol Integer]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Unipol (F p) -> m [(Unipol (F p), Natural)]
forall (m :: * -> *) k.
(MonadRandom m, CoeffRing k, FiniteField k, Random k) =>
Unipol k -> m [(Unipol k, Natural)]
factorise Unipol (F p)
f0
    let gs :: Vector (Unipol Integer)
gs = [Unipol Integer] -> Vector (Unipol Integer)
forall a. [a] -> Vector a
V.fromList ([Unipol Integer] -> Vector (Unipol Integer))
-> [Unipol Integer] -> Vector (Unipol Integer)
forall a b. (a -> b) -> a -> b
$ a
-> Prime Integer
-> Unipol Integer
-> [Unipol Integer]
-> [Unipol Integer]
postProcess a
info Prime Integer
p Unipol Integer
f [Unipol Integer]
fps
        count :: Int
count = Vector (Unipol Integer) -> Int
forall a. Vector a -> Int
V.length Vector (Unipol Integer)
gs
        alives :: [Int]
alives = [Int
0 .. Int
count Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1]
    [Unipol Integer] -> m [Unipol Integer]
forall (m :: * -> *) a. Monad m => a -> m a
return ([Unipol Integer] -> m [Unipol Integer])
-> [Unipol Integer] -> m [Unipol Integer]
forall a b. (a -> b) -> a -> b
$ Integer
-> [Int]
-> Int
-> Int
-> Unipol Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
loop (a -> Prime Integer -> Unipol Integer -> Integer
modulus a
info Prime Integer
p Unipol Integer
f) [Int]
alives Int
count Int
1 Unipol Integer
f Vector (Unipol Integer)
gs []
  where
    lc :: Coefficient (Unipol Integer)
lc = Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f
    isGoodPrime :: Prime Integer -> Bool
isGoodPrime Prime Integer
p = Integer
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField (Prime Integer -> Integer
forall a. Prime a -> a
unPrime Prime Integer
p) ((forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool)
-> (forall (p :: Nat). KnownNat p => Proxy (F p) -> Bool) -> Bool
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
fp ->
      Integer
Coefficient (Unipol Integer)
lc Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Prime Integer -> Integer
forall a. Prime a -> a
unPrime Prime Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 Bool -> Bool -> Bool
&& Unipol (F p) -> Bool
forall poly.
(IsOrderedPolynomial poly, GCDDomain poly) =>
poly -> Bool
isSquareFree ((Integer -> F p) -> Unipol Integer -> Unipol (F p)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
fp) Unipol Integer
f)
    kA :: Norm (Coefficient (Unipol Integer))
kA = Unipol Integer -> Norm (Coefficient (Unipol Integer))
forall poly.
(IsPolynomial poly, Normed (Coefficient poly)) =>
poly -> Norm (Coefficient poly)
maxNorm Unipol Integer
f
    n :: Int
n = Unipol Integer -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol Integer
f
    kB :: Double
    kB :: Double
kB = Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall r. Additive r => r -> r -> r
+ Double
1) Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
* Double
2 Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
P.^ Int
n Double -> Double -> Double
forall r. Multiplicative r => r -> r -> r
* Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer
Norm (Coefficient (Unipol Integer))
kA Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Integer
Coefficient (Unipol Integer)
lc)
    loop :: Integer
-> [Int]
-> Int
-> Int
-> Unipol Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
loop !Integer
pk [Int]
alives !Int
count !Int
l !Unipol Integer
h Vector (Unipol Integer)
gs [Unipol Integer]
acc
      | Int
2 Int -> Int -> Int
forall r. Multiplicative r => r -> r -> r
* Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
count =
          if Unipol Integer
h Unipol Integer -> Unipol Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Unipol Integer
forall r. Unital r => r
one then [Unipol Integer]
acc else Unipol Integer
h Unipol Integer -> [Unipol Integer] -> [Unipol Integer]
forall a. a -> [a] -> [a]
: [Unipol Integer]
acc
      | Bool
otherwise = Integer
-> [Int]
-> Int
-> Int
-> [([Int], [Int])]
-> Unipol Integer
-> Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
go Integer
pk [Int]
alives Int
count Int
l
          (Int -> [Int] -> [([Int], [Int])]
combinations Int
l [Int]
alives) Unipol Integer
h (Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
h) Vector (Unipol Integer)
gs [Unipol Integer]
acc
    go :: Integer
-> [Int]
-> Int
-> Int
-> [([Int], [Int])]
-> Unipol Integer
-> Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
go Integer
pk [Int]
alives !Int
count Int
l [] Unipol Integer
f' Integer
_ Vector (Unipol Integer)
pols [Unipol Integer]
acc = Integer
-> [Int]
-> Int
-> Int
-> Unipol Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
loop Integer
pk [Int]
alives Int
count (Int
l Int -> Int -> Int
forall r. Additive r => r -> r -> r
+ Int
1) Unipol Integer
f' Vector (Unipol Integer)
pols [Unipol Integer]
acc
    go Integer
pk [Int]
alives !Int
count Int
l (([Int]
gs, [Int]
hs) : [([Int], [Int])]
cands) Unipol Integer
f' Integer
b Vector (Unipol Integer)
pols [Unipol Integer]
acc
      | Unipol Integer -> Norm (Coefficient (Unipol Integer))
forall poly.
(IsPolynomial poly, Normed (Coefficient poly),
 Monoidal (Norm (Coefficient poly))) =>
poly -> Norm (Coefficient poly)
oneNorm Unipol Integer
g Integer -> Integer -> Integer
forall r. Multiplicative r => r -> r -> r
* Unipol Integer -> Norm (Coefficient (Unipol Integer))
forall poly.
(IsPolynomial poly, Normed (Coefficient poly),
 Monoidal (Norm (Coefficient poly))) =>
poly -> Norm (Coefficient poly)
oneNorm Unipol Integer
h Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor Double
kB =
          Integer
-> [Int]
-> Int
-> Int
-> Unipol Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
loop Integer
pk [Int]
hs (Int
count Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
l) Int
l (Unipol Integer -> Unipol Integer
forall poly.
(Euclidean (Coefficient poly), IsPolynomial poly) =>
poly -> poly
pp Unipol Integer
h) Vector (Unipol Integer)
pols (Unipol Integer -> Unipol Integer
forall poly.
(Euclidean (Coefficient poly), IsPolynomial poly) =>
poly -> poly
pp Unipol Integer
g Unipol Integer -> [Unipol Integer] -> [Unipol Integer]
forall a. a -> [a] -> [a]
: [Unipol Integer]
acc)
      | Bool
otherwise = Integer
-> [Int]
-> Int
-> Int
-> [([Int], [Int])]
-> Unipol Integer
-> Integer
-> Vector (Unipol Integer)
-> [Unipol Integer]
-> [Unipol Integer]
go Integer
pk [Int]
alives Int
count Int
l [([Int], [Int])]
cands Unipol Integer
f' Integer
b Vector (Unipol Integer)
pols [Unipol Integer]
acc
      where
        g :: Unipol Integer
g = Integer -> Unipol Integer -> Unipol Integer
normalizeMod Integer
pk (Unipol Integer -> Unipol Integer)
-> Unipol Integer -> Unipol Integer
forall a b. (a -> b) -> a -> b
$ Integer
b Integer -> Unipol Integer -> Unipol Integer
forall r m. LeftModule r m => r -> m -> m
.* Mult (Unipol Integer) -> Unipol Integer
forall a. Mult a -> a
runMult ((Int -> Mult (Unipol Integer)) -> [Int] -> Mult (Unipol Integer)
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap (Unipol Integer -> Mult (Unipol Integer)
forall a. a -> Mult a
Mult (Unipol Integer -> Mult (Unipol Integer))
-> (Int -> Unipol Integer) -> Int -> Mult (Unipol Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Unipol Integer) -> Int -> Unipol Integer
forall a. Vector a -> Int -> a
V.unsafeIndex Vector (Unipol Integer)
pols) [Int]
gs)
        h :: Unipol Integer
h = Integer -> Unipol Integer -> Unipol Integer
normalizeMod Integer
pk (Unipol Integer -> Unipol Integer)
-> Unipol Integer -> Unipol Integer
forall a b. (a -> b) -> a -> b
$ Integer
b Integer -> Unipol Integer -> Unipol Integer
forall r m. LeftModule r m => r -> m -> m
.* Mult (Unipol Integer) -> Unipol Integer
forall a. Mult a -> a
runMult ((Int -> Mult (Unipol Integer)) -> [Int] -> Mult (Unipol Integer)
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap (Unipol Integer -> Mult (Unipol Integer)
forall a. a -> Mult a
Mult (Unipol Integer -> Mult (Unipol Integer))
-> (Int -> Unipol Integer) -> Int -> Mult (Unipol Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Unipol Integer) -> Int -> Unipol Integer
forall a. Vector a -> Int -> a
V.unsafeIndex Vector (Unipol Integer)
pols) [Int]
hs)

-- | Given that @f = gh (mod m)@ with @sg + th = 1 (mod m)@ and @leadingCoeff f@ isn't zero divisor mod m,
--   @henselStep m f g h s t@ calculates the unique (g', h', s', t') s.t.
--   @f = g' h' (mod m^2), g' = g (mod m), h' = h (mod m), s' = s (mod m), t' = t (mod m)@, @h'@ monic.
henselStep
  :: (Eq r, Euclidean r)
  => r        -- ^ modulus
  -> Unipol r
  -> Unipol r
  -> Unipol r
  -> Unipol r
  -> Unipol r
  -> (Unipol r, Unipol r, Unipol r, Unipol r)
henselStep :: r
-> Unipol r
-> Unipol r
-> Unipol r
-> Unipol r
-> Unipol r
-> (Unipol r, Unipol r, Unipol r, Unipol r)
henselStep r
m Unipol r
f Unipol r
g Unipol r
h Unipol r
s Unipol r
t =
  let modCoeff :: Unipol r -> Unipol r
modCoeff = (r -> r) -> Unipol r -> Unipol r
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (r -> r -> r
forall d. Euclidean d => d -> d -> d
`rem` r
mr -> Natural -> r
forall r. Unital r => r -> Natural -> r
^Natural
2)
      divModSq :: Unipol r -> Unipol r -> (Unipol r, Unipol r)
divModSq Unipol r
u Unipol r
v =
        (r -> Fraction r) -> Unipol r -> Unipol (Fraction r)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (r -> r -> Fraction r
forall d. GCDDomain d => d -> d -> Fraction d
F.% r
forall r. Unital r => r
one) Unipol r
u Unipol (Fraction r)
-> Unipol (Fraction r)
-> (Unipol (Fraction r), Unipol (Fraction r))
forall d. Euclidean d => d -> d -> (d, d)
`divide` (r -> Fraction r) -> Unipol r -> Unipol (Fraction r)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (r -> r -> Fraction r
forall d. GCDDomain d => d -> d -> Fraction d
F.% r
forall r. Unital r => r
one) Unipol r
v
        (Unipol (Fraction r), Unipol (Fraction r))
-> ((Unipol (Fraction r), Unipol (Fraction r))
    -> (Unipol r, Unipol r))
-> (Unipol r, Unipol r)
forall a b. a -> (a -> b) -> b
& (Unipol (Fraction r) -> Identity (Unipol r))
-> (Unipol (Fraction r), Unipol (Fraction r))
-> Identity (Unipol r, Unipol r)
forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both ((Unipol (Fraction r) -> Identity (Unipol r))
 -> (Unipol (Fraction r), Unipol (Fraction r))
 -> Identity (Unipol r, Unipol r))
-> (Unipol (Fraction r) -> Unipol r)
-> (Unipol (Fraction r), Unipol (Fraction r))
-> (Unipol r, Unipol r)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Fraction r -> r) -> Unipol (Fraction r) -> Unipol r
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Maybe r -> r
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe r -> r) -> (Fraction r -> Maybe r) -> Fraction r -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. r -> Fraction r -> Maybe r
forall s. (Euclidean s, Eq s) => s -> Fraction s -> Maybe s
modFraction (r
mr -> Natural -> r
forall r. Unital r => r -> Natural -> r
^Natural
2))
      e :: Unipol r
e = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
f Unipol r -> Unipol r -> Unipol r
forall r. Group r => r -> r -> r
- Unipol r
g Unipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
* Unipol r
h
      (Unipol r
q, Unipol r
r) = Unipol r -> Unipol r -> (Unipol r, Unipol r)
divModSq (Unipol r
sUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
e) Unipol r
h
      g' :: Unipol r
g' = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
g Unipol r -> Unipol r -> Unipol r
forall r. Additive r => r -> r -> r
+ Unipol r
t Unipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
* Unipol r
e Unipol r -> Unipol r -> Unipol r
forall r. Additive r => r -> r -> r
+ Unipol r
q Unipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
* Unipol r
g
      h' :: Unipol r
h' = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
h Unipol r -> Unipol r -> Unipol r
forall r. Additive r => r -> r -> r
+ Unipol r
r
      b :: Unipol r
b  = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
sUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
g' Unipol r -> Unipol r -> Unipol r
forall r. Additive r => r -> r -> r
+ Unipol r
tUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
h' Unipol r -> Unipol r -> Unipol r
forall r. Group r => r -> r -> r
- Unipol r
forall r. Unital r => r
one
      (Unipol r
c, Unipol r
d) = Unipol r -> Unipol r -> (Unipol r, Unipol r)
divModSq (Unipol r
sUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
b) Unipol r
h'
      s' :: Unipol r
s' = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
s Unipol r -> Unipol r -> Unipol r
forall r. Group r => r -> r -> r
- Unipol r
d
      t' :: Unipol r
t' = Unipol r -> Unipol r
modCoeff (Unipol r -> Unipol r) -> Unipol r -> Unipol r
forall a b. (a -> b) -> a -> b
$ Unipol r
t Unipol r -> Unipol r -> Unipol r
forall r. Group r => r -> r -> r
- Unipol r
tUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
b Unipol r -> Unipol r -> Unipol r
forall r. Group r => r -> r -> r
- Unipol r
cUnipol r -> Unipol r -> Unipol r
forall r. Multiplicative r => r -> r -> r
*Unipol r
g'
  in (Unipol r
g', Unipol r
h', Unipol r
s', Unipol r
t')

-- | Repeatedly applies hensel lifting for monics.
repeatHensel
  :: Integer -> Int
  -> Unipol Integer
  -> Unipol Integer -> Unipol Integer
  -> Unipol Integer -> Unipol Integer
  -> (Unipol Integer, Unipol Integer, Unipol Integer, Unipol Integer)
repeatHensel :: Integer
-> Int
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> (Unipol Integer, Unipol Integer, Unipol Integer, Unipol Integer)
repeatHensel !Integer
m Int
0 Unipol Integer
_ Unipol Integer
g Unipol Integer
h Unipol Integer
s Unipol Integer
t = (Integer -> Unipol Integer -> Unipol Integer
normalizeMod Integer
m Unipol Integer
g, Integer -> Unipol Integer -> Unipol Integer
normalizeMod Integer
m Unipol Integer
h, Unipol Integer
s, Unipol Integer
t)
repeatHensel !Integer
m Int
n Unipol Integer
f Unipol Integer
g Unipol Integer
h Unipol Integer
s Unipol Integer
t =
  let (Unipol Integer
g', Unipol Integer
h', Unipol Integer
s', Unipol Integer
t') = Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> (Unipol Integer, Unipol Integer, Unipol Integer, Unipol Integer)
forall r.
(Eq r, Euclidean r) =>
r
-> Unipol r
-> Unipol r
-> Unipol r
-> Unipol r
-> Unipol r
-> (Unipol r, Unipol r, Unipol r, Unipol r)
henselStep Integer
m Unipol Integer
f Unipol Integer
g Unipol Integer
h Unipol Integer
s Unipol Integer
t
  in Integer
-> Int
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> (Unipol Integer, Unipol Integer, Unipol Integer, Unipol Integer)
repeatHensel (Integer
m Integer -> Natural -> Integer
forall r. Unital r => r -> Natural -> r
^ Natural
2) (Int
n Int -> Int -> Int
forall r. Group r => r -> r -> r
- Int
1) Unipol Integer
f Unipol Integer
g' Unipol Integer
h' Unipol Integer
s' Unipol Integer
t'

-- | Monic hensel lifting for many factors.
multiHensel :: Natural          -- ^ prime @p@
            -> Int              -- ^ iteration count @k@.
            -> Unipol Integer   -- ^ original polynomial
            -> [Unipol Integer] -- ^ monic coprime factorisation mod @p@
            -> [Unipol Integer] -- ^ coprime factorisation mod @p^(2^k)@.
{-# INLINE multiHensel #-}
multiHensel :: Natural
-> Int -> Unipol Integer -> [Unipol Integer] -> [Unipol Integer]
multiHensel (Natural -> Integer
forall a. Integral a => a -> Integer
toInteger -> Integer
p) Int
n = \Unipol Integer
f -> DList (Unipol Integer) -> [Unipol Integer]
forall a. DList a -> [a]
DL.toList (DList (Unipol Integer) -> [Unipol Integer])
-> ([Unipol Integer] -> DList (Unipol Integer))
-> [Unipol Integer]
-> [Unipol Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Unipol Integer -> Vector (Unipol Integer) -> DList (Unipol Integer)
go Unipol Integer
f (Vector (Unipol Integer) -> DList (Unipol Integer))
-> ([Unipol Integer] -> Vector (Unipol Integer))
-> [Unipol Integer]
-> DList (Unipol Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Unipol Integer] -> Vector (Unipol Integer)
forall a. [a] -> Vector a
V.fromList
  where
    go :: Unipol Integer -> Vector (Unipol Integer) -> DList (Unipol Integer)
go Unipol Integer
f Vector (Unipol Integer)
gs
      | Int
len Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = DList (Unipol Integer)
forall a. DList a
DL.empty
      | Int
len Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = Unipol Integer -> DList (Unipol Integer)
forall a. a -> DList a
DL.singleton (Unipol Integer -> DList (Unipol Integer))
-> Unipol Integer -> DList (Unipol Integer)
forall a b. (a -> b) -> a -> b
$
          let q :: Integer
q = Integer
p Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
P.^ Int
n
              (Integer
u,Integer
_,Integer
bInv) = Integer -> Integer -> (Integer, Integer, Integer)
forall d. PID d => d -> d -> (d, d, d)
egcd Integer
q (Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f)
          in (Integer -> Integer) -> Unipol Integer -> Unipol Integer
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
q)
            (Unipol Integer -> Unipol Integer)
-> Unipol Integer -> Unipol Integer
forall a b. (a -> b) -> a -> b
$ (Integer
bInv Integer -> Integer -> Integer
forall d. Euclidean d => d -> d -> d
`quot` Integer
u) Integer -> Unipol Integer -> Unipol Integer
forall r m. Module (Scalar r) m => r -> m -> m
.*. Unipol Integer
f
      | Bool
otherwise  = Integer
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> DList (Unipol Integer))
-> DList (Unipol Integer)
forall a.
Integer -> (forall (p :: Nat). KnownNat p => Proxy (F p) -> a) -> a
reifyPrimeField Integer
p ((forall (p :: Nat).
  KnownNat p =>
  Proxy (F p) -> DList (Unipol Integer))
 -> DList (Unipol Integer))
-> (forall (p :: Nat).
    KnownNat p =>
    Proxy (F p) -> DList (Unipol Integer))
-> DList (Unipol Integer)
forall a b. (a -> b) -> a -> b
$ \Proxy (F p)
fp ->
        let k :: Int
k = Int
len Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
            d :: Int
d = Int -> Int
ceilingLogBase2 Int
n
            (Vector (Unipol Integer)
ls, Vector (Unipol Integer)
rs) = Int
-> Vector (Unipol Integer)
-> (Vector (Unipol Integer), Vector (Unipol Integer))
forall a. Int -> Vector a -> (Vector a, Vector a)
V.splitAt Int
k Vector (Unipol Integer)
gs
            (Unipol Integer
g0, Unipol Integer
h0) = (Unipol Integer -> Coefficient (Unipol Integer)
forall poly. IsOrderedPolynomial poly => poly -> Coefficient poly
leadingCoeff Unipol Integer
f Integer -> Unipol Integer -> Unipol Integer
forall r m. Module (Scalar r) m => r -> m -> m
.*. Vector (Unipol Integer) -> Unipol Integer
forall (f :: * -> *) r. (Foldable f, Unital r) => f r -> r
product Vector (Unipol Integer)
ls, Vector (Unipol Integer) -> Unipol Integer
forall (f :: * -> *) r. (Foldable f, Unital r) => f r -> r
product Vector (Unipol Integer)
rs)
              (Unipol Integer, Unipol Integer)
-> ((Unipol Integer, Unipol Integer)
    -> (Unipol Integer, Unipol Integer))
-> (Unipol Integer, Unipol Integer)
forall a b. a -> (a -> b) -> b
& (Unipol Integer -> Identity (Unipol Integer))
-> (Unipol Integer, Unipol Integer)
-> Identity (Unipol Integer, Unipol Integer)
forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both ((Unipol Integer -> Identity (Unipol Integer))
 -> (Unipol Integer, Unipol Integer)
 -> Identity (Unipol Integer, Unipol Integer))
-> (Unipol Integer -> Unipol Integer)
-> (Unipol Integer, Unipol Integer)
-> (Unipol Integer, Unipol Integer)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (Integer -> Integer) -> Unipol Integer -> Unipol Integer
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
p)
            (Unipol (F p)
_, Unipol (F p)
s0, Unipol (F p)
t0) = Unipol (F p)
-> Unipol (F p) -> (Unipol (F p), Unipol (F p), Unipol (F p))
forall d. PID d => d -> d -> (d, d, d)
egcd
              ((Integer -> F p) -> Unipol Integer -> Unipol (F p)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
fp) Unipol Integer
g0)
              ((Integer -> F p) -> Unipol Integer -> Unipol (F p)
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Proxy (F p) -> Integer -> F p
forall k (proxy :: * -> *) (p :: k).
Reifies p Integer =>
proxy (F p) -> Integer -> F p
modNat' Proxy (F p)
fp) Unipol Integer
h0)
            (Unipol Integer
s, Unipol Integer
t) = (Unipol (F p)
s0, Unipol (F p)
t0) (Unipol (F p), Unipol (F p))
-> ((Unipol (F p), Unipol (F p))
    -> (Unipol Integer, Unipol Integer))
-> (Unipol Integer, Unipol Integer)
forall a b. a -> (a -> b) -> b
& (Unipol (F p) -> Identity (Unipol Integer))
-> (Unipol (F p), Unipol (F p))
-> Identity (Unipol Integer, Unipol Integer)
forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both ((Unipol (F p) -> Identity (Unipol Integer))
 -> (Unipol (F p), Unipol (F p))
 -> Identity (Unipol Integer, Unipol Integer))
-> (Unipol (F p) -> Unipol Integer)
-> (Unipol (F p), Unipol (F p))
-> (Unipol Integer, Unipol Integer)
forall s t a b. ASetter s t a b -> (a -> b) -> s -> t
%~ (F p -> Integer) -> Unipol (F p) -> Unipol Integer
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol F p -> Integer
forall k (p :: k). F p -> Integer
naturalRepr
            (Unipol Integer
g, Unipol Integer
h, Unipol Integer
_, Unipol Integer
_) = Integer
-> Int
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> Unipol Integer
-> (Unipol Integer, Unipol Integer, Unipol Integer, Unipol Integer)
repeatHensel Integer
p Int
d Unipol Integer
f Unipol Integer
g0 Unipol Integer
h0 Unipol Integer
s Unipol Integer
t
        in Unipol Integer -> Vector (Unipol Integer) -> DList (Unipol Integer)
go Unipol Integer
g Vector (Unipol Integer)
ls DList (Unipol Integer)
-> DList (Unipol Integer) -> DList (Unipol Integer)
forall a. Semigroup a => a -> a -> a
<> Unipol Integer -> Vector (Unipol Integer) -> DList (Unipol Integer)
go Unipol Integer
h Vector (Unipol Integer)
rs
      where len :: Int
len = Vector (Unipol Integer) -> Int
forall a. Vector a -> Int
V.length Vector (Unipol Integer)
gs

recipMod :: (Euclidean a, Eq a) => a -> a -> Maybe a
recipMod :: a -> a -> Maybe a
recipMod a
m a
u =
  let (a
a, a
r, a
_) : [(a, a, a)]
_ = a -> a -> [(a, a, a)]
forall d. Euclidean d => d -> d -> [(d, d, d)]
euclid a
u a
m
  in if a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall r. Unital r => r
one
     then a -> Maybe a
forall a. a -> Maybe a
Just a
r else Maybe a
forall a. Maybe a
Nothing

modFraction :: (Euclidean s, Eq s) => s -> Fraction s -> Maybe s
modFraction :: s -> Fraction s -> Maybe s
modFraction s
m Fraction s
d = ((Fraction s -> s
forall t. Fraction t -> t
numerator Fraction s
d s -> s -> s
forall d. Euclidean d => d -> d -> d
`rem` s
m) s -> s -> s
forall r. Multiplicative r => r -> r -> r
*) (s -> s) -> Maybe s -> Maybe s
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> s -> s -> Maybe s
forall a. (Euclidean a, Eq a) => a -> a -> Maybe a
recipMod s
m (Fraction s -> s
forall t. Fraction t -> t
denominator Fraction s
d)

combinations :: Int -> [Int] -> [([Int], [Int])]
combinations :: Int -> [Int] -> [([Int], [Int])]
combinations Int
n = FMList ([Int], [Int]) -> [([Int], [Int])]
forall (t :: * -> *) a. Foldable t => t a -> [a]
FML.toList (FMList ([Int], [Int]) -> [([Int], [Int])])
-> ([Int] -> FMList ([Int], [Int])) -> [Int] -> [([Int], [Int])]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Int] -> FMList ([Int], [Int])
forall r a. (Eq r, Num r, Group r) => r -> [a] -> FMList ([a], [a])
go Int
n
  where
    go :: r -> [a] -> FMList ([a], [a])
go r
0 [a]
xs  = ([a], [a]) -> FMList ([a], [a])
forall a. a -> FMList a
FML.singleton ([], [a]
xs)
    go r
_ [] = FMList ([a], [a])
forall (f :: * -> *) a. Alternative f => f a
FML.empty
    go r
k (a
x:[a]
xs) =
      let present :: FMList ([a], [a])
present = r -> [a] -> FMList ([a], [a])
go (r
k r -> r -> r
forall r. Group r => r -> r -> r
- r
1) [a]
xs
          absent :: FMList ([a], [a])
absent = r -> [a] -> FMList ([a], [a])
go r
k [a]
xs
      in (([a], [a]) -> ([a], [a]))
-> FMList ([a], [a]) -> FMList ([a], [a])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first (([a] -> [a]) -> ([a], [a]) -> ([a], [a]))
-> ([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall a b. (a -> b) -> a -> b
$ (:) a
x) FMList ([a], [a])
present
        FMList ([a], [a]) -> FMList ([a], [a]) -> FMList ([a], [a])
forall a. Semigroup a => a -> a -> a
<> (([a], [a]) -> ([a], [a]))
-> FMList ([a], [a]) -> FMList ([a], [a])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second (([a] -> [a]) -> ([a], [a]) -> ([a], [a]))
-> ([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall a b. (a -> b) -> a -> b
$ (:) a
x) FMList ([a], [a])
absent

normalizeMod :: Integer -> Unipol Integer -> Unipol Integer
normalizeMod :: Integer -> Unipol Integer -> Unipol Integer
normalizeMod Integer
p = (Integer -> Integer) -> Unipol Integer -> Unipol Integer
forall b a. DecidableZero b => (a -> b) -> Unipol a -> Unipol b
mapCoeffUnipol (Integer -> Integer -> Integer
forall r. Group r => r -> r -> r
subtract Integer
half (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
p) (Integer -> Integer) -> (Integer -> Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall r. Additive r => r -> r -> r
+ Integer
half))
  where half :: Integer
half = Integer
p Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2

isSquareFree :: forall poly. (IsOrderedPolynomial poly, GCDDomain poly) => poly -> Bool
isSquareFree :: poly -> Bool
isSquareFree poly
f = (poly
f poly -> poly -> poly
forall d. GCDDomain d => d -> d -> d
`gcd` Ordinal (Arity poly) -> poly -> poly
forall poly.
IsOrderedPolynomial poly =>
Ordinal (Arity poly) -> poly -> poly
diff Ordinal (Arity poly)
0 poly
f) poly -> poly -> Bool
forall a. Eq a => a -> a -> Bool
== poly
forall r. Unital r => r
one

-- | Reducibility test for univariate polynomials over finite fields.
isReducible
  :: forall k. (FiniteField k, Eq k)
  => Unipol k -> Bool
isReducible :: Unipol k -> Bool
isReducible Unipol k
f
  | Natural
n Natural -> Natural -> Bool
forall a. Eq a => a -> a -> Bool
== Natural
0 = Bool
True
  | Natural
n Natural -> Natural -> Bool
forall a. Eq a => a -> a -> Bool
== Natural
1 = Bool
False
  | Bool
otherwise =
    let divisors :: [Natural]
divisors = ((Prime Natural, Word) -> Natural)
-> [(Prime Natural, Word)] -> [Natural]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
map (Prime Natural -> Natural
forall a. Prime a -> a
PRIMES.unPrime (Prime Natural -> Natural)
-> ((Prime Natural, Word) -> Prime Natural)
-> (Prime Natural, Word)
-> Natural
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Prime Natural, Word) -> Prime Natural
forall a b. (a, b) -> a
fst) ([(Prime Natural, Word)] -> [Natural])
-> [(Prime Natural, Word)] -> [Natural]
forall a b. (a -> b) -> a -> b
$ Natural -> [(Prime Natural, Word)]
forall a. UniqueFactorisation a => a -> [(Prime a, Word)]
PRIMES.factorise Natural
n
        q :: Natural
q = Proxy k -> Natural
forall k (proxy :: * -> *). FiniteField k => proxy k -> Natural
order (Proxy k -> Natural) -> Proxy k -> Natural
forall a b. (a -> b) -> a -> b
$ Proxy k
forall k (t :: k). Proxy t
Proxy @k
        xqn :: Unipol k
xqn = Unipol k -> Natural -> Unipol k -> Unipol k
forall poly. Euclidean poly => poly -> Natural -> poly -> poly
modPow IsLabel "x" (Unipol k)
Unipol k
#x (Natural
qNatural -> Natural -> Natural
forall r. Unital r => r -> Natural -> r
^Natural
n) Unipol k
f
              -- FIXME: Use the fast modular composition algorithm
              -- using Horner's rule
        test :: Natural -> Bool
test Natural
t =
          let b :: Unipol k
b = Unipol k -> Natural -> Unipol k -> Unipol k
forall poly. Euclidean poly => poly -> Natural -> poly -> poly
modPow IsLabel "x" (Unipol k)
Unipol k
#x (Natural
q Natural -> Natural -> Natural
forall r. Unital r => r -> Natural -> r
^ (Natural
n Natural -> Natural -> Natural
forall a. Integral a => a -> a -> a
`div` Natural
t)) Unipol k
f
              -- FIXME: Use the fast modular composition algorithm
              -- using Horner's rule
          in Unipol k -> Unipol k -> Unipol k
forall d. GCDDomain d => d -> d -> d
gcd (Unipol k
b Unipol k -> Unipol k -> Unipol k
forall r. Group r => r -> r -> r
- IsLabel "x" (Unipol k)
Unipol k
#x) Unipol k
f Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= Unipol k
forall r. Unital r => r
one
    in Unipol k
xqn Unipol k -> Unipol k -> Bool
forall a. Eq a => a -> a -> Bool
/= IsLabel "x" (Unipol k)
Unipol k
#x Bool -> Bool -> Bool
|| (Natural -> Bool) -> [Natural] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any Natural -> Bool
test [Natural]
divisors
  where
    n :: Natural
n = (Integral Int, Num Natural) => Int -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral @_ @Natural (Int -> Natural) -> Int -> Natural
forall a b. (a -> b) -> a -> b
$ Unipol k -> Int
forall poly. IsPolynomial poly => poly -> Int
totalDegree' Unipol k
f