Constraining Types with Regular Expressions

Structures with Constraints
Here's a picture of an element of a binary tree type:

The leaves correspond to elements and the letters indicate the type of those elements. If we read the leaves from left to right they match the regular expression A^*B^*. We can define a binary type whose leaves always match this string by making use of dissections. Similarly we can construct all kinds of other structures, such as lists and more general trees, whose leaves match this expression.

But can we make structures that match other regular expressions? Here are some examples:

1. Structures that match the expression Aⁿ for some fixed n. These are structures, like trees, that have a fixed size. For the case n=5 we need a type that includes elements like:

2. Structures that match (AA)^*. These are structures like trees with an even number of elements. This is easy for lists. We can just use a list of pairs. But it's more complex for trees because subtrees may have odd or even size even though the overall structure is even in size. We could also generalise to ((AA)ⁿ)^* for some fixed n.

3. Structures that match A^*1A^*. Here's an example:

The element marked 1 is of unit type (). It's a 'hole'. So this regular expression is a derivative type. We can also construct a variety of types like those matching A^*BA^*CB^* by using a mixture of dissection and differentiation:

4. A bytestring- or rope-like type used to represent a string that is statically guaranteed to match a given regular expression.

5. Many kinds of other constraints you could imagine on a datastructure. Like trees which are guaranteed not to have two neighbouring leaves of the same type, or whose sequence of leaves never contain a certain subsequence.

So the challenge is this: can we implement a uniform way of building any container type to match any regular expression we want?

We'll need some extensions:

> {-# OPTIONS_GHC -fwarn-incomplete-patterns #-}
> {-# LANGUAGE TypeFamilies, EmptyDataDecls, UndecidableInstances, TypeOperators #-}

First Example: Differentiating Lists
Let l(x) be the type of lists of elements of type x. We can define a list type through

data List x = Nil | Cons x (List x)

Algebraically we can write this as

l(x) = 1+x l(x)

Let's try approaching the derivative of a list from the perspective of regular expressions. We know from Kleene's Theorem that the set of strings (ie. the language) matching a regular expression is precisely the language accepted by a finite state automaton. Let's consider the case x^*1x^*. This is the language we get by considering all paths from state 0 to state 1 in the following automaton:

Note that I have overloaded the numeral '1' to mean both the type with one element (when labelling an edge) and to mean the state numbered 1 (when labelling a vertex).

Define the L_ij(x) to be the type of lists whose sequence of elements correspond to all possible paths from state i to state j. A list in L_ij(x) is either an empty list, or constructed from a first element and a list. If we combine an element and a list, the combination has to match a possible path from i to j. There are a number of ways we could do this. The element could correspond to a transition from i to k. But if this is the case, then the remainder of the list must correspond to a path from k to j. So we must replace Cons with something that for its first argument takes a type corresponding to a single automaton step from i to k. For its second argument it must take an element of L_kj(x). The set of paths from i to j that correspond to a single step is the (i,j)-th element of the transition matrix for the automaton. This is X_ij(x) where

X(x) = (Xij(x)) =	(x	1)
	(0	x)

On the other hand, we also need to replace Nil with a version that respects transitions too. As Nil takes no arguments, it must correspond to paths of length zero in the automaton. The only such paths are zero length paths from a state to itself. So the matrix for such paths is:

I =	(1	0)
	(0	1)

Let's also define the matrix L(x)= (L_ij(x)).

The words above boil down to:

L_ij(x) = I_ij+Σ_k X_ik(x)L_kj(x)

where the sum is over all the places we might visit on the first step of the journey from i to j.

We can rewrite this using standard matrix notation:

L(x) = I + X(x)L(x)

Compare with the definition of ordinary lists given above. We get the type of constrained lists by taking the original definition of a list and replacing everything with matrices. We replace 1 with I. We replace x with the transition matrix of the automaton. And we replace the structure we're trying to define with a family of types - one for each pair of start and end states for the automaton. We can describe this replacement more formally: it's a homomorphism from the set of types to the set of matrices of types. (Actually, it's a bit more subtle than that. This isn't quite the usual semiring of types. For one thing, the order of multiplication matters.) And it doesn't just apply to lists. We can apply this rule to any container type. For example, suppose we wish to repeat the above for trees. Then we know that for ordinary binary trees, t(x), we have

t(x) = x+t(x)²

We replace this with

T(x) = X(x)+T(x)²

We started this section by considering the specific case of the pattern x^*1x^* with a corresponding matrix X(x). Because X₁₀=0 and X₀₀=X₁₁ it's not hard to see that any type T we constrain using this regular expression will also have similar 'shape', ie. T₁₀=0 and T₀₀=T₁₁. So we can write

T =	(t(x)	t'(x))
	(0	t(x))

where by definition, t(x)=T₀₀(x) and t'(x)=T₀₁(x). Suppose we have two such collections of types, (S_ij) and (T_ij). Now consider the types of pairs where the first element is of type S_ij and the second of T_jk. Then the leaves of the pair structure correspond to a path from i to k. So we have

(st(x)	(st)'(x))	=	(s(x)	s'(x))	(t(x)	t'(x))
(0	st(x))		(0	s(x))	(0	t(x))

Multiply out and we find that

(st)'(x) = s(x)t'(x)+s'(x)t(x)

In other words, the usual Leibniz rule for differentiation is nothing more than a statement about transitions for the automaton I drew above. To get a transition 0→1 you either go 0→0→1 or 0→1→1.

Although I talked specifically about differentiation, much of what I said above applies for any finite state automaton whose edges are labelled by types. The best thing now is probably to put together some code to see how this all looks.

A Specification
If you haven't checked Brent Yorgey's solution to my problem last week, now is a good opportunity. My code is a generalisation of that but it may be helpful to look at Brent's specialisation first.

The goal is to be able to define a transition matrix like this. ('K' is an abbreviation for 'constant' matrix.)

> type D a = K22 (a,    ())
>                (Void, a )

And then define a functor like this:

> type ListF = I :+: (X :*: Y)

Think of ListF being a bifunctor taking arguments X and Y.

We'd like then to able to form the matrix of fixed points of Y = ListF X Y. In this case, ordinary lists should appear as the element at position (0,0) in the matrix:

> type List x = Fix (I0, I0) (D x) ListF

I'm using In to represent the integer n at the type level.

Derivatives of lists should appear at (0,1) so we want

> type List' x = Fix (I0, I1) (D x) ListF

But Fix is intended to be completely generic. So it needs to be defined in a way that also works for trees:

> type TreeF = X :+: (Y :*: Y)
> type Tree  = Fix (I0, I0) (D Int) TreeF
> type Tree' = Fix (I0, I1) (D Int) TreeF

And of course it needs to work with other transition matrices. For example x^*1y^*1z^* has the following transition diagram and matrix:

> type E x y z = K33 (x,    (),   Void)
>                    (Void, y,    ()  )
>                    (Void, Void, z   )

So we'd expect

> type DDTree x y z = Fix (I0, I1) (E x y z) TreeF

to define the second divided difference of trees.

The Implementation

We'll need a type with no elements:

> data Void

And some type level integers:

> data Zero
> data S a
> type I0 = Zero
> type I1 = S I0
> type I2 = S I1
> type I3 = S I2

Now we'll need some type-level matrices. For any square matrix, we need a type-level function to give its dimension and another to access its (i, j)-th element:

> type family Dim x :: *
> type family (:!) x ij :: *

(Thanks to Brent's tutorial that code is much better than how it used to look.)

We can now define matrix addition through pointwise addition of the elements:

> data (:+) m n
> type instance Dim (m :+ n) = Dim m
> type instance (m :+ n) :! ij = Either (m :! ij) (n :! ij)

And similarly we can define multiplication. I'm using the type-level function Product' to perform the loop required in the definition of matrix multiplication:

> data (:*) m n
> type instance Dim (m :* n) = Dim m
> type instance (m :* n) :! ij = Product' I0 (Dim m) m n :! ij

> data Product' i k m n
> type instance Product' p I1 m n :! (i, j) = (m :! (i, p), n :! (p, j))
> type instance Product' p (S (S c)) m n :! (i, j) = Either
>    (m :! (i, p), n :! (p, j))
>    (Product' (S p) (S c) m n :! (i, j))
> 

(Weird seeing all that familiar matrix multiplication code at the type level.)

Now we need some types to represent our functors:

> data I
> data X
> data Y
> data K n
> data (f :+: g)
> data (f :*: g)
> data F m f

(I think phantom empty data types should be called ethereal types.)

To turn these into usable types we need to implement the homomorphism I described above. So here are the rules laid out formally:

> type Id = K22 ((),   Void)
>               (Void, ())

> type family   Hom self m f :: *
> type instance Hom self m I = Id
> type instance Hom self m X = m
> type instance Hom self m Y = self
> type instance Hom self m (K n) = n
> type instance Hom self m (f :+: g) = Hom self m f :+ Hom self m g
> type instance Hom self m (f :*: g) = Hom self m f :* Hom self m g
> type instance Dim (F m f) = Dim m

> data Fix ij m f = Fix (Hom (F m f) m f :! ij)
> type instance (:!) (F m f) ij = Fix ij m f

That's more or less it. We can now go ahead and try to construct some elements. We could (as Brent suggests) write some smart constructors to make our life easier. But for now I'm writing everything explicitly so you can see what's going on:

> x0 = Fix (Left ())                         :: List Int
> x1 = Fix (Right (Left (1, Fix (Left ())))) :: List Int

x0 is the empty list. x1 is the list [1]. The Left and Right get a bit tedious to write. But this is intended as a proof that the concept works rather than a user-friendly API.

We can explicitly implement the isomorphism with the more familiar list type:

> iso1 :: [x] -> List x
> iso1 []     = Fix (Left ())
> iso1 (a:as) = Fix (Right (Left (a, iso1 as)))

> iso1' :: List x -> [x]
> iso1' (Fix (Left ()))              = []
> iso1' (Fix (Right (Left (a, as)))) = a : iso1' as
> iso1' (Fix (Right (Right a)))      = error "Can't be called as a is void"

So that's it! If we can write our container as the fixed point of a polynmomial functor, and if we can convert our regular expression to a finite state automaton, then Fix completely automatically builds the constrained container type.

What have we learnt?

We haven't just solved the original problem. We've shown that derivatives and dissections are special cases of a more general operation. Take a look at the definition of D x again. we can think of it as xI+Delta where Delta is the matrix

> type Delta = K22 (Void, ()  )
>                  (Void, Void)

This matrix has the property that its square is zero. It's the 'infinitesimal type' I described here. In other words, this is type-level automatic differentiation. We've also been doing type-level automatic divided differencing.

We can now go back and look at the matrix form of divided differences on wikipedia. I hope you can now see that the matrix T_id(x₀,...,x_n-1) defined there is nothing other than a transition matrix for this automaton:

In fact, we can use what we've learnt about regular expressions here to solve some numerical problems. But I won't write about that until the next article.

By the way, I think what I've described here can be viewed as an application of what Backhouse talks about in these slides.

I think that for any automaton we have a 2-category. The 0-cells are states, the 1-cells are the types associated with paths from one state to another, and the 2-cells are functions between types that respect the constraint. I haven't worked out the details however. The 2-category structure is probably important. As things stand, I've just shown how to make the types. But we don't yet have an easy way to write functions that respect these constraints. I suspect 2-categories give a language to talk about these things. But that's just speculation right now.

By the way, I couldn't write a working Show instance for Fix. Can you write one? And an implementation of arbitrary for QuickCheck?

And I hope you can now solve my problem from last week.

Leftover bits of code

K22 and K33 are constructors for 2×2 and 3×3 matrices. It would probably have been better to have used lists like Brent did.

> data K22 row0 row1
> type instance Dim (K22 row0 row1) = I2
> 
> type instance (:!) (K22 (m00, m01) row1) (I0, I0) = m00
> type instance (:!) (K22 (m00, m01) row1) (I0, I1) = m01
> type instance (:!) (K22 row0 (m10, m11)) (I1, I0) = m10
> type instance (:!) (K22 row0 (m10, m11)) (I1, I1) = m11
> 
> data K33 row0 row1 row2
> type instance Dim (K33 row0 row1 row2) = I3
> 
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I0) = m00
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I1) = m01
> type instance (:!) (K33 (m00, m01, m02) row1 row2) (I0, I2) = m02
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I0) = m10
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I1) = m11
> type instance (:!) (K33 row0 (m10, m11, m12) row2) (I1, I2) = m12
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I0) = m20
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I1) = m21
> type instance (:!) (K33 row0 row1 (m20, m21, m22)) (I2, I2) = m22

