t distribution and GSL random number generator

t distribution with large degrees of freedom is the same as Normal (Gaussian) distribution with mean 0 and standard devitation 1 approximately.

• normal[x_] := PDF[NormalDistribution[0, 1], x]
• t[x_, nu_] := PDF[StudentTDistribution[nu]], x]
• Plot[{normal[x], t[x, 3], t[x, 6], t[x, 20]}, {x, -6, 6},  PlotStyle -> {Black, Gray, Blue, Red}]

• Five real random numbers x0 from the Gaussian distribution (0,1):

• Five real random numbers x1 from the Gaussian distribution (mu,sigma) :

• There is the following relation : x1 = x0*sigma + mu

Therefore, from a real random number from t distribution random number generator with nu, I can get the random magnetic field strength B_t as following:

,where t_random is generated by gsl_ran_tdist function as

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>   
double gsl_ran_tdist (const gsl_rng * r, double nu) ,where nu is the degrees of freedom. 

double random_t_distribution(double ndf)
{
  unsigned long random_seed   = 0;
  double        random_number = 0.0;
 
  const gsl_rng_type *T;
  gsl_rng *r;


  T = gsl_rng_default; // gsl_rng+mt19937l
  r = gsl_rng_alloc (T);
  random_seed = random_seed_number();
  gsl_rng_set(r, random_seed);

  random_number = gsl_ran_tdist(r, ndf);
  gsl_rng_free (r);
  return random_number;
}

In the random_t_distribution function, I used a sort of random number generator in order to get a random seed. The original random number generator, which is based on /dev/random, is shown in Advanced Linux Programming Chapter 6. I modified it as following:

#include <stdio.h>
#include <assert.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <fcntl.h>
#include <unistd.h>
unsigned long random_seed_number()
{
  /* Store a file descriptor opened to /dev/random in a static
     variable. That way, we don't need to open the file every time
     this function is called. */
  static int dev_random_fd = -1;
  char* next_random_byte;
  int bytes_to_read;
  unsigned long random_value;

  if (dev_random_fd == -1)
    {
      dev_random_fd = open ("/dev/random", O_RDONLY);
      assert (dev_random_fd != -1);
    }

  /* Read enough random bytes to fill an integer variable. */
  next_random_byte = (char*) &random_value;
  bytes_to_read    = sizeof (random_value);

  /* Loop until we've read enough bytes. Because /dev/random is filled
     from user-generated actions, the read may block and may only
     return a single random byte at a time. */
  do
    {
      int bytes_read;
      bytes_read = read (dev_random_fd, next_random_byte, bytes_to_read);
      bytes_to_read -= bytes_read;
      next_random_byte += bytes_read;
    } while (bytes_to_read > 0);
   return random_value;
}

Secure Programming Cookbook for C and C++ (11.3 Using the Standard Unix Randomness Infrastructure) is also useful reference book on the web.