```{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}
{- |
Copyright   :  (c) Henning Thielemann 2004-2005

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

For a multi-set of numbers,
we describe a sequence of the sums of powers of the numbers in the set.
These can be easily converted to polynomials and back.
Thus they provide an easy way for computations on the roots of a polynomial.
-}
module MathObj.PowerSum where

import qualified MathObj.Polynomial  as Poly
import qualified MathObj.PowerSeries as PS

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

import Algebra.Module((*>))

import Control.Monad(liftM2)
import qualified Data.List as List
import NumericPrelude.List (shearTranspose, sieve)

import PreludeBase as P hiding (const)
import NumericPrelude as NP

newtype T a = Cons {sums :: [a]}

{- * Conversions -}

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

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

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

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

{- Newton-Girard formulas,  cf. Modula-3: arithmetic/RootBasic.mg
s'/s = p -}

{-
s[k] - the elementary symmetric polynomial of degree k
p[k] - sum of the k-th power

s(x0,x1,x2) = 1
s(x0,x1,x2) = x0+x1+x2
s(x0,x1,x2) = x0*x1+x1*x2+x2*x0
s(x0,x1,x2) = x0*x1*x2
s(x0,x1,x2) = 0

p(x0,x1,x2) =  1   +  1   +  1
p(x0,x1,x2) = x0   + x1   + x2
p(x0,x1,x2) = x0^2 + x1^2 + x2^2
p(x0,x1,x2) = x0^3 + x1^3 + x2^3
p(x0,x1,x2) = x0^4 + x1^4 + x2^4

s(t) := s + s*t + s*t^2 + ...
p(t) :=        p*t + p*t^2 + ...

Then it holds
t*s'(t) + p(-t)*s(t) = 0
This can be proven by considering p as sum of geometric series
and differentiating s in the root-wise factored form.

Note that we index the coefficients the other way round
and that the coefficients of the polynomial
are not pure elementary symmetric polynomials of the roots
but have alternating signs, too.
-}
fromElemSym :: (Eq a, Ring.C a) => [a] -> [a]
fromElemSym s =
fromIntegral (length s - 1) :
Poly.alternate (divOneFlip s (Poly.differentiate s))

divOneFlip :: (Eq a, Ring.C a) => [a] -> [a] -> [a]
divOneFlip (1:xs) =
let aux (y:ys) = y : aux (ys - Poly.scale y xs)
aux [] = []
in  aux
divOneFlip _ =
error "divOneFlip: first element must be one"

fromElemSymDenormalized :: (Field.C a, ZeroTestable.C a) => [a] -> [a]
fromElemSymDenormalized s =
fromIntegral (length s - 1) :
Poly.alternate (PS.derivedLog s)

toElemSym :: (Field.C a, ZeroTestable.C a) => [a] -> [a]
toElemSym p =
let s' = Poly.mul (Poly.alternate (tail p)) s
s  = Poly.integrate 1 s'
in  s

toElemSymInt :: (Integral.C a, ZeroTestable.C a) => [a] -> [a]
toElemSymInt p =
let s' = Poly.mul (Poly.alternate (tail p)) s
s  = Poly.integrateInt 1 s'
in  s

fromPolynomial :: (Field.C a, ZeroTestable.C a) => Poly.T a -> [a]
fromPolynomial =
let aux s =
fromIntegral (length s - 1) :
Poly.negate (PS.derivedLog s)
in  aux . reverse . Poly.coeffs

elemSymFromPolynomial :: Additive.C a => Poly.T a -> [a]
elemSymFromPolynomial = Poly.alternate . reverse . Poly.coeffs

{- toPolynomial is not possible because this had to consume the whole sum sequence. -}

binomials :: Ring.C a => [[a]]
binomials =  : binomials + map (0:) binomials

{- * Show -}

appPrec :: Int
appPrec  = 10

instance (Show a) => Show (T a) where
showsPrec p (Cons xs) =
showParen (p >= appPrec)
(showString "PowerSum.Cons " . shows xs)

{- * Additive -}

{- Use binomial expansion of (x+y)^n -}
add :: (Ring.C a) => [a] -> [a] -> [a]
add xs ys =
let powers = shearTranspose (Poly.tensorProduct xs ys)
in  zipWith Ring.scalarProduct binomials powers

instance (Ring.C a) => Additive.C (T a) where
zero   = const zero
(+)    = lift2 add
negate = lift1 Poly.alternate

{- * Ring -}

mul :: (Ring.C a) => [a] -> [a] -> [a]
mul xs ys = zipWith (*) xs ys

pow :: Integer -> [a] -> [a]
pow n =
if n<0
then error "PowerSum.pow: negative exponent"
else sieve (fromInteger n)
-- map head . iterate (List.genericDrop (toInteger n))

instance (Ring.C a) => Ring.C (T a) where
one           = const one
fromInteger n = const (fromInteger n)
(*)           = lift2 mul
x^n           = lift1 (pow n) x

{- * Module -}

instance (Module.C a v, Ring.C v) => Module.C a (T v) where
x *> y = lift1 (zipWith (*>) (iterate (x*) one)) y

instance (VectorSpace.C a v, Ring.C v) => VectorSpace.C a (T v)

{- * Field.C -}

instance (Field.C a, ZeroTestable.C a) => Field.C (T a) where
recip = lift1 (fromElemSymDenormalized . reverse . toElemSym)

{- * Algebra -}

root :: (Ring.C a) => Integer -> [a] -> [a]
root n xs =
let upsample m ys =
concat (List.intersperse
(List.genericReplicate (m - 1) zero)
(map (:[]) ys))
in  case compare n 0 of
LT -> upsample (-n) (reverse xs)
GT -> upsample n xs
EQ -> 

instance (Field.C a, ZeroTestable.C a) => Algebraic.C (T a) where
root n = lift1 (fromElemSymDenormalized . root n . toElemSym)

{- given the list of power sums @x1^j + ... + xn^j@
and a power series for the function @f@,
compute the series approximations of @f(x1) + ... + f(xn)@. -}
approxSeries :: Module.C a b => [b] -> [a] -> [b]
approxSeries y x =
scanl (+) zero (zipWith (*>) x y)

{- input lists contain roots -}
propOp :: (Eq a, Field.C a, ZeroTestable.C a) =>
([a] -> [a] -> [a]) -> (a -> a -> a) -> [a] -> [a] -> [Bool]
propOp powerOp op xs ys =
let zs = liftM2 op xs ys
xp = fromPolynomial (Poly.fromRoots xs)
yp = fromPolynomial (Poly.fromRoots ys)
ze = elemSymFromPolynomial (Poly.fromRoots zs)
in  zipWith (==) (toElemSym (powerOp xp yp)) ze
-- Poly.equal (toElemSym (powerOp xp yp)) ze
```