gsl_ran_gaussian函数改写

/*
int seed=1000; // 种子
srand(seed); // 初始化随机数生成器

( (double)rand() / ((double)(RAND_MAX)+(double)(1)) )
产生一个在(0, 1)区间上的双精度随机数;0和1都不包括

Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */

double gsl_ran_gaussian (const double sigma)
{
  double x, y, r2;

  // r 是双精度浮点随机数,范围0.0到1.0 
  double r = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) ); 

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
      x = -1 + 2 * ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );//
      y = -1 + 2 * ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );

      /* see if it is in the unit circle */
      r2 = x * x + y * y;
    }
  while (r2 > 1.0 || r2 == 0);

  /* Box-Muller transform */
  return sigma * y * sqrt (-2.0 * log (r2) / r2);
}