-----------------------------------------------------------------------------
-- |
-- 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
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 forall a. [a] -> Int -> a
!! (Int
v forall a. Num a => a -> a -> a
- Int
1)
  where
    xs :: [[NthRoot]]
xs = forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> UPolynomial Rational -> [NthRoot]
g) forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..] (forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate UPolynomial Rational -> UPolynomial Rational
f forall a b. (a -> b) -> a -> b
$ 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 = 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 .. forall a. Num a => Integer -> a
fromInteger Integer
n]
      let yi :: Rational
yi = if Int
i forall a. Eq a => a -> a -> Bool
== Int
1 then - (Int -> Rational
b Int
i) else - (Int -> Rational
b Int
i forall a. Fractional a => a -> a -> a
/ Int -> Rational
b (Int
iforall a. Num a => a -> a -> a
-Int
1))
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Integer -> Rational -> NthRoot
NthRoot (Integer
2 forall a b. (Num a, Integral b) => a -> b -> a
^ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
v) Rational
yi
      where
        bs :: IntMap Rational
bs = forall a. [(Int, a)] -> IntMap a
IM.fromList [(forall a. Num a => Integer -> a
fromInteger Integer
i, Rational
b) | (Rational
b,Monomial X
ys) <- forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial Rational
p, let i :: Integer
i = Integer
n forall a. Num a => a -> a -> a
- forall t. Degree t => t -> Integer
P.deg Monomial X
ys, Integer
i forall a. Eq a => a -> a -> Bool
/= Integer
0]
        b :: Int -> Rational
b Int
i = 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) forall a b. (Num a, Integral b) => a -> b -> a
^ (forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p) forall a. Num a => a -> a -> a
*
      forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [ (Rational
c, forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall t. Degree t => t -> Integer
P.deg Monomial X
xs forall a. Integral a => a -> a -> a
`mod` Integer
2 forall a. Eq a => a -> a -> Bool
== Integer
0) (forall a v. Var a v => v -> a
P.var X
X forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (forall t. Degree t => t -> Integer
P.deg Monomial X
xs forall a. Integral a => a -> a -> a
`div` Integer
2)))
                  | (Rational
c, Monomial X
xs) <- forall k v. Polynomial k v -> [Term k v]
P.terms (UPolynomial Rational
p forall a. Num 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 Rational
p (\X
X -> - 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 = forall k v. (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v
P.fromTerms [(Integer -> Rational
b Integer
k, forall a v. Var a v => v -> a
P.var X
X forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (Integer
n forall a. Num a => a -> a -> a
- Integer
k)) | Integer
k <- [Integer
0..Integer
n]]
  where
    n :: Integer
n = forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p

    a :: Integer -> Rational
    a :: Integer -> Rational
a Integer
k
      | Integer
n forall a. Ord a => a -> a -> Bool
>= Integer
k    = forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (forall a v. Var a v => v -> a
P.var X
X forall v. Ord v => Monomial v -> Integer -> Monomial v
`P.mpow` (Integer
n 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)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
k forall a. Num a => a -> a -> a
* (Integer -> Rational
a Integer
k)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
+ Rational
2 forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(-Rational
1)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
j forall a. Num a => a -> a -> a
* (Integer -> Rational
a Integer
j) forall a. Num a => a -> a -> a
* (Integer -> Rational
a (Integer
2forall a. Num a => a -> a -> a
*Integer
kforall a. Num a => a -> a -> a
-Integer
j)) | Integer
j <- [Integer
0..Integer
kforall 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 = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 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 = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
5 forall a. Num a => a -> a -> a
- UPolynomial Rational
3forall a. Num a => a -> a -> a
*UPolynomial Rational
x forall a. Num a => a -> a -> a
- UPolynomial Rational
1