Example programs for PHYS2020

Michael Ashley

The simplest possible C program

Here it is (it doesn't do anything useful):

main() {}

Hello world in C

Here is the generic "Hello world'' program in C:

#include <stdio.h>
int main(void) {
printf("hello world\n");
return 0;
}
Notes:

Note also the strange character sequence "\n" in the string constant "hello world\n". "\n" is called a newline character, and is regarded as a single character in C (i.e., it occupies one byte of storage). "\n" is known as an escape sequence, and the back-slash character is called the escape character.

C escape sequences

/* A program to print out all the special C escape sequences  */
/* Michael Ashley / UNSW / 03-Mar-2003 */

#include <stdio.h> // for printf definition

int main(void) {
printf("audible alert (bell) BEL \\a %d\n" , '\a');
printf("backspace BS \\b %d\n" , '\b');
printf("horizontal tab HT \\t %d\n" , '\t');
printf("newline LF \\n %d\n" , '\n');
printf("vertical tab VT \\v %d\n" , '\v');
printf("formfeed FF \\f %d\n" , '\f');
printf("carriage return CR \\r %d\n" , '\r');
printf("double quote \" \\\" %d\n", '\"');
printf("single quote \' \\\' %d\n", '\'');
printf("question mark ? \\? %d\n" , '\?');
printf("backslash \\ \\\\ %d\n", '\\');
return 0;
}
And here is the output it produces, when compiled with gcc on a GNU/Linux computer:
audible alert (bell) BEL  \a   7
backspace BS \b 8
horizontal tab HT \t 9
newline LF \n 10
vertical tab VT \v 11
formfeed FF \f 12
carriage return CR \r 13
double quote " \" 34
single quote ' \' 39
question mark ? \? 63
backslash \ \\ 92

Mandelbrot program

The Mandelbrot set is a fascinating example of a fractal complexity that can be generated from a very simple equation: z = z*z + c. Here is a program to generate an image of the Mandelbrot set:

/*
A program to generate an image of the Mandelbrot set.

Usage: ./mandelbrot > output
where "output" will be a binary image, 1 byte per pixel
The program will print instructions on stderr as to how to
process the output to produce a JPG file.

Michael Ashley / UNSW / 13-Mar-2003
*/

// Define the range in x and y here:

const double yMin = -1.0;
const double yMax = +1.0;
const double xMin = -2.0;
const double xMax = +0.5;

// And here is the resolution:

const double dxy = 0.005;

#include <stdio.h>
#include <limits.h>

int main(void) {

double cx, cy;
double zx, zy, new_zx;
unsigned char n;
int nx, ny;

// The Mandelbrot calculation is to iterate the equation
// z = z*z + c, where z and c are complex numbers, z is initially
// zero, and c is the coordinate of the point being tested. If
// the magnitude of z remains less than 2 for ever, then the point
// c is in the Mandelbrot set. We write out the number of iterations
// before the magnitude of z exceeds 2, or UCHAR_MAX, whichever is
// smaller.

for (cy = yMin; cy < yMax; cy += dxy) {
for (cx = xMin; cx < xMax; cx += dxy) {
zx = 0.0;
zy = 0.0;
n = 0;
while ((zx*zx + zy*zy < 4.0) & (n != UCHAR_MAX)) {
new_zx = zx*zx - zy*zy + cx;
zy = 2.0*zx*zy + cy;
zx = new_zx;
n++;
}
write(1, &n, sizeof(n)); // Write the result to stdout
}
}

// Now calculate the image dimensions. We use exactly the same
// for loops as above, to guard against any potential rounding errors.

nx = 0;
ny = 0;
for (cx = xMin; cx < xMax; cx += dxy) {
nx++;
}
for (cy = yMin; cy < yMax; cy += dxy) {
ny++;
}

fprintf(stderr, "To process the image: convert -depth 8 -size %dx%d gray:output out.jpg\n",
nx, ny);
return 0;
}

If you enjoyed that, have a look for the program "fracint" for Windows machines, and "mxp" for GNU/Linux. Here are some examples from "mxp", showing the complete Mandelbrot set first, and then enlargements of various sections of it:

The following is an improved version of the above program, using integer arithmetic in the "for" loop, to guarantee the size of the image.

/*
A program to generate an image of the Mandelbrot set.

Usage: ./mandelbrot > output
where "output" will be a binary image, 1 byte per pixel
The program will print instructions on stderr as to how to
process the output to produce a JPG file.

Michael Ashley / UNSW / 13-Mar-2003
*/

const double xCentre = -0.75;
const double yCentre = +0.0;
const int nx = 400;
const int ny = 400;

const double dxy = 0.005;

#include <stdio.h>
#include <limits.h>

int main(void) {

double cx, cy;
double zx, zy, new_zx;
unsigned char n;
int i, j;

// The Mandelbrot calculation is to iterate the equation
// z = z*z + c, where z and c are complex numbers, z is initially
// zero, and c is the coordinate of the point being tested. If
// the magnitude of z remains less than 2 for ever, then the point
// c is in the Mandelbrot set. We write out the number of iterations
// before the magnitude of z exceeds 2, or UCHAR_MAX, whichever is
// smaller.

for (j = 0; j < ny; j++) {
cy = yCentre + (j - ny/2)*dxy;
for (i = 0; i < nx; i++) {
cx = xCentre + (i - nx/2)*dxy;
zx = 0.0;
zy = 0.0;
n = 0;
while ((zx*zx + zy*zy < 4.0) & (n != UCHAR_MAX)) {
new_zx = zx*zx - zy*zy + cx;
zy = 2.0*zx*zy + cy;
zx = new_zx;
n++;
}
write(1, &n, sizeof(n)); // Write the result to stdout
}
}

fprintf(stderr, "To process the image: convert -depth 8 -size %dx%d gray:output out.jpg\n",
nx, ny);
return 0;
}

A circle printing program

/* 
A program to draw a circular area using ASCII characters.

Michael Ashley / UNSW / 12-Mar-2003
*/

#include <stdio.h>

const int xRange = 40;
const int yRange = 20;
const int radius = 15;

int main(void) {
int x, y;

for (y = -yRange; y < yRange; y++) {
for (x = -xRange; x < xRange; x++) {
if (x*x + y*y > radius * radius) {
printf (".");
} else {
printf ("#");
}
}
printf("\n");
}
return 0;
}

Note the use of const int definitions to describe dimensions of the drawing. It is good practice to define any special numbers in your program in this way, since it makes the logic easier to understand, and also makes it easier to change a number in the future.

And here is what the output looks like:

................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
........................................#.......................................
...................................###########..................................
.................................###############................................
...............................###################..............................
..............................#####################.............................
.............................#######################............................
............................#########################...........................
............................#########################...........................
...........................###########################..........................
...........................###########################..........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
.........................###############################........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
..........................#############################.........................
...........................###########################..........................
...........................###########################..........................
............................#########################...........................
............................#########################...........................
.............................#######################............................
..............................#####################.............................
...............................###################..............................
.................................###############................................
...................................###########..................................
........................................#.......................................
................................................................................
................................................................................
................................................................................
................................................................................

The GNU readline package

// A program to demonstrate the use of the GNU readline package.
// Michael Ashley / UNSW / 04 April 2003

#include <stdio.h>
#include <math.h>

// The following two lines are necessary for GNU readline:

#include <readline/readline.h>
#include <readline/history.h>

int main(void) {
char *cp;
float x;
int ret;

printf("This program calculates square roots\n");

while (1) {

// Grab an input line from the user. readline automatically
// allocates space for the line, and returns a pointer to it.

cp = readline ("please input a number (ENTER to exit) > ");

// If nothing was entered, exit the program.

if (cp == NULL) {
printf ("\n");
break;
}

// Else add the line to the history buffer, so that the user
// can use the up/down arrow keys to navigate through previous
// responses.

add_history(cp);

// Try to read a floating point number from the input line.

ret = sscanf(cp, "%f", &x);

// Free the space used by the input line.

free(cp);

// Now check the results from sscanf.

if (ret == 1) {
printf("The square root of %f is %f\n", x, sqrt(x));
} else if (ret == EOF) {
break;
} else {
printf("that wasn't a valid number...\n");
}
}

return 0;
}

To compile the above program, do, e.g., "gcc example.c -lm -lreadline -lcurses".

Random numbers

// An example showing how to generate random floating
// point numbers between 0 and 1.

// Michael Ashley / UNSW / 02-Apr-2003

#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

int main(void) {
struct timeval t;
int i;

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

printf("The time is %ld.%06ld seconds since the Unix epoch\n",
t.tv_sec, t.tv_usec);

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Generate and print ten random integers between 0 and RAND_MAX.

for (i = 0; i < 10; i++) {
printf("%ld\n", random());
}

// Generate and print ten random numbers between 0 and 1.0.

for (i = 0; i < 10; i++) {
printf("%.15f\n", random()/((double) RAND_MAX));
}

return 0;
}

Sorting things using qsort

/* 
An example showing the use of the "qsort" function for
sorting 100 random numbers.

Michael Ashley / UNSW / 11-Apr-2003
*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>

int myComparisonFunction(const void *x, const void *y) {

// x and y are pointers to doubles.

// Returns -1 if x < y
// 0 if x == y
// +1 if x > y

double dx, dy;

dx = *(double *)x;
dy = *(double *)y;

if (dx < dy) {
return -1;
} else if (dx > dy) {
return +1;
}
return 0;
}

int main(void) {
struct timeval t;
double data[100];
int i;

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Generate 100 random numbers between 0.0 and 1.0.

for (i = 0; i < ((sizeof data) / (sizeof data[0])); i++) {
data[i] = random()/((double) RAND_MAX);
}

// Now sort them into ascending order.

qsort(data, (sizeof data) / (sizeof data[0]), sizeof data[0],
&myComparisonFunction);

// And print them out.

for (i = 0; i < ((sizeof data) / (sizeof data[0])); i++) {
printf ("data[%d] = %f\n", i, data[i]);
}

return 0;
}

And here is a really tricky example showing how to sort character strings.

/* 
An example showing the use of the "qsort" function for
sorting up to 100 character strings.

Michael Ashley / UNSW / 11-Apr-2003
*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>
#include <string.h>

int myComparisonFunction(const void *x, const void *y) {

// This is a wrapper for the "strcmp" function.

// x and y are pointers to pointers to character strings

// Returns -1 if x < y
// 0 if x == y
// +1 if x > y

return (strcmp(*(const char **)x, *(const char **)y));

}

int main(void) {
struct timeval t;
char *data[100];
char *p;
int i, numStrings;

// Read up to 100 character strings from standard input.

for (i = 0; i < ((sizeof data) / (sizeof data[0])); i++) {
if (1 != scanf("%as", &data[i])) {
break;
}
}

numStrings = i;

printf("%d strings found\n", numStrings);

// Sort them alphabetically.

qsort(data, numStrings, sizeof data[0],
&myComparisonFunction);

// Print them out, and free the memory malloc'ed by scanf.

for (i = 0; i < numStrings; i++) {
printf("data[%d] = \"%s\"\n", i, data[i]);
free(data[i]);
}

return 0;
}

Numerical integration - gravity

// A program to illustrate numerical integration for solving gravitational
// interactions between two bodies.

#include <stdio.h>
#include <math.h>

// Physical constants.

const double bigG = 6.6726E-11; // Nm^2/kg^2

// Useful typedefs, i.e., our own data types for the problem at hand.

// A 3-D vector.

typedef struct {
double x, y, z;
} vector;

// All the information that we need to know about a test particle.

typedef struct {
double mass; // The mass of the particle, in kg.
double epoch; // The time in seconds at which the following
// two items are specified.
vector pos; // The particle's position.
vector vel; // The particle's velocity.
vector force; // The force acting on the particle.
double potential; // The gravitational potential that the particle is in.
double ke; // The particle's kinetic energy.
double pe; // The particle's potential energy.
vector I; // The particle's angular momentum.
} particle;

// Now we define some functions that we will need.

void gravity(particle *p1, particle *p2) {

// Calculates the force and potential acting on particle p1 due to
// particle p2.

double forceMagnitude;
double distance;
vector d;

// Calculate the position offset between the two particles.

d.x = p1->pos.x - p2->pos.x;
d.y = p1->pos.y - p2->pos.y;
d.z = p1->pos.z - p2->pos.z;

// The magnitude of the position offset vector is the distance between
// the particles.

distance = sqrt(d.x*d.x + d.y*d.y + d.z*d.z);

// Now calculate the magnitude of the force acting on p1 due to p2.

forceMagnitude = bigG * p1->mass * p2->mass / (distance * distance);

// And orient it correctly in 3 dimensions by multiplying by a unit
// vector from p2 to p1.

p1->force.x = - forceMagnitude * d.x / distance;
p1->force.y = - forceMagnitude * d.y / distance;
p1->force.z = - forceMagnitude * d.z / distance;

// Finally, calculate the potential at p1 due to p2.

p1->potential = - bigG * p2->mass / distance;
}

void update(particle *p, double deltaT) {

// Updates the properties of particle p after a time interval
// deltaT [seconds].

vector acc;

// Calculate the acceleration due to the force acting on the particle.

acc.x = p->force.x / p->mass;
acc.y = p->force.y / p->mass;
acc.z = p->force.z / p->mass;

// Update the position vector, assuming a constant acceleration
// over the time interval.

p->pos.x += (p->vel.x + 0.5 * acc.x * deltaT) * deltaT;
p->pos.y += (p->vel.y + 0.5 * acc.y * deltaT) * deltaT;
p->pos.z += (p->vel.z + 0.5 * acc.z * deltaT) * deltaT;

// Update the velocity vector.

p->vel.x += acc.x * deltaT;
p->vel.y += acc.y * deltaT;
p->vel.z += acc.z * deltaT;

// Update the epoch.

p->epoch += deltaT;

// Update the kinetic energy.

p->ke = 0.5 * p->mass * (p->vel.x * p->vel.x +
p->vel.y * p->vel.y +
p->vel.z * p->vel.z);

// Update the potential energy.

p->pe = p->mass * p->potential;

// Update the angular momentum.

p->I.x = p->mass * (p->pos.y * p->vel.z - p->pos.z * p->vel.y);
p->I.y = p->mass * (p->pos.z * p->vel.x - p->pos.x * p->vel.z);
p->I.z = p->mass * (p->pos.x * p->vel.y - p->pos.y * p->vel.x);
}

void displayParticle(particle *p) {

// Displays information about particle p.

printf("m %10.2e; x (%10.2e,%10.2e,%10.2e); v (%10.2e,%10.2e,%10.2e)",
p->mass,
p->pos.x, p->pos.y, p->pos.z,
p->vel.x, p->vel.y, p->vel.z);
}

void displayEnergy(particle *p) {

printf(" E (%10.2e =%10.2e +%10.2e)", p->ke + p->pe, p->ke, p->pe);
}

void displayI(particle *p) {

printf(" I (%10.2e,%10.2e, %10.2e)", p->I.x, p->I.y, p->I.z);
}

void displayCM(particle *p1, particle *p2) {
vector cm;

cm.x = (p1->mass * p1->pos.x + p2->mass * p2->pos.x) /
(p1->mass + p2->mass);
cm.y = (p1->mass * p1->pos.y + p2->mass * p2->pos.y) /
(p1->mass + p2->mass);
cm.z = (p1->mass * p1->pos.z + p2->mass * p2->pos.z) /
(p1->mass + p2->mass);
printf(" cm (%10.2e,%10.2e, %10.2e)", cm.x, cm.y, cm.z);
}


// Here are two representative particles:

particle sun = {1.989E30, // Mass in kg
0.0,
{0.0, 0.0, 0.0}, // The initial position
{0.0, 0.0, 0.0}}; // The initial velocity

particle earth = {5.976E24, // Mass in kg
0.0,
{149.6E9, 0.0, 0.0}, // The initial position
{0.0, 27.786E3, 0.0}}; // The initial velocity

const double deltaT = 3600.0; // The time interval in seconds

int main(void) {

double time;

// Choose the position of the Sun to put the centre of mass at (0,0,0).

sun.pos.x = - earth.pos.x * earth.mass/sun.mass;
sun.pos.y = - earth.pos.y * earth.mass/sun.mass;
sun.pos.z = - earth.pos.z * earth.mass/sun.mass;

// Choose the velocity of the Sun to keep the centre of mass
// at a constant position.

sun.vel.x = 0.0;
sun.vel.y = earth.vel.y * sun.pos.x / earth.pos.x;
sun.vel.z = 0.0;

// Now follow the motion for one year...

for (time = 0.0; time < 365.2422 * 24 * 3600; time += deltaT) {
gravity(&sun, &earth);
gravity(&earth, &sun);

update(&sun, deltaT);
update(&earth, deltaT);

// displayParticle(&sun);
// printf("\n");
#if 0
displayParticle(&sun);
displayEnergy(&sun);
displayI(&sun);
// displayCM(&earth, &sun);
printf("\n");
#endif
displayParticle(&earth);
displayEnergy(&earth);
displayI(&earth);

printf("\n");

}
return 0;
}

1D cellular automata - text output

// Cellular automata
// Wolfram, "A New Kind of Science", Chapter 2 and 3.

// Michael Ashley / UNSW / 11-May-2003

#define displayWidth 80
#define maxSteps 24

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

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the values of the neighbouring two cells.

The calculation of the new value therefore depends on the values of
3 single-bit quantities, giving 8 possible alternative configurations.

The calculation can therefore be represented using a array giving the
8 possibilities. We call this array "rule".
*/

int rule[8];

typedef struct {
unsigned char cell[displayWidth + 2 * maxSteps];
} state;

void initialise(state * s) {
int i;
for (i = 0; i < (sizeof s->cell) / (sizeof s->cell[0]); i++) {
s->cell[i] = 0;
if (i == maxSteps + displayWidth / 2)
s->cell[i] = 1;
}
}

void createRule(int r[8], unsigned char n) {
int i;
for (i = 0; i < 8; i++) {
r[i] = 0x01 & (n >> i);
}
}

void applyRule(state * prev, state * next, int i) {

/* This function takes an existing state, prev, and applies the
evolution rule to the i'th element, returning the result in next. */

next->cell[i] = rule[prev->cell[i - 1] * 4 +
prev->cell[i] * 2 +
prev->cell[i + 1]];
}

void evolve(state * prev, state * next) {

/* Applies the rule to all elements (apart from the end-points) of
an existing state, prev, and returns the result in next. */

int i;
for (i = 1; i < (sizeof prev->cell) / (sizeof prev->cell[0]) - 1; i++) {
applyRule(prev, next, i);
}
}

void displayState(state * s) {
int i;
for (i = 0; i < displayWidth; i++) {
if (s->cell[i + maxSteps]) {
printf("*");
} else {
printf(" ");
}
}
printf("\n");
}

void usage(char *prog) {
printf("Usage: %s rule
where rule is an integer between 0 and 255 inclusive\n", prog);
exit(1);
}

int main(int argc, char **argv) {
int n, i;
state s0, s1;

assert(maxSteps % 2 == 0); // Ensure we have an even number of steps.

initialise(&s0);
initialise(&s1);

if (argc != 2 ||
1 != sscanf(argv[1], "%d", &n) ||
n < 0 ||
n > 255)
usage(argv[0]);

createRule(rule, n);

for (i = 0; i < maxSteps / 2; i++) {
displayState(&s0);
evolve(&s0, &s1);
displayState(&s1);
evolve(&s1, &s0);
}
return 0;
}

Interesting features of this program include:

And here are some examples of the output from this program:

Rule 2, a class 1 cellular automata

                                        *                                       
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*

Rule 22, a class 2 cellular automata

                                        *                                       
***
* *
*** ***
* *
*** ***
* * * *
*** *** *** ***
* *
*** ***
* * * *
*** *** *** ***
* * * *
*** *** *** ***
* * * * * * * *
*** *** *** *** *** *** *** ***
* *
*** ***
* * * *
*** *** *** ***
* * * *
*** *** *** ***
* * * * * * * *
*** *** *** *** *** *** *** ***

Rule 30, a class 3 cellular automata

                                        *                                       
***
** *
** ****
** * *
** **** ***
** * * *
** **** ******
** * *** *
** **** ** * ***
** * * **** ** *
** **** ** * * ****
** * *** ** ** * *
** **** ** *** *** ** ***
** * * *** * *** * *
** **** ** * * ***** *******
** * *** **** * *** *
** **** ** *** ** ** * ***
** * * *** * ** *** **** ** *
** **** ** * ****** * * *** ****
** * *** **** **** *** ** * *
** **** ** *** * ** * * * *** ***
** * * *** * *** ** * *** ** * * * *
** **** ** * *** * * **** * * ** ******

Rule 110, a class 4 cellular automata

                                        *                                       
**
***
** *
*****
** *
*** **
** * ***
******* *
** ***
*** ** *
** * *****
***** ** *
** * *** **
*** **** * ***
** * ** ***** *
******** ** ***
** **** ** *
*** ** * *****
** * *** **** *
***** ** *** * **
** * ***** * ** ***
*** ** ** ******** *
** * ****** ** ***

1D cellular automata - image output

// Cellular automata
// Wolfram, "A New Kind of Science", Chapter 2 and 3.

// Michael Ashley / UNSW / 11-May-2003

#define displayWidth 640
#define maxSteps 480

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

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the values of the neighbouring two cells.

The calculation of the new value therefore depends on the values of
3 single-bit quantities, giving 8 possible alternative configurations.

The calculation can therefore be represented using a array giving the
8 possibilities. We call this array "rule".
*/

int rule[8];

typedef struct {
unsigned char cell[displayWidth + 2 * maxSteps];
} state;

void initialise(state * s) {
int i;
for (i = 0; i < (sizeof s->cell) / (sizeof s->cell[0]); i++) {
s->cell[i] = 0;
if (i == maxSteps + displayWidth / 2)
s->cell[i] = 1;
}
}

void createRule(int r[8], unsigned char n) {
int i;
for (i = 0; i < 8; i++) {
r[i] = 0x01 & (n >> i);
}
}

void applyRule(state * prev, state * next, int i) {

/* This function takes an existing state, prev, and applies the
evolution rule to the i'th element, returning the result in next. */

next->cell[i] = rule[prev->cell[i - 1] * 4 +
prev->cell[i] * 2 +
prev->cell[i + 1]];
}

void evolve(state * prev, state * next) {

/* Applies the rule to all elements (apart from the end-points) of
an existing state, prev, and returns the result in next. */

int i;
for (i = 1; i < (sizeof prev->cell) / (sizeof prev->cell[0]) - 1; i++) {
applyRule(prev, next, i);
}
}

void displayState(state * s) {
assert(displayWidth == write(1, &s->cell[maxSteps], displayWidth));
}

void usage(char *prog) {
printf("Usage: %s rule | display -size %dx%d -depth 1 gray:-
where rule is an integer between 0 and 255 inclusive\n",
prog, displayWidth, maxSteps);
exit(1);
}

int main(int argc, char **argv) {
int n, i;
state s0, s1;

assert(maxSteps % 2 == 0); // Ensure we have an even number of steps.

initialise(&s0);
initialise(&s1);

if (argc != 2 ||
1 != sscanf(argv[1], "%d", &n) ||
n < 0 ||
n > 255)
usage(argv[0]);

createRule(rule, n);

for (i = 0; i < maxSteps / 2; i++) {
displayState(&s0);
evolve(&s0, &s1);
displayState(&s1);
evolve(&s1, &s0);
}
return 0;
}

Interesting features of this program include:

And here are some examples of the output from this program:

Rule 2, a class 1 cellular automata

Rule 22, a class 2 cellular automata

Rule 30, a class 3 cellular automata

Rule 110, a class 4 cellular automata

John Horton Conway's Game of Life - simple version

// John Horton Conway's game of Life.

// Michael Ashley / UNSW / 23-May-2003

#define displayWidth 80
#define displayHeight 24

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <unistd.h>
#include <sys/time.h>

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the sum of the values of the neighbouring 8 cells.

The rules are:

Death: If an occupied cell has 0, 1, 4, 5, 6, 7, or 8 occupied
neighbours, the organism dies (0, 1 neighbours: of loneliness;
4 thru 8: of overcrowding).

Survival: If an occupied cell has two or three neighbours, the organism
survives to the next generation.

Birth: If an unoccupied cell has three occupied neighbours, it becomes
occupied.

These rules can be written in terms of a 2D array, where the first index
is either 0 or 1 depending on the value of cell under consideration, and
the second index ranges from 0 to 8 inclusive, and is the number of
occupied nearest neighbours. The value of the array gives the state of
the cell in the next generation.
*/

int rule[2][9] = {{0,0,0,1,0,0,0,0,0},
{0,0,1,1,0,0,0,0,0}};

/*
Now we create a type that contains all the information we need to
know about the state of the system.
*/

typedef struct {
unsigned char cell[displayHeight][displayWidth];
} state;

void initialise(state * s) {

// Initialises the state pointed to by s. This is where we put our
// initial conditions.

int i, j;
struct timeval t;

#if 0
// Zero the state.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = 0;
}
}

// A "glider" pattern.

s->cell[40][10] = 1;
s->cell[41][10] = 1;
s->cell[42][10] = 1;
s->cell[42][11] = 1;
s->cell[41][12] = 1;
#endif

// A random pattern.

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Here we randomly choose 1/8th of the cells to be alive.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = random() > 7*(RAND_MAX/8);
}
}

}

