Wednesday, April 11, 2007

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
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

Labels: , sigfpe said...

Nice!

I think Haskell would be a great language for implementing a full symbolic algebra system. I'd love to see a nice strongly typed algebra system like Aldor with the laziness of Haskell.

Friday, April 13, 2007 at 5:24:00 PM GMT+1 Jon Harrop said...

Very cool!

I'd like to see an expression type that statically conveyed the fact that it was simplified (to some extent).

Wednesday, August 15, 2007 at 5:15:00 AM GMT+1 Anonymous said...

It seems that Haskell needs `Ord` for `Real` typeclasses. So deriving `Ord` requires:

No explicit implementation for