PowerSeries2.hs

{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}

{- |
Two-variate power series.
-}

module MathObj.PowerSeries2 where

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

import qualified Algebra.Differential   as Differential
import qualified Algebra.Vector         as Vector
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 qualified NumericPrelude as NP
import qualified PreludeBase as P

import Data.List (isPrefixOf)
import NumericPrelude.List (compareLength)

import PreludeBase    hiding (const)
import NumericPrelude hiding (negate, stdUnit,
                              sqrt, exp, log,
                              sin, cos, tan, asin, acos, atan)

{- |
In order to handle both variables equivalently
we maintain a list of coefficients for terms of the same total degree.
That is

> eval [[a], [b,c], [d,e,f]] (x,y) ==
>    a + b*x+c*y + d*x^2+e*x*y+f*y^2

Although the sub-lists are always finite and thus are more like polynomials than power series,
division and square root computation are easier to implement for power series.
-}
newtype T a = Cons {coeffs :: Core a} deriving (Ord)

type Core a = [[a]]

isValid :: [[a]] -> Bool
isValid = flip isPrefixOf [1..] . map length

check :: [[a]] -> [[a]]
check xs =
   zipWith (\n x ->
      if compareLength n x == EQ
        then x
        else error "PowerSeries2.check: invalid length of sub-list")
     (iterate (():) [()]) xs


fromCoeffs :: [[a]] -> T a
fromCoeffs  =  Cons . check

fromPowerSeries0 :: Ring.C a => PS.T a -> T a
fromPowerSeries0 x =
   fromCoeffs $
   zipWith (:) (PS.coeffs x) $
   iterate (0:) []

fromPowerSeries1 :: Ring.C a => PS.T a -> T a
fromPowerSeries1 x =
   fromCoeffs $
   zipWith (++) (iterate (0:) []) $
   map (:[]) (PS.coeffs x)


lift0 :: Core a -> T a
lift0 = Cons

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

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


lift0fromPowerSeries :: [PS.T a] -> Core a
lift0fromPowerSeries = map PS.coeffs

lift1fromPowerSeries :: ([PS.T a] -> [PS.T a]) -> (Core a -> Core a)
lift1fromPowerSeries f x0 = map PS.coeffs (f (map PS.fromCoeffs x0))

lift2fromPowerSeries :: ([PS.T a] -> [PS.T a] -> [PS.T a]) -> (Core a -> Core a -> Core a)
lift2fromPowerSeries f x0 x1 = map PS.coeffs (f (map PS.fromCoeffs x0) (map PS.fromCoeffs x1))


const :: a -> T a
const x = lift0 [[x]]


instance Functor T where
  fmap f (Cons xs) = Cons (map (map f) xs)

appPrec :: Int
appPrec  = 10

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


{- * Series arithmetic -}

add, sub :: (Additive.C a) => Core a -> Core a -> Core a
add = PS.add
sub = PS.sub

negate :: (Additive.C a) => Core a -> Core a
negate = PS.negate


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 PS.negate
    (+)    = lift2 PS.add
    (-)    = lift2 PS.sub
    zero   = lift0 []


scale :: Ring.C a => a -> Core a -> Core a
scale = map . (Vector.*>)

mul :: Ring.C a => Core a -> Core a -> Core a
mul = lift2fromPowerSeries PS.mul

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

instance Vector.C T where
   zero  = zero
   (<+>) = (+)
   (*>)  = Vector.functorScale


divide :: (Field.C a) =>
   Core a -> Core a -> Core a
divide = lift2fromPowerSeries PS.divide


instance (Field.C a) => Field.C (T a) where
  (/) = lift2 divide


sqrt :: (Field.C a) =>
   (a -> a) -> Core a -> Core a
sqrt fSqRt = lift1fromPowerSeries $ PS.sqrt (PS.const . (\[x] -> fSqRt x) . PS.coeffs)


instance (Algebraic.C a) => Algebraic.C (T a) where
   sqrt   = lift1 (sqrt Algebraic.sqrt)
--   x ^/ y = lift1 (pow (Algebraic.^/ y)
--                       (fromRational' y)) x


swapVariables :: Core a -> Core a
swapVariables = map reverse


differentiate0 :: (Ring.C a) => Core a -> Core a
differentiate0 =
   swapVariables . differentiate1 . swapVariables

differentiate1 :: (Ring.C a) => Core a -> Core a
differentiate1 = lift1fromPowerSeries $ map Differential.differentiate

integrate0 :: (Field.C a) => [a] -> Core a -> Core a
integrate0 cs =
   swapVariables . integrate1 cs . swapVariables

integrate1 :: (Field.C a) => [a] -> Core a -> Core a
integrate1 = zipWith PS.integrate




{- |
Since the inner series must start with a zero,
the first term is omitted in y.
-}
comp :: (Ring.C a) => [a] -> Core a -> Core a
comp = lift1fromPowerSeries . PS.comp . map PS.const