On Removing Singularities from Rational Functions

Introduction
Suppose we have the function

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

Some basic algebraic manipulation shows that in the limit as x→0, f(x)→1. But we can't simply compute f 0 because this computation involves division by zero at intermediate stages. How can we automate the process of computing the limit without implementing symbolic algebra?

I've already described one way to remove singularities from a function. But that approach is very limited in its applicability.

This article is about a variation on the approach to formal power series that nicely showcases some advantages of lazy lists. It will allow us to form Laurent series of functions so we can keep track of the singularities.

The usual Haskell approach to power series allows you to examine the coefficients of any term in the power series of the functions you can form. These series can't be used, however, to evaluate the function. Doing so requires summing an infinite series, but we can't do so reliably because no matter how many terms in a power series we add, we can never be sure that there aren't more large terms further downstream that we haven't reached yet. And if we want to perform computations completely over the rationals, say, we don't want to be dealing with infinite sums.

I'd like to look at a way of working with power series that allows us to perform exact computations making it possible to answer questions like "what is the sum of all the terms in this power series starting with the x^n term?" By extending to Laurent series, and implementing the ability to selectively sum over just the terms with non-negative powers, we can compute functions like f above at 0 and simply skip over the troublesome http://en.wikipedia.org/wiki/Pole_%28complex_analysis%29.

Power Series
When I previously discussed power series I used code that worked with the coefficients of the power series. This time we want to work with values of the function so it makes sense to store, not the coefficients a_i but the terms themselves, a_ixⁱ. So instead of a list of coefficients, Num a => [a] we need a representation that looks a little like:

> data Power a = Power (a -> [a])

where we pass x in as an argument to the function contained in a Power. But we also want to allow Laurent series so we need to also store an offset to say which (possibly negative) term our series starts with:

> data Laurent a = Laurent (a -> (Int, [a]))

But this fails us for at least two reasons:

1. We have the individual terms, but to evaluate the function requires summing all of the terms in an infinite list.
2. If we have a Laurent series, then we need to store values of a_ixⁱ for x=0 and i<0. We'll end up with division by zero errors.

Partial Sum Series
So here's what we'll do instead. Suppose our power series is Σ_i=n^∞a_ixⁱ. We'll store the terms s_j=Σ_i=j^∞a_ix^i-j. Our type will look like:

> data Partial a = Partial (a -> (Int, [a]))

> instance Eq (Partial a)
> instance Show (Partial a)

It's straightforward to add two functions in this form. We just add them term by term after first aligning them so that the xⁱ term in one is lined up with the xⁱ term in the other:

> instance Num a => Num (Partial a) where

>  Partial f + Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         pad 0 _ ys = ys
>         pad n x ys = let z:zs = pad (n-1) x ys
>                      in x*z : z : zs
>         l = min m n
>     in (l, zipWith (+) (pad (m-l) x xs) (pad (n-l) x ys))

Notice the slight subtlety in the alignment routine pad. By the definition above, the jth term has a factor of x^j built into it. So we need to multiply by x each time we pad our list on the left.

Now we need to multiply series. We know from ordinary power series that we need some sort of convolution. But it looks like for this case we have an extra complication. We appear to need to difference our representation to get back the original terms, convolve, and then resum. Amazingly, we don't need to do this at all. We can convolve 'in place' so to speak.

Here's what an ordinary convolution looks like when we want to multiply the sequence of terms (a_i) by (b_i):

In this example, the blue diagonal corresponds to the terms that are summed to get the 4th term in the result.

However, we wish to work with partial sums s_j=Σ_i=j^∞a_ix^i-j and t_j=Σ_i=j^∞t_ix^i-j, constructing the partial sums of the convolution of a and b from s and t. The partial sums of the convolution can be derived from the partial sums by tweaking the convolution so it looks like this:

The blue terms work just like before and need to be summed. But we also need to subtract off the red terms, weighted by a factor of x. That's it! (I'll leave that as an exercise to prove. The inclusion-exclusion principle helps.)

The neat thing is that the red terms for each sum are a subset of the blue terms needed for the next element. We don't need to perform two separate sums. We can share much of the computation between the red and blue terms. All we need to do is write an ordinary convolution routine that additionally returns not just the blue terms, but a pair containing the blue sum and the red sum.

>  Partial f * Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         (outer, inner) = convolve xs ys
>         f' a b = a-x*b -- (the subtraction I mentioned above)
>     in (m+n, zipWith f' outer inner)
>  fromInteger n = let n' = fromInteger n in Partial $ \_ -> (0, n' : repeat 0)

>  negate (Partial f) = Partial $ \x -> let (m, xs) = f x
>                                 in (m, map negate xs)
>  signum = error "signum not implemented"
>  abs    = error "signum not implemented"

This is an ordinary convolution routine tweaked to return the partial sum inner.

> convolve (a0:ar@(a1:as)) ~(b0:br@(b1:bs)) =
>  let (inner, _) = convolve ar br
>      ab = map (a0 *) bs
>      ba = map (* b0) as
>  in (a0*b0 : a0*b1+a1*b0
>            : zipWith3 (\a b c -> a+b+c) inner ab ba, 0 : inner)

The code is very similar to the usual power series multiplication routine. We can also use the same method described by McIlroy to divide our series.

As our series are a munged up version of the usual power series it's pretty surprising that it's possible to divide with so little code:

> instance Fractional a => Fractional (Partial a) where
>   fromRational n = let n' = fromRational n in Partial $ \_ -> (0, n' : repeat 0)
>   recip (Partial f) = Partial $ \x ->
>      let nibble (n, x:xs) | x==0      = nibble (n+1, xs)
>                           | otherwise = (n, (x:xs))
>          (n, xs) = nibble (f x)
>      in (-n, rconvolve x xs)