int nearestNeighbours(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j]. We just sum up the neighbouring 8 cells, with
// careful allowance for hitting the boundary.

return
(i > 0 && j > 0 && s->cell[i-1][j-1]) +
(i > 0 && s->cell[i-1][j] ) +
(i > 0 && j < displayWidth-1 && s->cell[i-1][j+1]) +
( j > 0 && s->cell[i] [j-1]) +
( j < displayWidth-1 && s->cell[i] [j+1]) +
(i < displayHeight-1 && j > 0 && s->cell[i+1][j-1]) +
(i < displayHeight-1 && s->cell[i+1][j] ) +
(i < displayHeight-1 && j < displayWidth-1 && s->cell[i+1][j+1]);
}

void evolve(state * prev, state * next) {

// Evolves state *prev by one generation, returning the result in *next.

int i, j;

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];
}
}
}

void displayState(state * s) {

// Displays state *s using ASCII characters.

int i, j;

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
if (s->cell[i][j]) {
printf("*");
} else {
printf(" ");
}
}
printf("\n");
}
}

int main(int argc, char **argv) {
state s0, s1;
int i;

initialise(&s0);

// To display the state at each generation, uncomment the "displayState"
// lines. To slow down the display, uncomment the "usleep" lines.

for (i = 0; i < 50000; i++) {
// displayState(&s0);
evolve(&s0, &s1);
// usleep(100000);
// displayState(&s1);
evolve(&s1, &s0);
// usleep(100000);
}
return 0;
}

