{- |
module: Factor.Zx
description: The polynomial ring Z[x]
license: MIT

maintainer: Joe Leslie-Hurd <joe@gilith.com>
stability: provisional
portability: portable
-}
module Factor.Zx
where

import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap

import qualified Factor.Prime as Prime
import Factor.Term (Term(..),Var)
import qualified Factor.Term as Term
import Factor.Util

-------------------------------------------------------------------------------
-- Monomials in Z[x]
-------------------------------------------------------------------------------

data Monomial =
    Monomial
      {Monomial -> Int
degreeMonomial :: Int,
       Monomial -> Integer
coeffMonomial :: Integer}
  deriving (Monomial -> Monomial -> Bool
(Monomial -> Monomial -> Bool)
-> (Monomial -> Monomial -> Bool) -> Eq Monomial
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Monomial -> Monomial -> Bool
$c/= :: Monomial -> Monomial -> Bool
== :: Monomial -> Monomial -> Bool
$c== :: Monomial -> Monomial -> Bool
Eq,Eq Monomial
Eq Monomial
-> (Monomial -> Monomial -> Ordering)
-> (Monomial -> Monomial -> Bool)
-> (Monomial -> Monomial -> Bool)
-> (Monomial -> Monomial -> Bool)
-> (Monomial -> Monomial -> Bool)
-> (Monomial -> Monomial -> Monomial)
-> (Monomial -> Monomial -> Monomial)
-> Ord Monomial
Monomial -> Monomial -> Bool
Monomial -> Monomial -> Ordering
Monomial -> Monomial -> Monomial
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Monomial -> Monomial -> Monomial
$cmin :: Monomial -> Monomial -> Monomial
max :: Monomial -> Monomial -> Monomial
$cmax :: Monomial -> Monomial -> Monomial
>= :: Monomial -> Monomial -> Bool
$c>= :: Monomial -> Monomial -> Bool
> :: Monomial -> Monomial -> Bool
$c> :: Monomial -> Monomial -> Bool
<= :: Monomial -> Monomial -> Bool
$c<= :: Monomial -> Monomial -> Bool
< :: Monomial -> Monomial -> Bool
$c< :: Monomial -> Monomial -> Bool
compare :: Monomial -> Monomial -> Ordering
$ccompare :: Monomial -> Monomial -> Ordering
$cp1Ord :: Eq Monomial
Ord)

instance Show Monomial where
  show :: Monomial -> String
show Monomial
m = if Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Integer -> String
forall a. Show a => a -> String
show Integer
n else String
showCoeff String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
showPower
    where
      Monomial {degreeMonomial :: Monomial -> Int
degreeMonomial = Int
d, coeffMonomial :: Monomial -> Integer
coeffMonomial = Integer
n} = Monomial
m

      showCoeff :: String
showCoeff = case Integer
n of
                    Integer
1 -> String
""
                    -1 -> String
"-"
                    Integer
_ -> Integer -> String
forall a. Show a => a -> String
show Integer
n

      showPower :: String
showPower = String
"x" String -> ShowS
forall a. [a] -> [a] -> [a]
++ (if Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then String
"" else String
"^" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
d)

isZeroMonomial :: Monomial -> Bool
isZeroMonomial :: Monomial -> Bool
isZeroMonomial Monomial
m = Monomial -> Integer
coeffMonomial Monomial
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0

constantMonomial :: Integer -> Monomial
constantMonomial :: Integer -> Monomial
constantMonomial Integer
n = Monomial :: Int -> Integer -> Monomial
Monomial {degreeMonomial :: Int
degreeMonomial = Int
0, coeffMonomial :: Integer
coeffMonomial = Integer
n}

negateMonomial :: Monomial -> Monomial
negateMonomial :: Monomial -> Monomial
negateMonomial Monomial
m = Monomial
m {coeffMonomial :: Integer
coeffMonomial = Integer -> Integer
forall a. Num a => a -> a
Prelude.negate (Monomial -> Integer
coeffMonomial Monomial
m)}

toTermMonomial :: Var -> Monomial -> Term
toTermMonomial :: String -> Monomial -> Term
toTermMonomial String
v (Monomial {degreeMonomial :: Monomial -> Int
degreeMonomial = Int
d, coeffMonomial :: Monomial -> Integer
coeffMonomial = Integer
n}) =
    Term -> Term -> Term
MultiplyTerm (Integer -> Term
IntegerTerm Integer
n)
      (Term -> Term -> Term
ExpTerm (String -> Term
VarTerm String
v) (Integer -> Term
IntegerTerm (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
d)))

-------------------------------------------------------------------------------
-- The polynomial ring Z[x]
-------------------------------------------------------------------------------

data Zx =
    Zx
      {Zx -> Int
degree :: Int,
       Zx -> IntMap Integer
coeffMap :: IntMap Integer}
  deriving (Zx -> Zx -> Bool
(Zx -> Zx -> Bool) -> (Zx -> Zx -> Bool) -> Eq Zx
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Zx -> Zx -> Bool
$c/= :: Zx -> Zx -> Bool
== :: Zx -> Zx -> Bool
$c== :: Zx -> Zx -> Bool
Eq,Eq Zx
Eq Zx
-> (Zx -> Zx -> Ordering)
-> (Zx -> Zx -> Bool)
-> (Zx -> Zx -> Bool)
-> (Zx -> Zx -> Bool)
-> (Zx -> Zx -> Bool)
-> (Zx -> Zx -> Zx)
-> (Zx -> Zx -> Zx)
-> Ord Zx
Zx -> Zx -> Bool
Zx -> Zx -> Ordering
Zx -> Zx -> Zx
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Zx -> Zx -> Zx
$cmin :: Zx -> Zx -> Zx
max :: Zx -> Zx -> Zx
$cmax :: Zx -> Zx -> Zx
>= :: Zx -> Zx -> Bool
$c>= :: Zx -> Zx -> Bool
> :: Zx -> Zx -> Bool
$c> :: Zx -> Zx -> Bool
<= :: Zx -> Zx -> Bool
$c<= :: Zx -> Zx -> Bool
< :: Zx -> Zx -> Bool
$c< :: Zx -> Zx -> Bool
compare :: Zx -> Zx -> Ordering
$ccompare :: Zx -> Zx -> Ordering
$cp1Ord :: Eq Zx
Ord)

