SICP: The 8 queens puzzle

The Google AI Contest is over and now I continue to read the SICP book. Now we are solving the famous 8 queens puzzle. One way to solve the puzzle is to work across the board, placing a queen in each column. Once we have placed k - 1 queens, we must place the k^th queen in a position where it does not check any of the queens already on the board. We can formulate this approach recursively: Assume that we have already generated the sequence of all possible ways to place k - 1 queens in the first k - 1 columns of the board. For each of these ways, generate an extended set of positions by placing a queen in each row of the kth column. Now filter these, keeping only the positions for which the queen in the kth column is safe with respect to the other queens. This produces the sequence of all ways to place k queens in the first k columns. By continuing this process, we will produce not only one solution, but all solutions to the puzzle. We implement this solution as a procedure queens, which returns a sequence of all solutions to the problem of placing n queens on an n×n chessboard.

In the beginning we are given this procedure.

(define (queens board-size)
  (define (queen-cols k)
    (if (= k 0)
        (list empty-board)
        (filter
         (lambda (positions) (safe? k positions))
         (flatmap
          (lambda (rest-of-queens)
            (map (lambda (new-row)
                   (adjoin-position new-row k rest-of-queens))
                 (enumerate-interval 1 board-size)))
          (queen-cols (- k 1))))))
  (queen-cols board-size))

and we need to write all the sub-procedures that are used by the main one. Let's start from the simplest one.

(define empty-board '())

Next task is to write the safe? function. This should determine for a set of positions, whether the queen in the k^th column is safe with respect to the others. I decided to split it into three procedures. 1 checks the horizontal line, one diagonal up and the other diagonal down. We don't actually need the k argument, as we should always check the last queen only. I use reversed-positions as we are starting from the last queen and going to the first one.

(define (safe? k positions)
  (let ((reversed-positions (reverse positions))
        (last (car (reverse positions))))
  (and (horizontal-safe? last (cdr reversed-positions))
       (diagonal-up-safe? last (cdr reversed-positions))
       (diagonal-down-safe? last (cdr reversed-positions)))))

(define (horizontal-safe? q positions)
  (if (null? positions)
      true
      (and (not (= (car positions) q))
           (horizontal-safe? q (cdr positions)))))

(define (diagonal-up-safe? q reversed-positions)
  (if (null? reversed-positions)
      true
      (and (not(= (car reversed-positions) (+ q 1)))
           (diagonal-up-safe? (+ q 1) (cdr reversed-positions)))))

(define (diagonal-down-safe? q reversed-positions)
  (if (null? reversed-positions)
      true
      (and (not(= (car reversed-positions) (- q 1)))
           (diagonal-down-safe? (- q 1) (cdr reversed-positions)))))

The remaining procedures are much easier to write. flatmap should join all the lists, that represent queens on board, into one big list. enumerate-interval just returns a list with values between low and high. adjoin-position just adds new queen to the existing ones.

(define (flatmap proc lst)
  (foldr append '() (map proc lst)))

(define (enumerate-interval low high)
  (if (> low high)
      '()
      (append (list low) (enumerate-interval (+ low 1) high))))

(define (adjoin-position new-row k rest-of-queens)
  (append rest-of-queens (list new-row)))

As you can see, there are some parameters in procedures that are not used. They can be safely removed.

(queens 8)

Returns 92 solutions - which is right.

Determine If Number Is Prime

The SICP book is still being read...

Now I am going to introduce some algorithms that allow us to say whether the number is prime or not. Again we will start from the simplest one.

Take numbers 1 by 1 starting from 2 to √n and check, if it divides n. The complexity of the algorithm is Θ(√n).

(define (prime? n)
  (define (divides? k)
    (= (remainder n k) 0))
  (define (find-divisor n test)
    (cond ((> (sqr test) n) n)
          ((divides? test) test)
          (else (find-divisor n (+ test 1)))))

  (= n (find-divisor n 2)))

Fermat Primality Test

There is a test for primality with Θ(log(n)) complexity. It is based on Fermat little theorem.

an ≡ a (mod n), where 1 ≤ a < n

Which means, that if n is prime, then for any integer a, an − a will be evenly divisible by n.

;this procedure calculates baseexp mod m without using big numbers
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (remainder (sqr (expmod base (/ exp 2) m)) m))
        (else (remainder (* base (expmod base (- exp 1) m)) m))))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

Unfortunately there are some issues with this test - it is probabilistic. If the test fails, n is certainly composite. But if it passes, it does not guarantee that n is prime, although it shows high probability for n to be. So we can run as many tests as we need to make the mistake probability as low as possible.

