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