In my previous post on the β function, I mentioned that I would like to rewrite it using Java's BigInteger class in order to escape the overflow errors that arise on almost any list with more than 3 elements.  Although I never did this (I did mess around with redoing the OCaml implementation using their big_int data type), I decided to redo the implementation in Haskell.  Not only have I wanted to learn Haskell for assorted reasons, but I realized that it has an Integer data type which allows arbitrary precision integers.  So now I could kill two birds (learning Haskell and implementing the beta function with arbitrary precision integers) with one stone.

Below you will find the code of the Haskell implementation.  I will follow the code with some comments about it, including the features of Haskell that this project helped me uncover and that I like very much.  Note, however, that this post is NOT a general introduction to Haskell or a general comparison of Haskell with any other language.  Links to such resources will be provided at the end of the post. First, note that the Haskell syntax for type declarations is e::t instead of OCaml's e:t. (In Haskell, ':' is the list cons operation.)  I have also set up a git repository on GitHub containing all of the implementations.  The Haskell one can be found there as well.

Although I mentioned that 'e::t' are typing declarations, you may have noticed the notation 'beta::Integral a => [a] -> a -> a'. In this case, 'a' is a type variable. Thus the expression 'beta' has a polymorphic type, similar to 'a in OCaml.  Now, everything before the '=>' is a type constraint, that is, constraints that the type variables must satisfy.  There can be multiple such constraints even though all the examples in this case use only one.

