{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-} {- | Power series, either finite or unbounded. (zipWith does exactly the right thing to make it work almost transparently.) -} module MathObj.PowerSeries where import qualified MathObj.PowerSeries.Core as Core import qualified MathObj.Polynomial.Core as Poly import qualified Algebra.Differential as Differential import qualified Algebra.IntegralDomain as Integral import qualified Algebra.VectorSpace as VectorSpace import qualified Algebra.Module as Module import qualified Algebra.Vector as Vector import qualified Algebra.Transcendental as Transcendental import qualified Algebra.Algebraic as Algebraic 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 Algebra.Module((*>)) import NumericPrelude.Base hiding (const) import NumericPrelude.Numeric newtype T a = Cons {coeffs :: [a]} deriving (Ord) {-# INLINE fromCoeffs #-} fromCoeffs :: [a] -> T a fromCoeffs = lift0 {-# INLINE lift0 #-} lift0 :: [a] -> T a lift0 = Cons {-# INLINE lift1 #-} lift1 :: ([a] -> [a]) -> (T a -> T a) lift1 f (Cons x0) = Cons (f x0) {-# INLINE lift2 #-} lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a) lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1) {-# INLINE const #-} const :: a -> T a const x = lift0 [x] {- Functor instance is e.g. useful for showing power series in residue rings. @fmap (ResidueClass.concrete 7) (powerSeries [1,4,4::ResidueClass.T Integer] * powerSeries [1,5,6])@ -} instance Functor T where fmap f (Cons xs) = Cons (map f xs) {-# INLINE appPrec #-} appPrec :: Int appPrec = 10 instance (Show a) => Show (T a) where showsPrec p (Cons xs) = showParen (p >= appPrec) (showString "PowerSeries.fromCoeffs " . shows xs) {-# INLINE truncate #-} truncate :: Int -> T a -> T a truncate n = lift1 (take n) {- | Evaluate (truncated) power series. -} {-# INLINE evaluate #-} evaluate :: Ring.C a => T a -> a -> a evaluate (Cons y) = Core.evaluate y {- | Evaluate (truncated) power series. -} {-# INLINE evaluateCoeffVector #-} evaluateCoeffVector :: Module.C a v => T v -> a -> v evaluateCoeffVector (Cons y) = Core.evaluateCoeffVector y {-# INLINE evaluateArgVector #-} evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v evaluateArgVector (Cons y) = Core.evaluateArgVector y {- | Evaluate approximations that is evaluate all truncations of the series. -} {-# INLINE approximate #-} approximate :: Ring.C a => T a -> a -> [a] approximate (Cons y) = Core.approximate y {- | Evaluate approximations that is evaluate all truncations of the series. -} {-# INLINE approximateCoeffVector #-} approximateCoeffVector :: Module.C a v => T v -> a -> [v] approximateCoeffVector (Cons y) = Core.approximateCoeffVector y {- | Evaluate approximations that is evaluate all truncations of the series. -} {-# INLINE approximateArgVector #-} approximateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> [v] approximateArgVector (Cons y) = Core.approximateArgVector y {- Note that the derived instances only make sense for finite series. -} instance (Eq a, ZeroTestable.C a) => Eq (T a) where (Cons x) == (Cons y) = Poly.equal x y instance (Additive.C a) => Additive.C (T a) where negate = lift1 Poly.negate (+) = lift2 Poly.add (-) = lift2 Poly.sub zero = lift0 [] instance (Ring.C a) => Ring.C (T a) where one = const one fromInteger n = const (fromInteger n) (*) = lift2 Core.mul instance Vector.C T where zero = zero (<+>) = (+) (*>) = Vector.functorScale instance (Module.C a b) => Module.C a (T b) where (*>) x = lift1 (x *>) instance (Field.C a, Module.C a b) => VectorSpace.C a (T b) instance (Field.C a) => Field.C (T a) where (/) = lift2 Core.divide instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where divMod (Cons x) (Cons y) = let (d,m) = Core.divMod x y in (Cons d, Cons m) instance (Ring.C a) => Differential.C (T a) where differentiate = lift1 Core.differentiate instance (Algebraic.C a) => Algebraic.C (T a) where sqrt = lift1 (Core.sqrt Algebraic.sqrt) x ^/ y = lift1 (Core.pow (Algebraic.^/ y) (fromRational' y)) x instance (Transcendental.C a) => Transcendental.C (T a) where pi = const Transcendental.pi exp = lift1 (Core.exp Transcendental.exp) sin = lift1 (Core.sin Core.sinCosScalar) cos = lift1 (Core.cos Core.sinCosScalar) tan = lift1 (Core.tan Core.sinCosScalar) x ** y = Transcendental.exp (Transcendental.log x * y) {- This order of multiplication is especially fast when y is a singleton. -} log = lift1 (Core.log Transcendental.log) asin = lift1 (Core.asin Algebraic.sqrt Transcendental.asin) acos = lift1 (Core.acos Algebraic.sqrt Transcendental.acos) atan = lift1 (Core.atan Transcendental.atan) {- | It fulfills @ evaluate x . evaluate y == evaluate (compose x y) @ -} compose :: (Ring.C a, ZeroTestable.C a) => T a -> T a -> T a compose (Cons []) (Cons []) = Cons [] compose (Cons (x:_)) (Cons []) = Cons [x] compose (Cons x) (Cons (y:ys)) = if isZero y then Cons (Core.compose x ys) else error "PowerSeries.compose: inner series must not have an absolute term."