Programming I, 2018-9
Lecture 13 – Notes

sieve of Eratosthenes

We've seen in previous lectures that we can determine whether an integer n is prime using trial division, in which we attempt to divide n by successive integers. Because we must only check integers up to `sqrt(n)`, this primality test runs in time O(`sqrt n`).

Sometimes we may wish to generate all prime numbers up to some limit N. If we use trial division on each candidate, then we can find all these primes in time `O(N sqrt N`). But there is a faster way, using a classic algorithm called the Sieve of Eratosthenes.

It's not hard to carry out the Sieve of Eratosthenes using pencil and paper. It works as follows. First we write down all the integers from 2 to N in succession. We then mark the first integer on the left (2) as prime. Then we cross out all the multiples of 2. Now we mark the next unmarked integer on the left (3) as prime and cross out all the multiples of 3. We can now mark the next unmarked integer (5) as prime and cross out its multiples, and so on. Just as with trial division, we may stop when we reach an integer that is as large as `sqrt N`.

Here is the result of using the Sieve of Eratosthenes to generate all primes between 2 and 30:

23456789101112131415161718192021222324252627282930

We can easily implement the Sieve of Eratosthenes in Pascal:

{$mode delphi}

var
  n, i, j: integer;
  prime: array of boolean;

begin
  write('Enter n: ');
  readln(n);
  setLength(prime, n + 1);
  for i := 2 to n do
    prime[i] := true;
  
  for i := 2 to trunc(sqrt(n)) do
    if prime[i] then
      begin
        j := 2 * i;
        while j <= n do
          begin
            prime[j] := false;
            j := j + i;
          end;
      end;
  
  for i := 2 to n do
    if prime[i] then write(i, ' ');
  writeln;
end.

How long does the Sieve of Eratosthenes take to run?

The inner loop, in which we set prime[j] := false, runs N/2 times when we cross out all the multiples of 2, then N/3 times when we cross out multiples of 3, and so on. So its total number of iterations will be

`N(1/2 + 1/3 + 1/5 + ... + 1/p)`

where p <= `sqrt N`.

The series

`sum_(p prime) 1/p = 1/2 + 1/3 + 1/5 + …`

is called the prime harmonic series. How can we approximate its sum through a given element 1/p?

First, we may note that the elements of the prime harmonic series are a subset of elements of the ordinary harmonic series, which sums the reciprocals of all positive integers:

`1/1 + 1/2 + 1/3 + … `

The ordinary harmonic series diverges, and the sum of all its terms through 1/n is close to ln n. In fact the difference between this sum and ln n converges to a constant:

`lim_(n->oo) [sum_(k=1)^n 1/k - ln n] = 0.577...`

So since p <= `sqrt N` we have

`N(1/2 + 1/3 + 1/5 + ... + 1/p) < N(1/2 + 1/3 + 1/4 + .. + 1/p) = N * O(log sqrt N) = N * O(log N) = O(N log N)`

We can derive a tighter bound using properties of the prime harmonic series itself. Surprisingly, the prime harmonic series diverges (as was first demonstrated by Euler). It grows extremely slowly, and its partial sum through 1/p is close to ln (ln n):

`lim_(n->oo) [sum_(p<=n) 1/p - ln (ln n)] = 0.261...`

This shows in turn that

`N(1/2 + 1/3 + 1/5 + ... + 1/p) = N * O(log log (sqrt N)) = N * O(log log N) = O(N log log N)`

This is very close to O(N). (In fact more advanced algorithms can generate all primes through N in time O(N).)

priority queues

In recent lectures we learned about stacks and queues, which are abstract data types that we can implement in various ways, such as using an array or a linked list.

A priority queue is another abstract data type. Like a stack and an ordinary queue, it has one function for adding an element to a priority queue, and another function for removing an element. Here is an interface for a priority queue that holds integers:

type
  pqueue = …   // priority queue

procedure init(var q: pqueue);
procedure insert(var q: pqueue; i: integer);
function removeLargest(var q: pqueue): integer;

function isEmpty(q: pqueue): boolean;

