{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
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.Polynomial.Core as PolyCore
import qualified MathObj.PowerSeries.Core 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 Control.Monad(liftM2)
import qualified Data.List as List
import Data.List.HT (shearTranspose, sieve)

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


newtype T a = Cons {T a -> [a]
sums :: [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]

{- 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[0](x0,x1,x2) = 1
  s[1](x0,x1,x2) = x0+x1+x2
  s[2](x0,x1,x2) = x0*x1+x1*x2+x2*x0
  s[3](x0,x1,x2) = x0*x1*x2
  s[4](x0,x1,x2) = 0

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

  s(t) := s[0] + s[1]*t + s[2]*t^2 + ...
  p(t) :=        p[1]*t + p[2]*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 :: [a] -> [a]
fromElemSym [a]
s =
   Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
s Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
1) a -> [a] -> [a]
forall a. a -> [a] -> [a]
:
      [a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.alternate ([a] -> [a] -> [a]
forall a. (Eq a, C a) => [a] -> [a] -> [a]
divOneFlip [a]
s ([a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.differentiate [a]
s))

divOneFlip :: (Eq a, Ring.C a) => [a] -> [a] -> [a]
divOneFlip :: [a] -> [a] -> [a]
divOneFlip (a
1:[a]
xs) =
   let aux :: [a] -> [a]
aux (a
y:[a]
ys) = a
y a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
aux ([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
y [a]
xs)
       aux [] = []
   in  [a] -> [a]
aux
divOneFlip [a]
_ =
   [Char] -> [a] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"divOneFlip: first element must be one"

fromElemSymDenormalized :: (Field.C a, ZeroTestable.C a) => [a] -> [a]
fromElemSymDenormalized :: [a] -> [a]
fromElemSymDenormalized [a]
s =
   Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
s Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
1) a -> [a] -> [a]
forall a. a -> [a] -> [a]
:
      [a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.alternate ([a] -> [a]
forall a. C a => [a] -> [a]
PS.derivedLog [a]
s)


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

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



fromPolynomial :: (Field.C a, ZeroTestable.C a) => Poly.T a -> [a]
fromPolynomial :: T a -> [a]
fromPolynomial =
   let aux :: [a] -> [a]
aux [a]
s =
          Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
s Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
1) a -> [a] -> [a]
forall a. a -> [a] -> [a]
:
             [a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.negate ([a] -> [a]
forall a. C a => [a] -> [a]
PS.derivedLog [a]
s)
   in  [a] -> [a]
forall a. C a => [a] -> [a]
aux ([a] -> [a]) -> (T a -> [a]) -> T a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> (T a -> [a]) -> T a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> [a]
forall a. T a -> [a]
Poly.coeffs

elemSymFromPolynomial :: Additive.C a => Poly.T a -> [a]
elemSymFromPolynomial :: T a -> [a]
elemSymFromPolynomial = [a] -> [a]
forall a. C a => [a] -> [a]
PolyCore.alternate ([a] -> [a]) -> (T a -> [a]) -> T a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> (T a -> [a]) -> T a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> [a]
forall a. T a -> [a]
Poly.coeffs

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



binomials :: Ring.C a => [[a]]
binomials :: [[a]]
binomials = [a
1] [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
forall a. C a => [[a]]
binomials [[a]] -> [[a]] -> [[a]]
forall a. C a => a -> a -> a
+ ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) [[a]]
forall a. C a => [[a]]
binomials

{- * 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]
"PowerSum.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 :: (Ring.C a) => [a] -> [a] -> [a]
add :: [a] -> [a] -> [a]
add [a]
xs [a]
ys =
   let powers :: [[a]]
powers = [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
shearTranspose ([a] -> [a] -> [[a]]
forall a. C a => [a] -> [a] -> [[a]]
PolyCore.tensorProduct [a]
xs [a]
ys)
   in  ([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
Ring.scalarProduct [[a]]
forall a. C a => [[a]]
binomials [[a]]
powers

instance (Ring.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 => [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 :: (Ring.C a) => [a] -> [a] -> [a]
mul :: [a] -> [a] -> [a]
mul [a]
xs [a]
ys = (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]
xs [a]
ys

pow :: Integer -> [a] -> [a]
pow :: Integer -> [a] -> [a]
pow Integer
n =
   if Integer
nInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<Integer
0
     then [Char] -> [a] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"PowerSum.pow: negative exponent"
     else Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
sieve (Integer -> Int
forall a. C a => Integer -> a
fromInteger Integer
n)
       -- map head . iterate (List.genericDrop (toInteger n))

instance (Ring.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 => [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. Integer -> [a] -> [a]
pow Integer
n) T a
x


{- * Module -}

instance (Module.C a v, Ring.C v) => Module.C a (T v) where
   a
x *> :: a -> T v -> T v
*> T v
y = ([v] -> [v]) -> T v -> T v
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> v -> v) -> [a] -> [v] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> v -> v
forall a v. C a v => a -> v -> v
(*>) ((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)) T v
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 :: T a -> T a
recip = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ([a] -> [a]
forall a. (C a, C a) => [a] -> [a]
fromElemSymDenormalized ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [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 b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
toElemSym)


{- * Algebra -}

root :: (Ring.C a) => Integer -> [a] -> [a]
root :: Integer -> [a] -> [a]
root Integer
n [a]
xs =
   let upsample :: i -> [a] -> [a]
upsample i
m [a]
ys =
          [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
List.intersperse
             (i -> a -> [a]
forall i a. Integral i => i -> a -> [a]
List.genericReplicate (i
m i -> i -> i
forall a. C a => a -> a -> a
- i
1) a
forall a. C a => a
zero)
             ((a -> [a]) -> [a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[]) [a]
ys))
   in  case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Integer
n Integer
0 of
          Ordering
LT -> Integer -> [a] -> [a]
forall i a. (Integral i, C i, C a) => i -> [a] -> [a]
upsample (-Integer
n) ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
xs)
          Ordering
GT -> Integer -> [a] -> [a]
forall i a. (Integral i, C i, C a) => i -> [a] -> [a]
upsample Integer
n [a]
xs
          Ordering
EQ -> [a
1]

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 ([a] -> [a]
forall a. (C a, C a) => [a] -> [a]
fromElemSymDenormalized ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [a] -> [a]
forall a. C a => Integer -> [a] -> [a]
root Integer
n ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. (C a, C a) => [a] -> [a]
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 :: [b] -> [a] -> [b]
approxSeries [b]
y [a]
x =
   (b -> b -> b) -> b -> [b] -> [b]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl b -> b -> b
forall a. C a => a -> a -> a
(+) b
forall a. C a => a
zero ((a -> b -> b) -> [a] -> [b] -> [b]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> b -> b
forall a v. C a v => a -> v -> v
(*>) [a]
x [b]
y)


{- input lists contain roots -}
propOp :: (Eq a, Field.C a, ZeroTestable.C a) =>
   ([a] -> [a] -> [a]) -> (a -> a -> a) -> [a] -> [a] -> [Bool]
propOp :: ([a] -> [a] -> [a]) -> (a -> a -> a) -> [a] -> [a] -> [Bool]
propOp [a] -> [a] -> [a]
powerOp a -> a -> a
op [a]
xs [a]
ys =
   let zs :: [a]
zs = (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 -> a
op [a]
xs [a]
ys
       xp :: [a]
xp = T a -> [a]
forall a. (C a, C a) => T a -> [a]
fromPolynomial ([a] -> T a
forall a. C a => [a] -> T a
Poly.fromRoots [a]
xs)
       yp :: [a]
yp = T a -> [a]
forall a. (C a, C a) => T a -> [a]
fromPolynomial ([a] -> T a
forall a. C a => [a] -> T a
Poly.fromRoots [a]
ys)
       ze :: [a]
ze = T a -> [a]
forall a. C a => T a -> [a]
elemSymFromPolynomial ([a] -> T a
forall a. C a => [a] -> T a
Poly.fromRoots [a]
zs)
   in  (a -> a -> Bool) -> [a] -> [a] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(==) ([a] -> [a]
forall a. (C a, C a) => [a] -> [a]
toElemSym ([a] -> [a] -> [a]
powerOp [a]
xp [a]
yp)) [a]
ze
       -- PolyCore.equal (toElemSym (powerOp xp yp)) ze