And here is an example of a random initial condition, with 1/8th of the cells alive:

*     *  *  *    *    *        *       *  *                     *    *     * *  
*** ** * * *** * * ** * *
* * * * ** * ** *
* * * * * * * * *
** * * * * * **
** * * * * * ** * * *
* * * ** * * * *
** ** ** * * * * * ** *
* ** * * * * * * *
* * * * ** *** ** *
* * * * * *
* * * * ** * * *
* * * * * * * * * ** *
* * * * * * * * * * *
** * * * * * * ** * * ** * * *
* ** * * * * * * * * *
* * * * ** ** * * * * * * * * * * * **
* **** * * * * * * **
* * * * * * * * * * * * *
* * * * ** * * * * ** *
* * * * * * * ** * *
* * *** * * * * *
* * * * * * * * * *
* * * * ** * * **
** ** ** ** **
** ** ** ** **

And this is the result after 100 generations:

                                                                                
**
*
*
* *
* **
**
*
*
*


***
**
** * * * * * **
** * * * * **
* * * * *
**
***


John Horton Conway's Game of Life - inline function version

// John Horton Conway's game of Life.

// Michael Ashley / UNSW / 23-May-2003

#define displayWidth 80
#define displayHeight 24

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <unistd.h>
#include <sys/time.h>

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the sum of the values of the neighbouring 8 cells.