A priority queue differs from a stack and an ordinary queue in the order in which elements are removed. A stack is last in first out: the pop function removes the element that was added most recently. An ordinary queue is first in first out: the dequeue function removes the element that was added least recently. In a priority queue, the removeLargest function removes the element with the largest value.

The interface above describes a max-queue, in which we can efficiently remove the largest value. Alternatively we can build a min-queue, which has a removeSmallest element that removes the smallest value; this is more convenient for some applications. (Generally any data structure that implements a max-queue can be trivially modified to produce a min-queue, by changing the direction of element comparisons.)

In theory we could implement a priority queue using a binary search tree. If we did so and the tree was balanced, then insert and removeLargest would run in time O(log N), where N is the number of elements in the tree. But there are more efficient data structures for implementing priority queues, such as binary heaps, to be discussed next.

binary heaps

A binary heap is a binary tree that satisfies two properties.

If a node with value p has a child with value c, then p >= c.
All levels of the tree are complete except possibly the last level, which may be missing some nodes on the right side.

For example, here is a binary heap:

tree

This heap looks like a complete binary tree with three nodes missing on the right in the last level.

The height of a binary heap with N nodes is floor(log₂(N)) = O(log N). In other words, a binary heap is always balanced.

Typically we don’t store a binary heap using dynamically allocated nodes and pointers. Instead, we use an array, which is possible because of the shape of the binary heap tree structure. The binary heap above can be stored in an array a like this:

Notice that the array values are in the same order in which they appear in the tree, reading across tree levels from top to bottom and left to right.

Here is a tree showing the indices at which heap nodes are stored:

tree
From this tree we can see the following patterns. If a heap node N has index i, then

N’s left child (if any) has index i * 2 + 1
N’s right child (if any) has index i * 2 + 2
N’s parent (if any) has index (i - 1) div 2

So we can easily move between related tree nodes by performing index arithmetic.

In Pascal, we can store a heap like this:

type
  heap = record
    a: array of integer;  // elements a[0 .. n - 1] store the heap
    n: integer;
  end;

procedure init(var h: heap);
begin
  setLength(h.a, 1);
  h.n := 0;
end;

We store the heap size in an integer field n, which will sometimes be less than the size of the underlying array a.

Here are functions for moving through a heap using index arithmetic:

function parent(i: integer): integer;
begin
  parent := (i - 1) div 2;
end;

function leftChild(i: integer): integer;
begin
  leftChild := 2 * i + 1;
end;

function rightChild(i: integer): integer;
begin
  rightChild := 2 * i + 2;
end;

heap operations

We will now describe operations that will let us use a heap as a priority queue.

Suppose that a heap structure satisfies the heap properties, except that one node N has a value v which is smaller than one or both of its children. We can restore the heap properties by performing a down heap operation, in which the value v moves downward in the tree to an acceptable position. Let v₁ be the value of the largest of N’s children. We swap v with v₁. Now N’s value is v₁, which restores the heap property for this node, since v₁ > v and v₁ is also greater than or equal to the other child node (if there is one). We then continue this process, swapping v downward as many times as necessary until v reaches a point where it has no smaller children. The process is guaranteed to terminate successfully, since if v eventually descends to a leaf node there will be no children and the condition will be satisfied.

Here is Pascal code that implements the down heap operation:

procedure swap(var i, j: integer);
var
  k: integer;
begin
  k := i;
  i := j;
  j := k;
end;

// Move the value at index i downward until its children are smaller than it
procedure down(var h: heap; i: integer);
var
  l, r, largest: integer;
begin
  l := leftChild(i);
  r := rightChild(i);
  largest := i;
  if (l < h.n) and (h.a[l] > h.a[i]) then
    largest := l;
  if (r < h.n) and (h.a[r] > h.a[largest]) then
    largest := r;
  if largest <> i then  // some child is larger
    begin
      swap(h.a[i], h.a[largest]);
      down(h, largest);
    end;
end;

