{-# LANGUAGE RebindableSyntax #-}
{- |
Copyright   :  (c) Henning Thielemann 2004-2005

Maintainer  :  numericprelude@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

Computations on the set of roots of a polynomial.
These are represented as the list of their elementar symmetric terms.
The difference between a polynomial and the list of elementar symmetric terms
is the reversed order and the alternated signs.

Cf. /MathObj.PowerSum/ .
-}
module MathObj.RootSet where

import qualified MathObj.Polynomial      as Poly
import qualified MathObj.Polynomial.Core as PolyCore
import qualified MathObj.PowerSum        as PowerSum

import qualified Algebra.Algebraic    as Algebraic
import qualified Algebra.IntegralDomain as Integral
import qualified Algebra.Field        as Field
import qualified Algebra.Ring         as Ring
import qualified Algebra.Additive     as Additive
import qualified Algebra.ZeroTestable as ZeroTestable

import qualified Algebra.RealRing     as RealRing

import qualified Data.List.Match as Match
import qualified Data.List.Key as Key
import Control.Monad (liftM2, replicateM, )

import NumericPrelude.Base as P hiding (const)
import NumericPrelude.Numeric as NP


newtype T a = Cons {T a -> [a]
coeffs :: [a]}


{- * Conversions -}

lift0 :: [a] -> T a
lift0 :: [a] -> T a
lift0 = [a] -> T a
forall a. [a] -> T a
Cons

lift1 :: ([a] -> [a]) -> (T a -> T a)
lift1 :: ([a] -> [a]) -> T a -> T a
lift1 [a] -> [a]
f (Cons [a]
x0) = [a] -> T a
forall a. [a] -> T a
Cons ([a] -> [a]
f [a]
x0)

lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a)
lift2 :: ([a] -> [a] -> [a]) -> T a -> T a -> T a
lift2 [a] -> [a] -> [a]
f (Cons [a]
x0) (Cons [a]
x1) = [a] -> T a
forall a. [a] -> T a
Cons ([a] -> [a] -> [a]
f [a]
x0 [a]
x1)


const :: (Ring.C a) => a -> T a
const :: a -> T a
const a
x = [a] -> T a
forall a. [a] -> T a
Cons [a
1,a
x]


toPolynomial :: T a -> Poly.T a
toPolynomial :: T a -> T a
toPolynomial (Cons [a]
xs) = [a] -> T a
forall a. [a] -> T a
Poly.fromCoeffs ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
xs)

fromPolynomial :: Poly.T a -> T a
fromPolynomial :: T a -> T a
fromPolynomial T a
xs = [a] -> T a
forall a. [a] -> T a
Cons ([a] -> [a]
forall a. [a] -> [a]
reverse (T a -> [a]
forall a. T a -> [a]
Poly.coeffs T a
xs))