The rules are:

Death: If an occupied cell has 0, 1, 4, 5, 6, 7, or 8 occupied
neighbours, the organism dies (0, 1 neighbours: of loneliness;
4 thru 8: of overcrowding).

Survival: If an occupied cell has two or three neighbours, the organism
survives to the next generation.

Birth: If an unoccupied cell has three occupied neighbours, it becomes
occupied.

These rules can be written in terms of a 2D array, where the first index
is either 0 or 1 depending on the value of cell under consideration, and
the second index ranges from 0 to 8 inclusive, and is the number of
occupied nearest neighbours. The value of the array gives the state of
the cell in the next generation.
*/

int rule[2][9] = {{0,0,0,1,0,0,0,0,0},
{0,0,1,1,0,0,0,0,0}};

/*
Now we create a type that contains all the information we need to
know about the state of the system.
*/

typedef struct {
unsigned char cell[displayHeight][displayWidth];
} state;

void initialise(state * s) {

// Initialises the state pointed to by s. This is where we put our
// initial conditions.

int i, j;
struct timeval t;

#if 0
// Zero the state.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = 0;
}
}

// A "glider" pattern.

s->cell[40][10] = 1;
s->cell[41][10] = 1;
s->cell[42][10] = 1;
s->cell[42][11] = 1;
s->cell[41][12] = 1;
#endif