instance Show Zx where
  show :: Zx -> String
show = Term -> String
Term.toString (Term -> String) -> (Zx -> Term) -> Zx -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> Zx -> Term
toTerm String
"x"

valid :: Zx -> Bool
valid :: Zx -> Bool
valid Zx
f =
    (Monomial -> Bool) -> [Monomial] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Bool -> Bool
not (Bool -> Bool) -> (Monomial -> Bool) -> Monomial -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Monomial -> Bool
isZeroMonomial) [Monomial]
ms Bool -> Bool -> Bool
&&
    (if [Monomial] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Monomial]
ms then Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
1 else Int
0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Int] -> Int
forall a. [a] -> a
head [Int]
ds Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Int] -> Int
forall a. [a] -> a
last [Int]
ds)
  where
    d :: Int
d = Zx -> Int
degree Zx
f
    ms :: [Monomial]
ms = Zx -> [Monomial]
toMonomials Zx
f
    ds :: [Int]
ds = (Monomial -> Int) -> [Monomial] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Monomial -> Int
degreeMonomial [Monomial]
ms

fromNormCoeffMap :: IntMap Integer -> Zx
fromNormCoeffMap :: IntMap Integer -> Zx
fromNormCoeffMap IntMap Integer
c | IntMap Integer -> Bool
forall a. IntMap a -> Bool
IntMap.null IntMap Integer
c = Zx
zero
fromNormCoeffMap IntMap Integer
c | Bool
otherwise = Zx :: Int -> IntMap Integer -> Zx
Zx {degree :: Int
degree = Int
d, coeffMap :: IntMap Integer
coeffMap = IntMap Integer
c}
  where (Int
d,Integer
_) = IntMap Integer -> (Int, Integer)
forall a. IntMap a -> (Int, a)
IntMap.findMax IntMap Integer
c

-------------------------------------------------------------------------------
-- Polynomial operations
-------------------------------------------------------------------------------

isZero :: Zx -> Bool
isZero :: Zx -> Bool
isZero Zx
f = Zx -> Int
degree Zx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0

isOne :: Zx -> Bool
isOne :: Zx -> Bool
isOne Zx
f = Zx
f Zx -> Zx -> Bool
forall a. Eq a => a -> a -> Bool
== Zx
one

isConstant :: Zx -> Bool
isConstant :: Zx -> Bool
isConstant Zx
f = Zx -> Int
degree Zx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0

isLinear :: Zx -> Bool
isLinear :: Zx -> Bool
isLinear Zx
f = Zx -> Int
degree Zx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1

isMonic :: Zx -> Bool
isMonic :: Zx -> Bool
isMonic Zx
f = Zx -> Integer
leadingCoeff Zx
f Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1

powerCoeff :: Zx -> Int -> Integer
powerCoeff :: Zx -> Int -> Integer
powerCoeff Zx
f Int
i = Integer -> Int -> IntMap Integer -> Integer
forall a. a -> Int -> IntMap a -> a
IntMap.findWithDefault Integer
0 Int
i (Zx -> IntMap Integer
coeffMap Zx
f)

constantCoeff :: Zx -> Integer
constantCoeff :: Zx -> Integer
constantCoeff Zx
f = Zx -> Int -> Integer
powerCoeff Zx
f Int
0

linearCoeff :: Zx -> Integer
linearCoeff :: Zx -> Integer
linearCoeff Zx
f = Zx -> Int -> Integer
powerCoeff Zx
f Int
1

leadingCoeff :: Zx -> Integer
leadingCoeff :: Zx -> Integer
leadingCoeff Zx
f = Zx -> Int -> Integer
powerCoeff Zx
f (Zx -> Int
degree Zx
f)

trailingCoeff :: Zx -> (Int,Integer)
trailingCoeff :: Zx -> (Int, Integer)
trailingCoeff Zx
f =
    case Zx -> [(Int, Integer)]
monomials Zx
f of
      [] -> String -> (Int, Integer)
forall a. HasCallStack => String -> a
error String
"trailing coefficient is undefined for zero polynomial"
      (Int, Integer)
m : [(Int, Integer)]
_ -> (Int, Integer)
m

