# Dual numbers

26 Feb 2022By horsing around on YouTube these days I found
this video
about *dual numbers*,
which I found actually quite interesting. Since I have a computer engineering
background, I immediately asked myself what can these numbers used for. It turns
out, for some pretty nifty stuff.

### What dual numbers are

Essentialy, they are expressions of the form^{1}

$$ a + b\varepsilon \quad\quad a\in\R, b\in\R, \varepsilon^2 = 0 $$

Addition and multiplication are straightforward

$$ \begin{align*} (a + b\varepsilon) + (c + d\varepsilon) &= (a + c) + (b + d)\varepsilon \\ (a + b\varepsilon)(c + d\varepsilon) &= ac + (ad + bc)\varepsilon \end{align*} $$

If we define the conjugate of a dual number $z$ as $z^* = a - b\varepsilon$, we can define its absolute value as

$$ |z|^2 = zz^* = a^2 $$

And division is very similar to the complex numbers case

$$ \begin{align*} \frac{a+b\varepsilon}{c+d\varepsilon} &= \frac{a+b\varepsilon}{c+d\varepsilon}\frac{c-d\varepsilon}{c-d\varepsilon} \ &= \cdots \ &= \frac{a}{c} + \frac{bc - ad}{c^2}\varepsilon \end{align*} $$

### What are they useful for

There are various uses, but one example could be *automatic differentiation*.

As highlighted both in the previous video and on Wikipedia, given any real polynomial $ P(x) $ it's straightforward to extend its to domain to the dual numbers, and then you can obtain the following

$$ P(x + \varepsilon) = P(x) + \varepsilon P'(x) $$

Where $P'(x)$ is the first derivative of $P$ in $x$. Essentialy, you can in one pass compute both the value of a (polynomial) function and the value of its derivative in a specific point.

### Trying to put everything into Haskell

For this I've used the `GADTs`

language extension, so that I can use all numeric types to build the `Dual`

datatype.
The definition is pretty simple

```
{-# LANGUAGE GADTs #-}
data Dual a where
(:+) :: Num a => a -> a -> Dual a
eps :: Num a => Dual a
eps = 0 :+ 1 -- A little utility for later
```

Here it's necessary to define `Dual a`

and not only `Dual`

, otherwise we'll run into some problems when defining
instances. A downside of using a GADT is having to define manually (or standalone-derive) the `Show`

instance.

```
instance Show a => Show (Dual a) where
show (a :+ b) = show a <> " :+ " <> show b
```

So, dual numbers are *numbers* after all, so a `Num`

and a `Fractional`

instances won't hurt^{2}, and they are pretty much a translation
of the rules we seen before

```
instance Num a => Num (Dual a) where
(a :+ b) + (c :+ d) = (a + c) :+ (b + d)
(a :+ b) * (c :+ d) = (a * c) :+ (a*d + b*c)
fromInteger n = fromInteger n :+ 0
negate x = fromInteger (-1) * x
abs (a :+ _) = abs a :+ 0
signum = undefined
instance Fractional a => Fractional (Dual a) where
fromRational n = fromRational n :+ 0
(a :+ b) / (c :+ d) = (a / c) :+ ((b*c - a*d) / c ^ 2)
```

Finally, we can test if it all works by using a simple polynomial function and its derivative

$$ \begin{aligned} f(x) &= 3x^2 - 4x \ f'(x) &= 6x - 4 \end{aligned} $$

Let's try some random inputs

$$x$$ | $$f(x)$$ | $$f'(x)$$ |
---|---|---|

$$0$$ | $$0$$ | $$-4$$ |

$$1$$ | $$-1$$ | $$2$$ |

$$-1$$ | $$7$$ | $$-10$$ |

$$7$$ | $$119$$ | $$38$$ |

$$3$$ | $$15$$ | $$14$$ |

$$-2$$ | $$20$$ | $$-16$$ |

Finally, we can define all the necessary functions in Haskell^{3}

```
xs = [0, 1, -1, 7, 3, -2]
diff :: (Num a, Num b) => (Dual a -> b) -> Dual a -> b
diff f x = f (fromInteger x + eps)
xs' = map (diff f) xs
```

Now we get

```
ghci> mapM_ print xs'
0 :+ -4
-1 :+ 2
7 :+ -10
119 :+ 38
15 :+ 14
20 :+ -16
```

Which are all correct! We can see that the "real" part of the numbers is the value of $f$ in the specific point, while the "epsilon" part is the value of $f'$ in the same point.

I think it's really nice how these things arise from well defined and surprisingly simple rules and definitions.

The full code:

^{1}

Notice the similarity with complex numbers

^{2}

I haven't found if it's reasonable to talk about a `signum`

of a dual
number, so for now I'll leave it undefined, maybe I'll find something in the
future

^{3}

I used `Integer`

here, but it's totally possible to use any other numeric type