C/C++ create random noise (gaussian noise/white noise)

I was amazed when use randn command at Matlab. randn command will generate random data every we call that command. After I search at google, I found how to make this happen. I get this code at seismic unix source code. This code will generate random noise or white noise with Gaussian method. Code for main.c is :

#include <stdio.h>
#include <stdlib.h>
#include "frannor.h"

int main(int argc, char **argv)
{
	double *noise=NULL;
	int i;
	int ndata;
	unsigned int seed;

	if(argc!=2)
	{
		fprintf(stderr, "usage : random 5");
		exit(0);
	}

	ndata = atoi(argv[1]);

	if (-1 == (seed = (unsigned int) time((time_t *) NULL)))
	{
		fprintf(stderr, "time() failed to set seed");
		exit(0);
	}

	srannor(seed);  //seed random number generator
	noise = (double*) calloc (ndata, sizeof(double));
	for(i=0; i<ndata; i++)  //create random data
	{
		/* Compute noise vector elements in [-1, 1] */
		/* GAUSS METHOD. frannor gives elements in N(0,1)--ie. pos & negs */
		noise[i] = (double) frannor();
	}

	for(i=0; i<ndata; i++)
		printf("%3.2f   ", noise[i]);
	printf("\n");

	free(noise);  //free allocated memory
	return(1);
}

This is code for frannor.h :

#ifndef FRANNOR_H_
#define FRANNOR_H_

#include <math.h>
#define ABS(x) ((x) < 0 ? -(x) : (x))
#define AA 12.37586
#define Bsu 0.4878992
#define Csu 12.67706
#define C1 0.9689279
#define C2 1.301198
#define PC 0.01958303
#define XN 2.776994
#define OXN 0.3601016
#define NBITS 24

/* macro to generate a random number uni uniform on [0,1) */
#define UNI(uni) \
	uni = u[i]-u[j]; if (uni<0.0) uni += 1.0; u[i] = uni; \
	if (--i<0) i = 16;  if (--j<0) j = 16

float frannor(void);
void srannor (int seed);

#endif /* FRANNOR_H_ */

This is code for frannor.c :

/* Copyright (c) Colorado School of Mines, 2010.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
FRANNOR - functions to generate a pseudo-random float normally distributed
		with N(0,1); i.e., with zero mean and unit variance.

frannor		return a normally distributed random float
srannor		seed random number generator for normal distribution

******************************************************************************
Function Prototypes:
float frannor (void);
void srannor (int seed);

******************************************************************************
frannor:
Input:		(none)
Returned:	normally distributed random float

srannor:
Input:
seed		different seeds yield different sequences of random numbers.

******************************************************************************
Notes:
Adapted from subroutine rnor in Kahaner,  Moler and Nash (1988)
which in turn was based on  an algorithm by
Marsaglia and Tsang (1984).

******************************************************************************
References:
"Numerical Methods and Software", D.
Kahaner, C. Moler, S. Nash, Prentice Hall, 1988.

Marsaglia G. and Tsang, W. W., 1984,
A fast, easily implemented method for sampling from decreasing or symmetric
unimodal density functions:  SIAM J. Sci. Stat. Comput., v. 5, no. 2,
p. 349-359.

*****************************************************************************
Author:  Dave Hale, Colorado School of Mines, 01/21/89
*****************************************************************************/
/**************** end self doc ********************************/

#include "frannor.h"

static float v[]={
	0.3409450, 0.4573146, 0.5397793, 0.6062427, 0.6631691,
	0.7136975, 0.7596125, 0.8020356, 0.8417227, 0.8792102, 0.9148948,
	0.9490791, 0.9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917,
	1.1332738, 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635,
	1.2982650, 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929,
	1.4569915, 1.4834526, 1.5100121, 1.5367061, 1.5635712, 1.5906454,
	1.6179680, 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422,
	1.7896223, 1.8200099, 1.8510770, 1.8829044, 1.9155830, 1.9492166,
	1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713,
	2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117,
	2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943
};

/* internal state variables for uniform random number generator */
static  int i=16,j=4;
static  float u[]={
	0.8668672834288,  0.3697986366357,  0.8008968294805,
	0.4173889774680,  0.8254561579836,  0.9640965269077,
	0.4508667414265,  0.6451309529668,  0.1645456024730,
	0.2787901807898,  0.06761531340295, 0.9663226330820,
	0.01963343943798, 0.02947398211399, 0.1636231515294,
	0.3976343250467,  0.2631008574685
};

float frannor(void)
/**************************************************************************
return a normally distributed random float
 ***************************************************************************
Returned:	normally distributed random float
 **************************************************************************/
{
	int k;
	float uni,vni,rnor,x,y,s,bmbx,xnmx;

	/* uni is uniform on [0,1) */
	UNI(uni);

	/* vni is uniform on [-1,1) */
	vni = uni+uni-1.0;

	/* k is in range [0,63] */
	k = ((int)(u[i]*128))%64;

	/* fast part */
	rnor = vni*v[k+1];
	if (ABS(rnor)<=v[k]) return rnor;

	/* slow part */
	x = (ABS(rnor)-v[k])/(v[k+1]-v[k]);
	UNI(y);
	s = x+y;
	if (s<=C2) {
		if (s<=C1) return rnor;
		bmbx = Bsu-Bsu*x;
		if (y<=Csu-AA*exp(-0.5*bmbx*bmbx)) {
			if (exp(-0.5*v[k+1]*v[k+1])+y*PC/v[k+1] <=
					exp(-0.5*rnor*rnor)) return rnor;
			do {
				UNI(y);
				x = OXN*log(y);
				UNI(y);
			} while (-2.0*log(y)<=x*x);
			xnmx = XN-x;
			return (rnor>=0.0 ? ABS(xnmx) : -ABS(xnmx));
		}
	}
	bmbx = Bsu-Bsu*x;
	return (rnor>=0.0 ? ABS(bmbx) : -ABS(bmbx));
}

void srannor (int seed)
/*****************************************************************************
seed random number generator
 ******************************************************************************
Input:
seed		different seeds yield different sequences of random numbers.
 *****************************************************************************/
{
	int ii,jj,ia,ib,ic,id;
	float s,t;

	i = 16;
	j = 4;
	ia=ABS(seed)%32707;
	ib=1111;
	ic=1947;
	for (ii=0; ii<17; ii++) {
		s = 0.0;
		t = 0.5;
		for (jj=0; jj<64; jj++) {
			id = ic-ia;
			if (id<0) {
				id += 32707;
				s += t;
			}
			ia = ib;
			ib = ic;
			ic = id;
			t *= 0.5;
		}
		u[ii] = s;
	}
}

To compile this source code, use command :

gcc main.c frannor.c -lm -o random

This is output from program :

toto@toto-laptop:/home/toto/$ ./random 6
0.08   1.78   -1.35   -1.03   1.62   0.07
toto@toto-laptop:/home/toto/$ ./random 6
1.12   1.32   0.73   -0.34   0.25   -0.58

Source : http://www.cwp.mines.edu/cwpcodes/index.html