## Usage:  X = compare_rngs ([ N])
##
## Time the computation of a bunch of N(0,1) random numbers using the
## default RNG from RANDLIB, and the Mersenne Twister, and plot the result
##
## If not argument specified, compute at N = logspace(3,6,50);

function X = compare_rngs ( N)

  if (nargin < 1)
    N = logspace(3,6,50);
    printf("Computing at N = logspace(3,6,50) which might take a while ...\n");
    fflush(stdout);
  endif

  for i = 1:length(N)

    for n = 1:5

      t = cputime();
      E = randn(N(i),1);
      Y(n,1) = cputime() - t;

      t = cputime();
      E = randmtn(N(i),1);
      Y(n,2) = cputime() - t;

    endfor

    X(i, :) = mean(Y(3:5,:));

  endfor

  gset key top left
  xlabel("random numbers generated")
  ylabel("time in seconds")
  title("Cost of random number generation")    
  plot(N,X(:,1),";default generator;",N,X(:,2),";mersenne twister;")

endfunction