Recursion: Functions that call themselves

Comments and questions to John Rowe.

God is subtle but he is not malicious.
Einstein

Recursion is an extremely simple concept: a function simply calls itself.

Recursion refers to a function that calls itself either directly or indirectly.

The test to see if a function must call itself

Of course, a function can't always call itself, that would create an infinite loop, so all recursive functions have some sort of test before they call themselves.

void myfun(some_args) {

 /* Do some stuff */
  if (some_test)
     myfun(some_other_args);
}

A slighty more spohisticated alternative is to have a loop for the function to call itself as many times as possible and for that loop to be empty when the function does not need to call itself.

There must always be a test (or a possibly-empty loop) to see if a function must call itself.

A better test

A nice variation on the above is for the function to always call itself and for the (next generation of ) the function to see if it actually needed to be called and to immediately return if not:

void myfun(some_args) {
  if (I_didnt_need_to_be_called)
    return;

    // Do some stuff
    myfun(some_other_args);
}

Where it works this is neater (as we do the test only once).

A function that examines its arguments to see if it doesn't need to do anything is always easier to use than one where the person calling it has to do so.

A function may itself contain the test and immediately return if nothing is required..

Try googling the word "recursion"

The test is the key

Recursion can lead to subtle errors: either something gets treated twice or not at all, or our recursion goes into an infinite loop.

The key to recursion is to understand the "test" inside the recursive function that decides whether or not it calls itself again.

Recursively-called functions have their own local variables

When functions call themselves each instance has its own local variables, unless they are static.

Example: a factorial function

The simplest example of recursion is the following somewhat pretentious factorial function:

long factorial(int i) {
  if ( i > 0 )
    return  i * factorial(i - 1);

  return 1;
}

Step through this code

NB: in the above stepped example, the initial call to factorial() from main() is labelled just "factorial", recursive calls are labelled "factorial-2", "factorial-3", etc.

This illustrates the concepts of recursion:

1. A function called factorial() that calls itself.
2. A test to make sure the sequence of function calls stops at some point
3. A local variable, in this case called "i", which local to each invocation of factorial(): when factorial() calls itself each instance has its own copy of i.

It also illustrates one possible pitfall: it's possible to use recursion just to show off. The above example would have been better written as a simple loop.

Error checking and wrappers

It's always nice if a function does some sanity checking on its arguments and we could add some simple argument error checking as follows:

long factorial(int i) {
  if ( i > 0 )
    return  i * factorial(i - 1);
  else if (i < 0 ) {
    fprintf(stderr, "Factorial must be positive\n");
    return -1;
  }

  return 1;
}

Another option would be to make the publicly-callable function a wrapper around the real function. This avoids the test being done many times over when it only needs to be done once:

// The actual factorial function
long factorial_real(int i) {
  if ( i <= 1 )
    return 1;
  return  i * factorial_real(i - 1);
}

// Publically-callable front-end to the 
// factorial function with error checking
long factorial(int i) {
  if (i < 0 ) { // Error check
    fprintf(stderr, "Factorial must be positive\n");
    return -1;
  }
 return factorial_real(i);
}

Step through this code

This is slightly over the top in this case but is a useful illustration.

When recursion goes wrong

Recursion can be hard to debug. In non-recursive programs we can get a long way by knowing what function we are in and what function called it and with what arguments. But with problems with recursion we know what the function is and it's usually called by itself! The biggest problem tends to be the internal test: either the recursion never stops or some possibilities get missed out.

The biggest source of problems with recursion is the internal test.

Imagine we had made a mistake in our previous rather over-the-top way of calculating a factorial:

long factorial(int i) {

  if ( 1 > 0 ) // Bug!
    return  i * factorial(i - 1);
  return 1;
}

You can see what it's meant to do but we've written '1' instead of 'i' so it will go on for ever. The first stage of debugging is just to print out its arguments:

long factorial(int i) {

  fprintf(stderr, "i: %d\n", i);
  if ( 1 > 0 ) // Bug!
    return  i * factorial(i - 1);
  return 1;
}

}

Try printing out the function's arguments.

Try adding a depth parameter

We can get even more details by adding an extra debugging argument to tell us the depth:

#include <stdlib.h>
#include <assert.h>

