## 1. Approximating pi in C. 
## 2. dnorm in C. 
## 3. Sum of squared differences between observations in C. 

Reminder -- hw2 is due Wed Oct23. 

Reminder -- no lecture Thu Oct31. 

1. Approximate pi in C.

In mypi.c,
    #include <R.h>
    #include <Rmath.h>
    void pi2 (int *n, double *y){
      int i;
      double x[*n];
      x[0] = 1.0;
      y[0] = sqrt(6.0);
      for(i = 1; i < *n; i++) {
          x[i] = x[i-1] + 1.0 / ((i+1.0)*(i+1.0));
          /* or x[i] = x[i-1] + 1.0 / pow(i+1.0,2.0); */ 
          y[i] = sqrt(6.0 * x[i]);
      }
    }
      

In R,
    ## set working directory to the one containing mypi.c. 
    system("R CMD SHLIB mypi.c")
    dyn.load("mypi.so")
    pi3 = function(n){
		.C("pi2",as.integer(n), y = double(n))
	}
    b = pi3(1000000)
    b$y[1000000]

Note that you have to be incredibly careful in C when doing arithmetic between 
integers and non-integers.
If instead of 1.0 / ((i+1.0)*(i+1.0));
you do         1  / ((i+1.0)*(i+1.0));
or            1.0 / (i+1.0)^2;
crazy stuff might happen.
^ is a bitwise operator meaning "XOR", i.e. X^Y = 1 if X=1 or Y=1 but not both.
    
2. dnorm in C.
You can access C versions of many basic R functions, including for 
instance dnorm(), rnorm(), etc.
The syntax in C of dnorm is 
double dnorm(double x, double mu, double sigma, int give_log)
 
in mydn.c,
    #include <R.h>
    #include <Rmath.h>
    void norm2 (int *n, double *upper, double *bw, double *y){
      int i;
      double x, inc;
      x = -1.0 * *upper;
      inc = 2.0 * *upper / *n;
      for(i = 0; i < *n; i++) {
      y[i] = dnorm(x / *bw, 0,1,0);
      x += inc;
      }
    }
    
In R, 
    system("R CMD SHLIB mydn.c")
    dyn.load("mydn.so")
    norm3 = function(n, u, b){
		d = .C("norm2", as.integer(n), as.double(u), 
		as.double(b), y = double(n))
	d$y
	}
    b = 12.4
    n = 100000
    u = 5*b
    a = norm3(n,u,b)
    title2 = paste("normal density with sd ", as.character(b))
    plot(seq(-u,u,length=n), a, type="l", main=title2,xlab="x", ylab="f(x)")

3. Sum of squared differences between observations in C.

In sumsq.c,

	#include <R.h>
	#include <Rmath.h>

	void ss2 (double *x, int *n, double *y)
	/* x will be the vector of data of length n,
	   and y will be a vector of squared differences from obs i
	   to the other n-1 observations.
	*/
	{
		int i,j;
		double a;
		for(i = 0; i < *n; i++){
			a = 0.0;
			for(j=0; j < *n; j++){
				a += pow(x[i] - x[j], 2);
			}
			y[i] = a;
		}
	}

in R,
	system("R CMD SHLIB sumsq.c")
	dyn.load("sumsq.so")
	sum3 = function(data2){
		n = length(data2)
		a = .C("ss2", as.double(data2), as.integer(n),
			y=double(n))
		a$y  ## or equivalently a[[3]]
	}

	b = c(1,3,4)
        sum3(b)

	n = c(100, 1000, 2000, 3000, 5000, 7000, 8000, 10000)
	t2 = rep(0,8)
	for(i in 1:8){
	  b = runif(n[i])
	  timea = Sys.time()
	  d = sum3(b)
	  timeb = Sys.time()
	  t2[i] = timeb-timea
	  cat(n[i]," ")
	}
	par(mfrow=c(1,2))
	plot(n,t2,ylab="time (sec)",main="C")
	
	## Now try the same thing in R, without C.
	
	sum4 = function(data2){
	  n = length(data2)
	  x = rep(0,n)
	  for(i in 1:n){
	    for(j in 1:n){
	      x[i] = x[i] + (data2[i] - data2[j])^2
	      }
	    }
	  x
	  }

	b = c(1,3,4)
        sum4(b)

	n = c(1:8)*100
	t3 = rep(0,8)
	for(i in 1:8){
	  b = runif(n[i])
	  timea = Sys.time()
	  d = sum4(b)
	  timeb = Sys.time()
	  t3[i] = timeb-timea
	cat(n[i]," ")
	}
	plot(n,t3,ylab="time (sec)",main="R")