{-# 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 :: forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly = 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 :: forall k.
(Fractional k, Ord k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
rootAdd UPolynomial k
p1 UPolynomial k
p2 = 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 :: forall k.
(Fractional k, Ord k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
rootMul UPolynomial k
p1 UPolynomial k
p2 = 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 :: forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootShift k
0 Polynomial k X
p = Polynomial k X
p
rootShift k
r Polynomial k X
p = forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly forall a b. (a -> b) -> a -> b
$ forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst Polynomial k X
p (\X
X -> forall a v. Var a v => v -> a
P.var X
X forall a. Num a => a -> a -> a
- 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 :: forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootScale k
0 UPolynomial k
p = forall a v. Var a v => v -> a
P.var X
X
rootScale k
r UPolynomial k
p = forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly forall a b. (a -> b) -> a -> b
$ 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 -> forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant (forall a. Fractional a => a -> a
recip k
r) forall a. Num a => a -> a -> a
* forall a v. Var a v => v -> a
P.var X
X)

rootRecip :: (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
rootRecip :: forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
rootRecip UPolynomial k
p = forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly forall a b. (a -> b) -> a -> b
$ forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [(k
c, forall a v. Var a v => v -> a
P.var X
X forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (Integer
d forall a. Num a => a -> a -> a
- forall t. Degree t => t -> Integer
P.deg Monomial X
xs)) | (k
c, Monomial X
xs) <- forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial k
p]
  where
    d :: Integer
d = 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 :: forall k l.
(Fractional k, Ord k) =>
(l -> UPolynomial k) -> UPolynomial l -> UPolynomial k
rootSimpPoly l -> UPolynomial k
f UPolynomial l
p = forall k.
(Fractional k, Ord k) =>
Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly (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 = forall a b. [a] -> [b] -> [(a, b)]
zip (forall a. Set a -> [a]
Set.toAscList forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList [l -> UPolynomial k
f l
c | (l
c, 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 = forall k a. [(k, a)] -> Map k a
Map.fromDistinctAscList [(UPolynomial k, Var)]
ys

    p' :: Polynomial k Var
    p' :: Polynomial k Var
p' = forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> forall a v. Var a v => v -> a
P.var Var
0) (forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff (\l
c -> forall a v. Var a v => v -> a
P.var (Map (UPolynomial k) Var
m 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' forall a. a -> [a] -> [a]
: [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 -> 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 :: forall k.
(Fractional k, Ord k) =>
Integer -> UPolynomial k -> UPolynomial k
rootNthRoot Integer
n Polynomial k X
p = forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst Polynomial k X
p (\X
X -> (forall a v. Var a v => v -> a
P.var X
X)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 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
f UPolynomial k
p1 UPolynomial k
p2 = 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 = forall a. Num a => a -> a -> a
f (forall a v. Var a v => v -> a
P.var Var
a) (forall a v. Var a v => v -> a
P.var Var
b)

    gbase :: [Polynomial k Var]
    gbase :: [Polynomial k Var]
gbase = [ 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 -> forall a v. Var a v => v -> a
P.var Var
a), 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 -> 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 :: forall k.
(Fractional k, Ord k) =>
Polynomial k Var -> [Polynomial k Var] -> UPolynomial k
findPoly Polynomial k Var
c [Polynomial k Var]
ps = forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant k
coeff forall a. Num a => a -> a -> a
* (forall a v. Var a v => v -> a
P.var X
X) forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n | (Integer
n,k
coeff) <- forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] [k]
coeffs]
  where
    vn :: Var
    vn :: Var
vn = if forall a. Set a -> Bool
Set.null Set Var
vs then Var
0 else forall a. Set a -> a
Set.findMax Set Var
vs forall a. Num a => a -> a -> a
+ Var
1
      where
        vs :: Set Var
vs = forall a. Ord a => [a] -> Set a
Set.fromList [Var
x | Polynomial k Var
p <- (Polynomial k Var
cforall a. a -> [a] -> [a]
:[Polynomial k Var]
ps), (k
_,Monomial Var
xs) <- forall k v. Polynomial k v -> [Term k v]
P.terms Polynomial k Var
p, (Var
x,Integer
_) <- forall v. Monomial v -> [(v, Integer)]
P.mindices Monomial Var
xs]

    coeffs :: [k]
    coeffs :: [k]
coeffs = forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. [Maybe a] -> [a]
catMaybes forall a b. (a -> b) -> a -> b
$ [[Polynomial k Var] -> Maybe [k]
isLinearlyDependent [Polynomial k Var]
cs2 | [Polynomial k Var]
cs2 <- forall a. [a] -> [[a]]
inits [Polynomial k Var]
cs]
      where
        cmp :: MonomialOrder Var
cmp = forall v. Ord v => MonomialOrder v
P.grlex
        ps' :: [Polynomial k Var]
ps' = 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  = forall a. (a -> a) -> a -> [a]
iterate (\Polynomial k Var
p -> 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 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 forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (k
0forall a. Eq a => a -> a -> Bool
/=) [k]
sol then forall a. a -> Maybe a
Just [k]
sol else forall a. Maybe a
Nothing
      where
        cs2 :: [(Var, Polynomial k Var)]
cs2 = forall a b. [a] -> [b] -> [(a, b)]
zip [Var
vn..] [Polynomial k Var]
cs
        sol :: [k]
sol = forall a b. (a -> b) -> [a] -> [b]
map (\(Var
l,Polynomial k Var
_) -> forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\Var
_ -> k
1) forall a b. (a -> b) -> a -> b
$ 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 (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   = forall v. Ord v => MonomialOrder v
P.grlex
        gbase2 :: [Polynomial k Var]
gbase2 = 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 = forall k a. Map k a -> [a]
Map.elems forall a b. (a -> b) -> a -> b
$ forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith forall a. Num a => a -> a -> a
(+) forall a b. (a -> b) -> a -> b
$ do
          (k
n,Monomial Var
xs) <- forall k v. Polynomial k v -> [Term k v]
P.terms forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall a v. Var a v => v -> a
P.var Var
ln 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' = forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
xs
              ys :: Monomial Var
ys = forall v. Ord v => Map v Integer -> Monomial v
P.mfromIndicesMap forall a b. (a -> b) -> a -> b
$ forall k a. (k -> a -> Bool) -> Map k a -> Map k a
Map.filterWithKey (\Var
x Integer
_ -> Var
x forall a. Ord a => a -> a -> Bool
<  Var
vn) Map Var Integer
xs'
              zs :: Monomial Var
zs = forall v. Ord v => Map v Integer -> Monomial v
P.mfromIndicesMap forall a b. (a -> b) -> a -> b
$ forall k a. (k -> a -> Bool) -> Map k a -> Map k a
Map.filterWithKey (\Var
x Integer
_ -> Var
x forall a. Ord a => a -> a -> Bool
>= Var
vn) Map Var Integer
xs'
          forall (m :: * -> *) a. Monad m => a -> m a
return (Monomial Var
ys, forall k v. (Eq k, Num k) => Term k v -> Polynomial k v
P.fromTerm (k
n,Monomial Var
zs))