monomials :: Zx -> [(Int,Integer)]
monomials :: Zx -> [(Int, Integer)]
monomials = IntMap Integer -> [(Int, Integer)]
forall a. IntMap a -> [(Int, a)]
IntMap.toAscList (IntMap Integer -> [(Int, Integer)])
-> (Zx -> IntMap Integer) -> Zx -> [(Int, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> IntMap Integer
coeffMap

lengthMonomials :: Zx -> Int
lengthMonomials :: Zx -> Int
lengthMonomials = IntMap Integer -> Int
forall a. IntMap a -> Int
IntMap.size (IntMap Integer -> Int) -> (Zx -> IntMap Integer) -> Zx -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> IntMap Integer
coeffMap

filterMonomials :: (Int -> Integer -> Bool) -> Zx -> Zx
filterMonomials :: (Int -> Integer -> Bool) -> Zx -> Zx
filterMonomials Int -> Integer -> Bool
p = IntMap Integer -> Zx
fromNormCoeffMap (IntMap Integer -> Zx) -> (Zx -> IntMap Integer) -> Zx -> Zx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer -> Bool) -> IntMap Integer -> IntMap Integer
forall a. (Int -> a -> Bool) -> IntMap a -> IntMap a
IntMap.filterWithKey Int -> Integer -> Bool
p (IntMap Integer -> IntMap Integer)
-> (Zx -> IntMap Integer) -> Zx -> IntMap Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> IntMap Integer
coeffMap

constant :: Integer -> Zx
constant :: Integer -> Zx
constant = Int -> Integer -> Zx
monomial Int
0

variable :: Zx
variable :: Zx
variable = Int -> Integer -> Zx
monomial Int
1 Integer
1

monomial :: Int -> Integer -> Zx
monomial :: Int -> Integer -> Zx
monomial Int
_ Integer
0 = Zx
zero
monomial Int
d Integer
n = Zx :: Int -> IntMap Integer -> Zx
Zx {degree :: Int
degree = Int
d, coeffMap :: IntMap Integer
coeffMap = Int -> Integer -> IntMap Integer
forall a. Int -> a -> IntMap a
IntMap.singleton Int
d Integer
n}

simpleRoot :: Integer -> Zx
simpleRoot :: Integer -> Zx
simpleRoot Integer
r = Zx -> Zx -> Zx
Factor.Zx.subtract Zx
variable (Integer -> Zx
constant Integer
r)

evaluate :: Zx -> Integer -> Integer
evaluate :: Zx -> Integer -> Integer
evaluate Zx
f Integer
x = Int -> (Int, Integer) -> Integer
forall b. Integral b => b -> (b, Integer) -> Integer
align Int
0 ((Int, Integer) -> Integer) -> (Int, Integer) -> Integer
forall a b. (a -> b) -> a -> b
$ (Int -> Integer -> (Int, Integer) -> (Int, Integer))
-> (Int, Integer) -> IntMap Integer -> (Int, Integer)
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Integer -> (Int, Integer) -> (Int, Integer)
forall b.
Integral b =>
b -> Integer -> (b, Integer) -> (b, Integer)
fma (Int
0,Integer
0) (IntMap Integer -> (Int, Integer))
-> IntMap Integer -> (Int, Integer)
forall a b. (a -> b) -> a -> b
$ Zx -> IntMap Integer
coeffMap Zx
f
  where
    fma :: b -> Integer -> (b, Integer) -> (b, Integer)
fma b
i Integer
n (b, Integer)
z = (b
i, b -> (b, Integer) -> Integer
forall b. Integral b => b -> (b, Integer) -> Integer
align b
i (b, Integer)
z Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n)
    align :: b -> (b, Integer) -> Integer
align b
i (b
d,Integer
n) = let k :: b
k = b
d b -> b -> b
forall a. Num a => a -> a -> a
- b
i in if b
k b -> b -> Bool
forall a. Ord a => a -> a -> Bool
<= b
0 then Integer
n else Integer
xInteger -> b -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^b
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n

derivative :: Zx -> Zx
derivative :: Zx -> Zx
derivative (Zx {degree :: Zx -> Int
degree = Int
d, coeffMap :: Zx -> IntMap Integer
coeffMap = IntMap Integer
c}) =
    if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 then Zx
zero
    else Int -> Zx -> Zx
multiplyPower (-Int
1) (Zx -> Zx) -> Zx -> Zx
forall a b. (a -> b) -> a -> b
$ Zx :: Int -> IntMap Integer -> Zx
Zx {degree :: Int
degree = Int
d, coeffMap :: IntMap Integer
coeffMap = IntMap Integer -> IntMap Integer
deriv IntMap Integer
c}
  where
    deriv :: IntMap Integer -> IntMap Integer
deriv = (Int -> Integer -> Integer) -> IntMap Integer -> IntMap Integer
forall a b. (Int -> a -> b) -> IntMap a -> IntMap b
IntMap.mapWithKey Int -> Integer -> Integer
forall a. Integral a => a -> Integer -> Integer
multDeg (IntMap Integer -> IntMap Integer)
-> (IntMap Integer -> IntMap Integer)
-> IntMap Integer
-> IntMap Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> IntMap Integer -> IntMap Integer
forall a. Int -> IntMap a -> IntMap a
IntMap.delete Int
0
    multDeg :: a -> Integer -> Integer
multDeg a
i Integer
n = a -> Integer
forall a. Integral a => a -> Integer
toInteger a
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n

fromMonomial :: Monomial -> Zx
fromMonomial :: Monomial -> Zx
fromMonomial (Monomial {degreeMonomial :: Monomial -> Int
degreeMonomial = Int
d, coeffMonomial :: Monomial -> Integer
coeffMonomial = Integer
n}) = Int -> Integer -> Zx
monomial Int
d Integer
n

fromMonomials :: [Monomial] -> Zx
fromMonomials :: [Monomial] -> Zx
fromMonomials = [Zx] -> Zx
Factor.Zx.sum ([Zx] -> Zx) -> ([Monomial] -> [Zx]) -> [Monomial] -> Zx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Monomial -> Zx) -> [Monomial] -> [Zx]
forall a b. (a -> b) -> [a] -> [b]
map Monomial -> Zx
fromMonomial

toMonomials :: Zx -> [Monomial]
toMonomials :: Zx -> [Monomial]
toMonomials = ((Int, Integer) -> Monomial) -> [(Int, Integer)] -> [Monomial]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Integer) -> Monomial
mk ([(Int, Integer)] -> [Monomial])
-> (Zx -> [(Int, Integer)]) -> Zx -> [Monomial]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> [(Int, Integer)]
monomials
  where mk :: (Int, Integer) -> Monomial
mk (Int
d,Integer
n) = Monomial :: Int -> Integer -> Monomial
Monomial {degreeMonomial :: Int
degreeMonomial = Int
d, coeffMonomial :: Integer
coeffMonomial = Integer
n}

fromCoeff :: [Integer] -> Zx
fromCoeff :: [Integer] -> Zx
fromCoeff = [Zx] -> Zx
Factor.Zx.sum ([Zx] -> Zx) -> ([Integer] -> [Zx]) -> [Integer] -> Zx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer -> Zx) -> [Int] -> [Integer] -> [Zx]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Integer -> Zx
monomial [Int
0..]

toCoeff :: Zx -> [Integer]
toCoeff :: Zx -> [Integer]
toCoeff Zx
f = (Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Zx -> Int -> Integer
powerCoeff Zx
f) [Int
0 .. Zx -> Int
degree Zx
f]

