dual.hs

{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}

module Number.Dual where

import qualified Algebra.Absolute as Absolute
import qualified Algebra.Additive as Additive
import qualified Algebra.Algebraic as Algebraic
import qualified Algebra.Differential as Differential
import qualified Algebra.Field as Field
import qualified Algebra.Module as Module
import qualified Algebra.Ring as Ring
import qualified Algebra.Transcendental as Trans
import qualified Algebra.Vector as Vector
import qualified Algebra.VectorSpace as VectorSpace
import qualified Algebra.ZeroTestable as ZeroTestable

import Text.Show.HT (showsInfixPrec)
import Text.Read.HT (readsInfixPrec)

import NumericPrelude

infix 6 +:, `Cons`

-- Can I restrict this to Ring.C a?
data T a = Cons { real :: !a, nonreal :: !a } deriving Eq

{-# INLINE epsilon #-}
epsilon :: Ring.C a => T a
epsilon = zero +: one

instance Show a => Show (T a) where
	showsPrec prec (Cons x y) = showsInfixPrec "+:" 6 prec x y

instance Read a => Read (T a) where
	readsPrec prec = readsInfixPrec "+:" 6 prec (+:)

{-# INLINE (+:) #-}
(+:) :: a -> a -> T a
(+:) = Cons

{-# INLINE (-:) #-}
(-:) :: (Additive.C a) => a -> a -> T a
a -: b = Cons a (negate b)

-- These next two functions seem like they should belong in a Functor/Monad,
-- but I don't know how that would work
unwrap :: Ring.C a => (T a -> T b) -> a -> b
unwrap f = real . f . todual

todual :: Additive.C a => a -> T a
todual x = x +: zero

-- The type should be something like (T a -> T a) -> T a -> T a because deriv is a transformation
-- In fact, we should (logically) be able to do (a -> b) -> T a -> T b for all functions once we fix Absolute
deriv :: Additive.C a => (T (T a) -> T (T b)) -> T a -> T b
deriv f x@(Cons a b) = (nonreal . f) (x +: (todual b))

instance Additive.C a => Additive.C (T a) where
	zero = Cons zero zero
	Cons a b + Cons c d = (a + c) +: (b + d)
	Cons a b - Cons c d = (a - c) +: (b - d)

instance Ring.C a => Ring.C (T a) where
	one = Cons one zero
	Cons a b * Cons c d = (a * c) +: (a * d + b * c)

instance Field.C a => Field.C (T a) where
	Cons a b / Cons c d = (a / c) +: ((b * c - a * d)/(c ^ 2))

instance Module.C a b => Module.C a (T b) where
	n *> Cons a b = (n *> a) +: (n *> b)

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

instance Algebraic.C a => Algebraic.C (T a) where
	Cons a b ^/ c = (a ^/ c) +: (fromRational' c * b * a ^/ (c - 1))

instance Trans.C a => Trans.C (T a) where
	pi = Cons pi zero
	exp (Cons a b) = (exp a) +: (b * exp a)
	log (Cons a b) = (log a) +: (b / a)
	sin (Cons a b) = (sin a) +: (b * cos a)
	cos (Cons a b) = (cos a) +: (b * negate (sin a))
	atan (Cons a b) = (atan a) +: (b / (one + a ^ 2))

-- Can I remove the dependency on Field?
-- We need a way to return an indeterminate form in signum when a == zero
instance Absolute.C a => Absolute.C (T a) where
	abs (Cons a b) = (Absolute.abs a) +: (b * (signum a))
	signum (Cons a b) = (signum a) +: zero

instance ZeroTestable.C a => ZeroTestable.C (T a) where
	isZero = isZero . real

instance Ring.C a => Differential.C (T (T a) -> T (T b)) where
	differentiate f = deriv f