// A random pattern.

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Here we randomly choose 1/8th of the cells to be alive.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = random() > 7*(RAND_MAX/8);
}
}

}

inline int nearestNeighbours(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j]. We just sum up the neighbouring 8 cells, with
// careful allowance for hitting the boundary.

return
(i > 0 && j > 0 && s->cell[i-1][j-1]) +
(i > 0 && s->cell[i-1][j] ) +
(i > 0 && j < displayWidth-1 && s->cell[i-1][j+1]) +
( j > 0 && s->cell[i] [j-1]) +
( j < displayWidth-1 && s->cell[i] [j+1]) +
(i < displayHeight-1 && j > 0 && s->cell[i+1][j-1]) +
(i < displayHeight-1 && s->cell[i+1][j] ) +
(i < displayHeight-1 && j < displayWidth-1 && s->cell[i+1][j+1]);
}

inline void evolve(state * prev, state * next) {

// Evolves state *prev by one generation, returning the result in *next.

int i, j;

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];
}
}
}

void displayState(state * s) {

// Displays state *s using ASCII characters.

int i, j;

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
if (s->cell[i][j]) {
printf("*");
} else {
printf(" ");
}
}
printf("\n");
}
}

int main(int argc, char **argv) {
state s0, s1;
int i;

initialise(&s0);

// To display the state at each generation, uncomment the "displayState"
// lines. To slow down the display, uncomment the "usleep" lines.

for (i = 0; i < 50000; i++) {
// displayState(&s0);
evolve(&s0, &s1);
// usleep(100000);
// displayState(&s1);
evolve(&s1, &s0);
// usleep(100000);
}
return 0;
}