toTerm :: Var -> Zx -> Term
toTerm :: String -> Zx -> Term
toTerm String
v =
    Term -> Term
Term.simplify (Term -> Term) -> (Zx -> Term) -> Zx -> Term
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
    Term -> Term
Term.nnf (Term -> Term) -> (Zx -> Term) -> Zx -> Term
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
    [Term] -> Term
Term.mkSum ([Term] -> Term) -> (Zx -> [Term]) -> Zx -> Term
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
    (Monomial -> Term) -> [Monomial] -> [Term]
forall a b. (a -> b) -> [a] -> [b]
map (String -> Monomial -> Term
toTermMonomial String
v) ([Monomial] -> [Term]) -> (Zx -> [Monomial]) -> Zx -> [Term]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
    [Monomial] -> [Monomial]
forall a. [a] -> [a]
reverse ([Monomial] -> [Monomial])
-> (Zx -> [Monomial]) -> Zx -> [Monomial]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
    Zx -> [Monomial]
toMonomials

-------------------------------------------------------------------------------
-- Ring operations
-------------------------------------------------------------------------------

zero :: Zx
zero :: Zx
zero = Zx :: Int -> IntMap Integer -> Zx
Zx {degree :: Int
degree = -Int
1, coeffMap :: IntMap Integer
coeffMap = IntMap Integer
forall a. IntMap a
IntMap.empty}

one :: Zx
one :: Zx
one = Integer -> Zx
constant Integer
1

negate :: Zx -> Zx
negate :: Zx -> Zx
negate Zx
f | Zx -> Bool
isZero Zx
f = Zx
zero
negate Zx
f = Zx
f {coeffMap :: IntMap Integer
coeffMap = (Integer -> Integer) -> IntMap Integer -> IntMap Integer
forall a b. (a -> b) -> IntMap a -> IntMap b
IntMap.map Integer -> Integer
forall a. Num a => a -> a
Prelude.negate (Zx -> IntMap Integer
coeffMap Zx
f)}

add :: Zx -> Zx -> Zx
add :: Zx -> Zx -> Zx
add Zx
f Zx
g | Zx -> Bool
isZero Zx
f = Zx
g
add Zx
f Zx
g | Zx -> Bool
isZero Zx
g = Zx
f
add Zx
f Zx
g | Bool
otherwise = IntMap Integer -> Zx
fromNormCoeffMap IntMap Integer
c
  where
    c :: IntMap Integer
c = (Int -> Integer -> Integer -> Maybe Integer)
-> (IntMap Integer -> IntMap Integer)
-> (IntMap Integer -> IntMap Integer)
-> IntMap Integer
-> IntMap Integer
-> IntMap Integer
forall a b c.
(Int -> a -> b -> Maybe c)
-> (IntMap a -> IntMap c)
-> (IntMap b -> IntMap c)
-> IntMap a
-> IntMap b
-> IntMap c
IntMap.mergeWithKey Int -> Integer -> Integer -> Maybe Integer
forall a p. (Eq a, Num a) => p -> a -> a -> Maybe a
addCoeff IntMap Integer -> IntMap Integer
forall a. a -> a
id IntMap Integer -> IntMap Integer
forall a. a -> a
id (Zx -> IntMap Integer
coeffMap Zx
f) (Zx -> IntMap Integer
coeffMap Zx
g)
    addCoeff :: p -> a -> a -> Maybe a
addCoeff p
_ a
fn a
gn = let n :: a
n = a
fn a -> a -> a
forall a. Num a => a -> a -> a
+ a
gn in if a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then Maybe a
forall a. Maybe a
Nothing else a -> Maybe a
forall a. a -> Maybe a
Just a
n

sum :: [Zx] -> Zx
sum :: [Zx] -> Zx
sum = (Zx -> Zx -> Zx) -> Zx -> [Zx] -> Zx
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Zx -> Zx -> Zx
add Zx
zero

subtract :: Zx -> Zx -> Zx
subtract :: Zx -> Zx -> Zx
subtract Zx
f Zx
g = Zx -> Zx -> Zx
add Zx
f (Zx -> Zx
Factor.Zx.negate Zx
g)

multiply :: Zx -> Zx -> Zx
multiply :: Zx -> Zx -> Zx
multiply Zx
f Zx
g | Zx -> Int
lengthMonomials Zx
g Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Zx -> Int
lengthMonomials Zx
f = Zx -> Zx -> Zx
multiply Zx
g Zx
f
multiply Zx
f Zx
_ | Zx -> Bool
isZero Zx
f = Zx
zero
multiply Zx
f Zx
g | Bool
otherwise = (Int -> Integer -> Zx -> Zx) -> Zx -> IntMap Integer -> Zx
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Integer -> Zx -> Zx
fma Zx
zero (Zx -> IntMap Integer
coeffMap Zx
f)
  where fma :: Int -> Integer -> Zx -> Zx
fma Int
i Integer
n Zx
z = Zx -> Zx -> Zx
add (Int -> Zx -> Zx
multiplyPower Int
i (Integer -> Zx -> Zx
multiplyConstant Integer
n Zx
g)) Zx
z

square :: Zx -> Zx
square :: Zx -> Zx
square Zx
f = Zx -> Zx -> Zx
multiply Zx
f Zx
f

product :: [Zx] -> Zx
product :: [Zx] -> Zx
product = (Zx -> Zx -> Zx) -> Zx -> [Zx] -> Zx
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Zx -> Zx -> Zx
multiply Zx
one

multiplyConstant :: Integer -> Zx -> Zx
multiplyConstant :: Integer -> Zx -> Zx
multiplyConstant Integer
0 Zx
_ = Zx
zero
multiplyConstant Integer
1 Zx
f = Zx
f
multiplyConstant Integer
n Zx
f = Zx
f {coeffMap :: IntMap Integer
coeffMap = (Integer -> Integer) -> IntMap Integer -> IntMap Integer
forall a b. (a -> b) -> IntMap a -> IntMap b
IntMap.map (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) Integer
n) (Zx -> IntMap Integer
coeffMap Zx
f)}

