C code discussed in class

This page summarizes C code discussed in class.

Our first C program tried to answer question 2 from the first assignment, which asked us to find all integers 0 \leq w \leq x \leq y < z < 10^6 such that w^4 + x^4 + y^4 = z^4.

The first difficulty we faced was that even with the unsigned long int data type, the maximum integer value we can store exactly is 2^{64} - 1, whereas we need to store up to 10^{6 \times 4} = 10^{24} \approx 2^{80}. (Here is a short description of data types. Note the lack of an exact specification of the amount of storage allocated for each type.)

As a first attempt, we tried to solve the slightly easier problem of finding all integers 0 \leq w \leq x \leq y < z < N such that w^m + x^m + y^m = z^m for smaller values of N and m. We came up with the following program.

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

/* This is a macro */
#define npower(x) ((x) * (x) * (x)) 	/* m = 3 */
#define M 10000

/* This is a function declaration.  The actual function is defined below main() */
double invpower(double x);

/* This is a global variable */
double p;

/* The main() function gets executed when a C program is run */
int main()
{
    unsigned long int x, y, z, w, k;
    double zroot;
    double p = 1.0 / 3.0;	/* assignment to global variable */

    for (x = 1; x <= M; x++)
	for (y = x; y<= M; y++)
	    for (w = y; w <= M; w++) {
		k = (npower(x) + npower(y) + npower(w));
		zroot = invpower((double) k);
		z = (int) zroot;
		if (npower(z) == k) {
		    printf("[A] %lu,%lu,%lu : %lu, %lu\n", 
			   x, y, w, k, z);
		}
		else if (npower(z+1) == k) {
		    printf("[B] %lu,%lu,%lu : %lu, %lu\n", 
			   x, y, w, k, z+1);
		}
	    }
    return 0;
}

double invpower(double x)
{
    return pow(x, p);
}

  

We also discussed a variant which defined invpower() as a nested function within main(). Nested functions are not part of the official C standard(s), but is supported by the GCC compiler as a language extension. The following works if we compile it with

$ gcc -lm -o euler euler.c

but we get a warning if we compile it with

$ gcc -lm -o euler -ansi -pedantic euler.c

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

#define npower(x) ((x) * (x) * (x)) 	/* m = 2 */
#define M 10000

int main()
{
    unsigned long int x, y, z, w, k;
    double zroot;

    double p = 1.0 / 3.0;

    double invpower(double x)
    {
	return pow(x, p);
    }

    for(x = 1; x <= M; x++)
	for(y = x; y<= M; y++)
	    for(w = y; w <= M; w++){
		k = (npower(x) + npower(y) + npower(w));
		zroot = invpower((double) k);
		z = (int) zroot;
		if (npower(z) == k) {
		    printf("[A] %lu,%lu,%lu : %lu, %lu\n", 
			   x, y, w, k, z);
		}
		else if (npower(z+1) == k) {
		    printf("[B] %lu,%lu,%lu : %lu, %lu\n", 
			   x, y, w, k, z+1);
		}
	    }
    return 0;
}

  

The second example we discussed was question 3 in the first assignment, which asked us to generate exponential random variates (with and without truncation). Our first attempt using drand48() failed, but the next attempt using rand() worked (after making a minor modification to the rng_unif() function defined in class).

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

#define N 100

#define one_minus_einv 0.632120558828558

double rng_unif() 
{
    return ((double) rand()) / RAND_MAX;
}

double rng1(double u) 
{ 
    return -log(u);
}

double rng2(double u) 
{ 
    return 1 - log(u);
}

double rng3(double u) 
{ 
    return -log(1 - u * one_minus_einv);
}

int main()
{
    int i;
    double a1[N], a2[N], a3[N];

    /* unsigned short rs[3]; */
    /* rs[0] = 234; rs[1] = 100; rs[2] = 15; */
    /* seed48(rs); */

    srand(22123);

    for (i = 0; i < N; i++) {
	a1[i] = rng1(rng_unif());
	a2[i] = rng2(rng_unif());
	a3[i] = rng3(rng_unif());
    }
    for (i = 0; i < N; i++) {
	printf("%f,%f,%f\n", a1[i], a2[i], a3[i]);
    }
    return 0;
}

---

Here is the R code we used in class to draw Q-Q plots.

expdat <- read.table("randexp2.csv", header = FALSE, sep = ",")