John Horton Conway's Game of Life - boundary version

// John Horton Conway's game of Life.

// Michael Ashley / UNSW / 23-May-2003

#define displayWidth 80
#define displayHeight 24

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <unistd.h>
#include <sys/time.h>

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the sum of the values of the neighbouring 8 cells.

The rules are:

Death: If an occupied cell has 0, 1, 4, 5, 6, 7, or 8 occupied
neighbours, the organism dies (0, 1 neighbours: of loneliness;
4 thru 8: of overcrowding).

Survival: If an occupied cell has two or three neighbours, the organism
survives to the next generation.

Birth: If an unoccupied cell has three occupied neighbours, it becomes
occupied.

These rules can be written in terms of a 2D array, where the first index
is either 0 or 1 depending on the value of cell under consideration, and
the second index ranges from 0 to 8 inclusive, and is the number of
occupied nearest neighbours. The value of the array gives the state of
the cell in the next generation.
*/

int rule[2][9] = {{0,0,0,1,0,0,0,0,0},
{0,0,1,1,0,0,0,0,0}};

/*
Now we create a type that contains all the information we need to
know about the state of the system.
*/

typedef struct {
unsigned char cell[displayHeight][displayWidth];
} state;

void initialise(state * s) {

// Initialises the state pointed to by s. This is where we put our
// initial conditions.

int i, j;
struct timeval t;

#if 0
// Zero the state.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = 0;
}
}

// A "glider" pattern.

s->cell[40][10] = 1;
s->cell[41][10] = 1;
s->cell[42][10] = 1;
s->cell[42][11] = 1;
s->cell[41][12] = 1;
#endif

// A random pattern.

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Here we randomly choose 1/8th of the cells to be alive.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = random() > 7*(RAND_MAX/8);
}
}

}

inline int nearestNeighbours(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j]. We just sum up the neighbouring 8 cells, with
// careful allowance for hitting the boundary.

return
(i > 0 && j > 0 && s->cell[i-1][j-1]) +
(i > 0 && s->cell[i-1][j] ) +
(i > 0 && j < displayWidth-1 && s->cell[i-1][j+1]) +
( j > 0 && s->cell[i] [j-1]) +
( j < displayWidth-1 && s->cell[i] [j+1]) +
(i < displayHeight-1 && j > 0 && s->cell[i+1][j-1]) +
(i < displayHeight-1 && s->cell[i+1][j] ) +
(i < displayHeight-1 && j < displayWidth-1 && s->cell[i+1][j+1]);
}

inline int nearestNeighbours2(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j], assuming that we are at least one cell away
// from the boundary.

return s->cell[i-1][j-1] +
s->cell[i-1][j] +
s->cell[i-1][j+1] +
s->cell[i] [j-1] +
s->cell[i] [j+1] +
s->cell[i+1][j-1] +
s->cell[i+1][j] +
s->cell[i+1][j+1];
}

inline void evolve(state * prev, state * next) {

// Evolves state *prev by one generation, returning the result in
// *next. This function makes special allowance for the boundary
// of the space.

int i, j;

// Process the first row.

i = 0;
for (j = 0; j < displayWidth; j++) {
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];
}

for (i = 1; i < displayHeight-1; i++) {

// Process the first column of each row.

j = 0;
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];


// Process cells that are not on the boundary. This is where
// almost all the computation takes place.

for (j = 1; j < displayWidth-1; j++) {
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours2(prev, i, j)];
}

// Process the last column of each row.

next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];
}

// Process the last row.

for (j = 0; j < displayWidth; j++) {
next->cell[i][j] =
rule[prev->cell[i][j]][nearestNeighbours(prev, i, j)];
}
}

void displayState(state * s) {

// Displays state *s using ASCII characters.

int i, j;

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
if (s->cell[i][j]) {
printf("*");
} else {
printf(" ");
}
}
printf("\n");
}
}

int main(int argc, char **argv) {
state s0, s1;
int i;

initialise(&s0);

// To display the state at each generation, uncomment the "displayState"
// lines. To slow down the display, uncomment the "usleep" lines.

for (i = 0; i < 50000; i++) {
// displayState(&s0);
evolve(&s0, &s1);
// usleep(100000);
// displayState(&s1);
evolve(&s1, &s0);
// usleep(100000);
}
return 0;
}

John Horton Conway's Game of Life - version to automatically find gliders