What does the constraint 'Integral a' mean though?  Integral is a type class in Haskell. A type class is akin to a Java interface in that it provides a list of functions (and their types) that must be applicable to a type in order for it to be an instance of the class. The Integral type class provides methods for integer division, modulus, and the like.  It inherits from the Num type class (Integral is actually a subclass of Real, not of Num, but that's irrelevant for our purposes) methods for addition, subtraction, and multiplication. Now, particular types are instances of a type class. The only types that are instances of Integral are Int and Integer. Int is akin to standard integer classes, while Integer (as mentioned at the outset) allows for arbitrary precision integers. For more on Haskell numbers, see this page.

Thus, in the case of 'beta::Integral a => [a] -> a-> a', this type says that beta is of type [a] -> a -> a ([a] is akin to 'a list in OCaml) subject to the constraint that 'a' must be an instance of the type class Integral (i.e. must be an Int or Integer). The beautiful thing is that integer constants in Haskell can be either Int or Integer; it will decide when it needs to use arbitrary precision and when not.  (Note: I do not know nearly enough about the inner workings of Haskell to know how it does so.)  So, in the file above, the factorial function can be called both on small numbers like 3 and 5, and as 'factorial 3813474071308470184' on numbers that are outside the range of Int.

A few other oddities that should be pointed out:

• The `mod` and `div` refer to the modulus and division functions of the Integral type class. Surrounding any Haskell function by backticks allows it to be used in infix notation. So 'mod 2 3' is equivalent to '2 `mod` 3' although the latter is more natural to read.
• The case pattern matching at the very end is necessary because the find function in Data.List has a return type of 'Maybe a'.  Thus, if the function finds no element satisfying the given predicate, it will come back with 'Nothing'.  Note that in such a case, OCaml would raise an exception.
• 'listeq' is a helper function that tests two lists for equality.  It also requires that the lists have the same elements in the same order; this requirement is actually critical for the β function to work as intended.
• '\x -> f' is Haskell's notation for lambda abstraction, i.e. for anonymous function declaration.  While I could have inlined the functions 'remlist' and 'beta_godel', I left them separately in order to resemble my previous implementations.  I left them separately in those implementations to remain similar to Enderton's exposition.
• '[0 ... n]' is a very nice way to denote the list of numbers from 0 to n.  Although I didn't use it in this implementation, Haskell has very nice list comprehensions.

To do

Although this implementation does work (I will post examples shortly), more remains to be done.

1. I still need to investigate . Two things directly affect the size of this number: and .
2. Haskell has a type inference engine so that the type declarations in this file are not necessary (although I do find them useful).  I'm curious to see whether Haskell will infer more general types than I've given the functions.
3. I need to add more detailed comments to the source files and clean up the git repository.
4. In due time, I will do real speed testing of the three implementations.

Resources

I refer the reader to my old post, linked at the beginning, for references about the logical details.

• Learn You a Haskell for Great Good! This is a great introductory book that can be read for free online.
• Haskell Wikibook
• Some notes on OCaml vs. Haskell.

Introduction

Gödel's β-function and the corresponding lemma perform a simple task: given a sequence of natural numbers, we can construct a function β that will extract the -th number in the sequence.

In essence, the β-function acts like OCaml's List.nth function:

let beta_ocaml (seq:'a list) (i:int) = List.nth seq i;;

Or generic array/list access:

a[i];

So what's all the fuss? Note that the above function takes a sequence of integers as its input. Gödel was of course interested in first-order arithmetic and so could work only with numbers, not sequences. While the construction of the β-function does not directly use the standard arithmetization of syntax, an appropriate number to represent the sequence must still be found.

Still, you might notice that we have -place functional symbols in arithmetic and so can define projection functions. The motivation for defining β is more subtle: using this function, one can represent exponentiation in Robinson's arithmetic (Peano minus induction) without exponentiation axioms. That is to say, using β allows exponentiation to be representable in which is a finitely axiomatizable fragment of first-order arithmetic. Because exponentiated can be represented in this fragment, it is a very weak theory that is still incomplete and undecidable by the same argument as Gödel's original.  The details of this representability will be omitted from this exposition, but resources will be mentioned at the end.

First,we'll take a look at the mathematics behind the function. This will be followed by implementations in OCaml and Python and then some remarks  on the requirements for the computation, which language was better suited, and how actual programming languages access lists.

The Mathematics

Gödel's β-function lemma. For any sequence of natural numbers , there exists a natural number and a recursive function such that for every .

The proof of this lemma generally uses an existence claim based on an ingenious use of the Chinese Remainder Theorem. Due to this existence claim and other number-theoretic trickery, it is hard to see from the proof (which I will give a sketch of shortly) how exactly this function works.  My roommate at Carnegie-Mellon's summer school in Logic and Formal Epistemology bemoaned this fact, claiming that he just wants to see it work for any sequence, even as simple as (more on this sequence later).

Due to this obtuseness, I decided to set out to write a simple program to implement the β-function. After all, it is a computable function.  Such a program would, in theory, shed some light on the inner workings of this function.  I will share the code and make some comments on language choices in a bit.  First, a sketch of the proof of the β-function lemma.

Proof sketch. First, we take a detour and show that and such that . We claim that (i.e. the remainder in  fits the bill. Let ; then it can be shown that are relatively prime in pairs. From this, the Chinese remainder theorem tells us that .

We now define in terms of . Let be the pairing function from . Because this is a bijection, we can numerically define projection functions such that and . Now, define . By construction, we have that , which completes the proof.

Implementations

First, the OCaml version. I wrote this implementation with the goal in mind of keeping as much of the spirit of the above proof in tact. Thus, a few helper functions are required:

(* This is Cantor's classic function from N*N -> N *)
let pairing ((a:int), (b:int)) = int_of_float (0.5 *. float_of_int ((a + b) * (a+b) + 3*a + b));;
 
(* Inverse of it; after all, Cantor showed that N*N and N are equinumerous,
i.e. the pairing function is 1-1 and onto *)
let inverse_pairing (z:int) =
	let w = int_of_float (floor (((sqrt (float_of_int (8*z + 1))) -. 1.0) /. 2.0))
	in let t = int_of_float (((float_of_int (w*w + w)) /. 2.0))
		in let x = z - t in
		(x,w-x);;
 
(* While I know that these projection functions will only be used on integers,
I'm defining them with polymorphic types for maximum generality. *)
let proj1 (x,y) = x;;
let proj2 (x,y) = y;;
 
(* Generic factorial algorithm *)
let rec factorial (n:int) =
	if n == 0 then 1 else n * factorial (n-1);;
 
(* List a ... b *)
let rec range a b =
	if a > b then []
	else a :: range (a+1) b;;
 
(* Calculates the maximum value of an integer list *)
let max (intlist:int list) = List.fold_left (function x -> function y -> if x > y then x else y) 0 intlist;;
 
(* product of a list of integers *)
let product (intlist:int list) = List.fold_left (function x -> function y -> x * y) 1 intlist;;
 
(* a = pairing (c, d), a la Enderton *)
let d (seq:int list) = let newseq = (List.length seq)::seq in
	let s = max newseq
	in (factorial s);;

Now with the preliminaries out of the way, we can get to the real β function. Even though OCaml has fantastic type inference, I've included many annotations to make it clear what the arguments to each function are. I've defined a general beta function to take an integer sequence and the index as its input (this is also where the bulk of the Chinese remainder theorem is implemented); it then calls a beta_godel function which resembles the actual β defined above.

let beta_godel ((a:int), (i:int)) = let m = (inverse_pairing a) in
	(proj1 m) mod ((i + 1) * (proj2 m) + 1);;
 
let beta (seq:int list) (i:int) =
	let d1 = (d seq) in
	let dlist = (List.map (function i -> 1+ (i+1)*d1) (range 0 ((List.length seq) - 1))) in
	let p = product dlist in
	let c = List.find (function x -> (remainder_list x dlist) = seq) (range 0 (p-1))
	in beta_godel ((pairing (c,d1)), i);;

In action:

# beta [1;2;3] 2;;
- : int = 3

There's a reason I chose a list as simple as [1; 2; 3], but more on that later. First, a quick implementation in Python. In some ways, this implementation is easier to read since it's only one function. In other ways, it's harder to read as it lacks typing annotations and makes frequent use of Python's lambda syntax; while this syntax is very appealing to logicians, it is rarely used by most programmers.

import math
 
def beta(seq, idx):
        seq.append(len(seq));
        d = math.factorial(max(seq));
        seq.pop(); #because of append in first line
        dlist = map(lambda i : 1 + (i+1)*d, range(0,len(seq)));
        p = reduce(lambda x,y : x*y, dlist);
        for c in range(0,p-1):
                remlist = map(lambda x : c % x, dlist);
                if remlist == seq:
                        return c % (dlist.pop(idx))

>>> godel_beta.beta([1,2,3],2)
3

Computational Limits

Recall that I omitted details of the proof of the Chinese remainder theorem and used it to prove an existence claim. To actually compute the relevant natural number--the one whose remainder when divided by --requires checking the range where . This product gets very large, very quickly. For the list [1, 2, 3], it is 1729. For [1,2,3,1], it is 8674225. Compare these results to the maximum integers in OCaml and Python (my machine being a Lenovo X60 tablet with 1.5Ghz Core 2 Duo):

# max_int;;
- : int = 1073741823

>>> import sys
>>> sys.maxint
2147483647

And if the numbers get any larger, the integer exceeds these relevant bounds. In the case of the sequence [1, 3, 2, 5]:

>>> godel_beta.beta([1,3,2,5],2)
5063545201
Traceback (most recent call last):
  File "", line 1, in
  File "godel_beta.py", line 11, in beta
    for c in range(0,p-1):
OverflowError: range() result has too many items

Most programming languages obviously do not do such arithmetic to access lists. They rather (as any C programmer will know) use pointer arithmetic to exploit the linear manner in which memory is stored on a chip. If you have the address for the location in memory of an element of a sequence, then the subsequent elements are located literally sequentially in memory. While this is a conjecture (which could be verified with more time), I imagine that implementations of functional languages like LISP and OCaml exploit the nature of memory in this way even if the syntax of the language does not reflect it.

What's Next

I have a few plans for how to extend this work. In particular, I could and/or will:

1. Do a more thorough analysis of the number to get fully clear on the computational limits.
2. Implement the β function in Java, using the BigInteger class, to circumvent these difficulties.
3. Write a note about some differences between OCaml and Python learned from this experience.
4. Write on representable relations, functions and show how exponentiation can be represented in a weak theory of arithmetic using β.

Resources

For more thorough expositions and other interest, I recommend:

1. Herbert B. Enderton, A Mathematical Introduction to Logic (2001), section 3.8 "Representing Exponentiation", pp. 276-281
2. Wikipedia: Gödel's β-function, Chinese Remainder Theorem, Computable function and related articles.
3. George Boolos, Richard Jeffrey and John Burgess, Computability and Logic.

Under supervision of Dr. J T Ratnanather, I have just been awarded a Provost Undergraduate Research Award for a project entitled "A New Approach to Biomedical Ontologies."

The basic idea is to develop a formal ontology in OWL that will properly encode metadata derived from the Large Deformation Diffeomorphic Metric Mapping algorithm in computational anatomy.  In this framework, an anatomy is the orbit of a diffeomorphism group acting on a subset of 3-dimensional space.  Two images (i.e. of the same anatomical structure in different patients) are then 'mapped onto' each other by computing the geodesic flow through the diffeomorphism group.  This also gives the notion of an initial momentum through the group and a metric distance between two anatomical structure.

By encoding this sort of `mechanical' metadata about anatomical structures in a semantically rich ontology, one hopes that new ways of classifying disease will arise. The core concept is that differences in the shape of certain structures can be correlated with disease.  While this project will test these notions on data already derived from the heart, the basic ideas line up with my other project on Planum Temporale Shape Analysis.

I'll post more information about this project on my research page as work begins at the end of this semester.

I've started doing some problems from Project Euler both as useful algorithms exercises and because I want to get familiar with Python.  I picked problem 22 because it hadn't been solved as much as some of the other problems.  That said, the problem itself is much more straightforward than I expected.

My solution is shown below; the corresponding names.txt file can be downloaded from the problem site linked above.

f = open("names.txt");
names = f.readline().split(',');
names.sort(key=str.lower);

# Loop 1
# This is an explicit, procedural way to do the problem
count = 0;
total = 0;
for name in names:
	name = name.strip('\"');
	word_count = 0;
	for letter in name:
		word_count += (ord(letter) - ord('A')) + 1
	count += 1
	total += word_count * count
print total

total = 0;
# Loop 2
# Same method as above, but slightly more concise
for index in range(1,len(names)+1):
	word_count = 0;
	for letter in names[index-1].strip('\"'):
		word_count += (ord(letter) - ord('A')) + 1
	total += word_count * index
print total