Update: I neglected to mention that there is a bit of subtlety with the issue of being able to create the same string by different walks through the automaton. I'll leave that as an exercise :-)

Divided Differences and the Tomography of Types

Health Warning
This article assumes knowledge of Conor McBride's work on the differentiation and dissection of types, even if the early parts look deceptively simple.

Tables of Divided Differences
Given a sequence of numbers like 2, 6, 12, 20, 30 we can use the well known method of finite differences to predict the next in the series. We write the numbers in a column. In the next column we write the differences between the numbers in the previous column and iterate like so:

2
	4
6		2
	6		0
12		2
	8		0
20		2
	10
30

On the assumption that the rightmost column is zero all the way down we can extend this table to:

2
	4
6		2
	6		0
12		2
	8		0
20		2
	10		0
30		2
	12
42

We can think of this as an exercise in polynomial fitting. We have some x-values, in this case: x₀=1, x₂=2,..., x₄=5, and some y-values: y₀=2, y₁=6,..., y₄=30. We hope to fit a polynomial f so that f(x_i)=y_i. The table of finite differences has the property that the (n+1)-th column becomes constant if a degree n polynomial fits the data.

But what happens if we want to fit data to x-values that aren't equally spaced? Then we can use the *divided* difference table instead. In this case we don't record the differences between the y-values, but the quotient between the y-differences and the x-differences. For x_0,1,2,3 = 1, 2, 4, 6 and y_0,1,2,3 = 3, 6, 18, 138 we get:

1	3
		(6-3)/(2-1)=3
2	6		(6-3)/(4-1)=1
		(18-6)/(4-2)=6		(1-1)/(6-1)=0
4	18		(10-6)/(6-2)=1
		(38-18)/(6-4)=10
6	38

Again we reach zero. This is because y_i = f(x_i) = x_i²+2 so we have a quadratic again.

Let's assume y_i = f(x_i) for some set of points. We'll leave x, y and f unknown and fill out the table:

x0	f(x0)
		(f(x1)-f(x0))/(x1-x0)= f[x0,x1]
x1	f(x1)		(f[x1,x2]-f[x0,x1])/(x2-x0)= f[x0,x1,x2]
		(f(x2)-f(x2))/(x2-x1)= f[x1,x2]		(f[x1,x2,x3]-f[x0,x1,x2])/(x3-x0)= f[x0,x1,x2,x3]
x2	f(x2)		(f[x2,x3]-f[x1,x2])/(x3-x1)= f[x1,x2,x3]
		(f(x3)-f(x2))/(x3-x2)= f[x2,x3]
x3	f(x3)

Note that for any f, this table also defines the generalised divided differences f[x,y], f[x,y,z] and so on. (Compare with the notation at Wikipedia).

Divided Differences of Types
Now suppose that f isn't a function of real numbers, but is a container type. Then as I noted here, the second column corresponds to Conor McBride's dissection types. f[x,y] is an f-container with a hole, with everything to the left of the hole containing an x, and everything to the right of the hole containing a y.

So what's the next column? Consider f(x)=x⁴. Then f[x,y] = (y⁴-x⁴)/(y-x) = x³+x²y+xy²+y³. This corresponds to the description I just gave involving a hole. Now consider f[x,y,z] = (f[y,z]-f[x,y])/(z-x) = z²+xz+x²+y²+yz+yx. This is the type consisting of a 4-tuple, with two holes, and with x left of the first hole, y between the first and second holes, and z to the right of the second hole. In other words, it's a trisection of a container. More generally, f[x₀,...,x_n] is the type corresponding to n holes and blocks of elements of x_i between them. I've only shown this for the type x⁴ but it works for all the same types in Conor's paper. But there's a catch: I've used a definition involving subtraction, and that makes no sense for types. Don't worry, I'll address that later.

We can give trisections and so on a computational interpretation like in Conor's paper. Dissections correspond to tail recursion elimination. They make explicit the state of a function traversing and transforming a recursive type. (I recommend Conor's description in the paper so I won't reproduce it here.) Trisections correspond to a pair of coroutines. The first is transforming a recursive type. The second is transformed the result of the first routine. The second one can (optionally) start consuming pieces of tree as soon as the first has started producing them. At any moment in time we have 3 pieces: (1) the part of tree that is as yet untouched, (2) the part that the first coroutine has produced, but which the second hasn't seen yet, and (3) the output from the second coroutine. Make that explicit, and you get a trisection. The deal's much the same for quadrisections and so on.

Rediscovering Newton
A consequence of this is that we can now give a computational interpretation to much of the description of divided differences at Wikipedia. Among other things, when Conor derives a type theoretical analogue of the Taylor series, he's actually rediscovering a form of the Newton polynomial.