multiplyPower :: Int -> Zx -> Zx
multiplyPower :: Int -> Zx -> Zx
multiplyPower Int
0 Zx
f = Zx
f
multiplyPower Int
i (Zx {degree :: Zx -> Int
degree = Int
d, coeffMap :: Zx -> IntMap Integer
coeffMap = IntMap Integer
c}) =
    if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then Zx
zero
    else Zx :: Int -> IntMap Integer -> Zx
Zx {degree :: Int
degree = Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i, coeffMap :: IntMap Integer
coeffMap = (Int -> Int) -> IntMap Integer -> IntMap Integer
forall a. (Int -> Int) -> IntMap a -> IntMap a
IntMap.mapKeysMonotonic (Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
i) IntMap Integer
c}

multiplyExp :: Zx -> Zx -> Integer -> Zx
multiplyExp :: Zx -> Zx -> Integer -> Zx
multiplyExp Zx
z Zx
_ Integer
_ | Zx -> Bool
isZero Zx
z = Zx
zero
multiplyExp Zx
z Zx
_ Integer
0 = Zx
z
multiplyExp Zx
_ Zx
f Integer
_ | Zx -> Bool
isZero Zx
f = Zx
zero
multiplyExp Zx
z Zx
f Integer
_ | Zx -> Bool
isOne Zx
f = Zx
z
multiplyExp Zx
z0 Zx
f0 Integer
k0 = Zx -> Zx -> Integer -> Zx
forall t. Integral t => Zx -> Zx -> t -> Zx
go Zx
z0 Zx
f0 Integer
k0
  where
    go :: Zx -> Zx -> t -> Zx
go Zx
z Zx
_ t
0 = Zx
z
    go Zx
z Zx
f t
k = Zx -> Zx -> t -> Zx
go Zx
z' Zx
f' t
k'
      where
        z' :: Zx
z' = if t -> Bool
forall a. Integral a => a -> Bool
even t
k then Zx
z else Zx -> Zx -> Zx
multiply Zx
z Zx
f
        f' :: Zx
f' = Zx -> Zx
square Zx
f
        k' :: t
k' = t
k t -> t -> t
forall a. Integral a => a -> a -> a
`div` t
2

exp :: Zx -> Integer -> Zx
exp :: Zx -> Integer -> Zx
exp = Zx -> Zx -> Integer -> Zx
multiplyExp Zx
one

-------------------------------------------------------------------------------
-- Polynomial composition
-------------------------------------------------------------------------------

compose :: Zx -> Zx -> Zx
compose :: Zx -> Zx -> Zx
compose Zx
f Zx
g = Int -> (Int, Zx) -> Zx
forall a. Integral a => a -> (a, Zx) -> Zx
align Int
0 ((Int, Zx) -> Zx) -> (Int, Zx) -> Zx
forall a b. (a -> b) -> a -> b
$ (Int -> Integer -> (Int, Zx) -> (Int, Zx))
-> (Int, Zx) -> IntMap Integer -> (Int, Zx)
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
IntMap.foldrWithKey Int -> Integer -> (Int, Zx) -> (Int, Zx)
forall a. Integral a => a -> Integer -> (a, Zx) -> (a, Zx)
fma (Int
0,Zx
zero) (IntMap Integer -> (Int, Zx)) -> IntMap Integer -> (Int, Zx)
forall a b. (a -> b) -> a -> b
$ Zx -> IntMap Integer
coeffMap Zx
f
  where
    fma :: a -> Integer -> (a, Zx) -> (a, Zx)
fma a
i Integer
c (a, Zx)
z = (a
i, Zx -> Zx -> Zx
add (a -> (a, Zx) -> Zx
forall a. Integral a => a -> (a, Zx) -> Zx
align a
i (a, Zx)
z) (Integer -> Zx
constant Integer
c))

    align :: a -> (a, Zx) -> Zx
align a
i (a
d,Zx
z) = if a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then Zx
z else Zx -> Zx -> Integer -> Zx
multiplyExp Zx
z Zx
g (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
k)
      where k :: a
k = a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
i

-------------------------------------------------------------------------------
-- Division
-------------------------------------------------------------------------------

division :: Zx -> Zx -> (Zx,Zx)
division :: Zx -> Zx -> (Zx, Zx)
division Zx
f0 Zx
g = if Int
gd Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then String -> (Zx, Zx)
forall a. HasCallStack => String -> a
error String
"Zx.division by zero" else Zx -> Zx -> (Zx, Zx)
go Zx
zero Zx
f0
  where
    go :: Zx -> Zx -> (Zx, Zx)
go Zx
q Zx
f =
        let fd :: Int
fd = Zx -> Int
degree Zx
f in
        let d :: Int
d = Int
fd Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
gd in
        if Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then (Zx
q,Zx
f)
        else case Integer -> Maybe Integer
gDiv (Zx -> Int -> Integer
powerCoeff Zx
f Int
fd) of
               Maybe Integer
Nothing -> (Zx
q,Zx
f)
               Just Integer
n -> Zx -> Zx -> (Zx, Zx)
go Zx
q' Zx
f'
                 where
                   -- f - n^d*g = q*g + r ==> f = (q + n^d)*g + r
                   q' :: Zx
q' = Zx -> Zx -> Zx
add Zx
q Zx
m
                   f' :: Zx
f' = Zx -> Zx -> Zx
Factor.Zx.subtract Zx
f (Zx -> Zx -> Zx
multiply Zx
m Zx
g)
                   m :: Zx
m = Int -> Integer -> Zx
monomial Int
d Integer
n
    gd :: Int
gd = Zx -> Int
degree Zx
g
    gDiv :: Integer -> Maybe Integer
gDiv = Integer -> Integer -> Maybe Integer
exactQuotient (Zx -> Int -> Integer
powerCoeff Zx
g Int
gd)

