-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.AlgebraicNumber.Graeffe
-- Copyright   :  (c) Masahiro Sakai 2012
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- Graeffe's Method
--
-- Reference:
--
-- * <http://mathworld.wolfram.com/GraeffesMethod.html>
--
-- * <http://en.wikipedia.org/wiki/Graeffe's_method>
--
-----------------------------------------------------------------------------

module ToySolver.Data.AlgebraicNumber.Graeffe
  ( NthRoot (..)
  , graeffesMethod
  ) where

import Control.Exception
import qualified Data.IntMap as IM
import ToySolver.Data.Polynomial (UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial as P

data NthRoot = NthRoot !Integer !Rational
  deriving (Int -> NthRoot -> ShowS
[NthRoot] -> ShowS
NthRoot -> String
(Int -> NthRoot -> ShowS)
-> (NthRoot -> String) -> ([NthRoot] -> ShowS) -> Show NthRoot
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [NthRoot] -> ShowS
$cshowList :: [NthRoot] -> ShowS
show :: NthRoot -> String
$cshow :: NthRoot -> String
showsPrec :: Int -> NthRoot -> ShowS
$cshowsPrec :: Int -> NthRoot -> ShowS
Show)

graeffesMethod :: UPolynomial Rational -> Int -> [NthRoot]
graeffesMethod :: UPolynomial Rational -> Int -> [NthRoot]
graeffesMethod UPolynomial Rational
p Int
v = [[NthRoot]]
xs [[NthRoot]] -> Int -> [NthRoot]
forall a. [a] -> Int -> a
!! (Int
v Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
  where
    xs :: [[NthRoot]]
xs = ((Int, UPolynomial Rational) -> [NthRoot])
-> [(Int, UPolynomial Rational)] -> [[NthRoot]]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> UPolynomial Rational -> [NthRoot])
-> (Int, UPolynomial Rational) -> [NthRoot]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> UPolynomial Rational -> [NthRoot]
g) ([(Int, UPolynomial Rational)] -> [[NthRoot]])
-> [(Int, UPolynomial Rational)] -> [[NthRoot]]
forall a b. (a -> b) -> a -> b
$ [Int] -> [UPolynomial Rational] -> [(Int, UPolynomial Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..] ([UPolynomial Rational] -> [UPolynomial Rational]
forall a. [a] -> [a]
tail ([UPolynomial Rational] -> [UPolynomial Rational])
-> [UPolynomial Rational] -> [UPolynomial Rational]
forall a b. (a -> b) -> a -> b
$ (UPolynomial Rational -> UPolynomial Rational)
-> UPolynomial Rational -> [UPolynomial Rational]
forall a. (a -> a) -> a -> [a]
iterate UPolynomial Rational -> UPolynomial Rational
f (UPolynomial Rational -> [UPolynomial Rational])
-> UPolynomial Rational -> [UPolynomial Rational]
forall a b. (a -> b) -> a -> b
$ MonomialOrder X -> UPolynomial Rational -> UPolynomial Rational
forall r v.
(Eq r, Fractional r) =>
MonomialOrder v -> Polynomial r v -> Polynomial r v
P.toMonic MonomialOrder X
P.nat UPolynomial Rational
p)

    n :: Integer
n = UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p

    g :: Int -> UPolynomial Rational -> [NthRoot]
    g :: Int -> UPolynomial Rational -> [NthRoot]
g Int
v UPolynomial Rational
p = do
      Int
i <- [Int
1::Int .. Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
n]
      let yi :: Rational
yi = if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then - (Int -> Rational
b Int
i) else - (Int -> Rational
b Int
i Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Int -> Rational
b (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
      NthRoot -> [NthRoot]
forall (m :: * -> *) a. Monad m => a -> m a
return (NthRoot -> [NthRoot]) -> NthRoot -> [NthRoot]
forall a b. (a -> b) -> a -> b
$ Integer -> Rational -> NthRoot
NthRoot (Integer
2 Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
v) Rational
yi
      where
        bs :: IntMap Rational
bs = [(Int, Rational)] -> IntMap Rational
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
i, Rational
b) | (Rational
b,Monomial X
ys) <- UPolynomial Rational -> [(Rational, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial Rational
p, let i :: Integer
i = Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Monomial X -> Integer
forall t. Degree t => t -> Integer
P.deg Monomial X
ys, Integer
i Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0]
        b :: Int -> Rational
b Int
i = Rational -> Int -> IntMap Rational -> Rational
forall a. a -> Int -> IntMap a -> a
IM.findWithDefault Rational
0 Int
i IntMap Rational
bs

f :: UPolynomial Rational -> UPolynomial Rational
f :: UPolynomial Rational -> UPolynomial Rational
f UPolynomial Rational
p = (-UPolynomial Rational
1) UPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^ (UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p) UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
*
      [(Rational, Monomial X)] -> UPolynomial Rational
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [ (Rational
c, Bool -> Monomial X -> Monomial X
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Monomial X -> Integer
forall t. Degree t => t -> Integer
P.deg Monomial X
xs Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
2 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0) (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` (Monomial X -> Integer
forall t. Degree t => t -> Integer
P.deg Monomial X
xs Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2)))
                  | (Rational
c, Monomial X
xs) <- UPolynomial Rational -> [(Rational, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms (UPolynomial Rational
p UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
* UPolynomial Rational
-> (X -> UPolynomial Rational) -> UPolynomial Rational
forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial Rational
p (\X
X -> - X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X)) ]

f' :: UPolynomial Rational -> UPolynomial Rational
f' :: UPolynomial Rational -> UPolynomial Rational
f' UPolynomial Rational
p = [(Rational, Monomial X)] -> UPolynomial Rational
forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [(Integer -> Rational
b Integer
k, 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
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
k)) | Integer
k <- [Integer
0..Integer
n]]
  where
    n :: Integer
n = UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p

    a :: Integer -> Rational
    a :: Integer -> Rational
a Integer
k
      | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
k    = Monomial X -> UPolynomial Rational -> Rational
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (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
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
k)) UPolynomial Rational
p
      | Bool
otherwise = Rational
0

    b :: Integer -> Rational
    b :: Integer -> Rational
b Integer
k = (-Rational
1)Rational -> Integer -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
k Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Integer -> Rational
a Integer
k)Rational -> Integer -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* [Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(-Rational
1)Rational -> Integer -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
j Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Integer -> Rational
a Integer
j) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Integer -> Rational
a (Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
j)) | Integer
j <- [Integer
0..Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1]]

test :: Int -> [NthRoot]
test Int
v = UPolynomial Rational -> Int -> [NthRoot]
graeffesMethod UPolynomial Rational
p Int
v
  where
    x :: UPolynomial Rational
x = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
- UPolynomial Rational
2

test2 :: Int -> [NthRoot]
test2 Int
v = UPolynomial Rational -> Int -> [NthRoot]
graeffesMethod UPolynomial Rational
p Int
v
 where
    x :: UPolynomial Rational
x = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
5 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
- UPolynomial Rational
3UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
*UPolynomial Rational
x UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
- UPolynomial Rational
1