(define (fermat-prime? n times)
  (cond ((= times 0) true)
        ((fermat-test n) (fermat-prime? n (- times 1)))
        (else false)))

Unfortunately that is not all. There are some numbers that pass the test (for every a!), without being prime. These are called Carmichael numbers. So to determine, if n is prime, we need to check if it passes the Fermat test and is not a Carmichael number, or use a similar Miller–Rabin primality test, that cannot be tricked.

Miller–Rabin Primality Test

It starts from an alternate form of Fermat's Little Theorem:

an-1 ≡ 1 (mod n), where 1 ≤ a < n

The difference is that whenever we perform the squaring step in expmod, we check to see if we have discovered a "nontrivial square root of 1 modulo p", that is, a number not equal to 1 or n - 1 whose square is equal to 1 mod n. It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing an-1 in this way will reveal a nontrivial square root of 1 mod n.

(define (miller-rabin-prime? n times)
  (cond ((= times 0) true)
        ((miller-rabin-test n) (miller-rabin-prime? n (- times 1)))
        (else false)))

(define (miller-rabin-test n)
  (define (try-it a)
    (= (expmod a (- n 1) n) 1))
  (try-it (+ 1 (random (- n 1)))))

(define (expmod base exp m)
  (define (!= x y) (not (= x y)))

  (cond ((= exp 0) 1)
        ((and (!= base 1) 
              (!= base (- m 1)) 
              (= (remainder( sqr base) m) 1)) 0)
        ((even? exp) (remainder (sqr (expmod base (/ exp 2) m)) m))
        (else (remainder (* base (expmod base (- exp 1) m)) m))))

Fibonacci Numbers

Continuing reading the SICP book...

In every book about algorithms I have read there is always a chapter about the Fibonacci numbers. Mainly on the different algorithms of calculating the n-th Fibonacci number and comparing their order of growth. First comes the algorithm, based on definition:

(define (fib n)
  (cond ((= n 0) 0)
        ((= n 1) 1)
        (else (+ (fib (- n 1))
                 (fib (- n 2))))))

This function is good at explaining tree recursion, but is awful for calculating Fibonacci numbers, as it makes too many unnecessary calculations. The order of growth is Θ(Fib(n)), which means exponential. The next comes iterative funcion, which solves the task much better, using less resources and time. Its order of growth is Θ(n) - linear.

(define (fib n)
  (fib-iter 1 0 n))

(define (fib-iter a b count)
  (if (= count 0)
     b
     (fib-iter (+ a b) a (- count 1))))

But there is also an algorithm that calculates Fibonacci number with Θ(log n) complexity. There, in the book they give you the function with gaps and some pointers ( no C here :) ), how to fill them.

(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   [??] calculate p'
                   [??] calculate q'
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

The authors remind us of the transformation of the state variables a and b in the previous (iterative) fib-iter procedure: a ← a + b and b ← a. Provided these state changes are labeled transformation T, applying T repeatedly for n iterations starting with a = 1 and b = 0 will produce the pair a = Fib(n + 1) and b = Fib(n). So the Fibonacci numbers are produced by the n^th power of the transformation T (Tⁿ), starting with the pair (1, 0).

Now consider the family of transformations T_pq which transforms the pair (a, b) according to the following rules:

a ← bq + aq + ap
b ← bp + aq

Where transformation T is just a special case of T_pq, where p = 0 and q = 1. Now it's up to you to find the transformation T_p'q' such so, if we apply T_pq twice, the effect is the same as using a single transformation T_p'q' of the same form. Compute p' and q' in terms of p and q and put into the function.

I will not provide the solution here, in case you want to handle it yourself. But if you are eager to see the result, visit some other blog that has the solution. :)

Structure and Interpretation of Computer Programs

Recently we decided to read the "Structure and Interpretation of Computer Programs" book through. And not just read it, but also to complete the exercises, given after each chapter. On Friday we had a first meeting, concerning the first chapters we read.

As a working environment the Racket was chosen. It turned out to be a very powerful instrument: not just the language and it's interpretor, but almost a complete IDE with autocomplete, unit tests and even debugging.

That is how errors are shown in Racket. Yes, these arrows are cute.

The meeting was nice and quite productive. We shared some solutions to exercises, and also some were solved just in place. We learned that

(A 2 6)

is a very big number :), where

(define (A x y)
  (cond ((= y 0) 0)
        ((= x 0) (* 2 y))
        ((= y 1) 2)
        (else (A (- x 1)
                 (A x (- y 1))))))

P.S. Thanks Anton for the picture :)