See the program hello.c in the notes.
See the program circle.c in the notes.
Most of the complexity in the solution below is in calculating the centroid. The main problems that people encountered with this assignment were (1) not initialising variables (such as those used to remember the min/max pixel values), (2) not using floating point numbers to calculate the centroid, (3) not using a 10x10 box, and (4) not worrying about the special case where the 10x10 box falls partially out of the image.
#include <stdio.h>
#include <unistd.h>
#include <assert.h>
// Michael Ashley / UNSW / 09-May-2004
const int xdim = 176;
const int ydim = 144;
int main(void) {
// A program to read the mass0.dat file, put the data into an
// unsigned char array, and calculate various things.
unsigned char image[ydim][xdim];
int i, j, n, imin, imax, jmin, jmax, min, max, minx, maxx, miny, maxy;
long sum, sumx, sumy;
double xcen, ycen;
// Read the data from standard input.
if (xdim*ydim != read(STDIN_FILENO, image, xdim*ydim)) {
fprintf(stderr, "read failure\n");
return 1;
}
// Find the min and max pixel coordinates and values.
// Begin by using the (0,0) pixel as both the min and max.
min = image[0][0];
max = image[0][0];
minx = 0;
miny = 0;
maxx = 0;
maxy = 0;
// Now see if we can do better...
for (j = 0; j < ydim; j++) {
for (i = 0; i < xdim; i++) {
n = image[j][i];
if (n < min) {
min = n;
minx = i;
miny = j;
}
if (n > max) {
max = n;
maxx = i;
maxy = j;
}
}
}
// Print the result. Note that if there are multiple min/max
// pixels, we will only find the first one.
printf("min of %d at (%d,%d)\n", min, minx, miny);
printf("max of %d at (%d,%d)\n", max, maxx, maxy);
// Now find the centroid in a 10x10 pixel box.
// Begin by finding the coodinates of the edges of the box, putting
// the maximum value close to the centre (which is tricky with a
// 10x10 box!).
imin = maxx - 5;
imax = maxx + 4;
jmin = maxy - 5;
jmax = maxy + 4;
// Now cope with the possibility that the box falls partially
// outside the image. We simply move the box so that one edge
// is on a boundary of the image, and leave the box as 10x10 in
// size. Note that this fails if the dimensions of the image are
// less than 10x10, so we assert that this isn't the case.
assert(xdim >= 10 && ydim >= 10);
if (imin < 0) {
imin = 0;
imax = imin + 9;
}
if (jmin < 0) {
jmin = 0;
jmax = jmin + 9;
}
if (imax >= xdim) {
imax = xdim;
imin = imax - 9;
}
if (jmax >= ydim) {
jmax = ydim;
jmin = jmax - 9;
}
// Initialise the sums for the centroid calculation. We use
// long data types to reduce the chance of overflows, while
// maintaining greater accuracy than a float or double.
sumx = 0;
sumy = 0;
sum = 0;
// Calculate the sums.
for (j = jmin; j <= jmax; j++) {
for (i = imin; i <= imax; i++) {
n = image[j][i];
sumx += i * n;
sumy += j * n;
sum += n;
}
}
// Calculate the centroid. Note that we allow for the pathological
// case of "sum" being zero, which can only happen if the entire image
// is zero.
if (sum != 0) {
xcen = ((double)sumx)/ sum;
ycen = ((double)sumy)/ sum;
} else {
xcen = maxx;
ycen = maxy;
}
printf("centroid at (%.2lf,%.2lf)\n", xcen, ycen);
return 0;
}
Points to note in the solution below: (1) the use of M_PI from <math.h> for the value of pi, and (2) the use of a logarithmically increasing number of iterations.
// A program to integrate sin(x) from 0 to pi/2 using Monte Carlo
// integration.
// Michael Ashley / UNSW / 12-Apr-2003
#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
int main(void) {
struct timeval t;
double x, y;
int i;
long in, total;
// 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);
in = 0;
total = 0;
while (1) {
// Generate a random number between 0 and pi/2.
x = (M_PI/2.0)*random()/((double) RAND_MAX);
// Generate a random number between 0 and 1.
y = random()/((double) RAND_MAX);
// Compare the random coordinate (x,y) with (x,sin(x)).
if (sin(x) > y) {
in++;
}
total++;
// Print out the approximation to the integral for various powers of
// ten.
if (total == 1 ||
total == 10 ||
total == 100 ||
total == 1000 ||
total == 10000 ||
total == 100000 ||
total == 1000000 ||
total == 10000000 ||
total == 100000000L) {
printf("%d, %f\n", total, ((M_PI/2.0)*in)/total);
}
if (total > 100000000L) break;
}
return 0;
}
and here is the output (the calculation takes about 45 seconds on a 850MHz Pentium III):
1, 1.570796
10, 0.942478
100, 0.973894
1000, 1.011593
10000, 0.990230
100000, 1.000409
1000000, 0.999669
10000000, 0.999943
100000000, 1.000090
Notes: the solution to this week's assignment was much simpler than most people thought. Note that the program does not need to know the value of the gravitational constant, or the mass of the sun or planet. This shows the value of thinking through the mathematics before starting to code.
// Michael Ashley 26-May-2004
#include <stdio.h>
#include <assert.h>
#include <math.h>
double maxVel(double a, double e, double P) {
// Returns the maximum velocity attained in a circular or elliptical orbit
// with semi-major axis a, eccentricity e, and period P. All quantities
// are in SI units.
assert(e >= 0.0 && e < 1.0); // Terminate if parabolic or hyperbolic.
assert(a > 0.0 && P > 0.0); // Sanity checks on a and P.
return sqrt(4*M_PI*M_PI*a*a*(1.0+e)/(P*P*(1.0-e)));
}
int main(void) {
printf("earth's maximum velocity = %e m/s\n",
maxVel(149.5985E9, 0.016718, 31556926.0));
return 0;
}
Notes: making a "bullet-proof", and accurate, quadratic equation solver is more difficult than you might think, as this example shows. Whenever you divide by a number, or take a sqrt, or logarithm, etc, ALWAYS check whether an illegal operation is being attempted (e.g., division by zero, square-root of a negative number). If you look at the following code, you will see that every division and square-root is safe.
A common error that people made was not to use exponential notation in the print statements, so that more digits of precision were visible.
#include <stdio.h>
#include <math.h>
#include <assert.h>
int solveQuadratic(double a, double b, double c, double *x1, double *x2) {
// Solves the quadratic equation ax^2 + bx + c = 0, and returns
// the number of roots (0, 1, or 2), or -1 if there are infinitely
// many roots (i.e., if a = b = c = 0).
double det, q;
det = b*b - 4*a*c;
if (det >= 0.0) {
if (b > 0.0) {
q = -0.5 * (b + sqrt(det));
} else {
q = -0.5 * (b - sqrt(det));
}
if (a != 0.0) {
*x1 = q/a;
if (det != 0.0 && q != 0.0) {
*x2 = c/q;
return 2;
}
return 1;
} else if (q != 0.0) {
*x1 = c/q;
return 1;
}
if (c == 0.0) {
return -1;
}
}
return 0;
}
// Here are 7 test polynomials, deliberately chosen to be pathological.
double a[] = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0};
double b[] = {123456789.0, 2.0, 2.00000000001, 1.99999999999, 2.0, 0.0, 0.0};
double c[] = {-0.000000123456789, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0};
int main(void) {
int i, n;
double x1, x2;
// We are guarding against mistakes in the initialisation of the coefficients.
// This is perhaps overkill, but it is good practice and will save you
// a lot of time one day.
assert((sizeof a) == (sizeof b) && (sizeof b) == (sizeof c));
// Note the use of "sizeof" to avoid having to hard-wire the number of
// polynomials.
for (i = 0; i < (sizeof a) / (sizeof a[0]); i++) {
n = solveQuadratic(a[i], b[i], c[i], &x1, &x2);
switch (n) {
case -1:
printf("infinitely many roots\n");
break;
case 0:
printf("no roots\n");
break;
case 1:
printf("one root %.18le\n", x1);
break;
case 2:
printf("two roots %.18le, %.18le\n", x1, x2);
break;
}
}
return 0;
}
Very few people got the right answers.
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
double test[] = {
10000000000001.23,
10000000000001.27,
10000000000001.02,
10000000000001.16,
10000000000001.19,
10000000000002.98};
void stats(double *data, int n, double *mean, double *avDev, double *stDev) {
int i;
double diff, sum0, sum1, sum2;
assert(n > 0);
sum0 = 0.0;
for (i = 0; i < n; i++) {
sum0 += data[i];
}
*mean = sum0 / n;
sum0 = sum1 = sum2 = 0.0;
for (i = 0; i < n; i++) {
diff = (data[i] - *mean);
sum0 += diff;
sum1 += fabs(diff);
sum2 += diff * diff;
}
*avDev = sum1 / n;
if (n > 1) {
*stDev = sqrt((sum2 - (sum0*sum0)/(n*n)) / (n-1));
} else {
*stDev = 0.0;
}
}
void statsClip(double *data, int n, double sclip,
double *mean, double *avDev, double *stDev) {
int i, j;
double *dp;
// Calculate the statistics using all the data.
stats(data, n, mean, avDev, stDev);
// Allocate enough memory to store a copy of the data array.
assert(NULL != (dp = (double *)malloc(n * sizeof(double))));
// Now copy across those elements of the array that are within
// the sigma-clipping bounds.
j = 0;
for (i = 0; i < n; i++) {
if (sclip * *stDev > fabs(data[i] - *mean)) {
dp[j++] = data[i];
}
}
// And redo the statistics.
stats(dp, j-1, mean, avDev, stDev);
free(dp);
}
int main(void) {
double mean, avDev, stDev;
stats(test, (sizeof test)/(sizeof test[0]), &mean, &avDev, &stDev);
printf("mean = %lf avDev = %lf stDev = %lf\n", mean, avDev, stDev);
statsClip(test, (sizeof test)/(sizeof test[0]), 1.5, &mean, &avDev, &stDev);
printf("mean = %lf avDev = %lf stDev = %lf\n", mean, avDev, stDev);
return 0;
}