{-# 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 NumericPrelude.Base    hiding (const)
import NumericPrelude.Numeric


{- $setup
>>> import qualified MathObj.PowerSeries.Core as PS
>>> import qualified MathObj.PowerSeries as PST
>>> import qualified Test.QuickCheck as QC
>>> import Test.NumericPrelude.Utility (equalTrunc, (/\))
>>> import NumericPrelude.Numeric as NP
>>> import NumericPrelude.Base as P
>>> import Prelude ()
-}


newtype T a = Cons {T a -> [a]
coeffs :: [a]} deriving (Eq (T a)
Eq (T a)
-> (T a -> T a -> Ordering)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> T a)
-> (T a -> T a -> T a)
-> Ord (T a)
T a -> T a -> Bool
T a -> T a -> Ordering
T a -> T a -> T a
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall a. (C a, Ord a) => Eq (T a)
forall a. (C a, Ord a) => T a -> T a -> Bool
forall a. (C a, Ord a) => T a -> T a -> Ordering
forall a. (C a, Ord a) => T a -> T a -> T a
min :: T a -> T a -> T a
$cmin :: forall a. (C a, Ord a) => T a -> T a -> T a
max :: T a -> T a -> T a
$cmax :: forall a. (C a, Ord a) => T a -> T a -> T a
>= :: T a -> T a -> Bool
$c>= :: forall a. (C a, Ord a) => T a -> T a -> Bool
> :: T a -> T a -> Bool
$c> :: forall a. (C a, Ord a) => T a -> T a -> Bool
<= :: T a -> T a -> Bool
$c<= :: forall a. (C a, Ord a) => T a -> T a -> Bool
< :: T a -> T a -> Bool
$c< :: forall a. (C a, Ord a) => T a -> T a -> Bool
compare :: T a -> T a -> Ordering
$ccompare :: forall a. (C a, Ord a) => T a -> T a -> Ordering
$cp1Ord :: forall a. (C a, Ord a) => Eq (T a)
Ord)

