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

## Saturday, April 28, 2007

## Friday, April 27, 2007

**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) epsSo 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.

## Wednesday, April 25, 2007

**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 -> RationalThe 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 _ = 1The

`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-50The 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 = undefinedThe 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 Prec50As 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 xIt 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 = undefinedAnd 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.00000000000000000000000000000000000000000000000000Naturally, 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/.

Labels: Haskell, overloading

## Friday, April 20, 2007

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 dThe Djinn source is at darcs.augustsson.net/Darcs/Djinn as usual. Thanks Cale!

Labels: Haskell

## Saturday, April 14, 2007

**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*and you want to know the derivative at

^{2}+ 1*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 xsSo 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*We want the derivative at

^{2}+ 1*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) 0A 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) 0We can now try the sqr example

Main> sqr (P 9 1) P 3.0 0.16666666666666666The 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) 0This 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-3And 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)) / 2And 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.4395397945007182It 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

_{1}*x*until we meet some convergence criterion. Using this, let's define a function that finds a zero of another function:

_{n+1}= x_{n}- f(x_{n})/f'(x_{n})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 xsThe 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.0Note 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.]

Labels: Haskell, overloading

## 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+3Look 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 yAnd 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

Labels: Haskell, overloading