Type Matrices
One interesting property of divided differences is that they have a description in terms of matrices. In particular, Wikipedia describes a homomorphism from the ring of functions of x to the ring of matrices whose elements are functions of n unknowns. We know that types don't form a ring, but they do form a semiring. So we can form an analogous semiring homomorphism from the set of container types to the set of matrices whose entries are *types*. Matrices of types aren't something we see every day. How can we put this notion to work?

Let's use the notation from Wikipedia and consider the matrices we get when n = 1. (At this point I recommend reading all of the Wikipedia article on divided differences.) Let's define ourselves a tree container:

> data F a = Leaf a | Form (F a) (F a)

We'll call the homomorphism T so that

T(f)(x0,x1)=	f(x0)	f[x0,x1]
	0	f(x1)

Our tree type satisfies the equation f(x)=x+f(x)². As we have a homomorphism, we also expect this to hold:

T(f)(x₀,x₁) = T(i)(x₀,x₁)+T(f)(x₀,x₁)² (Equation 1)

where i is the identity function i(x) = x.

We can easily compute T(i)(x₀,x₁) directly using (y-x)/(y-x)=1. We get

T(i)(x0,x1)=	x0	1
	0	x1

Multiplying out the matrices in equation 1 now gives us a bunch of equations:

f(x₀) = x₀ + f(x₀)²
f(x₁) = x₁ + f(x₁)²
f[x₀,x₁] = 1 + f(x₀)f[x₀,x₁]+f[x₀,x₁]f(x₁)

The first two equations are just the definition of our tree. The third line, in Haskell, is:

> data F' x0 x1 = Empty | ForkL (F x0) (F' x0 x1) | ForkR (F' x0 x1) (F x1)

This is essentially the application of Conor's version of the Leibniz law to trees. F' is the dissected tree. And note how we don't have any subtraction. By using the matrix formulation of divided differences we only have matrix multiplications and sums to deal with.

The image of f(x) = x+f(x)² under the homomorphism T yields the simultaneous definitions of f and its dissection. More generally, if we had chosen to work with larger n we'd get the simultaneous definition of trees, their dissections, their trisections and so on. And it'll work for any recursive container type.

So matrices of types are a meaningful concept. They give a way to organise the mutually recursive definitions of higher order divided differences. If you look at this hard enough you may also see that the matrix I defined above is playing the same role as the D type in my previous article. With a little template Haskell I think we could in fact implement automatic (higher order) divided differences at the type level.

But this is all just scratching the surface. The matrix and homomorphism defined above don't just apply to divided differences. Matrices of types have another deeper and surprising interpretation that will allow me to unify just about everything I've ever said on automatic differentiation, divided differences, and derivatives of types as well as solve a wide class of problems relating to building data types with certain constraints on them. I'll leave that for my next article.

Exercise
In preparation for the next installment, here's a problem to think about: consider the tree type above. We can easily build trees whose elements are of type A or of type B. We just need f(A+B). We can scan this tree from left to right building a list of elements of type A+B, ie. whose types are each either A or B. How can we redefine the tree so that the compiler enforces the constraint that at no point in the list, the types of four elements in a row spell the word BABA? Start with a simpler problem, like enforcing the constraint that AA never appears.

Apology
Sorry about the layout above. My literate Haskell->HTML program doesn't support tables yet so there was a lot of manual HTML. This meant that I didn't write stuff out as fully as I could have and it may be a bit sketchy in places. I may have to switch to PDF for my next post. (Or use Wordpress, but I had trouble there too. or I could use one of the Javascript TeX renderers, but I don't like the external dependency.)

Automatic Divided Differences

Divided Differences
I've previously talked about automatic differentiation here a few times. One of the standard arguments for using automatic differentiation is that it is more accurate than numeric differentiation implemented via divided differences. We can approximate f'(x) by using (f(x)-f(y))/(x-y) with a value of y near x. Accuracy requires y to be close to x, and that requires computing the difference between two numbers that are very close. But subtracting close numbers is itself a source of numerical error when working with finite precision. So you're doomed to error no matter how close you choose x and y to be.

However, the accuracy problem with computing divided differences can itself be fixed. In fact, we can adapt the methods behind automatic differentiation to work with divided differences too.

(This paragraph can be skipped. I just want to draw a parallel with what I said here. Firstly I need to correct the title of that article. I should have said it was about *divided differences*, not *finite differences*. The idea in that article was that the notion of a divided difference makes sense for types because for a large class of function you can define divided differences without using either differencing or division. You just need addition and multiplication. That's the same technique I'll be using here. I think it's neat to see the same trick being used in entirely different contexts.)

The Direct Approach
Firstly, here's a first attempt at divided differencing:

> diff0 f x y = (f x - f y)/(x - y)

We can try it on the function f:

> f x = (3*x+1/x)/(x-2/x)

diff0 f 1 1.000001 gives -14.0000350000029. Repeating the calculation with an arbitrary precision package (I used CReal) gives -14.000035000084000. We are getting nowhere near the precision we'd like when working with double precision floating point.

The Indirect Approach
Automatic differentiation used a bunch of properties of differentiation: linearity, the product rule and the chain rule. Similar rules hold for divided differences. First let me introduce some notation. If f is a function then I'll use f(x) for normal function application. But I'll use f[x,y] to mean the divided difference (f(x)-f(y))/(x-y). We have

(f+g)[x,y] = f[x,y]+g[x,y]
(fg)[x,y] = f(x)g[x,y]+f[x,y]g(y)
h[x,y] = f[g(x),g(y)]g[x,y] when h(x)=f(g(x))

We can modify the product rule to make it more symmetrical though it's not strictly necessary:

(fg)[x,y] = 0.5(f(x)+f(y))g[x,y]+0.5f[x,y] (g(x)+g(y))

(I got that from this paper by Kahan.)

In each case, given f evaluated at x and y, and its divided difference at [x, y], and the same for g, we can compute the corresponding quantities for the sum and product of f and g. So we can store f(x), f(y) and f[x,y] together in a single structure:

> data D a = D { fx :: a, fy :: a, fxy :: a } deriving (Eq, Show, Ord)

And now we can implement arithmetic on these structures using the rules above:

> instance Fractional a => Num (D a) where
>    fromInteger n = let m = fromInteger n in D m m 0
>    D fx fy fxy + D gx gy gxy = D (fx+gx) (fy+gy) (fxy+gxy)
>    D fx fy fxy * D gx gy gxy = D (fx*gx) (fy*gy) (0.5*(fxy*(gx+gy) + (fx+fy)*gxy))
>    negate (D fx fy fxy) = D (negate fx) (negate fy) (negate fxy)

I'll leave as an exercise the proof that this formula for division works:

> instance Fractional a => Fractional (D a) where
>    fromRational n = let m = fromRational n in D m m 0
>    D fx fy fxy / D gx gy gxy = D (fx/gx) (fy/gy) (0.5*(fxy*(gx+gy) - (fx+fy)*gxy)/(gx*gy))

For the identity function, i, we have i(x)=x, i(y)=y and i[x,y]=1. So for any x and y, the evaluation of the identity function at x, y and [x,y] is represented as D x y 1. To compute divided differences for any function f making use of addition, subtraction and division we need to simply apply f to D x y 1. We pick off the divided difference from the fxy element of the structure. Here's our replacement for diff0.

> diff1 f x y = fxy $ f (D x y 1)

This is all mimicking the construction for automatic differentiation.

Evaluating diff0 f 1 1.000001 gives -14.000035000083997. Much closer to the result derived using CReal. One neat thing about this is that we have a function that's well defined even in the limit as x tends to y. When we evaluate diff1 f 1 1 we get the derivative of f at 1.

I thought that this was a novel approach but I found it sketched at the end of this paper by Reps and Rall. (Though their sketch is a bit vague so it's not entirely clear what they intend.)

Both the Kahan paper and the Reps and Rall papers give some applications of computing divided diferences this way.

It's not clear how to deal with the standard transcendental functions. They have divided differences that are very complex compared to their derivatives.

Aside
There is a sense in which divided differences are uncomputable(!) and that what we've had to do is switch from an extensional description of functions to an intensional description to compute them. I'll write about this some day.

Note that the ideas here can be extended to higher order divided differences and that there are some really nice connections with type theory. I'll try to write about these too.

Update: I found another paper by Reps and Rall that uses precisely the method described here.

Death to Hydrae (or the operational semantics of ordinals)

Unprovable Propositions
Among other things, Godel's first incompleteness theorem allows us to construct a statement in the language of Peano arithmetic that can't be proved using the axioms of Peano arithmetic. Unfortunately, this statement is a highly contrived proposition whose sole purpose is to be unprovable. People who learn of Godel's theorems often ask if there are other more natural and uncontrived mathematical statements that can't be proved from the Peano axioms.

My goal in this post will be to describe one of these propositions. Not just uncontrived, but actually very useful. I only intend to tell half of the story here because I feel like there are many good treatments already out there that tell the rest. I'm just going to get to the point where I can state the unprovable proposition, and then sketch how it can be proved if you allow yourself a little Set Theory.

> {-# OPTIONS_GHC -fno-warn-missing-methods #-}
> import Prelude hiding ((^))
> infixr 8 ^
> type Natural = Integer

Termination
Suppose we implement a function to compute the Fibonacci numbers like so:

> fib 0 = 0
> fib 1 = 1
> fib n = fib (n-2) + fib (n-1)

How do we know that fib terminates for all natural number arguments? One approach is this: if we pass in the argument n it clearly never recurses more than n levels. Each time it recurses it calls itself at most twice. So it must terminate in O(2ⁿ) steps (assuming that the primitive operations such as addition take constant time). We can think of this code in a kind of imperative way. It's a bit like n nested loops, each loop going round up to two times.

Suppose instead that we have some kind of recursive function g that goes n levels deep but for which the number of calls of g to itself is no longer two. In fact, suppose the number of self-calls is very large. Even worse, suppose that each time g is called, it calls itself many more times than it did previously, maybe keeping track of this ever growing number through a global variable. Or instead of a global variable, maybe an evil demon decides how many times g calls itself at each stage. Can you still be sure of termination?

A Simple Machine
In order to look at this question, we'll strip a computer right down to the bare minimum. It will have an input (that the evil demon could use) for natural numbers and will output only one symbol. Here's a design for such a machine:

> data Machine = Done | Output Machine | Input (Natural -> Machine)

A value of type Machine represents the state of the machine. Done means it has finished running. Output s means output a symbol and continue in state s. Input f means stop to input a number from the demon (or elsewhere), call it i, and then continue from state f i. This is very much in the style discussed by apfelmus and I in recent blog posts.

Here's an interpreter for one of these machines:

> run1 Done       = return ()
> run1 (Output x) = print "*" >> run1 x
> run1 (Input f)  = readLn >>= (run1 . f)

For any n we can easily build a machine to output n stars. This is such a natural machine to want to build it seems only right to give it the name n. If we want to do this then we need to make Machine an instance of Num and define fromInteger for it:

> instance Num Machine where
>    fromInteger 0 = Done
>    fromInteger n = Output (fromInteger (n-1))

Typing run1 8, say, will output 8 stars.

Now given two of these machines there is a natural notion of adding them. a + b is the machine that does everything b does followed by everything a does. (Remember, that's b then a.) To do this we need to dig into b and replace every occurrence of Done in it with a. That way, instead of finishing like b, it leads directly into a. In the case of a + Input f, for each number i we need to dig into f i replacing each Done with a:

>    a + Done     = a
>    a + Output b = Output (a + b)
>    a + Input f  = Input (\i -> a + f i)

There's a natural way to multiply these machines too. The idea is that in a * b we run machine b. But each time the Output command is run, instead of printing a star it executes a. You can think of this as a control structure. If n is a natural number then a * n means running machine a n times. In the case of a * Input f, instead of multiplying by a fixed natural number, we get an input from the user and multiply by f i instead:

>    _ * Done     = Done
>    a * Output b = a*b + a
>    a * Input f  = Input (\i -> a * f i)

We can make a machine to input a number and then output that many stars. Here it is:

> w = Input fromInteger

Try running run1 w.

Can you guess what the machine w * w does? Your first guess might be that it inputs two numbers and outputs as many stars as the product of the two numbers. Try it. What actually happens is that we're computing w * Input fromInteger. Immediately from the definition of * we get Input (\i -> w*i). In other words, the first input gives us an input i, and then w is run i times. So if we initially input i, we are then asked for i more inputs and after each input, the corresponding number of stars is output. Although the original expression contains just two occurrences of w, we are required to enter i+1 numbers.

Given the definitions of + and * it seems natural to define the power operation too:

> (^) :: Machine -> Machine -> Machine
> a ^ Done     = Output Done
> a ^ Output b = a^b * a
> a ^ Input f  = Input (\i -> a ^ f i)

The power operation corresponds to the nesting of loops. So, for example, w ^ n can be thought of loops nested n deep.

Try working out what w ^ w does when executed with run1.

Consider the set M of all machines built using just a finite number of applications of the three operators +, * and ^ to w and the non-zero naturals. (The non-zero condition means we exclude machines like 0*w that accept an input and do nothing with it.) Any such expression can be written as f w, where the definition of f makes no mention of w.

Suppose we use run1 and we always enter the same natural n. Then each occurrence of w acts like n. So if we start with some expression in w, say f w, then always inputting n results in f n stars. We could test this with run1 (w^w^w^w), always entering 2, but it would require a lot of typing. Instead we can write another intepreter that consumes its inputs from a list rather from the user (or demon). And instead of printing stars it simply prints out the total number of stars at the end:

> run2 Done    _        = 0
> run2 (Output x)   as  = 1 + run2 x as
> run2 (Input f) (a:as) = run2 (f a) as

Now you can try run2 (w^w^w^w) [2,2..] and see that we (eventually) get 2²²².

Termination Again
If we run a machine in M there's a pattern that occurs again and again. We input a number, and then as a result we go into a loop requesting more numbers. These inputs may in turn request more inputs. Like the mythological hydra, every input we give may spawn many more requests for inputs. As the number of inputs required may depend on our previous inputs, and we may input numbers as large as we like, these machines may run for a long time. Suppose our machine terminates after requesting n inputs. Then there must be some highest number that we entered. Call it m. Then if the original machine was f w (with f defined in terms of the 3 operators and non-zero naturals), the machine must have terminated outputting no more than f m stars. So if our machine terminates, we can bound how many steps it took.

But do our machines always terminate? The input we give to the machine might not be bounded. If we run run2 (w^w) [4,5..], say, the inputs grow and grow. If these inputs grow faster than we can chop off the heads of our hydra, we might never reach termination.

Consider a program to input n and then output fib n. It accepts an input, recurses to a depth of at most n, and calls itself at most twice in each recursion. Compare with the machine 2 ^ w. This accepts an input n, recurses to a depth n, calling itself exactly twice each time. So if 2 ^ w terminates, so does fib. The more complex example above where I introduced the evil demon will terminate if w ^ w does, as long as the demon doesn't stop inputting numbers. So if we can show in one proof that every machine of type Machine terminates, then there are many programs whose termination we could easily prove.

Let's consider an example like run1 (w ^ w) with inputs 2, 3, 4, ...

We start with w ^ w. Examining the definition of the operator ^ we see that this proceeds by requesting an input. The first input is 2. Now we're left with w ^ 2. This is w * w. Again it accepts an input. This time 3. Now we go to state w * 3. This is w*2 + w. Again we accept an input. This time 4. We are led to w*2 + 4. This now outputs 4 stars and we are left with w * 2 which is w + w. We accept an input 5, output 5 stars and are left with w. After a further input of 6, it outputs 6 stars and terminates. Or we could just run run2 (w ^ w) [2,3..] and get 15(=4+5+6) as output.

The transitions are:
w ^ w
-> w ^ 2
-> w * w
-> w * 3
-> w*2 + w
-> w*2 + 4 -> ... -> w * 2
-> w + w
-> w + 5 -> ... w
-> 6 -> ... -> 0.

Now for some Set Theory. Rewrite the above sequence using the transfinite ordinal ω instead of w. The sequence becomes a sequence of ordinals. Any time we accept an input, the rightmost ω becomes a finite ordinal. So we have a descending sequence of ordinals. This is true whatever ordinal we start with. The execution of either Input or Output always strictly decreases our ordinal, and any descending sequence of ordinals must eventually terminate. Therefore every machine in M eventually terminates.

But here's the important fact: to show termination we used the ordinal ω, and this required the axiom of infinity and some Set Theory. Instead we could encode the termination question, via Godel numbering, as a proposition of Peano arithmetic. If we do this, then we hit against an amazing fact. It can't be proved using the axioms of Peano arithmetic. So we have here a useful fact, not a contrived self-referential one, that can't be proved with Peano arithmetic.

Why can't it be proved using just the Peano axioms?
A few years back, Jim Apple made a post about constructing (some) countable ordinals in Haskell. His construction nicely reflects the definitions a set theorist might make, but the code doesn't actually do anything. Later I learned from Hyland and Power how you can interpret algebraic structures as computational effects. apfelmus illustrates nicely how an abstract datatype can be made to do things with the help of an interpreter. Roughly speaking, doing this is what is known as operational semantics. So I thought, why not apply this approach to the algebraic rules for defining and combining ordinals. The result are the interpreters run1 and run2 above.

run1 gives an example of a Hydra game. In fact, its precisely the hydra game described in this paper because it always chops off the rightmost head. The Kirby-Paris theorem tells us we can't prove this game terminates using just the Peano axioms. A web search on Goodstein's theorem will reveal many great articles with the details.

A well-ordered quantity that you can keep decreasing as a program runs, and that can be used to prove termination, is an example of a loop variant. Loop variants are often natural numbers but the above shows that transfinite ordinals make fine loop variants. But in the interest of being fair and balanced, here's a dissenting view. The author has a point. If you are forced to use transfinite ordinals to show your program terminates, the age of the universe will probably be but the briefest flicker compared to your program's execution. On the other hand, if you don't want an actual bound on the execution time, ordinals can provide very short proofs of termination for useful programs.

(By the way, this article is a sequel to my own article.)

Exercise
Algebraic structures give rise to monads. Can you see how to generalise the definition of Machine to make it a monad? If you pick the right definition then the substitution of a for Done in the definition of + should give you a particularly simple definition of ordinal addition. (See this for a hint on how substitution works in monads.)

> instance Show Machine
> instance Eq Machine
> instance Ord Machine

Constructing Intermediate Values

Introduction
The Intermediate Value Theorem, as it is usually told, tells us that if we have a continuous real-valued function f on the closed interval [0,1], with f(0)<0 and f(1)>0 then there is a point x in (0,1) such that f(x)=0. Here's a sketch of a proof:

Define three sequences l_i, u_i, x_i recursively as follows:

x_i+1 = (l_i+u_i)/2
(l_i+1, u_i+1) = (l_i, x_i+1) if f(x_i+1)>0
(l_i+1, u_i+1) = (x_i+1, u_i) otherwise

Set l₀=0, u₀=1.
The x_i+1 are bracketed by the l_i and u_i, and |u_i-l_i|=2^-i so the x_i form a Cauchy sequence with some limit x. Because f is continuous, we get f(x)=0.

This proof not only shows that we have an intermediate value x where f crosses the x-axis, it also seems to describe a procedure for computing x. Let's use this strategy to write a numerical method.

A Numerical Method

> import Ratio
> infixl 6 .*

> ivt0 :: (Float -> Float) -> Float
> ivt0 f = ivt' 0 1 where
>    ivt' l u =
>          let z = 0.5*(l+u)
>          in if z==l || z==u || f z == 0
>               then z
>               else if f z < 0
>                   then ivt' z u
>                   else ivt' l z

Here are some simple functions we can use to test it:

> f0, g0 :: Float -> Float
> f0 x = 3*x-1
> g0 x = 2*x-1

ivt0 f0 and ivt0 f1 give the expected results.

Of course, it might not actually be the case that f (ivt f) == 0. We're dealing with the type Float here and there might not actually be a solution to f x == 0 to the precision of a Float. And it's not clear that the word continuous means anything here anyway because the limited precision of a Float precludes us from constructing the deltas and epsilons the definition of continuity demands.

Switching to Exact Real Arithmetic
So can we do better? Can we implement a function to compute intermediate values for exact real arithmetic? I recently talked about computable real numbers, and they are what I intend to use here. I'll represent the real number x using a sequence of rationals, x_i, such that |x_i-x|<2^-i. Here's a suitable type:

> data R = R { runR :: Integer -> Rational }

(I'll only be calling it for non-negative arguments, but alas, Haskell has no natural number type.)

We'll need to display our numbers. Here's something cheap and cheerful that's good enough for this article:

> instance Show R where
>    show (R f) = show (fromRational (f 30) :: Float)

It displays a Float that is formed from a rational approximation that is within 2^-30 of our value.

We can inject the rationals into the computable reals by using constant sequences:

> instance Fractional R where
>    fromRational x = R $ const $ fromRational x

An important point here is that two different sequences might define the same real. So when we construct functions mapping reals to reals we must be sure that they take different representations of the same real to representations of the same real. I'll call these 'well-defined'.

Now our problems begin. Haskell wants us to implement Eq before Num, but equality of computable reals is not decidable. Still, it is semi-decidable in that if two reals differ, we can detect this. Here's a partial equality function:

> instance Eq R where
>    R f == R g = eq' 2 0 f g where
>        eq' delta i f g = 
>               let fi = f i
>                   gi = g i
>               in if fi<gi-delta
>                   then False
>                   else if fi>gi+delta
>                       then False
>                       else eq' (delta/2) (i+1) f g

The idea is that if f i and g i can be separated by more than 2*2^-i then the reals that f and g represent can't be equal.

We have to implement Ord too. If one computable real is bigger than another then they must differ by more than some power of 1/2. This means we can eventually distinguish between them. As a result, compare will always terminate if we compare distinct numbers:

> instance Ord R where
>    compare (R f) (R g) = compare' 2 0 f g where
>        compare' delta i f g =
>               let fi = f i
>                   gi = g i
>               in if fi<gi-delta
>                   then LT
>                   else if fi>gi+delta
>                       then GT
>                       else compare' (delta/2) (i+1) f g

We can define min and max straightforwardly:

>    min (R f) (R g) = R $ \i -> min (f i) (g i)
>    max (R f) (R g) = R $ \i -> max (f i) (g i)

And now we can define addition. We use the fact that if |x-x'|<ε and |y-y'|<ε then |(x+y)-(x'+y')|<2ε:

> instance Num R where
>    fromInteger n = R $ const $ fromInteger n
>    R f + R g = R $ \i -> let j = i+1 in f j + g j
>    negate (R f) = R $ negate . f

I'm going to skip multiplication of two reals. But it's easy to define multiplication of a rational with a real. Here I'm using the fact that if |x-x'|<ε, then |ax-ax'|<|a|ε:

> x .* R f = R $ mul x x f where
>    mul x x' f i | abs x' < 1 = x * f i
>                 | otherwise  = mul x (x'/2) f (i+1)

Now we're ready! It's straightforward to translate our previous algorithm. Note that we don't need a function to take the limit of a sequence. By definition we use sequences of rational approximations that are accurate to within 2^-i so our ivt1 function merely needs to compute an approximation to our intermediate value that is accurate to within 2^-i.

> ivt1 :: (R -> R) -> R
> ivt1 f = R $ \i -> ivt' 0 1 i where
>    ivt' :: Rational -> Rational -> Integer -> Rational
>    ivt' x y i =
>           let z = (1%2)*(x+y)
>           in if i==0
>                then z
>                else if f (fromRational z) < 0
>                    then ivt' z y (i-1)
>                    else ivt' x z (i-1)

Some new test functions:

> f, g :: R -> R
> f x = 3.*x-1
> g x = 2.*x-1

Note that (amazingly!) any well-defined function R -> R we implement is continuous because of the argument I sketched here.

And you'll find that ivt1 f gives the result you expect. But sadly, ivt1 g never returns. What went wrong?

Clutching Defeat from the Jaws of Victory
After one iteration, ivt1 sets z to the value we want. This seems like a particularly easy case. But the algorithm immediately fails because it then proceeds to compare f z with zero. We already know that our equality test must fail for this case. How do we fix our algorithm?

Consider the actions of any higher order function that operates on a continuous function f type R -> R. It will evaluate f at various real values. The return values from f will then be sampled at various values of i. Each sample corresponds to a bracket around a value of f. Consider a higher order function acting on the function below:



If our higher order function terminates, it can only compute a finite number of these brackets. So I can choose an ε small enough that I could nudge the graph up or down by ε without invalidating the brackets. I could provide the higher order function with a different function R -> R that responds in exactly the same way to the finite number of samples.

Now consider the case of an ivt function. It will respond the same way to two functions differing only by ε. The function in the graph above has been designed so that it is flat from 0.4 to 0.6. Nudging the function above by ε>0 up and down causes the crossing of the x-axis to move from x<0.4 to x>0.6. In other words, no algorithm could possibly return a valid intermediate value for all possible inputs because I can always contrive a pair of functions, to which ivt would respond identically, and yet which have widely separated points where they cross the x-axis. This shows that the problem of determining a crossing point of a continuous function that starts off less than zero, and ends up greater than zero, is uncomputable. We can do it for some functions. But no function can solve this problem for all continuous inputs. What's really weird about this is that it fails precisely when the problem seems easiest - when there are lots of zeroes!

A Revised Approach
But we don't need to give up. If we could ensure that at each bisection we could pick a point that wasn't exactly a zero, then we'd be fine. But how can we do this? Well here's a cheating way: we don't bother. We let that problem be the caller's responsibility. Instead of bisecting at each stage, we'll trisect. We'll pick a point that isn't a zero in the middle third by using a function pased in by the caller. Call this function g. g is of type Rational -> Rational -> Rational and the specification is that g x y must be a rational number between x and y that isn't a zero. Here's our new ivt2 function:

> ivt2 :: (Rational -> Rational -> Rational) -> (R -> R) -> R
> ivt2 g f = R $ \i ->
>     let delta = (1%2)^i
>         ivt' :: Rational -> Rational -> Rational
>         ivt' x y =
>               let l = (1%3)*(2*x+y)
>                   u = (1%3)*(x+2*y)
>                   z = g l u
>               in if u-l < delta
>                    then z
>                    else if f (fromRational z) < 0
>                        then ivt' z y
>                        else ivt' x z
>     in ivt' 0 1

How can we make use of this? Well suppose we want to find where a polynomial, f, crosses the x-axis. If the polynomial has degree n, it can only cross at n points. So if we pick n+1 values, f must evaluate to a non-zero value at one of them. We can examine the n+1 values of f(x) in parallel to find one that is non-zero. This must terminate. Here's an implementation for the case n+1=2. If we pass two real values to pickNonZero it will True or False according to whether the first or second argument is non-zero. It must only be called with a pair of numbers where one is non-zero:

> pickNonZero (R f) (R g) = pickNonZero' 0 1 where
>     pickNonZero' i e = let fi = f i
>                            gi = g i
>                        in if fi > e || fi < -e
>                           then True
>                           else if gi > e || gi < -e
>                               then False
>                               else pickNonZero' (i+1) ((1%2)*e)

Note how it examines the two values in parallel. If it simply tested the two numbers in turn against zero it would fail if the first were zero. (Interesting side note: we could write this using par.)

Suppose the function f is strictly monotonically increasing. Then we know that if y>x, then f(y)>f(x) so we can't have that both f(x) and f(y) are zero. We can use this to write a suitable helper function to write an intermediate value function that is 100% guaranteed to work for all monotonically increasing functions (assuming we don't run out of memory):

> monotonicPicker :: (R -> R) -> Rational -> Rational -> Rational
> monotonicPicker f x y = if pickNonZero
>                               (f (fromRational x)) (f (fromRational y))
>                             then x
>                             else y

> ivt3 :: (R -> R) -> R
> ivt3 f = ivt2 (monotonicPicker f) f

By varying monotonicPicker we can find crossing points for all manner of continuous function. Try ivt3 g. But no matter how hard we try, we can never write one that works for all continuous functions.

What went wrong?
So why did the Intermediate Value Theorem, despite its appearance of constructing an intermediate value, fail to yield an algorithm? For an answer to that, you'll have to wait until my next article. The surprising thing is that it's not a problem of analysis, but a problem of logic.

Acknowledgements
I have borrowed liberally from Paul Taylor, especially from here. However, I've chosen to phrase everything in terms of classical reasoning about Haskell programs. Any errors are, of course, mine. It was probably Andrej Bauer's awesome blog that got me interested in this stuff.

Afterword
I want to stress that I have made absolutely no mention of intuitionism or constructivism in what I've said above. I have simply used ordinary classical mathematics to reason about some Haskell programs. I'm leaving that for the next article.

---

Optimising pointer subtraction with 2-adic integers.

Here is a simple C type and a function definition:

struct A
{
    char x[7];
};

int diff(struct A *a, struct A *b)
{
    return a-b;
}

It doesn't seem like there could be much to say about that. The A structure is 7 bytes long so the subtraction implicitly divides by 7. That's about it. But take a look at the assembly language generated when it's compiled with gcc:

movl 4(%esp), %eax
subl 8(%esp), %eax
imull $-1227133513, %eax, %eax
ret

Where is the division by 7? Instead we see multiplication by -1227133513. A good first guess is that maybe this strange constant is an approximate fixed point representation of 1/7. But this is a single multiplication with no shifting or bit field selection tricks. So how does this work? And what is -1227133513? Answering that question will lead us on a trip through some suprising and abstract mathematics. Among other things, we'll see how not only can you represent negative numbers as positive numbers in binary using twos complements, but that we can also represent fractions similarly in binary too.

But first, some history.

Some n-bit CPUs



That's an Intel 4004 microprocessor, the first microprocessor completely contained on a single chip. It was a 4 bit processor equipped with 4 bit registers. With 4 bits we can represent unsigned integers from 0 to 15 in a single register. But what happens if we want to represent larger integers?

Let's restrict ourselves to arithmetic operations using only addition, subtraction and multiplication and using one register per number. Then a curious thing happens. Take some numbers outside of the range 0 to 15 and store only the last 4 bits of each number in our registers. Now perform a sequence of additions, subtractions and multiplications. Obviously we usually expect to get the wrong result because if the final result is outside of our range we can't represent it in a single register. But the result we do get will have the last 4 bits of the correct result. This happens because in the three operations I listed, the value of a bit in a result doesn't depend on higher bits in the inputs. Information only propagates from low bit to high bit. We can think of a 4004 as allowing us to correctly resolve the last 4 bits of a result. From the perspective of a 4004, 1 and 17 and 33 all look like the same number. It doesn't have the power to distinguish them. But if we had a more powerful 8 bit processor like the 6502 in this machine, we could distinguish them:

Introduction
This is analogous to the situation we have with distances in physical world. With our eyes we can resolve details maybe down to 0.5mm. If we want to distinguish anything smaller we need more powerful equipment, like a magnifying class. When that fails we can get a microscope, or an electron microscope, or these days even an atomic force microscope. The more we pay, the smaller we can resolve. We can think of the cost of the equipment required to resolve two points as being a kind of measure of how close they are.

We have the same with computers. To resolve 1 and 17 we need an 8-bit machine. To resolve 1 and 65537 we need a 32-bit machine. And so on. So if we adopt a measure based on cost like in the previous paragraph, there is a sense in which 1 is close to 17, but 1 is even closer to 257, and it's closer still to 65537. We have this inverted notion of closeness where numbers separated by large (in the usual sense) powers of two are close in this new sense.

We have an interesting relationship between computing machines with different 'resolving' power. If we take an arithmetical computation on an N-bit machine, and then take the last M bits of the inputs and result, we get exactly what the M-bit machine would have computed. So an M-bit machine can be thought of as a kind of window onto the last M-bits onto an N-bit machine. Here's a sequence of machines:

Each machine provides a window onto the low bits of the previous machine in the sequence. But what happens at the "..." on the left? That suggests the bizarre idea that maybe all of these finite machines could be thought of as window to some infinite bit machine. Does that idea make any kind of sense?

I'll try to convince you that's a sensible idea by pointing out that it's something familiar to anyone who's taken a rigorous analysis course. (And I'll mention in passing that the above diagram illustrates a limit in an appropriate category!)

Mathematicians (often) build the real numbers from the rational numbers by a process known as completion. Consider a sequence like 1, 14/10, 141/100, 1414/1000, ... . The nth term is the largest fraction, with 10ⁿ in the denominator, such that its square is less than 2. It's well known that there is no rational number whose square is 2. And yet it feels like this sequence ought to be converging to something. It feels this way because successive terms in the sequence get as close to each other as you like. If you pick any ε there will be a term in the series, say x, with the property that later terms never deviate from x by more than ε. Such a sequence is called a Cauchy sequence. But these sequences don't all converge to rational numbers. A number like √2 is a gap. What are we to do?

Mathematicians fill the gap by defining a new type of number, the real number. These are by definition Cauchy sequences. Now every Cauchy sequence converges to a real number because, by definition, the real number it converges to is the sequence. For this to be anything more than sleight of hand we need to prove that we can do arithmetic with these sequences. But that's just a technical detail that can be found in any analysis book. So, for example, we can think of the sequence I gave above as actually being the square root of two. In fact, the decimal notation we use to write √2, 1.414213..., can be thought of as shorthand for the sequence (1, 14/10, 141/100, ...).

The notion of completeness depends on an an idea of closeness. I've described an alternative to the usual notion of closeness and so we can define an alternative notion of Cauchy sequence. We'll say that the sequence x₁, x₂, ... is a Cauchy sequence in the new sense if all the numbers from x_n onwards agree on their last n bits. (This isn't quite the usual definition but it'll do for here.) For example, 1, 3, 7, 15, 31, ... define a Cauchy sequence. We consider a Cauchy sequence equal to zero if x_n always has zero for its n lowest bits. So 2, 4, 8, 16, 32, ... is a representation of zero. We can add, subtract and multiply Cauchy sequences pointwise, so, for example, the product and sum of x_n and y_n has terms x_ny_n. Two Cauchy sequences are considered equal if their difference is zero. These numbers are called 2-adic integers.

Exercise: prove that if x is a 2-adic integer then x+0=x and that 0x=0.

There's another way of looking at 2-adic integers. They are infinite strings of binary digits, extending to the left. The last n digits are simply given by the last n digits of x_n. For example we can write 1, 3, 7, 31, ... as ...1111111. Amazingly we can add subtract and multiply these numbers using the obvious extensions of the usaul algorithms. Let's add ...1111111 to 1:

...11111111
...00000001
-----------
...00000000

We get a carry of 1 that ripples off to infinity and gives us zeroes all the way.

We can try doing long multiplication of ...111111 with itself. We get:

...11111111
...1111111 
...111111  
...11111   
...
-----------
...00000001

It's important to notice that even though there are an infinite number of rows and columns in that multiplication you only need to multiply and add a finite number of numbers to get any digit of the result. If you don't like that infinite arrangement you can instead compute the last n digits of the product by multiplying 11...n digits...111 by itself and taking the last n digits. The infinite long multiplication is really the same as doing this for all n and organising it in one big table.

So ...1111111 has many of the properties we expect of -1. Added to 1 we get zero and squaring it gives 1. It is -1 in the 2-adic integers. This gives us a new insight into twos complement arithmetic. The negative twos-complements are the truncated last n digits of the 2-adic representations of the negative integers. We should properly be thinking of twos-complement numbers as extending out to infinity on the left.

The field of analysis makes essential use of the notion of closeness with its δ and ε proofs. Many theorems from analysis carry over to the 2-adic integers. We find ourselves in a strange alternative number universe which is a sort of mix of analysis and number theory. In fact, people have even tried studying physics in p-adic universes. (p-adics are what you get when you repeat the above for base p numbers, but I don't want to talk about that now.) One consequence of analysis carrying over is that some of our intuitions about real numbers carry over to the 2-adics, even though some of our intuitive geometric pictures seem like they don't really apply. I'm going to concentrate on one example.

The Newton-Raphson Method
I hope everyone is familiar with the Newton-Raphson method for solving equations. If we wish to solve f(x)=0 we start with an estimate x_n. We find the tangent to y=f(x) at x=x_n. The tangent line is an approximation to the curve y=f(x) so we solve the easy problem of finding where the tangent line crosses the x-axis to get a new estimate x_n+1. This gives the formula

x_n+1 = x_n-f(x_n)/f'(x_n).

With luck the new estimate will be closer than the old one. We can do some analysis to get some sufficient conditions for convergence.

The surprise is this: the Newton-Raphson method often works very well for the 2-adic integers even though the geometric picture of lines crossing axes doesn't quite make sense. In fact, it often works much better than with real numbers allowing us to state very precise and easy to satisfy conditions for convergence.

Now let's consider the computation of reciprocals of real numbers. To find 1/a we wish to solve f(x)=0 where f(x)=1/x-a. Newton's method gives the iteration x_n+1 = x_n(2-ax_n). This is a well know iteration that is used internally by CPUs to compute reciprocals. But for it to work we need to start with a good estimate. The famous Pentium divide bug was a result of it using an incorrect lookup table to provide the first estimate. So let's say we want to find 1/7. We might start with an estimate like 0.1 and quickly get estimates 0.13, 0.1417, 0.142848, ... . It's converging to the familiar 0.142857...

But what happens if we start with a bad estimate like 1. We get the sequence:

1
-5
-185
-239945
-403015701065
-1136951587135200126341705

It's diverging badly. But now let's look at the binary:

00000000000000000000000000000000000000000000000000000000000000000000000000000000001
11111111111111111111111111111111111111111111111111111111111111111111111111111111011
11111111111111111111111111111111111111111111111111111111111111111111111111101000111
11111111111111111111111111111111111111111111111111111111111111111000101011010110111
11111111111111111111111111111111111111111111010001000101010011001000110110110110111
11100001111001111011011101100100000010000000011011010110110110110110110110110110111

Our series may be diverging rapidly in the usual sense, but amazingly it's converging rapidly in our new 2-adic sense!

If it's really converging to a meaningful reciprocal we'd expect that if we multiplied the last n digits of these numbers by 7 then we'd get something that agreed with the number 1 in the last 7 digits. Let's take the last 32 digits:

10110110110110110110110110110111

and multiply by 7:

10100000000000000000000000000000001

The last 32 bits are

00000000000000000000000000000001.

So if we're using 32 bit registers, and we're performing multiplication, addition and subtraction, then this number is, to all intents and purposes, a representation of 1/7. If we interpret as a twos complements number, then in decimal it's -1227133513. And that's the mysterious number gcc generated.

Epilogue
There are many things to follow up with. I'll try to be brief.

Try compiling C code with a struct of size 14. You'll notice some extra bit shifting going on. So far I've only defined the 2-adic integers. But to get he reciprocal of every non-zero number we need numbers whose digits don't just extend leftwards to infinity but also extend a finite number of steps to the right of the "binary point". These are the full 2-adic numbers as opposed to merely the 2-adic integers. That's how the extra shifts can be interpreted. Or more simply, if you need to divide by 14 you can divide by 2 first and then use the above method to divide by 7.

I don't know how gcc generates its approximate 2-adic reciprocals. Possibly it uses something based on the Euclidean GCD algorithm. I wasn't able to find the precise line of source in a reasonable time.

An example of a precise version of the Newton-Raphson method for the p-adics is the Hensel lemma.

The last thing I want to say is that all of the above is intended purely to whet your appetite and point out that a curious abstraction from number theory has an application to compiler writing. It's all non-rigourous and hand-wavey. Recommend reading further at Wikipedia. I learnt most of what I know on the subject from the first few chapters of Koblitz's book p-adic Numbers, p-adic Analysis, and Zeta functions. The proof of the von Staudt–Clausen theorem in that book is mindblowing. It reveals that the real numbers and the p-adic numbers are equally valid ways to approximately get a handle on rational numbers and that there are whole alternative p-adic universes out there inhabited by weird versions of familiar things like the Riemann zeta function.

(Oh, and please don't take the talk of CPUs too literally. I'm fully aware that you can represent big numbers even on a 4 bit CPU. But what what I say about a model of computation restricted to multiplication, addition and subtraction in single n-bit registers holds true.)

Some Challenges
1. Prove from first principles that the iteration for 1/7 converges. Can you prove how many digits it generates at a time?
2. Can you find a 32 bit square root of 7? Using the Newton-Raphson method? Any other number? Any problems?

Acknowledgements
Update: I have replaced the images with images that are, to the best of my knowledge, public domain, or composited from public domain images. (Thanks jisakujien. My bad.)

Update
If you want to play with some p-adics yourself, there is some code to be found here. That also has code for transcendental functions applied to p-adics.

Here's some C code to compute inverses of odd numbers modulo 2³² (assuming 32 bit ints). Like the real valued Newton method, it doubles the number of correct digits at each step so we only need 5 iterations to get 32 bits. (C code as I think it's traditional to twiddle one's bits in C rather than Haskell.)

#include <stdio.h>
#include <assert.h>

typedef unsigned int uint;

uint inverse(uint x)
{
   uint y = 2-x;
   y = y*(2-x*y);
   y = y*(2-x*y);
   y = y*(2-x*y);
   y = y*(2-x*y);

   return y;
}

int main()
{
   uint i;
   for (i = 1; i<0xfffffffe; i += 2)
   {
       assert (i*inverse(i) == 1);
   }
}

---

On representing some real numbers exactly

There are uncountably many real numbers, and only countably many finite strings of symbols, so we know for sure that there is no scheme to represent all real numbers as finite strings of symbols in such a way that different reals get different representations. However, there are useful subsets of the real numbers that *are* equipped with finite encoding schemes. I hope to give some examples of such sets, including one, the set of periods, that isn't currently well known outside the mathematical world.

It's traditional to think of the real numbers as divided up into three types of number: the rational numbers, the algebraic numbers, and the leftovers known collectively as the transcendental numbers. If you're not familiar with those terms I'll be explaining them below. We can represent this as follows:

ℚ⊂A⊂ℝ

I'm using unicode for blackboard bold R and Q there. I hope your browser supports them. ℝ is the set of real numbers, ℚ is the set of rational numbers, and A is the algebraic numbers.

But there are many different kinds of transcendental number in existence, so the subject of this article will be three more levels of classification we could insert between the set of algebraic numbers and the set of real numbers in the sequence above.

But first, let's revisit ℚ and A.

The Rational Numbers, ℚ
The rational numbers are fairly straightforward to deal with. They are nothing more than the ratios of integers. Given two such ratios we can easily compare them to see whether they are equal. For example, we're taught at a young age how to tell if 2/3 is the same as 4/6.

However, it didn't take long for the Ancient Greeks to realise that they had some equations they could solve approximately using rationals, but couldn't solve exactly. The best known example is solving x²=2.

Assume x is rational, so that x=p/q, with p and q integers with highest factors cancelled out. We have p²=2q², and so p must be divisible by 2. But then p² is divisible by 4, and so q is divisible by 2. This means that we can cancel 2 from the top and bottom of p/q. That contradicts our assumption about cancellation. So the equation has no rational solution.

If we expect x²=2 to have a solution we must work with numbers that can't be represented as fractions.

Notice how the rational numbers come equipped with a notation to describe them. We can just write one integer over another, with a horizontal line between them. Given such a description we can immediately tell if it's valid (it's only invalid if the denominator is 0) and whether it equals some other rational number.

One interesting property of the rational numbers is that they are *dense* in the real numbers. This is just another way of saying that any real number can be approximated as well as we like by a rational number. For example, we can approximate √2 to within 1/1000 by the rational number 141421/10000. This is an important property that will come up again later.

We also know that there are only countably many rational numbers because there are only countably many pairs of integers.

The Real Algebraic Numbers, A
The real algebraic numbers are all the real numbers that we can obtain by finding roots of polynomial equations whose coefficients are rational numbers. For example, the real solutions to x⁵-2x²+½=0 are all algebraic numbers. Among other things, this includes the solution to x²=2 that wasn't contained in the rationals. The algebraic real numbers also include the solutions to "insoluble" algebraic equations. (By 'algebraic' equations, I just mean polynomial equations.) Even though Galois proved we can't write expressions for the solutions to all algebraic equations using the four arithmetic operations and nth roots, these solutions still exist in the real numbers, and they are all, by definition, algebraic numbers.

We can describe an algebraic number by writing down the equation it solves, and additionally providing some description to say which root of the equation we're interested in. Unfortunately, given any algebraic number there is always an infinite number of algebraic equations that it satisfies. So like with the rationals there is some redundancy. The good news, however, is that given two such descriptions of algebraic numbers there is an algorithm to tell whether or not they describe equal numbers.

Even though there are vastly more algebraic numbers than rational numbers, they still only form a countable set. As mentioned in the previous paragraph, every algebraic number can be described by a finite string of symbols, and there are only countably many such strings.

The Computable Real Numbers, C
Now I want to insert a set between A and ℝ. These are the computable real numbers. But first think back to what I said about the rational numbers.

Pick a real number x. Given any ε we can always find a rational number that comes within ε of x. We can define a function f on the rationals with the property that f(ε) is rational and within ε of x for any rational ε. So f gives rational approximations to x to any desired accuracy. Such a function uniquely specifies x. In fact x is the limit of f(ε) as ε tends to zero. The computable real numbers are the real numbers specified in this way by *computable* functions. Notice how there are no infinities involved here. A computable real number is represented by nothing more than a finite string of symbols forming a computer program that processes finite integers. We can find out what the real number is to any desired accuracy in a finite time.

Because we can write computer programs to find roots of algebraic equations to any desired accuracy we know that the algebraic numbers are contained in the computables. So now we have:

ℚ⊂A ⊂C ⊂ℝ

Where I'm using C to represent the computables.

Almost any number we ever want to work with is computable. Whether it's √2, or π, or the value of some humongous integral that's appeared in your engineering problem, chances are there's an algorithm somewhere that can approximate it as accurately as you like, given enough CPU time and RAM.

Note that every computable is described by a finite string of symbols, and so the computables form a countable set.

There are some problems with computable numbers. Given two representations of computable numbers we have no way of telling whether or not they are equal. It's not just that this is hard to do. There simply is no algorithm to prove that the numbers are equal. The problem is that the two computer programs may give exactly the same approximations to every degree of accuracy down to a certain ε. But for higher accuracy they may give different results. We have no way of knowing in advance at what value of ε they'll start differing. And anyway, there simply is no algorithm for telling if two computer programs will always generate the same results. So paradoxically, testing the equality of two computable numbers is itself uncomputable.

Actually, the situation is much worse. We can't even tell if we have a valid representation of a computable real. In order to do that, we need to know that our computer program to generate approximations terminates. But solving that would solve the halting problem.

By the way, despite these issues this representation of computable real numbers as functions leads to practical ways to do arbitrary precision arithmetic on a computer in a way that gives us guaranteed bounds on the accuracy of our results.

From here we could go two ways. We could try to rein in our class of number to something more manageable. Or we could try to push further and try to find ways to represent real numbers that aren't even computable. Let's do the former first.

The Periods, P
The periods are an interesting class of real number that doesn't seem to be well known. The definition may be old, but it was only recently that they were being promoted in mathematical circles as an interesting thing to study. So I thought I'd contribute to their promotion.

Consider the number π. It arises straightforwardly in geometric problems. An example is computing the area of a circle. Back in the 18th century Lambert proved that it was irrational and in the 19th century Lindemann proved that it was transcendental. But in a sense the transcendental numbers are simply a rubbish heap into which the leftover numbers have been discarded. Can we get some kind of handle on at least some of these numbers in a way that puts π back on the map?

The real number π is the area of a unit circle. So it is the area of the region of the plane given by the equation x²+y²<1.



Here's a similar construction of another real number:



The area of that region is log(2). (Natural logarithm of course!) We can construct the logarithm of any positive rational number in a similar way. These are all transcendental numbers. But in a way they are nice transcendental numbers. They arise from considering areas of the plane described by straightforward algebraic inequalities. (As in the previous section, I'm using 'algebraic inequalities' to mean 'polynomial inequalities'.)

That suggests a class of real number: those numbers that can be represented as the volume of a region defined by a bunch of algebraic inequalities with rational coefficients. These are known as the periods. I'm talking about generalised volumes in n-dimensions, not just areas in the plane. Clearly π and log(2) are periods.

But does this representation solve the problem we had with computable numbers where we were unable to guarantee we could check the equality of two numbers? Now it gets interesting. The answer is: we don't know. It is conjectured that if we have the same number represented in two different ways as a period that we can transform one representation into another simply by using a small set of elementary operations. It is also conjectured that there is a terminating algorithm for finding such a sequence of operations, or if one doesn't exist, demonstrating this.

One obvious question now is this: are there any computable reals that aren't periods? A recent paper, Periods and elementary real numbers claims to exhibit one by means of a kind of diagonalisation argument. But in general it's hard to prove that a number isn't a period. I don't believe there is a proof yet that e = 2.718... isn't a period, but mathematicians expect that it isn't.

I think the name 'period' comes from the fact that the period of a pendulum of rational length, in a gravitational field of rational strength, is a period. Computing this requires elliptic functions, and these often result in periods when applied to rational numbers.

Although mathematics is about numbers, few mathematical publications talk much about specific real numbers. Curiously, when they do, they often talk about numbers that are periods. For example there have been many recent papers on the values of the Riemann zeta function for interesting arguments. These are often periods.

It has been suggested that the study of periods is actually the study of algebraic geometry in disguise. There is certainly no end to the interesting mathematics we can do using only periods.

We now have:

Q ⊂A ⊂P ⊂C ⊂R

The Definables
In the section on computable numbers I mentioned that we could go the other way and try to find a bigger class of real numbers that could be represented by finite strings of symbols. We could simply go "all the way" and consider those real numbers that can be defined, by any means possible, using the symbols of mathematics. We can try to pin this down a bit better. We'll work with the language of set theory. In this language we can write strings of symbols like S = "x>0 and x²=2". Such a string uniquely defines a real number if when we glue the string "there exists a unique x such that" onto the beginning of it we get a true proposition. We can now represent the number √2 as the string S.

It looks like we now have:

Q ⊂A ⊂P ⊂C ⊂D ⊂R

with D the set of definables.

Except there's a problem. The set D is not definable! The problem is this: we can talk about strings of symbols representing real numbers. To talk about these in set theory we'd encode these strings of symbols as sets so that we can apply the language of set theory. In order to use our attempted definition of definable we need some way to say when a string of symbols represents a true proposition. But to define a set of definables we need to talk about true propositions in the language of set theory. Godel showed us how to talk about the provability of a proposition within set theory. He did this by showing that provability is about mechanical operations we can perform on strings. But there's nothing analogous for talking about the truth of propositions. In fact, Tarski showed us this is impossible. So while we can talk about all kinds of individual numbers as being definable, we can't construct the set of definable numbers.

You may be interested to see an example of a definable number that isn't computable. Probably the most publicised example is Chaitin's constant. It represents the probability that a randomly generated string of symbols is a computer program that terminates in a finite time. We can't actually compute this number because it requires us to solve the halting problem. Nonetheless, it's perfectly well defined.

You can find the original paper on periods by Zagier and Kontsevich here. I first found out about periods from the book Mathematics Unlimited.

Update: Jared asked if I intended the real algebraic numbers. That's what I said in the text, but I confusingly wrote ℚ for this set. But that usually means the *complex* algebraic numbers. So I've changed notation and now use A for the real algebraic numbers.

Similarly, the periods are usually defined to be complex numbers whose real and imaginary parts are given by algebraically specified volumes. I was trying to avoid mention of the complex numbers to keep prerequisites to a minimum and the essential ideas here work without mention of the complex numbers.