In effect, rconvolve solves the equation convolve a b==1:

> rconvolve x (a0:ar@(a1:as)) =
>   let (outer, inner) = convolve ar result
>       f a b = x*b-a
>       r = -1/f a0 a1
>       result = recip a0 : (map (r *) $ zipWith f outer inner)
>   in result

Note one ugly quirk of this code. I need to 'nibble' off leading zeroes from the series. This requires our underlying type a to have computable equality. (In principle we can work around this using parallel or).

That's it. We can now write a function to compute the positive part of a rational function. (By positive part, I mean all of the terms using non-negative powers of x.)

> pos f z = let Partial g = f $ Partial $ \x -> (1, 1 : repeat 0)
>               (n, xs) = g z
>           in if n>0
>               then z^n*head xs
>               else xs!!(-n)

Here are some examples:

> test1 = let f x = (1+2*x)/(3-4*x*x)
>         in pos (\x -> 1/(f x-f 0)/x) (0::Rational)

> test2 = pos (\x -> 1/(1+4*x+3*x^2+x^3) - 1/(1+x)) (1::Rational)

The original example I started with:

> test3 = pos (\x -> 1/(x+x^2) - 1/(x+2*x^2)) (0::Rational)

No division by zero anywhere!

Conclusions
The code works. But it does have limitations. As written it only supports rational functions. It's not hard to extend to square roots. (Try writing the code - it makes a nice exercise.) Unfortunately, any implementation of square root will (I think) require a division by x. This means that you'll be able to compute the positive part away from zero, but not at zero.

This method can't be extended fully to transcendental functions. But it is possible to add partial support for them. In fact, So with a little work we can still compute the positive part of functions like 1/sqrt(cos x-1) away from x==0. But applying cos to an arbitrary rational function may need more complex methods. I encourage you to experiment.

Note that this code makes good use of laziness. If your function has no singularities then you might find it performs no computations beyond what is required to compute the ordinary numerical value.

Automatic even/odd splitting

Statement of the Problem
Suppose you have a real valued function on the reals, say f. We can split it into the sum of an even and odd part:

f(x) = f_odd(x)+f_even(x)

where

f_odd(x) = (f(x)-f(-x))/2, f_even(x) = (f(x)+f(-x))/2

If f_odd has a power series around zero, then all of its terms must have odd powers in x. So f_odd(x)/x must have all even powers and it becomes natural to 'compress' down the terms to form f'(x) = f_odd(sqrt(x))/sqrt(x). (Take that as a definition of f' for this article.) We can implement this operation as a higher order function:

> import Prelude hiding (odd)

> odd f x = let s = sqrt x in 0.5*(f s - f (-s))/s

Here's a simple example:

> f x = x*(1+x*(2+x*(3+7*x)))

> test0 = odd f 3

But there's something not quite right about this. If f is rational, then so is odd f. But the implementation of odd involves square roots. Among other things, square roots introduce inaccuracy. As square roots don't appear in the final result, can we eliminate them from the intermediate steps of the computation too?

A Better Solution
Let's use the results of last week's article to compute odd another way. We want a linear function that maps as follows:

1	->	0
x	->	1
x2	->	0
x3	->	x
x2n	->	0
x2n+1	->	xn

Here's an automaton:

If we start at 0, end at 1, and take exactly n steps, then the product of the factors we collect up along the way is given by the second column of that table. In n is even there is no such path so we collect up 0. As we're working with polynomials over the reals, rather than types, we have x1=1x and so on. We can construct a transition matrix:

0	1
x	0

We now do something like we did last time. Any time we have a function of some variable x, we replace x with the transition matrix. Our functions now take matrix values like

a	b
c	d

Any polynomial of our transition matrix always gives us equal elements along the diagonal. This is true even if we form the inverse of the transition matrix. So we don't need to store d. So now we implement a simple matrix type needing only to store three elements instead of four:

> data O a = O a a a deriving (Show, Eq)

> instance Num a => Num (O a) where
>    O a b c + O a' b' c' = O (a+a') (b+b') (c+c')
>    O a b c * O a' b' c' = O (a*a'+b*c') (a*b'+b*a') (c*a'+a*c')
>    fromInteger n = let i = fromInteger n in O i 0 0
>    negate (O a b c) = O (negate a) (negate b) (negate c)

Notice how similar this is to automatic differentiation. We can extend to reciprocals too:

> instance Fractional a => Fractional (O a) where
>    fromRational n = let i = fromRational n in O i 0 0
>    recip (O a b c) = let idet = recip (a*a-b*c) in O (idet*a) (idet*negate b) (idet*negate c)

And now we can implement a replacement for odd:

> transition x = O 0 1 x
> odd' f x = let O _ fx _ = f (transition x) in fx

> test1 = odd' f 3

Another example:

> test3 = odd' (\x -> (x+3*x*x-1/x)/(x*x)) 2
> test4 = odd  (\x -> (x+3*x*x-1/x)/(x*x)) 2

This new version has many advantages: it uses only rational functions, it's more accurate, and it's well defined at zero.

Conclusion
Automatic differentiation is just one of a family of methods that can be used to compute a wide variety of functions of real-valued functions. Essentially we're just working over real-valued matrices instead of real numbers. By using automata we can simplify the process of working out which matrices to use. (Though for the simple example above, you may have been able to guess the matrix without any other help).

(BTW I think there are hidden automata lurking in a few places in mathematics. For example, in Umbral calculus.)

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);
   }
}

---