quotient :: Zx -> Zx -> Zx
quotient :: Zx -> Zx -> Zx
quotient Zx
f Zx
g = (Zx, Zx) -> Zx
forall a b. (a, b) -> a
fst ((Zx, Zx) -> Zx) -> (Zx, Zx) -> Zx
forall a b. (a -> b) -> a -> b
$ Zx -> Zx -> (Zx, Zx)
Factor.Zx.division Zx
f Zx
g

remainder :: Zx -> Zx -> Zx
remainder :: Zx -> Zx -> Zx
remainder Zx
f Zx
g = (Zx, Zx) -> Zx
forall a b. (a, b) -> b
snd ((Zx, Zx) -> Zx) -> (Zx, Zx) -> Zx
forall a b. (a -> b) -> a -> b
$ Zx -> Zx -> (Zx, Zx)
Factor.Zx.division Zx
f Zx
g

divides :: Zx -> Zx -> Bool
divides :: Zx -> Zx -> Bool
divides Zx
f Zx
g = Zx -> Bool
isZero Zx
g Bool -> Bool -> Bool
|| (Bool -> Bool
not (Zx -> Bool
isZero Zx
f) Bool -> Bool -> Bool
&& Bool
trDiv Bool -> Bool -> Bool
&& Zx -> Bool
isZero (Zx -> Zx -> Zx
remainder Zx
g Zx
f))
  where
    trDiv :: Bool
trDiv = Int
ft Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
gt Bool -> Bool -> Bool
&& Integer -> Integer -> Bool
Factor.Util.divides Integer
fn Integer
gn  -- fast check
    (Int
ft,Integer
fn) = Zx -> (Int, Integer)
trailingCoeff Zx
f
    (Int
gt,Integer
gn) = Zx -> (Int, Integer)
trailingCoeff Zx
g

exactQuotientConstant :: Integer -> Zx -> Zx
exactQuotientConstant :: Integer -> Zx -> Zx
exactQuotientConstant Integer
0 Zx
_ = String -> Zx
forall a. HasCallStack => String -> a
error String
"Z[x] exact quotient by constant 0"
exactQuotientConstant Integer
1 Zx
f = Zx
f
exactQuotientConstant Integer
n Zx
f =
    Zx
f {coeffMap :: IntMap Integer
coeffMap = (Integer -> Integer) -> IntMap Integer -> IntMap Integer
forall a b. (a -> b) -> IntMap a -> IntMap b
IntMap.map (Maybe Integer -> Integer
forall p. Maybe p -> p
nDiv (Maybe Integer -> Integer)
-> (Integer -> Maybe Integer) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer -> Maybe Integer
exactQuotient Integer
n) (Zx -> IntMap Integer
coeffMap Zx
f)}
  where
    nDiv :: Maybe p -> p
nDiv Maybe p
Nothing = String -> p
forall a. HasCallStack => String -> a
error String
"Z[x] exact quotient does not divide all coefficents"
    nDiv (Just p
c) = p
c

-------------------------------------------------------------------------------
-- Primitive part and content
--
-- https://en.wikipedia.org/wiki/Primitive_part_and_content
-------------------------------------------------------------------------------

content :: Zx -> Integer
content :: Zx -> Integer
content Zx
f =
    case IntMap Integer -> Maybe (Integer, IntMap Integer)
forall a. IntMap a -> Maybe (a, IntMap a)
IntMap.minView (Zx -> IntMap Integer
coeffMap Zx
f) of
      Maybe (Integer, IntMap Integer)
Nothing -> Integer
1
      Just (Integer
h,IntMap Integer
t) -> (Integer -> Integer -> Integer)
-> Integer -> IntMap Integer -> Integer
forall a b. (a -> b -> a) -> a -> IntMap b -> a
IntMap.foldl Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
Prelude.gcd (Integer -> Integer
forall a. Num a => a -> a
abs Integer
h) IntMap Integer
t

isPrimitive :: Zx -> Bool
isPrimitive :: Zx -> Bool
isPrimitive = (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
(==) Integer
1) (Integer -> Bool) -> (Zx -> Integer) -> Zx -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> Integer
content

contentPrimitive :: Zx -> (Integer,Zx)
contentPrimitive :: Zx -> (Integer, Zx)
contentPrimitive Zx
f = (Integer
c, Integer -> Zx -> Zx
exactQuotientConstant Integer
c Zx
f)
  where c :: Integer
c = Zx -> Integer
content Zx
f

primitive :: Zx -> Zx
primitive :: Zx -> Zx
primitive = (Integer, Zx) -> Zx
forall a b. (a, b) -> b
snd ((Integer, Zx) -> Zx) -> (Zx -> (Integer, Zx)) -> Zx -> Zx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> (Integer, Zx)
contentPrimitive

-------------------------------------------------------------------------------
-- Greatest common divisor using the subresultant pseudo-remainder sequence
--
-- https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor
-------------------------------------------------------------------------------

gcd :: Zx -> Zx -> Zx
gcd :: Zx -> Zx -> Zx
gcd Zx
f Zx
g = Integer -> Zx -> Zx
multiplyConstant Integer
hk Zx
hp
  where
    (Integer
fc,Zx
fp) = Zx -> (Integer, Zx)
contentPrimitive Zx
f
    (Integer
gc,Zx
gp) = Zx -> (Integer, Zx)
contentPrimitive Zx
g
    hc :: Integer
hc = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
Prelude.gcd Integer
fc Integer
gc
    hp :: Zx
hp = Zx -> Zx -> Zx
gcdPrimitive Zx
fp Zx
gp
    hk :: Integer
hk = if Zx -> Integer
leadingCoeff Zx
hp Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 then Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
hc else Integer
hc