plot(sort(expdat$V1), qexp(seq(0, 1, length = 10000)))
abline(0, 1)


qexp.trunc <- function(q) 
{
    -log(1 - q * (1 - exp(-1)))
}

plot(sort(expdat$V3), qexp.trunc(seq(0, 1, length = 10000)))
abline(0, 1)

Built-in R functions that are useful in this context are quantile(), ppoints(), and qqmath() in the lattice package.

Finally, here is Kaustav's program to generate standard normal random numbers, modified to use static variables so that random numbers are not wasted.

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

#define N 100

#ifndef M_PI
#define M_PI 3.14159265358979
#endif

double normal();

int main(void){
    int i;
    srand(1234567891);
    for(i = 0; i < N; i++)
	printf("%f\n",normal());
    return 0;
}

double normal()
{
    double x, y;
    static double z1, z2;
    static char new = 0;
    if (new == 0) {
	x = ((double)rand())/RAND_MAX;
	y = ((double)rand())/RAND_MAX;
	z1 = sqrt(-2 * log(x)) * sin(2 * M_PI * y);
	z2 = sqrt(-2 * log(x)) * cos(2 * M_PI * y);
    }
    printf("Status: %d\t%f\t%f\n", new, z1, z2);
    if (new == 0) {
	new = 1;
	return z1;
    }
    else {
	new = 0;
	return z2;
    }
}

---

This program modifies the previous program to use the uniform random number generator provided by the "standalone R math library" instead of the one in stdlib.h. This example also illustrates the use of pointers and dynamic memory allocation.

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

#define MATHLIB_STANDALONE
#include <Rmath.h>

#ifndef M_PI
#define M_PI 3.14159265358979
#endif

double one_normal();
double my_runif();

void normal(int n, double *x);

int main(void) {
    int i, N = 10;
    double *z;
    z = (double*) malloc(N * sizeof(double));
    normal(N, z);
    for(i = 0; i < N; i++)
    	printf("%f\n", z[i]);
    free((void *)z);
    return 0;
}

double my_runif()
{
    /* return ((double)rand())/RAND_MAX; */
    return runif(0, 1);
}

void normal(int n, double *x)
{
    int i;
    for (i = 0; i < n; i++) {
	x[i] = one_normal();
    }
    return;
}


double one_normal()
{
    double x, y;
    static double z1, z2;
    static char new = 0;
    if (new == 0) {
	x = my_runif();
	y = my_runif();
	z1 = sqrt(-2 * log(x)) * sin(2 * M_PI * y);
	z2 = sqrt(-2 * log(x)) * cos(2 * M_PI * y);
    }
    if (new == 0) {
	new = 1;
	return z1;
    }
    else {
	new = 0;
	return z2;
    }
}

To compile this program, one would do something like

gcc -lRmath -lm -o rand_normal_rmathlib rand_normal_rmathlib.c

The next modification is aimed at getting direct access to our random numbers from inside R. The way this works is that we compile our code (which contains functions we can call from R) into a shared library (or DLL), which can then be loaded and used from within R.

#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265358979
#endif

double one_normal();
double my_runif();

SEXP normal(SEXP n) 
{
    int i, N = INTEGER_POINTER(n)[0];
    SEXP ans = PROTECT(NEW_NUMERIC(N));
    double *z = NUMERIC_POINTER(ans);
    for(i = 0; i < N; i++)
    	z[i] = one_normal();
    UNPROTECT(1);
    return ans;
}

double my_runif()
{
    return ((double)rand())/RAND_MAX;
}


double one_normal()
{
    double x, y;
    static double z1, z2;
    static char new = 0;
    if (new == 0) {
	x = my_runif();
	y = my_runif();
	z1 = sqrt(-2 * log(x)) * sin(2 * M_PI * y);
	z2 = sqrt(-2 * log(x)) * cos(2 * M_PI * y);
    }
    if (new == 0) {
	new = 1;
	return z1;
    }
    else {
	new = 0;
	return z2;
    }
}

To compile to a shared library, we use (assuming the code is saved as file r_ext.c)

R CMD SHLIB r_ext.c

which created the file r_ext.so. We can now load and use it from inside R; for example,

dyn.load("r_ext.so")
zz <- .Call("normal", 10)
qqnorm(zz)

See the "Writing R Extensions" manual for more details.