Now suppose that a heap structure satisfies the heap properties, except that one node has a value v which is larger than its parent. A complementary operation called up heap can move v upward to restore the heap properties. Suppose that v’s parent has value v₁. That parent may have a second child with value v₂. We begin by swapping v and its parent v₁. Now v’s children are v₁ and v₂ (if present). We know that v > v₁. If v₂ is present then v₁ > v₂, so v > v₂ by transitivity. Thus the heap properties have been restored for the node that now contains v. And now we continue this process, swapping v upward in the heap until it reaches a position where it is less than its parent, or until it reaches the root of the tree.

(You will implement the up heap operation for a homework exercise.)

We can use the up heap and down heap operations to implement priority queue operations.

We first consider removing the largest value from a heap. Since the largest value in a heap is always in the root node, we must place some other value there. So we take the value at the end of the heap array and place it in the root, decreasing the array size by 1 as we do so. We now perform a down heap operation on this value, which will lower it to a valid position in the tree. If there are N values in the heap, the tree height is floor(log₂(N)), so this process is guaranteed to run in O(log N) time.

Here is removeLargest in Pascal:

function removeLargest(var h: heap): integer;
begin
  removeLargest := h.a[0];
  h.a[0] := h.a[h.n - 1];
  h.n := h.n - 1;
  down(h, 0);
end;

Now we consider inserting a value v into a heap. To do this, we first add v to the end of the heap array, expanding the array size by 1. Now v is in a new leaf node at the end of the last heap level. We next perform an up heap operation on the value v, which will bring it upward to a valid position. Just like removeLargest, insert will always run in time O(log N).

heapsort

We’ve already seen several sorting algorithms in this course: bubble sort, insertion sort, merge sort and quicksort. We can use a binary heap to implement another sorting algorithm: heapsort. The essence of heapsort is simple: given a array of numbers to sort, we put them into a heap, and then we pull out the numbers one by one in sorted order.

Let’s examine these steps in a bit more detail. Given an array of N values, we first need to heapify it, i.e. turn it into a heap. One way to do this is to simply insert each number in turn using the insert operation described above. This will time O(N log N).

But there is actually a more efficient way to heapify the array: we walk backward through the array, running the down heap operation on each value in turn. Here is Pascal code to do that:

type
  intarray = array of integer;

procedure heapify(out h: heap; var a: intarray);
var
  i: integer;
begin
  h.a := a;
  h.n := length(a);
  
  for i := h.n - 1 downto 0 do
    down(h, i);
end;

How long will heapify take to run? Suppose that the heap has N nodes. There are approximately (N / 2) leaves. A down heap operation on a leaf does nothing since the leaf is at the bottom. There are approximately (N / 4) nodes immediately above the leaves. A down heap operation on each of these nodes may shift it downward 1 position in the tree. There are approximately (N / 8) nodes above those, and a down heap operation on those nodes may shift them downward up to 2 positions. The total number of downward shifts is at most approximately

(N / 4) + 2 (N / 8) + 3 (N / 16) + ... + (log₂ N)(1)

= (N / 4) + (N / 8 ) + (N / 16) + (N / 32) + ...
+ (N / 8) + (N / 16) + (N / 32) + …
+ (N / 16) + (N / 32) + …
+ (N / 32) + ...

= (N / 2) + (N / 4) + (N / 8) + ...

= N

This shows, informally, that heapify will run in O(N) (a somewhat surprising fact). We can turn an unordered array into a binary heap very efficiently.

Once we have heapified the array, sorting it is easy: we repeatedly call removeLargest to remove the largest element, which we place at the end of the array. Here is our Pascal implementation:

procedure heapsort(var a: intarray);
var
  h: heap;
  i: integer;
begin
  heapify(h, a);

  while h.n > 0 do
    begin
      i := h.n;
      h.a[i] := removeLargest(h);
    end;
end;

The total time for heapsort to sort N values is the time to heapify the array plus the time to remove each of the N values from the heap in turn. This is

O(N) + O(N log N) = O(N log N)

Thus heapsort has the same asymptotic running time as merge sort and quicksort. Some advantages of heapsort are that it runs on an array in place (unlike merge sort) and that it is guaranteed to run in time O(N log N) for any input (unlike quicksort).

grammars for arithmetic expressions

