Makefile

SRC=myRunifConcise.c randFortran.f90 randWrapper.c

all:
	R CMD COMPILE FCFLAGS='-O2' $(SRC)
	R CMD SHLIB $(SRC) -o myRand.so

clean: 
	rm *.o *.so

myRnorm.R

#myOtherRunif.R

myRnorm <-
function( n, mean = 0, sd = 1 ){

  rnorm <- numeric(n)
  normRandom <- .Fortran( 'myRnorm', rnorm = as.double(rnorm), 
    n = as.integer(n), mean = as.double(mean), sd = as.double(sd) )$rnorm

  return( normRandom )

}

myRunif.c

//myRunif.c 
#include <R.h>
#include <Rinternals.h>

SEXP myRunif( SEXP n, SEXP min, SEXP max ){

  SEXP statsPackage;
  PROTECT(
    statsPackage = eval( lang2( install("getNamespace"),
      ScalarString(mkChar("stats")) ),
      R_GlobalEnv
    )   
  );  

  SEXP RCallBack;
  PROTECT( RCallBack = allocVector(LANGSXP, 4 )); 
  SETCAR( RCallBack,
    findFun( install("runif"), statsPackage )
  );  

  SETCADR( RCallBack, n );
  SET_TAG( CDR( RCallBack ), install("n") );

  SETCADDR( RCallBack, min );
  SET_TAG( CDDR( RCallBack ), install("min") );

  SETCADDDR( RCallBack, max );
  SET_TAG( CDR(CDDR( RCallBack )), install("max"));

  SEXP randoms;
  PROTECT(
    randoms = eval( RCallBack, statsPackage )
  );  

  UNPROTECT(3);
  return( randoms );

}

myRunif.R

#myRunif.R
myRunif <-
function( n, min = 0, max = 1 ){

  unifRandom <- .Call( 'myRunif', n, min, max )

  return( unifRandom )

}

myOtherRunif <-
function( n, min = 0, max = 1 ){

  runif <- numeric(n)
  unifRandom <- .Fortran( 'myOtherRunif', runif = runif, n = n, 
    min = min, max = max )$runif

  return( unifRandom )

}

myRunifConcise.c

//Alternate myRunif.c 
#include <R.h>
#include <Rinternals.h>

SEXP myRunif( SEXP n, SEXP min, SEXP max ){

  SEXP statsPackage;
  PROTECT(
    statsPackage = eval( lang2( install("getNamespace"),
      ScalarString(mkChar("stats")) ),
      R_GlobalEnv
    )   
  );  

  SEXP randoms;
  PROTECT(
    randoms = eval( lang4( install("runif"),
     n,
     min,
     max),
     statsPackage
  ));

  UNPROTECT(2); 


  return( randoms );

}

randFortran.f90

subroutine myRnorm(rnorm, n, mean, sd)
    implicit none
    
    integer::n, i
    double precision::rnorm(n), mean, sd
    double precision:: normrnd
    
    call rndstart()
    
    do i=1,n
		rnorm(i) = mean + sd * normrnd()
    end do
    
    call rndend()
    
end subroutine

subroutine myOtherRunif(runif, n, min, max)
    implicit none
    
    integer::n, i
    double precision:: runif(n), min, max
    double precision:: unifrnd
    
    call rndstart()
    
    do i=1,n
		runif(i) = min + (max - min) * unifrnd()
    end do
	
    call rndend()
    
end subroutine

randWrapper.c

#include <R.h>

//Need to know the names of the R C routines for this to work 
void F77_SUB(rndstart)(void) { GetRNGstate(); }
void F77_SUB(rndend)(void) { PutRNGstate(); }
double F77_SUB(normrnd)(void) { return norm_rand(); }
double F77_SUB(unifrnd)(void) { return unif_rand(); }

testit.R

#!/usr/bin/env Rscript

dyn.load('myRand.so')

source('myRunif.R')
source('myRnorm.R')

cat('myRunif [0,1]:\n')
print(myRunif(5))
cat('myRunif [0,3]:\n')
print(myRunif(5,,3))

cat('myOtherRunif [0,1]:\n')
print(myRunif(10))
cat('myOtherRunif [0,3]:\n')
print(myRunif(10,,3))

cat('myRnorm N(0,1):\n')
print(myRnorm(10))

cat('myRnorm N(10,.5):\n')
print(myRnorm(10,10,.5))

zzz.R

#zzz.R
.First.lib <-
function( libname, pkgname ){

  library.dynam( pkgname, pkgname, libname )

}