long factorial(int i, int depth) {
  int d;

  // Print depth information for debugging
  for(d = 0; d <= depth; ++d)
    fprintf(stderr, "  ");
  fprintf(stderr, "i: %d (depth %d)\n", i, depth);

  // If we have gone too deep abort the program
  assert( depth < 12);

  if ( 1 > 0 ) // Bug!
    return  i * factorial(i - 1, depth + 1);
  return 1;
}

When debugging recursion the depth is always a useful number to know!

Here we have done several things:

• We have added an extra 'depth' argument which we increment by one every time.

• We have indented the output according to the depth. This isn't necessary but it is nice.

• Finally, and rather extremely, if factorial() goes too deep we kill the program using assert(). This is only useful when we are running under the debugger as we can look at the "layers" of recursive call.

assert(expression) is defined inside assert.h and kills the program if expression is false.

The output looks like this:

 
  i: 3 (depth 0)
    i: 2 (depth 1)
      i: 1 (depth 2)
        i: 0 (depth 3)
          i: -1 (depth 4)
            i: -2 (depth 5)
              i: -3 (depth 6)
                i: -4 (depth 7)
                  i: -5 (depth 8)
                    i: -6 (depth 9)
                      i: -7 (depth 10)
                        i: -8 (depth 11)
                          i: -9 (depth 12)
                            i: -10 (depth 13)