{-# INLINE fromCoeffs #-}
fromCoeffs :: [a] -> T a
fromCoeffs :: [a] -> T a
fromCoeffs = [a] -> T a
forall a. [a] -> T a
lift0

{-# INLINE lift0 #-}
lift0 :: [a] -> T a
lift0 :: [a] -> T a
lift0 = [a] -> T a
forall a. [a] -> T a
Cons

{-# INLINE lift1 #-}
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)

{-# INLINE lift2 #-}
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)

{-# INLINE const #-}
const :: a -> T a
const :: a -> T a
const a
x = [a] -> T a
forall a. [a] -> T a
lift0 [a
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 :: (a -> b) -> T a -> T b
fmap a -> b
f (Cons [a]
xs) = [b] -> T b
forall a. [a] -> T a
Cons ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
xs)

{-# INLINE appPrec #-}
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) (String -> ShowS
showString String
"PowerSeries.fromCoeffs " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> ShowS
forall a. Show a => a -> ShowS
shows [a]
xs)


{-# INLINE truncate #-}
truncate :: Int -> T a -> T a
truncate :: Int -> T a -> T a
truncate Int
n = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n)

{- |
Evaluate (truncated) power series.
-}
{-# INLINE evaluate #-}
evaluate :: Ring.C a => T a -> a -> a
evaluate :: T a -> a -> a
evaluate (Cons [a]
y) = [a] -> a -> a
forall a. C a => [a] -> a -> a
Core.evaluate [a]
y

{- |
Evaluate (truncated) power series.
-}
{-# INLINE evaluateCoeffVector #-}
evaluateCoeffVector :: Module.C a v => T v -> a -> v
evaluateCoeffVector :: T v -> a -> v
evaluateCoeffVector (Cons [v]
y) = [v] -> a -> v
forall a v. C a v => [v] -> a -> v
Core.evaluateCoeffVector [v]
y


{-# INLINE evaluateArgVector #-}
evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v
evaluateArgVector :: T a -> v -> v
evaluateArgVector (Cons [a]
y) = [a] -> v -> v
forall a v. (C a v, C v) => [a] -> v -> v
Core.evaluateArgVector [a]
y

{- |
Evaluate approximations that is evaluate all truncations of the series.
-}
{-# INLINE approximate #-}
approximate :: Ring.C a => T a -> a -> [a]
approximate :: T a -> a -> [a]
approximate (Cons [a]
y) = [a] -> a -> [a]
forall a. C a => [a] -> a -> [a]
Core.approximate [a]
y


{- |
Evaluate approximations that is evaluate all truncations of the series.
-}
{-# INLINE approximateCoeffVector #-}
approximateCoeffVector :: Module.C a v => T v -> a -> [v]
approximateCoeffVector :: T v -> a -> [v]
approximateCoeffVector (Cons [v]
y) = [v] -> a -> [v]
forall a v. C a v => [v] -> a -> [v]
Core.approximateCoeffVector [v]
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 :: T a -> v -> [v]
approximateArgVector (Cons [a]
y) = [a] -> v -> [v]
forall a v. (C a v, C v) => [a] -> v -> [v]
Core.approximateArgVector [a]
y


{-
Note that the derived instances only make sense for finite series.
-}

instance (Eq a, ZeroTestable.C a) => Eq (T a) where
   (Cons [a]
x) == :: T a -> T a -> Bool
== (Cons [a]
y) = [a] -> [a] -> Bool
forall a. (Eq a, C a) => [a] -> [a] -> Bool
Poly.equal [a]
x [a]
y

instance (Additive.C a) => Additive.C (T a) where
   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]
Poly.negate
   + :: 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]
Poly.add
   (-)    = ([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]
Poly.sub
   zero :: T a
zero   = [a] -> T a
forall a. [a] -> T a
lift0 []

{- |
prop> QC.choose (1,10) /\ \expon (QC.Positive x) xs -> let xt = x:xs in  equalTrunc 15 (PS.pow (const x) (1 % expon) (PST.coeffs (PST.fromCoeffs xt ^ expon)) ++ repeat zero) (xt ++ repeat zero)
-}
instance (Ring.C a) => Ring.C (T a) where
   one :: T a
one           = a -> T a
forall a. a -> T a
const a
forall a. C a => a
one
   fromInteger :: Integer -> T a
fromInteger Integer
n = a -> T a
forall 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]
Core.mul

instance Vector.C T where
   zero :: T a
zero  = T a
forall a. C a => a
zero
   <+> :: T a -> T a -> T a
(<+>) = T a -> T a -> T a
forall a. C a => a -> a -> a
(+)
   *> :: a -> T a -> T a
(*>)  = a -> T a -> T a
forall (v :: * -> *) a. (Functor v, C a) => a -> v a -> v a
Vector.functorScale

instance (Module.C a b) => Module.C a (T b) where
   *> :: a -> T b -> T b
(*>) a
x = ([b] -> [b]) -> T b -> T b
forall a. ([a] -> [a]) -> T a -> T a
lift1 (a
x a -> [b] -> [b]
forall a v. C a v => a -> v -> v
*>)

instance (Field.C a, Module.C a b) => VectorSpace.C a (T b)


instance (Field.C a) => Field.C (T a) where
   / :: 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]
Core.divide


instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where
   divMod :: T a -> T a -> (T a, T a)
divMod (Cons [a]
x) (Cons [a]
y) =
      let ([a]
d,[a]
m) = [a] -> [a] -> ([a], [a])
forall a. (C a, C a) => [a] -> [a] -> ([a], [a])
Core.divMod [a]
x [a]
y
      in  ([a] -> T a
forall a. [a] -> T a
Cons [a]
d, [a] -> T a
forall a. [a] -> T a
Cons [a]
m)


instance (Ring.C a) => Differential.C (T a) where
   differentiate :: T a -> T a
differentiate = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 [a] -> [a]
forall a. C a => [a] -> [a]
Core.differentiate


instance (Algebraic.C a) => Algebraic.C (T a) where
   sqrt :: T a -> T a
sqrt   = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> [a] -> [a]
Core.sqrt a -> a
forall a. C a => a -> a
Algebraic.sqrt)
   T a
x ^/ :: T a -> Rational -> T a
^/ Rational
y = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> a -> [a] -> [a]
forall a. C a => (a -> a) -> a -> [a] -> [a]
Core.pow (a -> Rational -> a
forall a. C a => a -> Rational -> a
Algebraic.^/ Rational
y)
                       (Rational -> a
forall a. C a => Rational -> a
fromRational' Rational
y)) T a
x


instance (Transcendental.C a) =>
             Transcendental.C (T a) where
   pi :: T a
pi = a -> T a
forall a. a -> T a
const a
forall a. C a => a
Transcendental.pi
   exp :: T a -> T a
exp = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> [a] -> [a]
Core.exp a -> a
forall a. C a => a -> a
Transcendental.exp)
   sin :: T a -> T a
sin = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> (a, a)) -> [a] -> [a]
forall a. C a => (a -> (a, a)) -> [a] -> [a]
Core.sin a -> (a, a)
forall a. C a => a -> (a, a)
Core.sinCosScalar)
   cos :: T a -> T a
