#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("%ft%fn", 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;
}