Frustration My trusty MacBook suddenly died, :( It can no longer boot, and the old way of getting an OpenFirmware prompt doesn't seem to work. Oh well.

Fixed precision, an update So I was a bit sloppy in my last post. When doing arithmetic it was performed exactly using Rational and not truncated according to the epsilon for the type. So, for instance, computing 4/10 + 4/10 with type Fixed Eps1 would give the answer 1. While this might seem nice, it's also wrong if every operation would be performed with epsilon 1, since 4/10 would be 0, and 0+0 is 0. So I'll amend my implementation a little.

instance (Epsilon e) => Num (Fixed e) where
    (+) = lift2 (+)
    (-) = lift2 (-)
    (*) = lift2 (*)
    negate (F x) = F (negate x)
    abs (F x) = F (abs x)
    signum (F x) = F (signum x)
    fromInteger = F . fromInteger

instance (Epsilon e) => Fractional (Fixed e) where
    (/) = lift2 (/)
    fromRational x = r
        where r = F $ approx x (getEps r)

lift2 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e
 -> Fixed e
lift2 op fx@(F x) (F y) = F $ approx (x `op` y) (getEps fx)

approx :: Rational -> Rational -> Rational
approx x eps = approxRational (x + eps/2) eps

So after each operation we add half an epsilon (so we get rounding instead of truncation) and call the magical approxRational to get the closest rational within an epsilon.

Overloading Haskell numbers, part 3, Fixed Precision. Fixed precision numbers can be handy to have sometimes. Haskell only provides integral, rational, and floating point in the Prelude, but I want more! I want to define a type of fixed precision numbers, i.e., numbers that have a certain number of decimals. Rational numbers are handy to represent these numbers so here we go.

newtype Fixed = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)

This uses an GHC extension that allows you to derive anything from the representation type so we don't have to define all those tedious instances. But wait, where is the fixed precision? What I just defined is just the Rational numbers with a different name. I really want a type which is parametrized over what precision it has. OK, so here we go again:

newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)

Here e is meant to be a type that specifies the precision. Let's capture those types that we allow as the precision as a type class:

class Epsilon e where
    eps :: e -> Rational

The idea here is that the actual argument to the eps function is not used, it's just a proxy and the type carries all the information. Here's an example:

data Eps1
instance Epsilon Eps1 where
    eps _ = 1

The eps returns 1 for the Eps1 type. So the intended meaning is that the type Fixed Eps1 are numbers within an epsilon of the right answer, and epsilon is 1. Similarly

data Prec50
instance Epsilon Prec50 where
    eps _ = 1e-50

The type Fixed Prec50 has an epsilon of 1e-50, i.e., 50 decimals. Btw, notice that what looks like a floating point literal is actually a rational number; one of the very clever decisions in the original Haskell design. We can also define type constructors like

data EpsDiv10 p
instance (Epsilon e) => Epsilon (EpsDiv10 e) where
    eps e = eps (un e) / 10
       where un :: EpsDiv10 e -> e
             un = undefined

The type Fixed (EpsDiv10 Prec50) has an epsilon of 1e-50/10 = 1e-51. Note how the function un is just a dummy to get the e. If we had used the scoped type variable extension we could have avoided this function. OK, so we can compute with these fixed precision numbers, but where does the epsilon enter? Nowhere, so far. All arithmetic is performed exactly by the rationals, but our first use of the epsilon is in the show function.

instance (Epsilon e) => Show (Fixed e) where
    showsPrec = showSigned showFixed
      where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r (eps $ un f)
              where q :: Integer
             (q, r) = properFraction x
      un :: Fixed e -> e
             un = undefined
             decimals a e | e >= 1 = ""
                          | otherwise = intToDigit b : decimals c (10 * e)
                 where (b, c) = properFraction (10 * a)

To print a fixed point number we first print the integer part. The properFraction method handily returns the integer part and some left over fraction that is <1. We then print decimal after decimal by multiplying the fraction by 10, converting the integer part to a digit, and then recursing. While generating these decimals we also multiply the epsilon by 10 each time and stop printing when epsilon is >= 1, because then we don't need any more decimals. All right, lets put it to the test:

Data.Number.Fixed> 1/3 :: Fixed Prec50
0.33333333333333333333333333333333333333333333333333
Data.Number.Fixed> 3/8 :: Fixed Prec10
0.3750000000
Data.Number.Fixed> (1 :: Fixed Prec10) + (0.5 :: Fixed Prec50)
    Couldn't match expected type `Prec10'
           against inferred type `Prec50'
      Expected type: Fixed Prec10
      Inferred type: Fixed Prec50

As we want, we can't mix different fixed point numbers with different precisions. Sometimes we do want to convert, so let's provide such a function

convertFixed :: Fixed e -> Fixed f
convertFixed (F x) = F x

It looks like this function does nothing, and it's true. It just converts the type, and since the epsilon type is not involved in the number itself it's trivial. Well, basic arithmetic isn't that much fun. What about the transcendental functions? Luckily Jan Skibinski has already implemented all of them. They compute the various function to with an epsilon using continued fractions. The instance declaration looks something like:

instance (Epsilon e) => Floating (Fixed e) where
    sqrt = toFixed1 F.sqrt
    exp = toFixed1 F.exp
    ...

toFixed1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e
toFixed1 f x@(F r) = F (f (eps (un x)) r)
  where un :: Fixed e -> e
        un = undefined

And all the clever code I stole from Jan is in the module F, see link below. Now it's getting more fun.

Data.Number.Fixed> pi :: Fixed Prec50
3.14159265358979323846264338327950288419716939937510
Data.Number.Fixed> exp 10 :: Fixed Prec10
22026.4657881843
Data.Number.Fixed> log 2 :: Fixed Prec50
0.69314718055994530941723212145817656807550013436026
Data.Number.Fixed> exp (sin (sqrt (54/42))) :: Fixed Prec10
2.4745696343
Data.Number.Fixed> sqrt 2 :: Fixed Prec500
1.41421356237309504880168872420969807856967187537694807317667973799073247
8462107038850387534327641572735013846230912297024924836055850737212644121
4970999358314132226659275055927557999505011527820605714701095599716059702
7453459686201472851741864088919860955232923048430871432145083976260362799
5251407989687253396546331808829640620615258352395054745750287759961729835
5752203375318570113543746034084988471603868999706990048150305440277903164
5424782306849293691862158057846311159666871301301561856898723724

Data.Number.Fixed Data.Complex> let i = 0 :+ 1 :: Complex (Fixed Prec50) 
Data.Number.Fixed Data.Complex> exp (i*pi)
(-0.99999999999999999999999999999999999999999999999999) :+ 0.00000000000000000000000000000000000000000000000000

Naturally, we need to package up all these definitions in a module that will keep the Fixed type abstract, and you can also keep the Epsilon class and the various Epsilon types abstract if you provide enough building blocks. Source at http://www.augustsson.net/Darcs/Data.Number/.

Following a suggestion from Cale Gibbard I added a convenient feature to Djinn. You can now define type classes and have contexts on functions. As with Djinn in general there is no polymorphic instantiation, so the methods of a class must not mention any type variables, but the class parameters. This makes the whole thing rather limited, but it was a quick hack. Sample session

Welcome to Djinn version 2007-04-20.
Type :h to get help.
Djinn> f ? (Eq a) => a -> a -> (b,c) -> Either b c
f :: (Eq a) => a -> a -> (b, c) -> Either b c
f a b =
    case a == b of
    False -> \ (c, _) -> Left c
    True -> \ (_, d) -> Right d

The Djinn source is at darcs.augustsson.net/Darcs/Djinn as usual. Thanks Cale!

Overloading Haskell numbers, part 2, Forward Automatic Differentiation. I will continue my overloading by some examples that have been nicely illustrated by an article by Jerzy Karczmarczuk. And blogged about by sigfpe. But at least I'll end this entry with a small twist I've not seen before. When computing the derivative of a function you normally do by either symbolic derivation, or by a numerical approximation. Say that you have a function     f(x) = x² + 1 and you want to know the derivative at x=5. Doing it symbolically you first get f'     f'(x) = 2x using high school calculus (maybe they don't teach it in high school anymore?), and then you plug in 5     f'(5) = 2*5 = 10 Computing it by numeric differentiation you compute     f'(x) = (f(x+h) - f(x)) / h for some small h. Let's pick h=1e-5, and we get f'(5) = 10.000009999444615. Close, but not that good. So why don't we always use the symbolic method? Well, some functions are not that easy to differentiate. Take this one     g x = if abs (x - 0.7) < 0.4 then x else g (cos x) What's the derivative? Well, it's tricky because this is not really a proper definition of g. It's an equation that if solved will yield a definition of g. And like equations in general, it could have zero, one, or many solutions. (If we happen to use CPOs there is always a unique smallest solution which is what programs compute, as if by magic.) If you think g is contrived, lets pick a different example: computing the square root with Newton-Raphson.

sqr x = convAbs $ iterate improve 1
  where improve r = (r + x/r) / 2
        convAbs (x1:x2:_) | abs (x1-x2) < 1e-10 = x2
 convAbs (_:xs) = convAbs xs

So symbolic is not so easy here, and numeric differentiation is not very accurate. But there is a third way! Automatic differentiation. The idea behind AD is that instead of computing with with just numbers, we instead compute with pairs of numbers. The first component is the normal number, and the second component is the derivative. What are the rules for these numbers? Let's look at addition     (x, x') + (y, y') = (x+y, x'+y') To add two numbers you just add the regular part and the derivatives. For multiplication you have to remember how to compute the derivative of a product:     (f(x)*g(x))' = f(x)*g'(x) + f'(x)*g(x) So for our pairs we get     (x, x') * (y, y') = (x*y, x*y' + x'*y) i.e., first the regular product, then the derivative according to the recipe above. Let's see how it works on     f(x) = x² + 1 We want the derivative at x=5. So what is the pair we use for x? It is (5, 1). Why? Well it has to be 5 for the regular part, and since this represents x and the derivative of x is 1, the pair is (5, 1). In the right hand side for f we need to replace 1 by (1,0), since the derivative of a constant is 0. So then we get     f (5,1) = (5,1)*(5,1) + (1,0) = (26,10) using the rules above. And look! There is the normal result, 26, as well as the derivative, 10. Let's turn this into Haskell, using the type PD to hold a pair of Doubles

data PD = P Double Double deriving (Eq, Ord, Show)
instance Num PD where
    P x x' + P y y' = P (x+y) (x'+y')
    P x x' - P y y' = P (x-y) (x'-y')
    P x x' * P y y' = P (x*y) (x*y' + y'*x)
    fromInteger i   = P (fromInteger i) 0

A first observation is that there is nothing Double specific in this definitions; it would work for any Num. So we can change it to

data PD a = P a a deriving (Eq, Ord, Show)
instance Num a => Num (PD a) where ...

Let's also add abs&signum and the Fractional instance

    ...
    abs (P x x') = P (abs x) (signum x * x')
    signum (P x x') = P (signum x) 0

instance Fractional a => Fractional (PD a) where
    P x x' / P y y' = P (x / y) ( (x'*y - x*y') / (y * y))
    fromRational r  = P (fromRational r) 0

We can now try the sqr example

Main> sqr (P 9 1)
P 3.0 0.16666666666666666

The derivative of x**0.5 is 0.5*x**(-0.5), i.e., 0.5*9**(-0.5) = 0.5/3 = 0.16666666666666666. So we got the right answer. BTW, if you want to be picky the derivative of signum is not 0. The signum function makes a jump from -1 to 1 at 0. So the "proper" value would be 2*dirac, if dirac is a Dirac pulse. But since we don't have numbers with Dirac pulses (yet), I'll just pretend the derivative is 0 everywhere. The very clever insight that Jerzy had was that when doing these numbers in Haskell there is no need to limit yourself to just the first derivative. Since Haskell is lazy we can easily keep an infinite list of all derivatives instead of just the first one. Let's look at how that definition looks. It's very similar to what we just did. But instead of the derivative being just a number, it's now one of our new numbers with a value, and all derivatives... Since we are now dealing with an infinite data structure we need to define our own show, (==), etc.

data Dif a = D a (Dif a)

val (D x _) = x

df (D _ x') = x'

dVar x = D x 1

instance (Show a) => Show (Dif a) where
    show x = show (val x) 

instance (Eq a) => Eq (Dif a) where
    x == y  =  val x == val y

instance (Ord a) => Ord (Dif a) where
    x `compare` y  =  val x `compare` val y

instance (Num a) => Num (Dif a) where
    D x x' + D y y'          =  D (x + y) (x' + y')
    D x x' - D y y'          =  D (x - y) (x' - y')
    p@(D x x') * q@(D y y')  =  D (x * y) (x' * q + p * y')
    fromInteger i            =  D (fromInteger i) 0
    abs p@(D x x')           =  D (abs x) (signum p * x')
    signum (D x _)           =  D (signum x) 0

instance (Fractional a) => Fractional (Dif a) where
    recip (D x x') = ip
 where ip = D (recip x) (-x' * ip * ip)
    fromRational r = D (fromRational r) 0

This looks simple, but it's rather subtle. For instance, take the 0 in the definition of fromInteger. It's actually of Dif type, so it's a recursive call to fromInteger. So let's try with our sqr function again, this time computing up to the third derivative. The dVar is used to create a value for "variable" where we want to differentiate.

Main> sqr $ dVar 9
3.0
Main> df $ sqr $ dVar 9
0.16666666666666669
Main> df $ df $ sqr $ dVar 9
-9.259259259259259e-3
Main> df $ df $ df $ sqr $ dVar 9
1.5432098765432098e-3

And the transcendentals in a similar way:

lift (f : f') p@(D x x') = D (f x) (x' * lift f' p)

instance (Floating a) => Floating (Dif a) where
    pi               = D pi 0
    exp (D x x')     = r where r = D (exp x) (x' * r)
    log p@(D x x')   = D (log x) (x' / p)
    sqrt (D x x')    = r where r = D (sqrt x) (x' / (2 * r))
    sin              = lift (cycle [sin, cos, negate . sin, negate . cos])
    cos              = lift (cycle [cos, negate . sin, negate . cos, sin])
    acos p@(D x x')  = D (acos x) (-x' / sqrt(1 - p*p))
    asin p@(D x x')  = D (asin x) ( x' / sqrt(1 - p*p))
    atan p@(D x x')  = D (atan x) ( x' / (p*p - 1))
    sinh x           = (exp x - exp (-x)) / 2
    cosh x           = (exp x + exp (-x)) / 2
    asinh x          = log (x + sqrt (x*x + 1))
    acosh x          = log (x + sqrt (x*x - 1))
    atanh x          = (log (1 + x) - log (1 - x)) / 2

And why not try the function g we defined above?

Main> g 10
0.6681539175313869
Main> g (dVar 10)
0.6681539175313869
Main> df $ g (dVar 10)
0.4047642621121782
Main> df $ df $ g (dVar 10)
0.4265424381635987
Main> df $ df $ df $ g (dVar 10)
-1.4395397945007182

It all works very nicely. So now when we can compute the derivative of a function, let's define something somewhat more interesting with it. Let's revisit the sqr function again. It uses Newton-Raphson to find the square root. How does Newton-Raphson actually work? Given a differentiable function, f(x), it finds a zero by starting with some x₁ and then iterating     x_n+1 = x_n - f(x_n)/f'(x_n) until we meet some convergence criterion. Using this, let's define a function that finds a zero of another function:

findZero f = convRel $ cut $ iterate step start
    where step x = x - val fx / val (df fx) where fx = f (dVar x)

   start = 1  -- just some value
   epsilon = 1e-10
   cut = (++ error "No convergence in 1000 steps") . take 1000
   convRel (x1:x2:_) | x1 == x2 || abs (x1+x2) / abs (x1-x2) > 1/epsilon = x2
   convRel (_:xs) = convRel xs

The only interesting part is the step function that does one iteration with Newton-Raphson. It computes f x and then divides the normal value with the derivative. We then produce the infinite list of approximations using step, then cut it of at some point (in case it doesn't converge), and then we look down the list for two values that are within some relative epsilon. And it even seems to work.

Main> findZero (\x -> x*x - 9)
3.0
Main> findZero (\x -> sin x - 0.5)
0.5235987755982989
Main> sin it
0.5
Main> findZero (\x -> x*x + 9)
*** Exception: No convergence in 1000 steps
Main> findZero (\x -> sqr x - 3)
9.0

Note how it finds a zero of the sqr function which is actually using recursion internally to compute the square root. So now we can compute numerical derivatives. But wait! We also have symbolic numbers. Can we combine them? Of course, that is the power of polymorphism. Let's load up both modules:

Data.Number.Symbolic Dif3> let x :: Num a => Dif (Sym a); x = dVar (var "x")
Data.Number.Symbolic Dif3> df $ x*x
x+x
Data.Number.Symbolic Dif3> df $ sin x
cos x
Data.Number.Symbolic Dif3> df $ sin (exp (x - 4) * x)
(exp (-4.0+x)*x+exp (-4.0+x))*cos (exp (-4.0+x)*x)
Data.Number.Symbolic Dif3> df $ df $ sin (exp (x - 4) * x)
(exp (-4.0+x)*x+exp (-4.0+x)+exp (-4.0+x))*cos (exp (-4.0+x)*x)+(exp (-4.0+x)*x+exp (-4.0+x))*(exp (-4.0+x)*x+exp (-4.0+x))*(-sin (exp (-4.0+x)*x))

We define x to be a differentiable number, "the variable", over symbolic numbers, over some numbers. And then we just happily use df to get the differentiated versions. So we set out to compute numeric derivatives, and we got these for free. Not too bad. One final note, the Dif type is defined above can be made more efficient by not keeping all the infinite tails with 0 derivatives around. In a real module for this, you'd want to make this optimization. [Edit: fixed typo.]

Overloading Haskell numbers, part 1, symbolic expressions. Haskell's overloaded numerical classes can be (ab)used to do some symbolic maths. This is in no way a new discovery, but I thought I'd write a few lines about it anyway since I've been playing with it the last few days. First we need a data type to represent expressions. We want constants, variables, and function applications. But we don't want to fix the type of the constants, so that will be a parameter to the type.

data Sym a = Con a | Var String | App String [Sym a]
    deriving (Eq, Show)

And we also take the opportunity to derive Eq and Show. Now we can actually claim that the type Sym N is a number if N is a number. Let do it:

instance (Num a) => Num (Sym a) where
    x + y         = App "+" [x, y]
    x - y         = App "-" [x, y]
    x * y         = App "*" [x, y]
    negate x      = App "negate" [x]
    abs    x      = App "abs"    [x]
    signum x      = App "signum" [x]
    fromInteger x = Con (fromInteger x)

A small interactive session shows that we are on the right track.

Sym1> let x = Var "x"
Sym1> x*x + 5
App "+" [App "*" [Var "x",Var "x"],Con 5]
Sym1> x*x*x
App "*" [App "*" [Var "x",Var "x"],Var "x"]
Sym1> 2 + 3 :: Sym Int
App "+" [Con 2,Con 3]

We can type in normal looking expressions, but when they are printed the Show instance is used so we get to see the raw syntax tree. That has its uses, but it gets old quickly. We want a pretty printer. To get the precedences right we need to define showsPrec and pass it the right arguments. It's a little tedious, but nothing strange.

instance (Show a) => Show (Sym a) where
    showsPrec p (Con c) = showsPrec p c
    showsPrec _ (Var s) = showString s
    showsPrec p (App op@(c:_) [x, y]) | not (isAlpha c) =
        showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y)
        where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [
                   ("**", (9,8,8)),
                   ("/",  (7,7,8)),
                   ("*",  (7,7,8)),
                   ("+",  (6,6,7)),
                   ("-",  (6,6,7))]
    showsPrec p (App "negate" [x]) =
        showParen (p>=6) (showString "-" . showsPrec 7 x)
    showsPrec p (App f xs) =
        showParen (p>10) (foldl (.) (showString f)
                                (map (\ x -> showChar ' ' . showsPrec 11 x) xs))

Let's try the same examples again:

Sym2> let x = var "x"
Sym2> x*x + 5
x*x+5
Sym2> x*x*x
x*x*x
Sym2> 2 + 3 :: Sym Int
2+3

Look we can type expressions and get them back again! The instance Num (Sym a) isn't too bad, the only fishy thing about it is the Eq superclass that is required for Num. We have Eq for Sym, but it doesn't really behave like it should. E.g., the expression 'x==1' would come out as False since the syntax trees are not equal. But this isn't really what we would like, ideally (==) would also turn into something symbol, but that is impossible with the standard Prelude. Let's make some more instances. A few of these definitions are just there to appease the Haskell numerical hierarchy and supply some operations it need.

instance (Fractional a) => Fractional (Sym a) where
    x / y          = App "/" [x, y]
    fromRational x = Con (fromRational x)

instance (Real a) => Real (Sym a) where
    toRational (Con c) = toRational c

instance (RealFrac a) => RealFrac (Sym a) where
    properFraction (Con c) = (i, Con c') where (i, c') = properFraction c

instance (Floating a) => Floating (Sym a) where
    pi = App "pi" []
    exp = app1 "exp"
    sqrt = app1 "sqrt"
    log = app1 "log"
    (**) = app2 "**"
    logBase = app2 "logBase"
    sin = app1 "sin"
    tan = app1 "tan"
    cos = app1 "cos"
    asin = app1 "asin"
    atan = app1 "atan"
    acos = app1 "acos"
    sinh = app1 "sinh"
    tanh = app1 "tanh"
    cosh = app1 "cosh"
    asinh = app1 "asinh"
    atanh = app1 "atanh"
    acosh = app1 "acosh"

instance (RealFloat a) => RealFloat (Sym a) where
    exponent _ = 0
    scaleFloat 0 x = x
    atan2 = app2 "atan2"

app1 :: String -> Sym a -> Sym a
app1 f x = App f [x]

app2 :: String -> Sym a -> Sym a -> Sym a
app2 f x y = App f [x, y]

Let's put this code to the test by bringing the Complex number module into scope.

Sym3> :m +Data.Complex
Sym3 Data.Complex> let x=Var "x"; y=Var "y"
Sym3 Data.Complex> sin (x:+y)
sin x*cosh y :+ cos x*sinh y

And by that last expression we have recovered the definition of complex sin as it is given in the Data.Complex module. Let's try another one.

Sym3 Data.Complex> asinh(x:+y)
log (sqrt ((x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+
 (x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0))*(x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x
 -y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0))+(y+abs (0.0+
 (x*y+y*x))/(sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+
 (x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0)*2.0))*(y+abs (0.0+(x*y+y*x))/(sqrt ((sqrt
 ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x
 -y*y)))/2.0)*2.0)))) :+ atan2 (y+abs (0.0+(x*y+y*x))/(sqrt ((sqrt ((1.0+(x*x-y*y))
 *(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+abs (1.0+(x*x-y*y)))/2.0)*2.0))
 (x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(0.0+(x*y+y*x))*(0.0+(x*y+y*x)))+
 abs (1.0+(x*x-y*y)))/2.0))

Hmmmm, that might be right, but it's rather ugly. There's also a lot of '0.0+...' in that expression. We need something that can simplify expressions. It would also be nice if all constant expressions were evaluated instead of stored. To achieve this we are going to change the representation a little. The App constructor will store the real function used to work on constants as well as the name of it. And while we are at it, we'll get rid of the Var constructor. We might as well use App with an empty argument list. Furthermore, since this is starting to look useful, we'll give the module a proper name and export only the interface we want to be visible. We will hide the details of the Sym type and just export some accessor functions. The simplification happens in the binOp and unOp functions. I have just put some algebraic laws there (assuming the underlying numeric type is a field). The list of rewrites performed by these functions is far from complete. It's just a few that I found useful. Note how the code in binOp pattern matches on constants like 0, 1, and -1 directly. This actually works because of the semantics of Haskell pattern matching against numeric literals. Also note that the constraint on `a' is just Num, even though we do some simplifications with (/) which belongs in Fractional. The instance declarations have been extended somewhat so that constant expressions in the Sym type will behave as the corresponding expressions in the underlying type. A small final run