cos = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> (a, a)) -> [a] -> [a]
forall a. C a => (a -> (a, a)) -> [a] -> [a]
Core.cos a -> (a, a)
forall a. C a => a -> (a, a)
Core.sinCosScalar)
   tan :: T a -> T a
tan = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> (a, a)) -> [a] -> [a]
forall a. C a => (a -> (a, a)) -> [a] -> [a]
Core.tan a -> (a, a)
forall a. C a => a -> (a, a)
Core.sinCosScalar)
   T a
x ** :: T a -> T a -> T a
** T a
y = T a -> T a
forall a. C a => a -> a
Transcendental.exp (T a -> T a
forall a. C a => a -> a
Transcendental.log T a
x T a -> T a -> T a
forall a. C a => a -> a -> a
* T a
y)
                {- This order of multiplication is especially fast
                   when y is a singleton. -}
   log :: T a -> T a
log  = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> [a] -> [a]
Core.log  a -> a
forall a. C a => a -> a
Transcendental.log)
   asin :: T a -> T a
asin = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> (a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> (a -> a) -> [a] -> [a]
Core.asin a -> a
forall a. C a => a -> a
Algebraic.sqrt a -> a
forall a. C a => a -> a
Transcendental.asin)
   acos :: T a -> T a
acos = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> (a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> (a -> a) -> [a] -> [a]
Core.acos a -> a
forall a. C a => a -> a
Algebraic.sqrt a -> a
forall a. C a => a -> a
Transcendental.acos)
   atan :: T a -> T a
atan = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 ((a -> a) -> [a] -> [a]
forall a. C a => (a -> a) -> [a] -> [a]
Core.atan a -> a
forall a. C a => a -> a
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 :: T a -> T a -> T a
compose (Cons [])    (Cons []) = [a] -> T a
forall a. [a] -> T a
Cons []
compose (Cons (a
x:[a]
_)) (Cons []) = [a] -> T a
forall a. [a] -> T a
Cons [a
x]
compose (Cons [a]
x) (Cons (a
y:[a]
ys)) =
   if a -> Bool
forall a. C a => a -> Bool
isZero a
y
     then [a] -> T a
forall a. [a] -> T a
Cons ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
Core.compose [a]
x [a]
ys)
     else String -> T a
forall a. HasCallStack => String -> a
error String
"PowerSeries.compose: inner series must not have an absolute term."

shrink :: Ring.C a => a -> T a -> T a
shrink :: a -> T a -> T a
shrink = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 (([a] -> [a]) -> T a -> T a)
-> (a -> [a] -> [a]) -> a -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
Poly.shrink

dilate :: Field.C a => a -> T a -> T a
dilate :: a -> T a -> T a
dilate = ([a] -> [a]) -> T a -> T a
forall a. ([a] -> [a]) -> T a -> T a
lift1 (([a] -> [a]) -> T a -> T a)
-> (a -> [a] -> [a]) -> a -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> [a] -> [a]
forall a. C a => a -> [a] -> [a]
Poly.dilate