recurse-err3.c:13: factorial: Assertion `depth < 12' failed.
Abort (core dumped)

Using a wrapper

long factorial_debug(int i, int depth) {
  // Real factorial code here...
}

long factorial(int i) {
  return  factorial_debug(i, 0);
}

We should be prepared to add some debugging to our recursive routine and also to add a "depth" parameter to the "real" recursive routine, perhaps with a simple publicly-callable wrapper routine.

The order of the actions

Often a recursive function has the choice of "do an action then call itself recursively" or "call itself recursively then do an action":

void myfun(some_args) {
  if (I_didnt_need_to_be_called)
    return;

  // Do some stuff here...
  myfun(some_other_args);
}

void myfun(some_args) {
  if (I_didnt_need_to_be_called)
    return;

  myfun(some_other_args);
  // Do some stuff here...
}

When a function performs an action and calls itself recursively it's usually important to do them in the right order.

Recursive data structures

One of the reasons why recursion is so useful is that we often have to deal with recursive things.

Example: a passive electrical circuit

A simple electrical circuit consists of resistors, capacitors and inductors in series or in parallel. But a more complicated cicuit may consist of sub-circuits in parallel. How do we find the impedance of (a resistor in series with a capacitor) in parallel with (another resistor in series with an inductor) ?

It turns out it's actually very easy. First, consider a simple data structure to represent a single passive electrical component. It has to tell us two things: the type of the component and its value.

Using named constants for the component type

Rather than using strings to represent the component type it's easier to use integer values, such as 0, 1. But we wish to avoid code like:

int type; // 0 is resistor, 1 capacitor, 2 inductor

as it becomes extremely easy to make a mistake. One solution is to use #define:

#define Resistor  0
#define Capacitor 1
#define Inductor  2

We can then use the words "Resistor" etc throughout the code rather than "0".

Enumerations

We cover enumerations more thoroughly in the next lecture.

C provides a more sophisticated way of doing this using enumerations. The following code:

typedef enum { Resistor, Capacitor, Inductor } Type;

The structure definition

We can now define a structure for a component as:

typedef enum { Resistor, Capacitor, Inductor } Type;
typedef struct component {
  double value;
  Type type;
} Component;

The impedance function

We can easily write a simple function to calculate the impedance of a single component at a specific frequency, taking a little care about division by zero for very low frequencies when dealing with capacitors.

Naive capacitor impedance

For a capacitor we might first write the code to return the impedance as:

 return 1.0/(2 * PI * I * c->value * freq);

But this has problems when the frequency is zero. We can get round this with the following code using a suitable value for INFINITY :

    // Handle DC and v low frequencies
    denom = 2 * PI * c->value * freq;
    if ( denom < 1.0/INFINITY )
     denom = 1.0/INFINITY;
    return 1.0/(I * denom);

The code

#include <complex.h>
#define INFINITY 1e100
#define PI 3.14159265358979

double complex impedance(Component *c, double freq) {
  double denom;
  
  switch (c->type) {
  case Resistor:
    return c->value;

  case Inductor:
    return 2 * PI * I * c->value * freq;

  case Capacitor:
    // Handle DC and v low frequencies
    denom = 2 * PI * c->value * freq;
    if ( denom < 1.0/INFINITY )
     denom = 1.0/INFINITY;
    return 1.0/(I * denom);
  }
}

The question is "how do we combine these to handle more-complicated circuits?".

Series and parallel circuits

Given an array of pointers to components arranged in series, the total impedance is just the sum of the individual impedances:

  // Series impedance
  double complex zs = 0;
  for (int i = 0; i < n; ++i)
    zs += impedance(component[i], freq);

The parallel case is almost as simple, it's the reciprocal of the sum of the reciprocals:

  // Parallel impedance
  double complex zp = 0;
  for (int i = 0; i < n; ++i)
    zp += 1.0/impedance(component[i], freq);
  
  zp = 1.0/zp;

This gives us the germ of an idea: given that we can find the impedance of a list of components arranged in series we can think of a series of components as a single, multipart component. The same applies to components arranged in parallel. So we can combine these multi=part components to build more-complicated circuits.

Representing series and parallel circuits

In this situation our instinctive reaction may be to have two types of structures, one for individual components and one for circuits consisting of several components.

typedef enum { Resistor, Capacitor, Inductor } Type;
typedef struct component {
  double value;
  Type type;
} Component;


 

But this becomes very complicated when we realise that a circuit may contain a series of sub-circuits; for example several parallel sub-circuits which are then arranged in series to form a larger circuit. How do we handle this?

But we can make life much simpler for ourselves by one simple, but perhaps slightly counterintuitive, trick: recognising that a "circuit" may consist of just one component. Thus we abolish the distinction between a component and a circuit. We say that a single component is just a circuit with only one thing in it and just merge the two structures together to look like this:

The value of NODE_MAX is a little tricky, but fortunately C gives us a neat way of handing this as described in this optional appendix. For the time being we shall just assume that NODE_MAX is large enough.

This is an example of what we might called the "Generalised Thing", or possible a "Compound Thing".

The "Generalised Thing"

Given a set of simple "Things" (in our case, electrical components) which can be combined together a useful concept is what we might call a "Generalised Thing". A Generalised Thing is either:

1. A single simple Thing.

   Or:

2. A collection of other Generalised Things, sometimes called nodes or children, plus enough information to know how they are combined together.

This simple, and indeed rather strange concept, is remarkably useful. Notice that it is recursive, each of the nodes may have other nodes, etc.

In our case the simple things are resistors, capacitors and inductors, and there are two ways of combining components, series or parallel. So a "generalised component" or circuit conists of either:

• A resistor or
• A capacitor or
• An inductor or
• Two or more sub-circuits arranged in series or
• Two or more sub-circuits arranged in parallel

Example: foods

If you've ever looked at the ingredients of a supermarket sandwich you may well know what we mean!

This is usually referred to as the list of ingredients. In programming the word "list" is often used for a "generalised thing" that is a simple thing or a collection of other things, a generalised thing that consists of a simple thing and possibly a collection of other things is called a tree. We don't get hung up on the difference, we just make sure we can handle either situation.

Another example of a "Generalised Thing" is a mathematical expression: a mathematical expression may contain a sub-expression in parentheses in which case the sub-expression is treated a a whole new top-level expression and may itself contain parentheses:

x * ( 1 + y * ( z - 2 * (a + b ) ) ) + a * ( x + y )

Example circuit

Consider the following circuit:

Ladder filter, image courtesy of Wikipedia (Fcorthay)

What impedance does this present at Uin if Uout is left open-circuit? Starting from the left, the circuit consists of:

R in series with ( C1 in parallel with ( L2 in series with ( C3 in parallel with ( L4 in series with R))))

The data structure

• How many nodes there are.
• Whether the nodes are in series or parallel
• An array of pointers, one for each node

The recursive impedance function

Our recursive impedance function is now just:

#define NODE_MAX 20
 
#define INFINITY 1e100
#define PI 3.14159265358979
double complex impedance(Component *c, double freq) {
  double complex z = 0;
  double denom;

  switch (c->type) {
  case Resistor:
    return c->value;

  case Inductor:
    return 2 * PI * I * c->value * freq;

  case Capacitor:
    // Handle DC and v low frequencies
    denom = 2 * PI * c->value * freq;
    if ( denom < 1.0/INFINITY )
     denom = 1.0/INFINITY;
    return 1.0/(I * denom);

  case Series:
    for (int n = 0; n < c->n; ++n)
      z += impedance(c->nodes[n], freq);
    return z;

  case Parallel:
    for (int n = 0; n < c->n; ++n)
      z += 1.0/impedance(c->nodes[n], freq);
    return 1.0/z;
  }
}

And we can now handle quite complicated circuits. (For an example of a cicuit this can't handle directly consider a pyramid with the corners as nodes and the vertices being the "wires" each with a component, as in this PDF courtesy of Charles Williams.)

See also the complete circuit program.

Directionality and trees

Our "generalised passive component" is an example of a tree. Each component, other than the top one, has precisely one parent and there are no loops. Trees are quite common in programming because the concept they represent, "things that can have other things to an unlimited depth", is quite common in real life. Recursion is an essential tool when dealing with them and they have the advantage that there is no need to deal with potential loops, as there are none, or things being handled multiple times as each node has precisely one parent.

This isn't always the case however: things may have multiple parents (as in the example of real-life family trees) or the concept of parenthood may not even exist, and/or we may have loops.

When data can contain loops we must make sure our recursive functions do not also loop, usually by "marking" visited nodes in some way.

The order in which this is done inside the function is critical and must be:

1. Check to see if node has been visited (or there is another reason to do nothing), and if so return.
2. Mark the node as visted.
3. Do stuff and Recurse.

The correct order: Check, Mark, Recurse

void myfun(some_args) {
  if (node_has_been_visited || I_didnt_need_to_be_called)
    return;                      //Check 

  // Mark this node as visited;  //Mark 

  // Do some stuff here...
  myfun(some_other_args);        //Recurse
}

We show some incorrect examples in the common mistakes with recursion page.

Flood fill

This example doesn't use structures at all, it's simply an array of characters (NB: in real life we would probably use integers, we are using characters here for clarity.)

A common siuation, with a number of variations, is when we have a rectangular grid of values, typically representing an image, and we wish to identify "islands", i.e. sets of neighbouring pixels (or cells) with the same value.

In-class-exercise

We can demonstrate the algorithm ourselves with people and empty seats.

The flood fill function

As an example of flood fill let's imaging we have an NY by NX character array where spaces represent nothing being there and a star, '*', represents something there. We wish to find and label the islands. In this case we assume diagonal cells are connected.

This way of storing the data totally non-recursive: just a two-dimensional array of chars. But the relationship is recursive: cells are connected to their neighbours, which are connected to their neighbours, which are connected to their neighbours, which are connected to their neighbours...

The following function looks for islands of stars ('*') in a sea of spaces or other characters and replaces then with the character area_id (typically 'A', 'B', 'C' etc). As a bonus it returns the size of the island.

int floodfill(char image[NY][NX], int x, int y, char area_id) {
  if ( x < 0 || x >= NX || y < 0 || y >= NY || image[y][x] != '*' ) {
    return 0;
  }

  image[y][x] = area_id; // BEFORE the recursive calls

  // Start at 12-o'clock and go round the cell clockwise
  // trying each of its neighbours in the following order:
  //
  // 8|1|2
  // 7| |3
  // 6|5|4
  //
  return 1 + floodfill(image, x, y-1, area_id)  // Call 1: 12 o'clock
     + floodfill(image, x+1, y-1, area_id)      // Call 2: 1  o'clock
     + floodfill(image, x+1, y, area_id)        // Call 3: 3  o'clock
     + floodfill(image, x+1, y+1, area_id)      // Call 4: 4  o'clock
     + floodfill(image, x, y+1, area_id)        // Call 5: 6  o'clock
     + floodfill(image, x-1, y+1, area_id)      // Call 6: 7  o'clock
     + floodfill(image, x-1, y, area_id)        // Call 7: 9  o'clock
     + floodfill(image, x-1, y-1, area_id) ;    // Call 8: 11 o'clock
}

Notes

1. The first thing we do is to see if the function needs to do anything, if not we immediately return.
   • We check the co-ordinates are in range before checking the pixel (image) value at that point.
2. Then the first thing we do is to mark the pixel as filled.
3. Then we try each of its neighbours, in this case starting with the pixel immediately above ("twelve o'clock") and working round clockwise.

Had we reversed the last two steps then the function would have been forever calling itself.

Showing the depth

Recall that the depth is always a useful number to know. We can add a little depth depth debugging by initially changing the grid character to '0', '1', '2', according to the depth, then doing the recursive function calls and only setting it to the proper character at the end:

The integer value of '1' is the integer value of '0' plus one, and so on up to the integer value of '9' being the integer value of '8' plus one.

int floodfill(int x, int y, char image[NY][NX], char area_id, int depth) {
  int count;
  if ( x < 0 || x >= NX || y < 0 || y >= NY || image[y][x] != '*' ) {
    return 0;
  }


  // When finished debugging do NOT remove this line, change it to area_id
  
  image[y][x] = '0' + (depth % 10);  // Was area_id

  count  = 1 + floodfill(x, y-1, image, area_id, depth + 1)
     + floodfill(x+1, y-1, image, area_id, depth + 1)
     + floodfill(x+1, y, image, area_id, depth + 1)
     + floodfill(x+1, y+1, image, area_id, depth + 1)
     + floodfill(x, y+1, image, area_id, depth + 1)
     + floodfill(x-1, y+1, image, area_id, depth + 1)
     + floodfill(x-1, y, image, area_id, depth + 1)
     + floodfill(x-1, y-1, image, area_id, depth + 1) ;
  
  
  image[y][x] = area_id;
  return count;
}

The whole program

All we need to do now is to go over each point in the image and call floodfill() for each point that needs it and update the area_id character if an island is found.

Notes:

1. To save time we are skippig over must of the steps of the program and only pausing on accesses to the image components and top-level calls to floodfill()
2. We are showing the image grid with non-blank elements colour-coded for clarity.

#include <stdio.h>
// Demonstrate a flood fill, using it to count and label non-empty
// contiguous regions in an image.  

#define NX 8
#define NY 6
void setup(int nx, int ny, char image[NY][NX]);

int floodfill(int x, int y, char image[NY][NX], char area_id, int depth) {
  int count;
  if ( x < 0 || x >= NX || y < 0 || y >= NY || image[y][x] != '*' ) {
    return 0;
  }


  // When finished debugging do NOT remove this line, change it to area_id
  
  image[y][x] = '0' + (depth % 10);  // Was area_id

  count  = 1 + floodfill(x, y-1, image, area_id, depth + 1)
     + floodfill(x+1, y-1, image, area_id, depth + 1)
     + floodfill(x+1, y, image, area_id, depth + 1)
     + floodfill(x+1, y+1, image, area_id, depth + 1)
     + floodfill(x, y+1, image, area_id, depth + 1)
     + floodfill(x-1, y+1, image, area_id, depth + 1)
     + floodfill(x-1, y, image, area_id, depth + 1)
     + floodfill(x-1, y-1, image, area_id, depth + 1) ;
  
  
  image[y][x] = area_id;
  return count;
}

#define IMAX 26
int main() {
  int found = 0;
  char image[NY][NX];
  int ix, iy, nf;
  char letters[32] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
  int sizes[IMAX];
  

  setup(NX, NY, image);
  
  for(iy = 0; iy < NY; ++iy)
    for(ix = 0; ix < NX; ++ix) {
      nf = floodfill(ix, iy, image, letters[found], 0);
      if ( nf ) {
        sizes[found] = nf;
        ++found;
      }
    }

  printf("I found %d islands\n", found);
  for (ix = 0; ix < found; ++ix)
    printf("Island %c size %d\n", letters[ix], sizes[ix]);
  return 0;
}

Step through this code

This will fill the first island/region with the value 'A', the second with the value 'B', etc. and print the total number of islands found together with their sizes.

If "mark as done" in the wrong place

In the following example we have disregarded our comment and are doing the "mark as done" after the recursive calls. We have put in a depth check to stop the program from entering an infinite recursive loop.

Summary

The text of each key point is a link to the place in the web page.