*Data.Number.Symbolic Data.Complex> 1+x+2
3+x
*Data.Number.Symbolic Data.Complex> 1+x*(y-y)-1
0
*Data.Number.Symbolic Data.Complex> sin(x:+1e-10)
sin x :+ 1.0e-10*cos x
*Data.Number.Symbolic Data.Complex> asinh(x:+y)
log (sqrt ((x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs
 (1.0+(x*x-y*y)))/2.0))*(x+sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*
 (x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))+(y+abs (x*y+y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x
 -y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0)))*(y+abs (x*y
 +y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs
 (1.0+(x*x-y*y)))/2.0))))) :+ atan2 (y+abs (x*y+y*x)/(2.0*sqrt ((sqrt ((1.0+(x*x
 -y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))) (x+sqrt
 ((sqrt ((1.0+(x*x-y*y))*(1.0+(x*x-y*y))+(x*y+y*x)*(x*y+y*x))+abs (1.0+(x*x-y*y)))/2.0))

As the final example shows, there is still a lot to do. Also note how the underlying numeric type has defaulted to Double, and we have a loss of precision in the second to last example. But an implementation of real numbers instead of floating point numbers will have to wait until a later posting.

module Data.Number.Symbolic(Sym, var, con, subst, unSym) where

import Data.Char(isAlpha)
import Data.Maybe(fromMaybe)
import Debug.Trace

data Sym a = Con a | App String ([a]->a) [Sym a]

instance (Eq a) => Eq (Sym a) where
    Con x      == Con x'        =  x == x'
    App f _ xs == App f' _ xs'  =  (f, xs) == (f', xs')
    _          == _             =  False

instance (Ord a) => Ord (Sym a) where
    Con x      `compare` Con x'        =  x `compare` x'
    Con _      `compare` App _ _ _     = LT
    App _ _ _  `compare` Con _         = GT
    App f _ xs `compare` App f' _ xs'  =  (f, xs) `compare` (f', xs')

var :: String -> Sym a
var s = App s undefined []

con :: a -> Sym a
con = Con

subst :: (Num a) => String -> Sym a -> Sym a -> Sym a
subst _ _ e@(Con _) = e
subst x v e@(App x' _ []) | x == x' = v
                         | otherwise = e
subst x v (App s f es) =
    case map (subst x v) es of
    [e] -> unOp (\ x -> f [x]) s e
    [e1,e2] -> binOp (\ x y -> f [x,y]) e1 s e2
    es' -> App s f es'

unSym :: (Show a) => Sym a -> a
unSym (Con c) = c
unSym e = error $ "unSym called: " ++ show e

instance (Show a) => Show (Sym a) where
    showsPrec p (Con c) = showsPrec p c
    showsPrec _ (App s _ []) = showString s
    showsPrec p (App op@(c:_) _ [x, y]) | not (isAlpha c) =
        showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y)
        where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [
                   ("**", (9,8,8)),
     ("/",  (7,7,8)),
     ("*",  (7,7,8)),
     ("+",  (6,6,7)),
     ("-",  (6,6,7))]
    showsPrec p (App "negate" _ [x]) =
        showParen (p>=6) (showString "-" . showsPrec 7 x)
    showsPrec p (App f _ xs) =
        showParen (p>10) (foldl (.) (showString f) (map (\ x -> showChar ' ' . showsPrec 11 x) xs))

instance (Num a) => Num (Sym a) where
    x + y         = binOp (+) x "+" y
    x - y         = binOp (-) x "-" y
    x * y         = binOp (*) x "*" y
    negate x      = unOp negate "negate" x
    abs    x      = unOp abs    "abs"    x
    signum x      = unOp signum "signum" x
    fromInteger x = Con (fromInteger x)

instance (Fractional a) => Fractional (Sym a) where
    x / y          = binOp (/) x "/" y
    fromRational x = Con (fromRational x)

-- Assume the numbers are a field and simplify a little
binOp :: (Num a) => (a->a->a) -> Sym a -> String -> Sym a -> Sym a
binOp f (Con x) _ (Con y) = Con (f x y)
binOp _ x "+" 0 = x
binOp _ 0 "+" x = x
binOp _ x "+" (App "+" _ [y, z]) = (x + y) + z
binOp _ x "+" y | isCon y && not (isCon x) = y + x
binOp _ x "+" (App "negate" _ [y]) = x - y
binOp _ x "-" 0 = x
binOp _ x "-" x' | x == x' = 0
binOp _ x "-" (Con y) | not (isCon x) = Con (-y) + x
binOp _ _ "*" 0 = 0
binOp _ x "*" 1 = x
binOp _ x "*" (-1) = -x
binOp _ 0 "*" _ = 0
binOp _ 1 "*" x = x
binOp _ (-1) "*" x = -x
binOp _ x "*" (App "*" _ [y, z]) = (x * y) * z
binOp _ x "*" y | isCon y && not (isCon x) = y * x
binOp _ x "*" (App "/" f [y, z]) = App "/" f [x*y, z]
{-
binOp _ x "*" (App "+" _ [y, z]) = x*y + x*z
binOp _ (App "+" _ [y, z]) "*" x = y*x + z*x
-}
binOp _ x "/" 1 = x
binOp _ x "/" (-1) = -x
binOp _ x "/" x' | x == x' = 1
binOp _ x "/" (App "/" f [y, z]) = App "/" f [x*z, y]
binOp f x op y = App op (\ [a,b] -> f a b) [x, y]

unOp :: (Num a) => (a->a) -> String -> Sym a -> Sym a
unOp f _ (Con c) = Con (f c)
unOp _ "negate" (App "negate" _ [x]) = x
unOp _ "abs" e@(App "abs" _ _) = e
unOp _ "signum" e@(App "signum" _ _) = e
unOp f op x = App op (\ [a] -> f a) [x]

isCon :: Sym a -> Bool
isCon (Con _) = True
isCon _ = False


instance (Real a) => Real (Sym a) where
    toRational (Con c) = toRational c

instance (RealFrac a) => RealFrac (Sym a) where
    properFraction (Con c) = (i, Con c') where (i, c') = properFraction c

instance (Floating a) => Floating (Sym a) where
    pi = var "pi"
    exp = unOp exp "exp"
    sqrt = unOp sqrt "sqrt"
    log = unOp log "log"
    x ** y = binOp (**) x "**" y
    logBase x y = binOp logBase x "logBase" y
    sin = unOp sin "sin"
    tan = unOp tan "tan"
    cos = unOp cos "cos"
    asin = unOp asin "asin"
    atan = unOp atan "atan"
    acos = unOp acos "acos"
    sinh = unOp sinh "sinh"
    tanh = unOp tanh "tanh"
    cosh = unOp cosh "cosh"
    asinh = unOp asinh "asinh"
    atanh = unOp atanh "atanh"
    acosh = unOp acosh "acosh"

instance (RealFloat a) => RealFloat (Sym a) where
    floatRadix = floatRadix . unSym
    floatDigits = floatDigits . unSym
    floatRange  = floatRange . unSym
    decodeFloat (Con c) = decodeFloat c
    encodeFloat m e = Con (encodeFloat m e)
    exponent (Con c) = exponent c
    exponent _ = 0
    significand (Con c) = Con (significand c)
    scaleFloat k (Con c) = Con (scaleFloat k c)
    scaleFloat _ x = x
    isNaN (Con c) = isNaN c
    isInfinite (Con c) = isInfinite c
    isDenormalized (Con c) = isDenormalized c
    isNegativeZero (Con c) = isNegativeZero c
    isIEEE = isIEEE . unSym
    atan2 x y = binOp atan2 x "atan2" y

Well, it was bound to happen. I have started a blog, just like everyone else.