{-# OPTIONS -fno-implicit-prelude #-} {- | 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.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 NumericPrelude.List (takeMatch) import Control.Monad (liftM2) import PreludeBase as P hiding (const) import NumericPrelude as NP newtype T a = Cons {coeffs :: [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] toPolynomial :: Poly.T a -> T a toPolynomial xs = Cons (reverse (Poly.coeffs xs)) fromPolynomial :: T a -> Poly.T a fromPolynomial (Cons xs) = Poly.fromCoeffs (reverse xs) toPowerSums :: (Field.C a, ZeroTestable.C a) => [a] -> [a] toPowerSums = PowerSum.fromElemSymDenormalized fromPowerSums :: (Field.C a, ZeroTestable.C a) => [a] -> [a] fromPowerSums = PowerSum.toElemSym {- | cf. 'MathObj.Polynomial.mulLinearFactor' -} addRoot :: Ring.C a => a -> [a] -> [a] addRoot x yt@(y:ys) = y : (ys + Poly.scale x yt) addRoot _ [] = error "addRoot: list of elementar symmetric terms must consist at least of a 1" fromRoots :: Ring.C a => [a] -> [a] fromRoots = foldl (flip addRoot) [1] liftPowerSum1Gen :: ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) liftPowerSum1Gen fromPS toPS op x = takeMatch x (fromPS (op (toPS x))) liftPowerSum2Gen :: ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a]) liftPowerSum2Gen fromPS toPS op x y = takeMatch (undefined : liftM2 (,) (tail x) (tail y)) (fromPS (op (toPS x) (toPS y))) liftPowerSum1 :: (Field.C a, ZeroTestable.C a) => ([a] -> [a]) -> ([a] -> [a]) liftPowerSum1 = liftPowerSum1Gen fromPowerSums toPowerSums liftPowerSum2 :: (Field.C a, ZeroTestable.C a) => ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a]) liftPowerSum2 = liftPowerSum2Gen fromPowerSums toPowerSums liftPowerSumInt1 :: (Integral.C a, Eq a, ZeroTestable.C a) => ([a] -> [a]) -> ([a] -> [a]) liftPowerSumInt1 = liftPowerSum1Gen PowerSum.toElemSymInt PowerSum.fromElemSym liftPowerSumInt2 :: (Integral.C a, Eq a, ZeroTestable.C a) => ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a]) liftPowerSumInt2 = liftPowerSum2Gen PowerSum.toElemSymInt PowerSum.fromElemSym {- * Show -} appPrec :: Int appPrec = 10 instance (Show a) => Show (T a) where showsPrec p (Cons xs) = showParen (p >= appPrec) (showString "RootSet.Cons " . shows xs) {- * Additive -} {- Use binomial expansion of (x+y)^n -} add :: (Field.C a, ZeroTestable.C a) => [a] -> [a] -> [a] add = liftPowerSum2 PowerSum.add addInt :: (Integral.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> [a] addInt = liftPowerSumInt2 PowerSum.add instance (Field.C a, ZeroTestable.C a) => Additive.C (T a) where zero = const zero (+) = lift2 add negate = lift1 Poly.alternate {- * Ring -} mul :: (Field.C a, ZeroTestable.C a) => [a] -> [a] -> [a] mul = liftPowerSum2 PowerSum.mul mulInt :: (Integral.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> [a] mulInt = liftPowerSumInt2 PowerSum.mul pow :: (Field.C a, ZeroTestable.C a) => Integer -> [a] -> [a] pow n = liftPowerSum1 (PowerSum.pow n) powInt :: (Integral.C a, Eq a, ZeroTestable.C a) => Integer -> [a] -> [a] powInt n = liftPowerSumInt1 (PowerSum.pow n) instance (Field.C a, ZeroTestable.C a) => Ring.C (T a) where one = const one fromInteger n = const (fromInteger n) (*) = lift2 mul x^n = lift1 (pow n) x {- * Field.C -} instance (Field.C a, ZeroTestable.C a) => Field.C (T a) where recip = lift1 reverse {- * Algebra -} instance (Field.C a, ZeroTestable.C a) => Algebraic.C (T a) where root n = lift1 (PowerSum.root n)