Physics and Astronomy 
Back to top
On this page
Contents RecursionFunctions that call themselvesTo iterate is human, to recurse divine.Computer Science saying
God is subtle but he is not malicious.
Recursion is an extremely simple concept: a function simply calls itself. Recursion refers to a function that calls itself either directly or indirectly. The "call myself again" testOf 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 slightly more sophisticated 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 possiblyempty loop) to see if a function must call itself. A better testA 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 and the calling function doesn't need to worry about it). 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 keyRecursion 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. Automatic and static variablesBy default variables and arrays inside functions are automatic, that is they automatically come into existence when a function starts and are automatically destroyed when it returns. This also means that if a function calls itself then each instance of the function has its own copies of variables. When functions call themselves each instance has its own local variables, unless they are static. Occasionally we may need to have a variable within a function keep its value between calls. The static keyword does this. Static variables retain their values between calls to functions. Static variables are permanent and initialised before the start of the program. Unlike ordinary variables, if no explicit initial value is specified a value of zero (or NULL for pointers) is assumed: Static variables are intialised once, when the program starts to run and default to zero if not explicitly initialised. Static variables are initialised to zero by default. Example of static variablesThe following, somewhat strange, function keeps track of the sum of its arguments. As a diagnostic, it also keeps track of how many times it has been called. // static demonstration #include <stdio.h> #include <math.h> double addup(double k) { static int called; static double sum; sum += k; printf("k is %g new sum is %f\n", k, sum); printf("The function has been called %d times\n", ++called); return sum; } int main() { int i, max=5; double mysum; for(i = 1; i <= max; ++i) mysum = addup(sqrt(i)); printf("The sum of the square roots of 1 to %d is %g\n", max, mysum); return 0; }
A simple example of recursionExample: a factorial functionThe 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; } NB: in the above stepped example, the initial call to factorial() from main() is labelled just "factorial", recursive calls are labelled "factorial2", "factorial3", etc. This illustrates the concepts of recursion:
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 wrappersIt'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 publiclycallable 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); } // Publicallycallable frontend 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); } This is slightly over the top in this case but is a useful illustration. When recursion goes wrongRecursion can be hard to debug. In nonrecursive 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 overthetop 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. Keeping track of the depthWe can get even more details by adding extra debugging argument to tell us the depth. The most obvious way to do this (we shall see another below) is just to add an extra parameter to the function for the function to increment this every time it calls itself: #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! We should be prepared to add some "depth" debugging to our recursive routine. Here we have done several things:
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) recurseerr3.c:13: factorial: Assertion `depth < 12' failed. Abort (core dumped)
Using a wrapperIt may be very inconvenient to add an extra argument to our routine because it is called elsewhere. One way around this is to rename our "real" function, say to factorial_debug(), and create a simple wrapper: long factorial_debug(int i, int depth) {
// Real factorial code here...
}
long factorial(int i) {
return factorial_debug(i, 0);
}
Using a static depth variableThe debugging codes above have the disadvantage of having to add an extra parameter to our function. An alternative way of keeping track of the depth is by having an infunction static variable. We can use a static variable to keep track of the depth as follows:
If a recursive function has more than one place it returns we must take care to decrement the value of depth at each place. This is a particular issue in our case where one of the return statements reads: return i * factorial(i  1); In order to decrement depth after the recursive call we have to declare a temporary variable: long fact1; fact1 = i * factorial(i  1); depth; return fact1; The complete code is shown below with the new parts highlighted in bold red. It can be seen that the code avoids having to declare an extra debugging paramater at the expense of some extra code inside the function. #include <stdlib.h> #include <assert.h> long factorial(int i) { long fact1; static int depth; ++depth; // Print depth information for debugging // for(int 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 ( i > 1 ) { fact1 = i * factorial(i  1); depth; return fact1; } depth; return 1; } Each method has its own advantages and disadvantages and as always we should choose the simplest and clearest option in any given situation. The order of the actionsOften a recursive function has the choice of "do an action then call itself recursively" or "call itself recursively then do an action":
When a function performs an action and calls itself recursively it's usually important to do them in the right order. Directionality and treesOur factorial is an example of a process with a clear order: The argument decreases by one each time. In this sense every recursive call has a "parent" and there is no danger of any loops. TreesSometimes a collection of data structures may also have an inbuilt hierarchy "direction". Consider the example of folders (or directories) on a computer. Folders may contain subfolders and/or normal files but any given folder may only be inside one other flder and there are no loops (a folder may not be its own decendent). We may imagine a recursive "list the contents of this folder and all its decendents" function as: listFolder(this folder):
This isn't always the case however: things may have multiple parents (as in the example of reallife 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:
The correct order: Check, Mark, Recursevoid 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 fillThis 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 situation, 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. InclassexerciseWe can demonstrate the algorithm ourselves with people and empty seats. Before you start: The aim of the exercise is to identify "islands" of people surrounded by empty seats. We will give each island a letter "A", "B" etc.
Note: we won't try all of these in the lecture Parallel, noncounting algorithmWhen a neighbour ("the original neighbour" ) tells you a letter (eg they say to you: "A"):
Expected result: a wave of people saying "A" etc spreading out from the original person. Also notice that if they are doing their job properly each of your neighbours should tell you a letter, even though you already have one. (This is a feature of all versions of the floodfill algorithm.) Completely broken parallel algorithmThe above algorithm is an example of the "If there is nothing to do then return, otherwise mark this node as visited, then do stiuff and recurse" model of recursion. So what happens if we "accidentally forget" to mark the node as visited? When a neighbour tells you a letter (eg they say to you: "A"):
Note the missing step and remember that your task is to see what would happen if something were to follow the instructions in a blind, unthinking manner. Expected result: an awful lot of noise! Serial algorithmThe parallel algorithn relies on the fact that each occupied cell has its own processor (a student!). It also relies on the initiator being able to listen for when the algorithm finishes. So it demonstrates the idea but we couldn't practically program it. We can fix this as follows: When a neighbour ("the original neighbour" ) gives you a letter(eg "A"):
Expected result: a much slower algorithm and a number of students who are waiting for people to reply. Counting serial algorithmWhen a neighbour ("the original neighbour" ) gives you a letter(eg "A"):
Expected result: a count of the size of each island. Version with depth debuggingSince we like to know the depth when looking at recursion we could add some depth debugging: the initial person is told "A zero", they then say "A one" to each of their neighbours, who then say "A two", etc. Expected result: we will audibly hear the depth of the recursion grow and shrink. The flood fill functionAs 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 nonrecursive: just a twodimensional 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. // See if point (x,y) is occupied and if so mark it as "aread_id" // and floodfill its neighbours. int floodfill(char image[NY][NX], int x, int y, char area_id) { // 1. Check if we need to do anything if ( x < 0  x >= NX  y < 0  y >= NY  image[y][x] != '*' ) { return 0; } // 2. Mark as already counted BEFORE the recursive calls image[y][x] = area_id; // 3. Start at 12o'clock and go round the cell clockwise // trying each of its neighbours in the following order: // // 812 // 7 3 // 654 // return 1 + floodfill(image, x, y1, area_id) // Call 1: 12 o'clock + floodfill(image, x+1, y1, 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, x1, y+1, area_id) // Call 6: 7 o'clock + floodfill(image, x1, y, area_id) // Call 7: 9 o'clock + floodfill(image, x1, y1, area_id) ; // Call 8: 11 o'clock } Notes
Had we reversed the last two steps then the function would have been forever calling itself. Showing the depthRecall 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, y1, image, area_id, depth + 1) + floodfill(x+1, y1, 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(x1, y+1, image, area_id, depth + 1) + floodfill(x1, y, image, area_id, depth + 1) + floodfill(x1, y1, image, area_id, depth + 1) ; image[y][x] = area_id; return count; } The whole programAll 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:
#include <stdio.h> // Demonstrate a flood fill, using it to count and label nonempty // 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, y1, image, area_id, depth + 1) + floodfill(x+1, y1, 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(x1, y+1, image, area_id, depth + 1) + floodfill(x1, y, image, area_id, depth + 1) + floodfill(x1, y1, image, area_id, depth + 1) ; image[y][x] = area_id; return count; } #define IMAX 26 int main() { char image[NY][NX]; int found = 0; 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; } 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 placeIn 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. Recursive data structuresNB: this is a fairly advanced topic. One of the reasons why recursion is so useful is that we often have to deal with recursive things. Example: a passive electrical circuitA simple electrical circuit consists of resistors, capacitors and inductors in series or in parallel. But a more complicated cicuit may consist of subcircuits 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 typeRather 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". EnumerationsWe 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;defines Resistor, Capacitor, Inductor as symbolic constants with the values 0, 1 and 2 respectively. In addition it defines "Type" as meaning an int which should only have the values, 0, 1, or 2. The structure definitionWe 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 functionWe 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 impedanceFor 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 morecomplicated circuits?". Series and parallel circuitsGiven 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 morecomplicated circuits. Representing series and parallel circuitsIn 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. 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: typedef enum type { Resistor, Capacitor, Inductor, Series, Parallel } Type; typedef struct component { double value; Type type; unsigned int n; // If parallel or serial struct component *nodes[]; // If parallel or serial } Component; Flexible arraysThe definition above apparently has a huge flaw: it doesn't say the size of the final array! In fact this is a handy feature, avalable only for arrays that are the final member of a dynamicallyallocated array. Imagine first that the array had an actual size, say 2. Then the following code would allocate an array with an array of two pointers: Component *c; c = malloc( sizeof *c); But the only thing a compiler cares about is where an array starts, not where it ends (that's up to us). So a neat trick is to write: Component *c; c = malloc( sizeof *c + n * sizeof *c>nodes); This would allocate us enough space for a structure with an array of n+2 pointers. So C allows us to do away with the dummy array size altogether and to declare an array of size zero and to then allocate as much memory as we like for the array to expand into. Of course, this can only work if the array is the final member of the structure and the memory is dynamically allocated. In this case sizeof *c is the size of the structure with a zerolength final array. 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:
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:
Example: foodsIf you've ever looked at the ingredients of a supermarket sandwich you may well know what we mean! Foods
may consist of several ingredients,
for example mayonnaise may consist of oil, eggs plus
some flavourings and preservatives. But some of the flavourings
such as "curry powder" may themselves have lists of ingredients
(and some of those ingredients may have..)
and the mayonnaise may then be an ingredient in a sandwich.
So a wellknown retailer's chicken tikka sandwich contains:
Question: how would we go about listing the ingredients in a retailer's sandwich platter? A Buffet containing a sandwich platter and some other foods?... 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. Given a simple "Thing", a "Generalised Thing" is
either:
Another example of a "Generalised Thing" is a mathematical expression: a mathematical expression may contain a subexpression in parentheses in which case the subexpression is treated a a whole new toplevel expression and may itself contain parentheses: x * ( 1 + y * ( z  2 * (a + b ) ) ) + a * ( x + y ) Example circuitConsider the following circuit:
What impedance does this present at Uin if Uout is left opencircuit? 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 structureTo represent our "generalised passive component" we need to know three new things:
The recursive impedance functionOur recursive impedance function is now just: #define NODE_MAX 20 typedef enum type { Resistor, Capacitor, Inductor, Series, Parallel } Type; typedef struct component { double value; Type type; unsigned int n; // If parallel or serial struct component *nodes[]; // If parallel or serial } Component; #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. SummaryThe text of each key point is a link to the place in the web page.
