Calculate a Gaussian distribution with zero mean and standard deviation I'm have to implement Gaussian Normal distribution in a programming language and now I'm trying to understand what is a Gaussian Normal distribution.
I have read its Wikipedia entry, but sorry, there are a lot of formula and a lot of things that I don't understand.
If I want to calculate a Gaussian distribution with zero mean and standard deviation $\sigma$, or $N(0,\sigma)$, do I need to implement the probability density function?
$$p(\Delta x_i) = \frac{1}{\sqrt{2\pi\sigma}}e^{-(\Delta x_i)^2/2\sigma^2}$$
NOTE:I have copied the above formula from the book Introduction to Evolutionary Computing which is the one I'm studying.
In other words: I need to calculate $x_i^{\prime} = x_i + N(0, \sigma)$ but I don't know how to calculate $N(0, \sigma)$.
I'm going to use C# to implement it and I don't have any other resources to do it, just Microsoft's .Net Framework.
 A: First I will sketch out a quick "Physicists Proof" of the method. This approach was taught to me by the Astrophysicist James Fry. 
To generate samples of a random variable $X$ which has a cumulative distribution function $F(x)$ you apply the inverse of $F$ to samples from a uniform random variable on the interval from $[0,1]$. 
This is simple enough, but the cumulative distribution function for a normal random variable doesn't have a simple closed form. If we want to stick with elementary functions, which we are likely to have well implemented in our libraries, we should avoid directly computing this inverse. What we do instead is re-express the problem in terms of a pair of random variables whose distribution functions are more easily inverted in terms of elementary functions. 
Let $X$ and $Y$ be independent normal random variables which obey $N(0,\sigma)$. These have a cumulative distribution function, 
$$ F(x,y) = \int_{-\infty}^{x} \int_{-\infty}^{y} dx'
 dy'
 \ \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{x'^2}{2\sigma^2} } \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{y'^2}{2\sigma^2} } $$
now we are going to change variables to polar coordinates, but I am not going to pay too much attention to the specific bounds. Recall that $x=r\cos(\theta)$ and $y=r\sin(\theta)$ and $dx\ dy = rdr\ d\theta$
$$ F(r,\theta) = \int_{0}^{r} \int_{0}^{\theta} \frac{r'dr'
 d\theta'}{2\pi\sigma^2} 
\  e^{-\frac{r'^2}{2\sigma^2} }$$
this integral can be factored, 
$$ F(r,\theta) = \int_{0}^{r}  \frac{r'dr'
 d\theta'}{\sigma^2} 
\  e^{-\frac{r'^2}{2\sigma^2} }\int_{0}^{\theta} \frac{d\theta'}{2\pi}$$
the $r$ integral can be modified by a change of variables $t=r^2/2\sigma^2$. 
$$ F(t,\theta) = \int_{0}^{t}  dt'
\  e^{-t'} \int_{0}^{\theta} \frac{d\theta'}{2\pi}$$
Each integral represents the cumulative distribution function of a random variable ($T$ and $\Theta$). These distribution functions are easy to compute.
$$F_T(t) = 1-e^{-t} \qquad F_\Theta(\theta) = \frac{\theta}{2\pi}$$  
Each of these functions can easily be inverted, 
$$F_T^{-1}(t) = -\ln(1-t) \qquad F_\Theta^{-1}(\theta) = 2\pi\theta,$$
so we can easily generate samples for the random variables $T$ and $\Theta$. These samples are then used to compute $r$ and subsequently $x$ and $y$. 

The procedure is as follows. 
1) By some means, generate a pair of random numbers which obey an uniform distribution on the interval [0,1]. Call these numbers $u_1$ and $u_2$. 
2) Compute numbers $t$ and $\theta$ by, $t=-\ln(1-u_1)$ and $\theta=2\pi u_2$. 
3) Compute a number $r$ by, $r=\sqrt{2\pi\sigma^2  t} $. 
4) Compute a pair of numbers $x$ and $y$ by, $x=r\cos(\theta)$ and $y=r\sin(\theta)$. 
5) $x$ and $y$ are samples from a $N(0,\sigma)$ distribution. If you want samples from a $N(\mu,\sigma)$ distribution then add $\mu$ to each of $x$ and $y$. 

Below is code I have in C for generating these numbers. To compile this code with a GNU C compile simply type "gcc name_of_file_with_code.c -lm".
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <complex.h>

double rUniform(double, double);
double complex rNormal(double, double);

const int ITMAX=1000;

int main(){
  srand(time(NULL));

  for( int i=0; i < ITMAX; i++){
    complex double z = rNormal(5,1); 
    printf("%f\t%f\n", creal(z),cimag(z));
  }
}

double rUniform(double a, double b){
  return rand()*(b-a)/RAND_MAX+a;
}

double complex rNormal(double mu, double sigma){
  // Generate a pair of random numbers which are
  // sampled from a uniform distribution on the
  // interval [0,1]. 
  double u1 = rUniform(0,1);
  double u2 = rUniform(0,1); 

  double pi = acos(-1); // Lazy way to get pi=3.14... from arccos(-1)

  // Generte samples for the variables Theta and T by
  // applying the inverses of the cumulative distribution
  // functions to uniform random variables. 
  double theta = 2*pi*u1;
  double t = -log(1-u2);


  // Generate r from t, recall that t=r^2/(2*sigma^2)
  // so we can get r from sqrt(2*sigma^2*t)
  double r = sqrt(2*sigma*sigma*t);

  // Generate the cartesian coordinates x and y from the polar
  // coordinates r and theta. Recall that these x and y are the
  // variables which actually obey N(mu,sigma).
  double x = r*cos(theta);
  double y = r*sin(theta);

  // We can return the pair in a complex variable z.
  double complex z =  (x+mu)+I*(y+mu);

  return z;  
}

