Computing derivatives of functions using dual numbers

DerivativesViaDualNumbers.hs

{-# LANGUAGE ViewPatterns #-}
module DerivativesViaDualNumbers where

{--

Dual numbers are an extension of the real numbers, like Complex
numbers, only with the identity that d^2 = 0. They are represented
similarly to complex numbers as well, i.e. a+b*d where 'a' and 'b' are
real numbers.

For any dual numbers d0 and d1, represented by a0+b0*d and a1+b1*d,
the following identities hold:

 Addition:       (a0+b0*d) + (a1+b1*d) == (a0+a1) + (b0+b1)*d
 Multiplication: (a0+b0*d) * (a1+b1*d) == (a0*a1) + (a0*b1 + a1*b0)*d
 Subtraction:    (a0+b0*d) - (a1+b1*d) == (a0-a1) + (b0-b1)*d
 Divison:        (a0+b0*d) / (a1+b1*d) == (a0/a1) + ((a1*b0 - a0*b1)/(a0^2))*d
--}

data Dual = Dual Float Float deriving (Eq, Show)

dual :: Float -> Dual
dual = flip Dual 0.0

infinitesimal :: Dual -> Float
infinitesimal (Dual _ x)  = x

instance Num Dual where
  (Dual a0 b0) + (Dual a1 b1) = Dual (a0 + a1) (b0 + b1)
  (Dual a0 b0) * (Dual a1 b1) = Dual (a0 * a1) ((a0 * b1) + (a1 * b0))
  (Dual a0 b0) - (Dual a1 b1) = Dual (a0 - a1) (b0 - b1)

  fromInteger = dual . fromIntegral

instance Fractional Dual where
  (Dual a  b)  / (Dual c  d)  = Dual p1 p2
    where p1 = a/c
          p2 = (c*b - a*d)/(c^2)

deriv :: (Dual -> Dual) -> Float -> Float
deriv f (dual -> x) = infinitesimal $ f (x + d)
  where d = Dual 0.0 1.0
  

f1 :: Num a => a -> a
f1 x = (x + 2) * (x + 1)

-- The derivative of f1, which is 2x+3
f1' :: Float -> Float
f1' = deriv f1

f2 :: Num a => a -> a
f2 x = x*x

-- The derivative of f2, i.e. 2*x
f2' :: Float -> Float
f2' = deriv f2