gcdPrimitive :: Zx -> Zx -> Zx
gcdPrimitive :: Zx -> Zx -> Zx
gcdPrimitive Zx
f Zx
g | Zx -> Int
degree Zx
f Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Zx -> Int
degree Zx
g = Zx -> Zx -> Zx
gcdPrimitive Zx
g Zx
f
gcdPrimitive Zx
f Zx
g | Zx -> Bool
isZero Zx
g = Zx
f
gcdPrimitive Zx
f Zx
g = Zx -> Zx
primitive (Zx -> Zx) -> Zx -> Zx
forall a b. (a -> b) -> a -> b
$ Zx -> Zx -> Int -> Integer -> Integer -> Integer -> Zx
go Zx
f Zx
g Int
d1 Integer
g1 Integer
s1 Integer
b1
  where
    go :: Zx -> Zx -> Int -> Integer -> Integer -> Integer -> Zx
go Zx
rip Zx
ri Int
di Integer
gi Integer
si Integer
bi =
        if Zx -> Bool
isZero Zx
r then Zx
ri else Zx -> Zx -> Int -> Integer -> Integer -> Integer -> Zx
go Zx
ri Zx
ri1 Int
di1 Integer
gi1 Integer
si1 Integer
bi1
      where
        r :: Zx
r = Zx -> Zx -> Zx
remainder (Integer -> Zx -> Zx
multiplyConstant (Integer
gi Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
diInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Zx
rip) Zx
ri
        ri1 :: Zx
ri1 = Integer -> Zx -> Zx
exactQuotientConstant Integer
bi Zx
r
        di1 :: Int
di1 = Zx -> Int
degree Zx
ri Int -> Int -> Int
forall a. Num a => a -> a -> a
- Zx -> Int
degree Zx
ri1
        gi1 :: Integer
gi1 = Zx -> Integer
leadingCoeff Zx
ri1
        (Integer
si1,Integer
_) = Integer -> Integer -> (Integer, Integer)
Factor.Util.division ((Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
gi) Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
di) (Integer -> Int -> Integer
forall b p. (Eq p, Num p, Integral b) => p -> b -> p
expm1 Integer
si Int
di)
        bi1 :: Integer
bi1 = Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
gi Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (Integer
si1 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
di1)

    expm1 :: p -> b -> p
expm1 p
si b
0 = if p
si p -> p -> Bool
forall a. Eq a => a -> a -> Bool
== -p
1 then -p
1 else String -> p
forall a. HasCallStack => String -> a
error String
"should have di>0 for all i>1"
    expm1 p
si b
di = p
si p -> b -> p
forall a b. (Num a, Integral b) => a -> b -> a
^ (b
di b -> b -> b
forall a. Num a => a -> a -> a
- b
1)

    d1 :: Int
d1 = Zx -> Int
degree Zx
f Int -> Int -> Int
forall a. Num a => a -> a -> a
- Zx -> Int
degree Zx
g
    g1 :: Integer
g1 = Zx -> Integer
leadingCoeff Zx
g
    s1 :: Integer
s1 = -Integer
1
    b1 :: Integer
b1 = if Int -> Bool
forall a. Integral a => a -> Bool
even Int
d1 then (-Integer
1) else Integer
1

-------------------------------------------------------------------------------
-- Square-free decomposition using Yun's algorithm
--
-- https://en.wikipedia.org/wiki/Square-free_polynomial#Yun's_algorithm
-------------------------------------------------------------------------------

squareFree :: Zx -> Bool
squareFree :: Zx -> Bool
squareFree Zx
f =
    if Zx -> Bool
isZero Zx
f then String -> Bool
forall a. HasCallStack => String -> a
error String
"Z[x] square-free property not defined for zero"
    else Zx -> Bool
isLinear Zx
f Bool -> Bool -> Bool
|| Zx -> Bool
isConstant (Zx -> Zx -> Zx
Factor.Zx.gcd Zx
f (Zx -> Zx
derivative Zx
f))

squareFreeDecomposition :: Zx -> [Zx]
squareFreeDecomposition :: Zx -> [Zx]
squareFreeDecomposition Zx
f =
    if Zx -> Bool
isZero Zx
f then
      String -> [Zx]
forall a. HasCallStack => String -> a
error String
"Z[x] square-free decomposition not defined for zero"
    else if Zx -> Bool
isLinear Zx
f Bool -> Bool -> Bool
|| Zx -> Bool
isConstant Zx
a0 then
      [Zx
f]
    else
      case Zx -> Zx -> Zx -> [Zx]
go Zx
fp Zx
fp' Zx
a0 of
        [] -> String -> [Zx]
forall a. HasCallStack => String -> a
error String
"Z[x] square-free decomposition returned an empty list"
        Zx
a1 : [Zx]
al -> Integer -> Zx -> Zx
multiplyConstant Integer
fk Zx
a1 Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: [Zx]
al
  where
    go :: Zx -> Zx -> Zx -> [Zx]
go Zx
bi Zx
di Zx
ai = if Zx -> Bool
isConstant Zx
bi1 then [] else Zx
ai1 Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: Zx -> Zx -> Zx -> [Zx]
go Zx
bi1 Zx
di1 Zx
ai1
      where
        bi1 :: Zx
bi1 = Zx -> Zx -> Zx
quotient Zx
bi Zx
ai
        di1 :: Zx
di1 = Zx -> Zx -> Zx
Factor.Zx.subtract (Zx -> Zx -> Zx
quotient Zx
di Zx
ai) (Zx -> Zx
derivative Zx
bi1)
        ai1 :: Zx
ai1 = Zx -> Zx -> Zx
Factor.Zx.gcd Zx
bi1 Zx
di1

    (Integer
fc,Zx
fp) = Zx -> (Integer, Zx)
contentPrimitive Zx
f
    fk :: Integer
fk = if Zx -> Integer
leadingCoeff Zx
f Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 then Integer -> Integer
forall a. Num a => a -> a
Prelude.negate Integer
fc else Integer
fc
    fp' :: Zx
fp' = Zx -> Zx
derivative Zx
fp
    a0 :: Zx
a0 = Zx -> Zx -> Zx
Factor.Zx.gcd Zx
fp Zx
fp'

squareFreeRecomposition :: [Zx] -> Zx
squareFreeRecomposition :: [Zx] -> Zx
squareFreeRecomposition [] = String -> Zx
forall a. HasCallStack => String -> a
error String
"Z[x] square-free recomposition of empty list"
squareFreeRecomposition [Zx]
al = (Zx, Zx) -> Zx
forall a b. (a, b) -> a
fst ((Zx, Zx) -> Zx) -> (Zx, Zx) -> Zx
forall a b. (a -> b) -> a -> b
$ (Zx -> (Zx, Zx) -> (Zx, Zx)) -> (Zx, Zx) -> [Zx] -> (Zx, Zx)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Zx -> (Zx, Zx) -> (Zx, Zx)
mult (Zx
one,Zx
one) [Zx]
al
  where mult :: Zx -> (Zx, Zx) -> (Zx, Zx)
mult Zx
a (Zx
f,Zx
g) = (Zx -> Zx -> Zx
multiply Zx
g' Zx
f, Zx
g') where g' :: Zx
g' = Zx -> Zx -> Zx
multiply Zx
a Zx
g

-------------------------------------------------------------------------------
-- Swinnerton-Dyer polynomials
--
-- https://mathworld.wolfram.com/Swinnerton-DyerPolynomial.html
-------------------------------------------------------------------------------

composePlusMinusSqrt :: Zx -> Integer -> Zx
composePlusMinusSqrt :: Zx -> Integer -> Zx
composePlusMinusSqrt = Zx -> Integer -> Zx
go
  where
    go :: Zx -> Integer -> Zx
go Zx
f Integer
0 = Zx -> Zx
square Zx
f
    go Zx
f Integer
a = [Zx] -> Zx
Factor.Zx.sum ([(Int, Integer)] -> [Zx]
mult (Zx -> [(Int, Integer)]
monomials Zx
f))
      where
        mult :: [(Int, Integer)] -> [Zx]
mult [] = []
        mult ((Int
n,Integer
c) : [(Int, Integer)]
ncs) = Zx -> Zx -> Zx
multiply Zx
xan ([Zx] -> Zx
Factor.Zx.sum [Zx]
row) Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: [(Int, Integer)] -> [Zx]
mult [(Int, Integer)]
ncs
          where
            row :: [Zx]
row = Integer -> Zx
constant Integer
c Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: ((Int, Integer) -> Zx) -> [(Int, Integer)] -> [Zx]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Integer) -> Zx
multm [(Int, Integer)]
ncs
            multm :: (Int, Integer) -> Zx
multm (Int
m,Integer
d) = Integer -> Zx -> Zx
multiplyConstant Integer
d ([Zx]
ds [Zx] -> Int -> Zx
forall a. [a] -> Int -> a
!! (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))
            xan :: Zx
xan = Integer -> Zx -> Zx
multiplyConstant Integer
c (Zx -> Integer -> Zx
Factor.Zx.exp Zx
xa (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n))
        xa :: Zx
xa = Zx -> Zx -> Zx
Factor.Zx.subtract (Zx -> Zx
square Zx
variable) (Integer -> Zx
constant Integer
a)
        ds :: [Zx]
ds = ([Zx] -> Zx) -> [[Zx]] -> [Zx]
forall a b. (a -> b) -> [a] -> [b]
map ((Zx -> Zx -> Zx) -> Zx -> [Zx] -> Zx
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Zx
m Zx
z -> Zx -> Zx -> Zx
add Zx
m (Integer -> Zx -> Zx
multiplyConstant Integer
a Zx
z)) Zx
zero) [[Zx]]
bs

    bs :: [[Zx]]
