{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.AlgebraicNumber.Root
-- Copyright   :  (c) Masahiro Sakai 2012-2016
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Manipulating polynomials for corresponding operations for algebraic numbers.
--
-- Reference:
--
-- * <http://www.dpmms.cam.ac.uk/~wtg10/galois.html>
--
-----------------------------------------------------------------------------
module ToySolver.Data.AlgebraicNumber.Root where

import Data.List
import Data.Maybe
import Data.Map (Map)
import qualified Data.Map as Map
import qualified Data.Set as Set

import ToySolver.Data.Polynomial (Polynomial, UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial as P
import qualified ToySolver.Data.Polynomial.GroebnerBasis as GB

type Var = Int

{--------------------------------------------------------------------
  Manipulation of polynomials
--------------------------------------------------------------------}

normalizePoly :: (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly :: UPolynomial k -> UPolynomial k
normalizePoly = MonomialOrder X -> UPolynomial k -> UPolynomial k
forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
P.toMonic MonomialOrder X
P.nat

rootAdd :: (Fractional k, Ord k) => UPolynomial k -> UPolynomial k -> UPolynomial k
rootAdd :: UPolynomial k -> UPolynomial k -> UPolynomial k
rootAdd UPolynomial k
p1 UPolynomial k
p2 = (forall a. Num a => a -> a -> a)
-> UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Fractional k, Ord k) =>
(forall a. Num a => a -> a -> a)
-> UPolynomial k -> UPolynomial k -> UPolynomial k
lift2 forall a. Num a => a -> a -> a
(+) UPolynomial k
p1 UPolynomial k
p2

rootMul :: (Fractional k, Ord k) => UPolynomial k -> UPolynomial k -> UPolynomial k
rootMul :: UPolynomial k -> UPolynomial k -> UPolynomial k
rootMul UPolynomial k
p1 UPolynomial k
p2 = (forall a. Num a => a -> a -> a)
-> UPolynomial k -> UPolynomial k -> UPolynomial k
forall k.
(Fractional k, Ord k) =>
(forall a. Num a => a -> a -> a)
-> UPolynomial k -> UPolynomial k -> UPolynomial k
lift2 forall a. Num a => a -> a -> a
(*) UPolynomial k
p1 UPolynomial k
p2

rootShift :: (Fractional k, Eq k) => k -> UPolynomial k -> UPolynomial k
rootShift :: k -> UPolynomial k -> UPolynomial k
rootShift k
0 UPolynomial k
p = UPolynomial k
p
rootShift k
r UPolynomial k
p = UPolynomial k -> UPolynomial k
forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly (UPolynomial k -> UPolynomial k) -> UPolynomial k -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ UPolynomial k -> (X -> UPolynomial k) -> UPolynomial k
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
p (\X
X -> X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
- k -> UPolynomial k
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
r)

rootScale :: (Fractional k, Eq k) => k -> UPolynomial k -> UPolynomial k
rootScale :: k -> UPolynomial k -> UPolynomial k
rootScale k
0 UPolynomial k
p = X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X
rootScale k
r UPolynomial k
p = UPolynomial k -> UPolynomial k
forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly (UPolynomial k -> UPolynomial k) -> UPolynomial k -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ UPolynomial k -> (X -> UPolynomial k) -> UPolynomial k
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
p (\X
X -> k -> UPolynomial k
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (k -> k
forall a. Fractional a => a -> a
recip k
r) UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X)

rootRecip :: (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
rootRecip :: UPolynomial k -> UPolynomial k
rootRecip UPolynomial k
p = UPolynomial k -> UPolynomial k
forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly (UPolynomial k -> UPolynomial k) -> UPolynomial k -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ [Term k X] -> UPolynomial k
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [(k
c, X -> Monomial X
forall a v. Var a v => v -> a
P.var X
X Monomial X -> Integer -> Monomial X
forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Monomial X -> Integer
forall t. Degree t => t -> Integer
P.deg Monomial X
xs)) | (k
c, Monomial X
xs) <- UPolynomial k -> [Term k X]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial k
p]
  where
    d :: Integer
d = UPolynomial k -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial k
p

-- | Given a polynomial P = Σ_i a_i x^i over L/K where each a_i is a root of polynomial @f a_i@ over K,
-- @rootSimpPoly f P@ computes a polynomial Q over K such that roots of P is a subset of roots of Q.
rootSimpPoly :: forall k l. (Fractional k, Ord k) => (l -> UPolynomial k) -> UPolynomial l -> UPolynomial k
rootSimpPoly :: (l -> UPolynomial k) -> UPolynomial l -> UPolynomial k
rootSimpPoly l -> UPolynomial k
f UPolynomial l
p = Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
forall k.
(Fractional k, Ord k) =>
Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly (Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
0) [Polynomial k Var]
ps
  where
    ys :: [(UPolynomial k, Var)]
    ys :: [(UPolynomial k, Var)]
ys = [UPolynomial k] -> [Var] -> [(UPolynomial k, Var)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Set (UPolynomial k) -> [UPolynomial k]
forall a. Set a -> [a]
Set.toAscList (Set (UPolynomial k) -> [UPolynomial k])
-> Set (UPolynomial k) -> [UPolynomial k]
forall a b. (a -> b) -> a -> b
$ [UPolynomial k] -> Set (UPolynomial k)
forall a. Ord a => [a] -> Set a
Set.fromList [l -> UPolynomial k
f l
c | (l
c, Monomial X
_) <- UPolynomial l -> [(l, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial l
p]) [Var
1..]

    m :: Map (UPolynomial k) Var
    m :: Map (UPolynomial k) Var
m = [(UPolynomial k, Var)] -> Map (UPolynomial k) Var
forall k a. [(k, a)] -> Map k a
Map.fromDistinctAscList [(UPolynomial k, Var)]
ys

    p' :: Polynomial k Var
    p' :: Polynomial k Var
p' = (X -> Polynomial k Var)
-> Polynomial (Polynomial k Var) X -> Polynomial k Var
forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
0) ((l -> Polynomial k Var)
-> UPolynomial l -> Polynomial (Polynomial k Var) X
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff (\l
c -> Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var (Map (UPolynomial k) Var
m Map (UPolynomial k) Var -> UPolynomial k -> Var
forall k a. Ord k => Map k a -> k -> a
Map.! (l -> UPolynomial k
f l
c))) UPolynomial l
p)

    ps :: [Polynomial k Var]
    ps :: [Polynomial k Var]
ps = Polynomial k Var
p' Polynomial k Var -> [Polynomial k Var] -> [Polynomial k Var]
forall a. a -> [a] -> [a]
: [UPolynomial k -> (X -> Polynomial k Var) -> Polynomial k Var
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
q (\X
X -> Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
x) | (UPolynomial k
q, Var
x) <- [(UPolynomial k, Var)]
ys]

rootNthRoot :: (Fractional k, Ord k) => Integer -> UPolynomial k -> UPolynomial k
rootNthRoot :: Integer -> UPolynomial k -> UPolynomial k
rootNthRoot Integer
n UPolynomial k
p = UPolynomial k -> (X -> UPolynomial k) -> UPolynomial k
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
p (\X
X -> (X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X)UPolynomial k -> Integer -> UPolynomial k
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
n)

lift2 :: forall k. (Fractional k, Ord k) => (forall a. Num a => a -> a -> a)
      -> UPolynomial k -> UPolynomial k -> UPolynomial k
lift2 :: (forall a. Num a => a -> a -> a)
-> UPolynomial k -> UPolynomial k -> UPolynomial k
lift2 forall a. Num a => a -> a -> a
f UPolynomial k
p1 UPolynomial k
p2 = Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
forall k.
(Fractional k, Ord k) =>
Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly Polynomial k Var
f_a_b [Polynomial k Var]
gbase
  where
    a, b :: Var
    a :: Var
a = Var
0
    b :: Var
b = Var
1

    f_a_b :: Polynomial k Var
    f_a_b :: Polynomial k Var
f_a_b = Polynomial k Var -> Polynomial k Var -> Polynomial k Var
forall a. Num a => a -> a -> a
f (Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
a) (Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
b)

    gbase :: [Polynomial k Var]
    gbase :: [Polynomial k Var]
gbase = [ UPolynomial k -> (X -> Polynomial k Var) -> Polynomial k Var
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
p1 (\X
X -> Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
a), UPolynomial k -> (X -> Polynomial k Var) -> Polynomial k Var
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial k
p2 (\X
X -> Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
b) ]

-- | Given a polynomial P and polynomials {P1,…,Pn} over K,
-- findPoly P [P1..Pn] computes a non-zero polynomial Q such that Q[P] = 0 modulo {P1,…,Pn}.
findPoly :: forall k. (Fractional k, Ord k) => Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly :: Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly Polynomial k Var
c [Polynomial k Var]
ps = UPolynomial k -> UPolynomial k
forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly (UPolynomial k -> UPolynomial k) -> UPolynomial k -> UPolynomial k
forall a b. (a -> b) -> a -> b
$ [UPolynomial k] -> UPolynomial k
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [k -> UPolynomial k
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
coeff UPolynomial k -> UPolynomial k -> UPolynomial k
forall a. Num a => a -> a -> a
* (X -> UPolynomial k
forall a v. Var a v => v -> a
P.var X
X) UPolynomial k -> Integer -> UPolynomial k
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n | (Integer
n,k
coeff) <- [Integer] -> [k] -> [(Integer, k)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] [k]
coeffs]
  where
    vn :: Var
    vn :: Var
vn = if Set Var -> Bool
forall a. Set a -> Bool
Set.null Set Var
vs then Var
0 else Set Var -> Var
forall a. Set a -> a
Set.findMax Set Var
vs Var -> Var -> Var
forall a. Num a => a -> a -> a
+ Var
1
      where
        vs :: Set Var
vs = [Var] -> Set Var
forall a. Ord a => [a] -> Set a
Set.fromList [Var
x | Polynomial k Var
p <- (Polynomial k Var
cPolynomial k Var -> [Polynomial k Var] -> [Polynomial k Var]
forall a. a -> [a] -> [a]
:[Polynomial k Var]
ps), (k
_,Monomial Var
xs) <- Polynomial k Var -> [(k, Monomial Var)]
forall k v. Polynomial k v -> [Term k v]
P.terms Polynomial k Var
p, (Var
x,Integer
_) <- Monomial Var -> [(Var, Integer)]
forall v. Monomial v -> [(v, Integer)]
P.mindices Monomial Var
xs]

    coeffs :: [k]
    coeffs :: [k]
coeffs = [[k]] -> [k]
forall a. [a] -> a
head ([[k]] -> [k]) -> [[k]] -> [k]
forall a b. (a -> b) -> a -> b
$ [Maybe [k]] -> [[k]]
forall a. [Maybe a] -> [a]
catMaybes ([Maybe [k]] -> [[k]]) -> [Maybe [k]] -> [[k]]
forall a b. (a -> b) -> a -> b
$ [[Polynomial k Var] -> Maybe [k]
isLinearlyDependent [Polynomial k Var]
cs2 | [Polynomial k Var]
cs2 <- [Polynomial k Var] -> [[Polynomial k Var]]
forall a. [a] -> [[a]]
inits [Polynomial k Var]
cs]
      where
        cmp :: MonomialOrder Var
cmp = MonomialOrder Var
forall v. Ord v => MonomialOrder v
P.grlex
        ps' :: [Polynomial k Var]
ps' = MonomialOrder Var -> [Polynomial k Var] -> [Polynomial k Var]
forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
GB.basis MonomialOrder Var
cmp [Polynomial k Var]
ps
        cs :: [Polynomial k Var]
cs  = (Polynomial k Var -> Polynomial k Var)
-> Polynomial k Var -> [Polynomial k Var]
forall a. (a -> a) -> a -> [a]
iterate (\Polynomial k Var
p -> MonomialOrder Var
-> Polynomial k Var -> [Polynomial k Var] -> Polynomial k Var
forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder Var
cmp (Polynomial k Var
c Polynomial k Var -> Polynomial k Var -> Polynomial k Var
forall a. Num a => a -> a -> a
* Polynomial k Var
p) [Polynomial k Var]
ps') Polynomial k Var
1

    isLinearlyDependent :: [Polynomial k Var] -> Maybe [k]
    isLinearlyDependent :: [Polynomial k Var] -> Maybe [k]
isLinearlyDependent [Polynomial k Var]
cs = if (k -> Bool) -> [k] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (k
0k -> k -> Bool
forall a. Eq a => a -> a -> Bool
/=) [k]
sol then [k] -> Maybe [k]
forall a. a -> Maybe a
Just [k]
sol else Maybe [k]
forall a. Maybe a
Nothing
      where
        cs2 :: [(Var, Polynomial k Var)]
cs2 = [Var] -> [Polynomial k Var] -> [(Var, Polynomial k Var)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Var
vn..] [Polynomial k Var]
cs
        sol :: [k]
sol = ((Var, Polynomial k Var) -> k) -> [(Var, Polynomial k Var)] -> [k]
forall a b. (a -> b) -> [a] -> [b]
map (\(Var
l,Polynomial k Var
_) -> (Var -> k) -> Polynomial k Var -> k
forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\Var
_ -> k
1) (Polynomial k Var -> k) -> Polynomial k Var -> k
forall a b. (a -> b) -> a -> b
$ MonomialOrder Var
-> Polynomial k Var -> [Polynomial k Var] -> Polynomial k Var
forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder Var
cmp2 (Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
l) [Polynomial k Var]
gbase2) [(Var, Polynomial k Var)]
cs2
        cmp2 :: MonomialOrder Var
cmp2   = MonomialOrder Var
forall v. Ord v => MonomialOrder v
P.grlex
        gbase2 :: [Polynomial k Var]
gbase2 = MonomialOrder Var -> [Polynomial k Var] -> [Polynomial k Var]
forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
GB.basis MonomialOrder Var
cmp2 [Polynomial k Var]
es
        es :: [Polynomial k Var]
es = Map (Monomial Var) (Polynomial k Var) -> [Polynomial k Var]
forall k a. Map k a -> [a]
Map.elems (Map (Monomial Var) (Polynomial k Var) -> [Polynomial k Var])
-> Map (Monomial Var) (Polynomial k Var) -> [Polynomial k Var]
forall a b. (a -> b) -> a -> b
$ (Polynomial k Var -> Polynomial k Var -> Polynomial k Var)
-> [(Monomial Var, Polynomial k Var)]
-> Map (Monomial Var) (Polynomial k Var)
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith Polynomial k Var -> Polynomial k Var -> Polynomial k Var
forall a. Num a => a -> a -> a
(+) ([(Monomial Var, Polynomial k Var)]
 -> Map (Monomial Var) (Polynomial k Var))
-> [(Monomial Var, Polynomial k Var)]
-> Map (Monomial Var) (Polynomial k Var)
forall a b. (a -> b) -> a -> b
$ do
          (k
n,Monomial Var
xs) <- Polynomial k Var -> [(k, Monomial Var)]
forall k v. Polynomial k v -> [Term k v]
P.terms (Polynomial k Var -> [(k, Monomial Var)])
-> Polynomial k Var -> [(k, Monomial Var)]
forall a b. (a -> b) -> a -> b
$ [Polynomial k Var] -> Polynomial k Var
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Var -> Polynomial k Var
forall a v. Var a v => v -> a
P.var Var
ln Polynomial k Var -> Polynomial k Var -> Polynomial k Var
forall a. Num a => a -> a -> a
* Polynomial k Var
cn | (Var
ln,Polynomial k Var
cn) <- [(Var, Polynomial k Var)]
cs2]
          let xs' :: Map Var Integer
xs' = Monomial Var -> Map Var Integer
forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
xs
              ys :: Monomial Var
ys = Map Var Integer -> Monomial Var
forall v. Ord v => Map v Integer -> Monomial v
P.mfromIndicesMap (Map Var Integer -> Monomial Var)
-> Map Var Integer -> Monomial Var
forall a b. (a -> b) -> a -> b
$ (Var -> Integer -> Bool) -> Map Var Integer -> Map Var Integer
forall k a. (k -> a -> Bool) -> Map k a -> Map k a
Map.filterWithKey (\Var
x Integer
_ -> Var
x Var -> Var -> Bool
forall a. Ord a => a -> a -> Bool
<  Var
vn) Map Var Integer
xs'
              zs :: Monomial Var
zs = Map Var Integer -> Monomial Var
forall v. Ord v => Map v Integer -> Monomial v
P.mfromIndicesMap (Map Var Integer -> Monomial Var)
-> Map Var Integer -> Monomial Var
forall a b. (a -> b) -> a -> b
$ (Var -> Integer -> Bool) -> Map Var Integer -> Map Var Integer
forall k a. (k -> a -> Bool) -> Map k a -> Map k a
Map.filterWithKey (\Var
x Integer
_ -> Var
x Var -> Var -> Bool
forall a. Ord a => a -> a -> Bool
>= Var
vn) Map Var Integer
xs'
          (Monomial Var, Polynomial k Var)
-> [(Monomial Var, Polynomial k Var)]
forall (m :: * -> *) a. Monad m => a -> m a
return (Monomial Var
ys, (k, Monomial Var) -> Polynomial k Var
forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
P.fromTerm (k
n,Monomial Var
zs))