toPowerSums :: (Field.C a, ZeroTestable.C a) => [a] -> [a]
toPowerSums :: [a] -> [a]
toPowerSums = [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
PowerSum.fromElemSymDenormalized

fromPowerSums :: (Field.C a, ZeroTestable.C a) => [a] -> [a]
fromPowerSums :: [a] -> [a]
fromPowerSums = [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
PowerSum.toElemSym


{- | cf. 'MathObj.Polynomial.mulLinearFactor' -}
addRoot :: Ring.C a => a -> [a] -> [a]
addRoot :: a -> [a] -> [a]
addRoot a
x yt :: [a]
yt@(a
y:[a]
ys) =
   a
y a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([a]
ys [a] -> [a] -> [a]
forall a. C a => a -> a -> a
+ a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
PolyCore.scale a
x [a]
yt)
addRoot a
_ [] =
   [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"addRoot: list of elementar symmetric terms must consist at least of a 1"

fromRoots :: Ring.C a => [a] -> [a]
fromRoots :: [a] -> [a]
fromRoots = ([a] -> a -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ((a -> [a] -> [a]) -> [a] -> a -> [a]
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
addRoot) [a
1]



liftPowerSum1Gen :: ([a] -> [a]) -> ([a] -> [a]) ->
   ([a] -> [a]) -> ([a] -> [a])
liftPowerSum1Gen :: ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
liftPowerSum1Gen [a] -> [a]
fromPS [a] -> [a]
toPS [a] -> [a]
op [a]
x =
   [a] -> [a] -> [a]
forall b a. [b] -> [a] -> [a]
Match.take [a]
x ([a] -> [a]
fromPS ([a] -> [a]
op ([a] -> [a]
toPS [a]
x)))

liftPowerSum2Gen :: ([a] -> [a]) -> ([a] -> [a]) ->
   ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a])
liftPowerSum2Gen :: ([a] -> [a])
-> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2Gen [a] -> [a]
fromPS [a] -> [a]
toPS [a] -> [a] -> [a]
op [a]
x [a]
y =
   [(a, a)] -> [a] -> [a]
forall b a. [b] -> [a] -> [a]
Match.take ((a, a)
forall a. HasCallStack => a
undefined (a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
: (a -> a -> (a, a)) -> [a] -> [a] -> [(a, a)]
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (,) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
x) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
y))
             ([a] -> [a]
fromPS ([a] -> [a] -> [a]
op ([a] -> [a]
toPS [a]
x) ([a] -> [a]
toPS [a]
y)))


liftPowerSum1 :: (Field.C a, ZeroTestable.C a) =>
   ([a] -> [a]) -> ([a] -> [a])
liftPowerSum1 :: ([a] -> [a]) -> [a] -> [a]
liftPowerSum1 = ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall a.
([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
liftPowerSum1Gen [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
fromPowerSums [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
toPowerSums

liftPowerSum2 :: (Field.C a, ZeroTestable.C a) =>
   ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a])
liftPowerSum2 :: ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2 = ([a] -> [a])
-> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a.
([a] -> [a])
-> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2Gen [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
fromPowerSums [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
toPowerSums

liftPowerSumInt1 :: (Integral.C a, Eq a, ZeroTestable.C a) =>
   ([a] -> [a]) -> ([a] -> [a])
liftPowerSumInt1 :: ([a] -> [a]) -> [a] -> [a]
liftPowerSumInt1 = ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall a.
([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
liftPowerSum1Gen [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
PowerSum.toElemSymInt [a] -> [a]
forall a. (Eq a, C a) => [a] -> [a]
PowerSum.fromElemSym

liftPowerSumInt2 :: (Integral.C a, Eq a, ZeroTestable.C a) =>
   ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a])
liftPowerSumInt2 :: ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSumInt2 = ([a] -> [a])
-> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a.
([a] -> [a])
-> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2Gen [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
PowerSum.toElemSymInt [a] -> [a]
forall a. (Eq a, C a) => [a] -> [a]
PowerSum.fromElemSym




{- * Show -}

appPrec :: Int
appPrec :: Int
appPrec  = Int
10

instance (Show a) => Show (T a) where
  showsPrec :: Int -> T a -> ShowS
showsPrec Int
p (Cons [a]
xs) =
    Bool -> ShowS -> ShowS
showParen (Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
appPrec)
       ([Char] -> ShowS
showString [Char]
"RootSet.Cons " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> ShowS
forall a. Show a => a -> ShowS
shows [a]
xs)


{- * Additive -}

{- Use binomial expansion of (x+y)^n -}
add :: (Field.C a, ZeroTestable.C a) => [a] -> [a] -> [a]
add :: [a] -> [a] -> [a]
add = ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a. (C a, C a) => ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2 [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PowerSum.add

addInt :: (Integral.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> [a]
addInt :: [a] -> [a] -> [a]
addInt = ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a.
(C a, Eq a, C a) =>
([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSumInt2 [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PowerSum.add

instance (Field.C a, ZeroTestable.C a) => Additive.C (T a) where
   zero :: T a
zero   = a -> T a
forall a. C a => a -> T a
const a
forall a. C a => a
zero
   + :: T a -> T a -> T a
(+)    = ([a] -> [a] -> [a]) -> T a -> T a -> T a
forall a. ([a] -> [a] -> [a]) -> T a -> T a -> T a
lift2 [a] -> [a] -> [a]
forall a. (C a, C a) => [a] -> [a] -> [a]
add
   negate :: T a -> T a
negate = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 [a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.alternate


{- * Ring -}

mul :: (Field.C a, ZeroTestable.C a) => [a] -> [a] -> [a]
mul :: [a] -> [a] -> [a]
mul = ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a. (C a, C a) => ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSum2 [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PowerSum.mul

mulInt :: (Integral.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> [a]
mulInt :: [a] -> [a] -> [a]
mulInt = ([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a.
(C a, Eq a, C a) =>
([a] -> [a] -> [a]) -> [a] -> [a] -> [a]
liftPowerSumInt2 [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PowerSum.mul


pow :: (Field.C a, ZeroTestable.C a) => Integer -> [a] -> [a]
pow :: Integer -> [a] -> [a]
pow Integer
n = ([a] -> [a]) -> [a] -> [a]
forall a. (C a, C a) => ([a] -> [a]) -> [a] -> [a]
liftPowerSum1 (Integer -> [a] -> [a]
forall a. Integer -> [a] -> [a]
PowerSum.pow Integer
n)

powInt :: (Integral.C a, Eq a, ZeroTestable.C a) => Integer -> [a] -> [a]
powInt :: Integer -> [a] -> [a]
powInt Integer
n = ([a] -> [a]) -> [a] -> [a]
forall a. (C a, Eq a, C a) => ([a] -> [a]) -> [a] -> [a]
liftPowerSumInt1 (Integer -> [a] -> [a]
forall a. Integer -> [a] -> [a]
PowerSum.pow Integer
n)


instance (Field.C a, ZeroTestable.C a) => Ring.C (T a) where
   one :: T a
one           = a -> T a
forall a. C a => a -> T a
const a
forall a. C a => a
one
   fromInteger :: Integer -> T a
fromInteger Integer
n = a -> T a
forall a. C a => a -> T a
const (Integer -> a
forall a. C a => Integer -> a
fromInteger Integer
n)
   * :: T a -> T a -> T a
(*)           = ([a] -> [a] -> [a]) -> T a -> T a -> T a
forall a. ([a] -> [a] -> [a]) -> T a -> T a -> T a
lift2 [a] -> [a] -> [a]
forall a. (C a, C a) => [a] -> [a] -> [a]
mul
   T a
x^ :: T a -> Integer -> T a
^Integer
n           = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 (Integer -> [a] -> [a]
forall a. (C a, C a) => Integer -> [a] -> [a]
pow Integer
n) T a
x


{- * Field.C -}

instance (Field.C a, ZeroTestable.C a) => Field.C (T a) where
   recip :: T a -> T a
recip = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 [a] -> [a]
forall a. [a] -> [a]
reverse


{- * Algebra -}

instance (Field.C a, ZeroTestable.C a) => Algebraic.C (T a) where
   root :: Integer -> T a -> T a
root Integer
n = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 (Integer -> [a] -> [a]
forall a. C a => Integer -> [a] -> [a]
PowerSum.root Integer
n)



{- |
Given an approximation of a root,
the degree of the polynomial and maximum value of coefficients,
find candidates of polynomials that have approximately this root
and show the actual value of the polynomial at the given root approximation.

This algorithm runs easily into a stack overflow, I do not know why.
We may also employ a more sophisticated integer relation algorithm,
like PSLQ and friends.
-}
{-# SPECIALISE approxPolynomial ::
       Int -> Integer -> Double -> (Double, Poly.T Double) #-}
{-# SPECIALISE approxPolynomial ::
       Int -> Integer -> Float -> (Float, Poly.T Float) #-}
approxPolynomial ::
   (RealRing.C a) =>
   Int -> Integer -> a -> (a, Poly.T a)
approxPolynomial :: Int -> Integer -> a -> (a, T a)
approxPolynomial Int
d Integer
maxCoeff a
x =
   let powers :: [a]
powers = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
dInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
xa -> a -> a
forall a. C a => a -> a -> a
*) a
forall a. C a => a
one
   in  -- List.minimumBy (\a b -> compare (abs (fst a)) (abs (fst b))) $
       ((a, T a) -> a) -> [(a, T a)] -> (a, T a)
forall b a. Ord b => (a -> b) -> [a] -> a
Key.minimum (a -> a
forall a. C a => a -> a
abs (a -> a) -> ((a, T a) -> a) -> (a, T a) -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, T a) -> a
forall a b. (a, b) -> a
fst) ([(a, T a)] -> (a, T a)) -> [(a, T a)] -> (a, T a)
forall a b. (a -> b) -> a -> b
$
       ([a] -> (a, T a)) -> [[a]] -> [(a, T a)]
forall a b. (a -> b) -> [a] -> [b]
map
          ((\[a]
cs -> ([a] -> a
forall a. C a => [a] -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(*) [a]
powers [a]
cs, [a] -> T a
forall a. [a] -> T a
Poly.fromCoeffs [a]
cs)) ([a] -> (a, T a)) -> ([a] -> [a]) -> [a] -> (a, T a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. [a] -> [a]
reverse)
          ((a -> [a] -> [a]) -> [a] -> [[a]] -> [[a]]
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 (:)
             ((Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> a
forall a. C a => Integer -> a
fromInteger [Integer
1 .. Integer
maxCoeff])
             (Int -> [a] -> [[a]]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
d ([a] -> [[a]]) -> [a] -> [[a]]
forall a b. (a -> b) -> a -> b
$ (Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> a
forall a. C a => Integer -> a
fromInteger [-Integer
maxCoeff .. Integer
maxCoeff]))