bs = (Zx -> [Zx]) -> [Zx] -> [[Zx]]
forall a b. (a -> b) -> [a] -> [b]
map (((Int, Integer) -> Zx) -> [(Int, Integer)] -> [Zx]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Integer) -> Zx
mon ([(Int, Integer)] -> [Zx])
-> (Zx -> [(Int, Integer)]) -> Zx -> [Zx]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Int, Integer)] -> [(Int, Integer)]
forall a. [a] -> [a]
reverse ([(Int, Integer)] -> [(Int, Integer)])
-> (Zx -> [(Int, Integer)]) -> Zx -> [(Int, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Integer) -> Bool) -> [(Int, Integer)] -> [(Int, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int, Integer) -> Bool
forall a a. (Ord a, Num a) => (a, a) -> Bool
pos ([(Int, Integer)] -> [(Int, Integer)])
-> (Zx -> [(Int, Integer)]) -> Zx -> [(Int, Integer)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Zx -> [(Int, Integer)]
monomials) ([Zx] -> [[Zx]]) -> [Zx] -> [[Zx]]
forall a b. (a -> b) -> a -> b
$ [Zx] -> [Zx]
forall a. [a] -> [a]
tail ([Zx] -> [Zx]) -> [Zx] -> [Zx]
forall a b. (a -> b) -> a -> b
$
         (Zx -> Zx) -> Zx -> [Zx]
forall a. (a -> a) -> a -> [a]
iterate (Zx -> Zx -> Zx
multiply (Zx -> Zx -> Zx
Factor.Zx.subtract Zx
variable Zx
one)) Zx
one
    pos :: (a, a) -> Bool
pos (a
_,a
c) = a
0 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
c
    mon :: (Int, Integer) -> Zx
mon (Int
n,Integer
c) = Int -> Integer -> Zx
monomial Int
n (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
c)

swinnertonDyer :: [Zx]
swinnertonDyer :: [Zx]
swinnertonDyer = Zx -> [Integer] -> [Zx]
go Zx
variable [Integer]
Prime.primes
  where
    go :: Zx -> [Integer] -> [Zx]
go Zx
f [Integer]
ps = Zx
g Zx -> [Zx] -> [Zx]
forall a. a -> [a] -> [a]
: Zx -> [Integer] -> [Zx]
go Zx
g ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
ps)
      where g :: Zx
g = Zx -> Integer -> Zx
composePlusMinusSqrt Zx
f ([Integer] -> Integer
forall a. [a] -> a
head [Integer]
ps)