Computable Numbers in Haskell

Numbers are all

In this post, we examine two examples set out in the classic paper by Turing. The first problem is to report the $n^{th}$ binary digit of 13. The second problem is to report the $n^{th}$ binary digit of $\sqrt{2}$. The 1936 paper On computable numbers, with an application to the Entscheidungsproble describes two machines to perform the tasks.

In this post, we give in Haskell two programs for the same. Both the numbers are computable as there are programs to compute the digit in any position. However, this computation (of reporting each of the numbers) is never ending though. What is more important is that each digit can be computed in a finite amount of time.

The program below computes the $n^{th}$ binary digit of a third. The program prints “010101…” as the binary digits of a third after the decimal point.

athird  n
     | mod n 2 == 1 = "0"
     | otherwise = "1"

The program above is correct because the number $S$ that it computes is given by

$$S = 1/{2^2} + 1/{2^4} + 1/{2^6} + \ldots$$ $$4*S = 1 + S$$

If the index is specified in binary as a string, then by examining the last bit we can decide whether the index is odd or even, which is a constant time operation. Above we use the mod function. The program is fairly efficient. In $O(1)$ time it prints the right digit as expected.

Is it Normal

Next, we tackle the problem of computing an arbitrary digit of $\sqrt{2}$. To make it a little bit of fun, we treat the binary numbers as strings, and perform their addition and multiplication recursively, without using any builtin arithmetic operators.

Function addc below implements single-digit addition with carry using pattern matching. The result is a pair where the first element is the result, and the second element is the resulting carry.

addc :: Char -> Char -> String -> (String, String)
addc '0' '0' "0" = ("0","0")
addc '0' '1' "0" = ("1","0")
addc '1' '0' "0" = ("1","0")
addc '1' '1' "0" = ("0","1")
addc '0' '0' "1" = ("1","0")
addc '0' '1' "1" = ("0","1")
addc '1' '0' "1" = ("0","1")
addc '1' '1' "1" = ("1","1")

The function add below implements the paper and pencil method of adding two numbers, bit by bit, from right to left. We assume that the smaller of the two numbers (length wise) is the first argument. The third argument is the carry forward. The sum is reported as a string.

add :: String -> String -> String -> String
add [] [] c = c
add a b c    -- length of a is at most the length of b    
   | length(a) > length(b) = add b a c
    | otherwise =     (add (take n an) (take n b) c1) ++ r
    where 
        (r,c1) = addc (last an) (last b) c
        an =  ['0' | x <- [1..d]]  ++ a
        d = (length b) - (length a)
        n = (length b) -1 

We perform a few unit checks.

add "111" "111" "0" 
add "111" "000" "0" 
add "101" "010" "0"
add "111" "11"

Next, we implement multiplication using the pencil and paper method again. The product is easy to compute if the multiplier is either “0” or “1”. These are the base cases listed below. The recursive step is to compute the product of a with the all but the last digits of b and add this result to a. The method shifts the rows to the right. We keep track of the number of zeros to pad using the third argument.

mult :: String -> String -> Integer -> String
mult a "0" d = ['0' | x <- [1..length(a)]]++['0' | x <- [1..d]]
mult a "1" d = a ++ ['0' | x <- [1..d]]
mult a b cnt = dropWhile (\x -> x == '0') $ 
  add (mult a [(last b)] (cnt)) (mult a (take ((length b)-1) b) (cnt+1)) "0"
mult "111" "11" 0            -- 7*3
mult "100" "10" 0           -- 4*2
mult "101" "10" 0           -- 5*2
mult "1" "1" 0           -- sqrt(2) = 1.01101010000010011110
mult "10" "10" 0
mult "101" "101" 0
mult "1011" "1011" 0
mult "10110" "10110" 0

The program works as follows: if “x” is the current guess then the program will compute “x1” * “x1”. If the product is less than $2$ then the guess is changed to “x1”, else the guess is updated to “x0”

guess :: String -> Int -> String    -- returns new guess
guess x n 
   | length(x) >= n = x
    | length(mult g g 0) < 2*length(g) = guess g n
    | otherwise = guess (x ++ "0") n
      where g = x ++ "1"

The first 200 binary digits of $\sqrt{2}$ can be computed as

*Main> guess "1" 200
"1011010100000100111100110011001111111001110111100110
 0100100001000101100101111101100010011011001101110101
 0100101010111110100111110001110101101111011000001011
 10101000100100111011101010000100110011101101"

For a moment assume that the work done is proportional to the number of single-digit multiplications. If the first $n-1$ digits are known then our recursive code takes $O(n^2)$ single-digit multiplications. Therefore, starting from a guess of “1”, it takes $O(n^3)$ single-digit multiplications. We can use a faster algorithm to multiply two $n$ digit binary numbers in $O(n \log{n})$ time. This means first $n$ binary digits of $\sqrt{2}$ can be computed in $O(n^2 \log{n})$ single-digit multiplications.

As an exercise try to speed up the method above using FFTs. An even faster algorithm would be great! As an aside, $\sqrt{2}$ is believed to be normal for any base, but there is no known proof for any base as of writing, so a faster method might be hard to find.

The code is available here.

Related