// John Horton Conway's game of Life.
// Modified to automatically find glider patterns.
// Michael Ashley / UNSW / 23-May-2003

#define displayWidth (8*80)
#define displayHeight (8*24)

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <unistd.h>
#include <sys/time.h>

/*
Each cell has a value of 0 or 1.

The value of a cell in the next generation depends on its current
value, and the sum of the values of the neighbouring 8 cells.

The rules are:

Death: If an occupied cell has 0, 1, 4, 5, 6, 7, or 8 occupied
neighbours, the organism dies (0, 1 neighbours: of loneliness;
4 thru 8: of overcrowding).

Survival: If an occupied cell has two or three neighbours, the organism
survives to the next generation.

Birth: If an unoccupied cell has three occupied neighbours, it becomes
occupied.

These rules can be written in terms of a 2D array, where the first index
is either 0 or 1 depending on the value of cell under consideration, and
the second index ranges from 0 to 8 inclusive, and is the number of
occupied nearest neighbours. The value of the array gives the state of
the cell in the next generation.
*/

int rule[2][9] = {{0,0,0,1,0,0,0,0,0},
{0,0,1,1,0,0,0,0,0}};

/*
Now we create a type that contains all the information we need to
know about the state of the system.
*/

typedef struct {
unsigned char cell[displayHeight][displayWidth];
} state;

void initialise(state * s) {

// Initialises the state pointed to by s. This is where we put our
// initial conditions.

int i, j;
struct timeval t;

// Zero the state.

for (i = 0; i < displayHeight; i++) {
for (j = 0; j < displayWidth; j++) {
s->cell[i][j] = 0;
}
}

#if 0
// A "glider" pattern.

s->cell[40][10] = 1;
s->cell[41][10] = 1;
s->cell[42][10] = 1;
s->cell[42][11] = 1;
s->cell[41][12] = 1;
#endif

// A random pattern.

// Obtain the time of day, to microsecond resolution.

assert(0 == gettimeofday(&t, NULL));

// Use the number of microseconds as the seed for the system
// random number generator.

srandom(t.tv_usec);

// Here we randomly choose 1/2th of the cells to be alive,
// within a central region of the state.

for (i = 3*displayHeight/8; i < 5*displayHeight/8; i++) {
for (j = 3*displayWidth/8; j < 5*displayWidth/8; j++) {
s->cell[i][j] = random() > 4*(RAND_MAX/8);
}
}
}

void displayState(state * s, int i, int j) {

// Displays up to a 20x20 pixel box including cell [i][j] of state
// *s using ASCII characters.

int iMin, iMax, jMin, jMax;

iMin = i - 10;
if (iMin < 0) iMin = 0;
jMin = j - 10;
if (jMin < 0) jMin = 0;
iMax = i + 10;
if (iMax >= displayHeight) iMax = displayHeight-1;
jMax = j + 10;
if (jMax >= displayWidth) jMax = displayWidth-1;

for (i = iMin; i < iMax; i++) {
for (j = jMin; j < jMax; j++) {
if (s->cell[i][j]) {
printf("*");
} else {
printf(" ");
}
}
printf("\n");
}
}

inline int nearestNeighbours(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j]. We just sum up the neighbouring 8 cells, with
// careful allowance for hitting the boundary.

return
(i > 0 && j > 0 && s->cell[i-1][j-1]) +
(i > 0 && s->cell[i-1][j] ) +
(i > 0 && j < displayWidth-1 && s->cell[i-1][j+1]) +
( j > 0 && s->cell[i] [j-1]) +
( j < displayWidth-1 && s->cell[i] [j+1]) +
(i < displayHeight-1 && j > 0 && s->cell[i+1][j-1]) +
(i < displayHeight-1 && s->cell[i+1][j] ) +
(i < displayHeight-1 && j < displayWidth-1 && s->cell[i+1][j+1]);
}

inline int nearestNeighbours2(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j], assuming that we are at least one cell away
// from the boundary.

return s->cell[i-1][j-1] +
s->cell[i-1][j] +
s->cell[i-1][j+1] +
s->cell[i] [j-1] +
s->cell[i] [j+1] +
s->cell[i+1][j-1] +
s->cell[i+1][j] +
s->cell[i+1][j+1];
}

inline int nearestNeighbours3(state *s, int i, int j) {

// Returns the number of nearest neighbours in the state *s at
// location [i][j], assuming that we are at least one cell away
// from the boundary, and that the previous cell (at [i][j-1])
// had zero nearest neighbours.

return s->cell[i] [j-1] +
s->cell[i-1][j+1] +
s->cell[i] [j+1] +
s->cell[i+1][j+1];
}

inline int processBoundary(state *prev, state *next, int i, int j) {

// Evolve the cell [i][j] on the boundary of the universe, from
// state *prev to state *next, returning the number of nearest
// neighbours. If a cell will become alive, display the
// neighbouring cells and exit.

int n;

n = nearestNeighbours(prev, i, j);
if (n != 0) {
displayState(prev, i, j);
exit(0);
}
next->cell[i][j] = rule[prev->cell[i][j]][n];
return n;
}

inline void evolve(state * prev, state * next) {

// Evolves state *prev by one generation, returning the result in
// *next. This function makes special allowance for the boundary
// of the space.

int i, j, n;

// Process the first row.

i = 0;
for (j = 0; j < displayWidth; j++) {
processBoundary(prev, next, i, j);
}

for (i = 1; i < displayHeight-1; i++) {

// Process the first column of each row.

j = 0;
n = processBoundary(prev, next, i, j);

// Process cells that are not on the boundary. This is where
// almost all the computation takes place.

for (j = 1; j < displayWidth-1; j++) {

// Handle the special case that the last cell had no nearest
// neighbours.

if (n == 0) {
n = nearestNeighbours3(prev, i, j);
next->cell[i][j] = rule[prev->cell[i][j]][n];
} else {
n = nearestNeighbours2(prev, i, j);
next->cell[i][j] = rule[prev->cell[i][j]][n];
}
}

// Process the last column of each row.

n = nearestNeighbours(prev, i, j);
if (n != 0) {
displayState(prev, i, j);
exit(0);
}
next->cell[i][j] = rule[prev->cell[i][j]][n];
}

// Process the last row.

for (j = 0; j < displayWidth; j++) {
processBoundary(prev, next, i, j);
}
}