In the last lecture we learned about prefix, infix and postfix syntaxes for arithmetic expressions. Continuing that discussion, we will now formally define these forms of arithmetic expressions. To keep things simple, we will work with expressions that contain only single-digit numbers.

We can define the syntax of infix arithmetic expressions using the following context-free grammar:

digit → '0' | '1' | '2' | '3' | '4' | '5' | '6' | '7' | '8' | '9'
op → '+' | '-' | '*' | '/'
expr → digit | '(' expr op expr ')'

Context-free grammars are commonly used to define programming languages. A grammar contains non-terminal symbols (such as digit, op, expr) and terminal symbols (such as '0', '1', '2', …, '+', '-', '*', '/'), which are characters. A grammar contains a set of production rules that define how each non-terminal symbol can be constructed from other symbols. These rules collectively define which strings of terminal symbols are syntactically valid.

For example, ((4 + 5) * 9) is a valid expression in the syntax defined by this grammar, because it can be constructed from the top-level non-terminal expr by successively replacing symbols using the production rules:

expr
→ ( expr op expr )
→ ( expr op digit )
→ (expr op 9)
→ (expr* 9)
→ ((expr op expr) * 9)
→ ((expr op digit) * 9)
→ ((expr op 5) * 9)
→ ((expr + 5) * 9)
→ ((digit + 5) * 9)
→ ((4 + 5) * 9)

You will learn much more about context-free grammars in more advanced courses, but it’s worth taking this first look at them now because they are so commonly used to specify languages.

To modify our language to use prefix operators, we need change only the last production rule above:

expr → op expr expr

Or, for postfix operators:

expr → expr expr op

Note again that the prefix and postfix expression languages have no parentheses.

evaluating arithmetic expressions

It is not difficult to write Pascal functions that evaluate arithmetic expressions, i.e. compute the numeric value that each expression represents. Because arithmetic expressions have a recursive structure, these functions will often be recursive as well.

First, here are some useful helper functions including evalOp, which applies an operator to two integer arguments:

uses character, sysutils;

procedure error(s: string);
begin
  writeln(s);
  halt;
end;

function readChar(): char;
var
  c: char;
begin
  if seekEoln then error('unexpected end of input');
  read(c);
  exit(c);
end;

function isOp(c: char): boolean;
begin

  exit(pos(c, '+-*/') > 0);
end;

function evalOp(c: char; i, j: integer): integer;
begin
  case c of
    '+': exit(i + j);
    '-': exit(i - j);
    '*': exit(i * j);
    '/': exit(i div j);
  end;
end;

We can evaluate a prefix expression using a recursive function. Its structure reflects the grammar for prefix expressions defined in the previous section.

function evalPrefix(): integer;
var
  c: char;
  i, j: integer;
begin
  c := readChar();
  if isDigit(c) then exit(strToInt(c));   // expr → digit
  if not isOp(c) then error('operator expected');
  i := evalPrefix();
  j := evalPrefix();
  exit(evalOp(c, i, j));    // expr → op expr expr
end;

We can evaluate an infix expression similarly:

function evalInfix(): integer;
var
  c: char;
  i, j: integer;
begin
  c := readChar();
  if isDigit(c) then exit(strToInt(c));   // expr → digit
  if c <> '(' then error('expected (');
  
  i := evalInfix();
  c := readChar();
  if not isOp(c) then error('operator expected');
  j := evalInfix();
  
  if readChar() <> ')' then error('expected )');
  exit(evalOp(c, i, j));  // expr → '(' expr op expr ')’
end;

The most straightforward way to evaluate a postfix expression is using a stack. As we read a postfix expression from left to right, each time we see a number we push it onto the stack. When we see an operator, we pop two numbers from the stack, apply the operator and then push the result back onto the stack. For example, in evaluating the postfix expression 4 5 + 2 1 + * we would perform the following operations:

push 4
push 5
pop 4 and 5, push 9
push 2
push 1
pop 2 and 1, push 3
pop 9 and 3, push 27

When we finish reading the expression, if it is valid then the stack will contain a single number, which is the result of our expression evaluation.

(You will write code to evaluate a postfix expression using a stack in an exercise this week.)

Programming I, 2018-9 Lecture 13 – Notes