int main(int argc, char **argv) {
state s0, s1;
int i;

initialise(&s0);

// To display the state at each generation, uncomment the "displayState"
// lines. To slow down the display, uncomment the "usleep" lines.

for (i = 0; i < 50000; i++) {
// displayState(&s0);
evolve(&s0, &s1);
// usleep(100000);
// displayState(&s1);
evolve(&s1, &s0);
// usleep(100000);
}
return 0;
}

And here is an example of its output, showing the classical Game of Life glider:

                    






* *
**
*


A program to fit straight lines to data

#include <math.h>
#include <stdio.h>
#include <assert.h>

/*
Michael Ashley / UNSW / 06-Jun-2004
*/

void linfit(double x[], double y[], double sy[], int n,
double *a, double *b, double *sa, double *sb) {

/*
Fits a straight line y = a + bx to the n data pairs (x,y).
sy is the standard deviation of y. The standard deviations
of a and b are also returned.
*/

double S=0.0, Sx=0.0, Sy=0.0, Stt=0.0, Sty=0.0;
double t;
int i;

assert(n > 1);

for (i = 0; i < n; i++) {
double ss;
ss = sy[i]*sy[i];
S += 1/ss;
Sx += x[i]/ss;
Sy += y[i]/ss;
}

for (i = 0; i < n; i++) {
t = (x[i]-Sx/S)/sy[i];
Stt += t*t;
Sty += t*y[i]/sy[i];
}

*b = Sty / Stt;
*a = (Sy - Sx* *b)/S;
*sa = sqrt((1+(Sx*Sx)/(S*Stt))/S);
*sb = sqrt(1/Stt);
}

int main(void) {
double x[3]={1.0, 2.0, 3.0};
double y[3]={1.1, 2.3, 3.5};
double sy[3]={0.01, 0.01, 0.01};

double a, b, sa, sb;

linfit(x, y, sy, (sizeof x)/(sizeof x[0]), &a, &b, &sa, &sb);

printf("y = (%lf+/-%lf) + (%lf+/-%lf) * x\n", a, sa, b, sb);
return 0;
}

A function for generating quasi-random Sobel sequences

#include <stdio.h>
#include <assert.h>

#define MAXBIT 30

void sobseq(float *x, float *y) {

// Returns two quasi-random numbers for a 2-dimensional Sobel
// sequence. Adapted from Numerical Recipies in C.

int j, k, l;
unsigned long i, im,ipp;

// The following variables are "static" since we want their
// values to remain stored after the function returns. These
// values represent the state of the quasi-random number generator.

static float fac;
static int init=0;
static unsigned long in, ix[3], *iu[2*MAXBIT+1];
static unsigned long mdeg[3]={0,1,2};
static unsigned long ip[3]={0,0,1};
static unsigned long iv[2*MAXBIT+1]=
{0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};

// Initialise the generator the first time the function is called.

if (!init) {
init = 1;
for (j = 1, k = 0; j <= MAXBIT; j++, k += 2) iu[j] = &iv[k];
for (k = 1; k <= 2; k++) {
for (j = 1; j <= mdeg[k]; j++) iu[j][k] <<= (MAXBIT-j);
for (j = mdeg[k]+1; j <= MAXBIT; j++) {
ipp = ip[k];
i = iu[j-mdeg[k]][k];
i ^= (i >> mdeg[k]);
for (l = mdeg[k]-1; l >= 1; l--) {
if (ipp & 1) i ^= iu[j-l][k];
ipp >>= 1;
}
iu[j][k] = i;
}
}
fac = 1.0/(1L << MAXBIT);
in = 0;
}

// Now calculate the next pair of numbers in the 2-dimensional
// Sobel sequence.

im = in;
for (j = 1; j <= MAXBIT; j++) {
if (!(im & 1)) break;
im >>= 1;
}
assert(j <= MAXBIT);
im = (j-1)*2;
ix[1] ^= iv[im+1];
ix[2] ^= iv[im+2];
*x = ix[1] * fac;
*y = ix[2] * fac;
in++;
}


int main(void) {

// Test the Sobel generator.

int i;
float x, y;

for (i = 0; i < 1000; i++) {
sobseq(&x, &y);
printf("%f %f\n", x, y);
}
return 0;
}

A program to fit straight lines to data

#include <math.h>
#include <stdio.h>
#include <assert.h>

/*
Michael Ashley / UNSW / 06-Jun-2004
*/

void linfit(double x[], double y[], double sy[], int n,
double *a, double *b, double *sa, double *sb) {

/*
Fits a straight line y = a + bx to the n data pairs (x,y).
sy is the standard deviation of y. The standard deviations
of a and b are also returned.
*/

double S=0.0, Sx=0.0, Sy=0.0, Stt=0.0, Sty=0.0;
double t;
int i;

assert(n > 1);
for (i = 0; i < n; i++) {
double ss;
ss = sy[i]*sy[i];
S += 1/ss;
Sx += x[i]/ss;
Sy += y[i]/ss;
}

for (i = 0; i < n; i++) {
t = (x[i]-Sx/S)/sy[i];
Stt += t*t;
Sty += t*y[i]/sy[i];
}

*b = Sty / Stt;
*a = (Sy - Sx* *b)/S;
*sa = sqrt((1+(Sx*Sx)/(S*Stt))/S);
*sb = sqrt(1/Stt);
}

int main(void) {
double x[3]={1.0, 2.0, 3.0};
double y[3]={1.1, 2.3, 3.5};
double sy[3]={0.01, 0.01, 0.01};

double a, b, sa, sb;

linfit(x, y, sy, (sizeof x)/(sizeof x[0]), &a, &b, &sa, &sb);

printf("y = (%lf+/-%lf) + (%lf+/-%lf) * x\n", a, sa, b, sb);
return 0;
}