ran.c

蒙特卡罗模拟光子成像C语言版,代码简洁专业

/*
 * This file is part of tMCimg.
 * 
 * tMCimg is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option)
 * any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "config.h"

#if (HAVE_ASSERT_H)
#include <assert.h>
#endif

#if (HAVE_GSL)
# include "gsl/gsl_rng.h"
# include <string.h>

/*
 * The period of "mt19937" is absolutely ludicrous (order 10^6000!),
 * and it's fast too.  "ranlxs0" is a good choice as well (with a period
 * of about 10^171)
 */

# define DEFAULT_RNG "mt19937"

static const gsl_rng_type *rng_type  = NULL;
static const gsl_rng      *rng_state = NULL;

static void seed_rng(unsigned long);

float ran(int *idum)
{
  static int seeded = 0;

  if (seeded == 0)
    {
      seed_rng((unsigned long)idum[0]);
      seeded = 1;
    }

  return (float)gsl_rng_uniform(rng_state);
}
   
static void seed_rng(unsigned long seed)
{
#if 0
  /* For the future */
  extern const gsl_rng_type *gsl_rng_default;
  extern unsigned long int gsl_rng_default_seed;

  if (getenv(GSL_RNG_TYPE) != NULL)
    {
      gsl_rng_env_setup();
      rng_type = gsl_rng_default;
    }
#endif

  if (rng_type == NULL)
    {
      const gsl_rng_type **t0, **t;

      t0 = gsl_rng_types_setup();

      for (t = t0; *t != NULL; t++)
	if (strcmp((*t)->name, DEFAULT_RNG) == 0)
	  {
	    rng_type = *t;
	    break;
	  }

      if (*t == NULL)
	{
	  fprintf(stderr,
		  "GSL random number generator %s not found\n", DEFAULT_RNG);
	  exit(1);
	}
    }

  /* Allocate storage for state information */

  if ((rng_state = gsl_rng_alloc(rng_type)) == NULL)
    exit(1);

#if (DEBUG)
  printf("Using GSL random number generator %s\n", gsl_rng_name(rng_state));
#endif

  /* Seed the RNG */

  gsl_rng_set(rng_state, seed);
}

#else                  /* HAVE_GSL */

/* ****************************************************************** */
/* GSL routines not available, use our own PRNG.  This one is pretty  */
/*  good actually, with a period of around 2^{191}, which is much     */
/*  longer than the stock PRNG's in Numerical Recipies.               */
/*                                                                    */
/* From "Pierre L'Ecuyer, `Good parameters and implementations for    */
/*   combined multiple recurxsive random number generators', 1998.    */
/* ****************************************************************** */

static double MRG32k3a(void);
static double s10, s11, s12, s20, s21, s22;

#define NORM 2.328306549295728e-10
#define M1   4294967087.0
#define M2   4294944443.0
#define A12     1403580.0
#define A13N     810728.0
#define A21      527612.0
#define A23N    1370589.0

float ran(int *idum)
{
  static int seeded = 0;
  double r;

  if (seeded == 0)
    {
      unsigned int seed;

      assert(idum != NULL);

      seed = ((unsigned int *)idum)[0];

      srand(seed);

      s10 = M1 * rand() / (double)RAND_MAX;
      s11 = M1 * rand() / (double)RAND_MAX;
      s12 = M1 * rand() / (double)RAND_MAX;

      s20 = M1 * rand() / (double)RAND_MAX;
      s21 = M1 * rand() / (double)RAND_MAX;
      s22 = M1 * rand() / (double)RAND_MAX;

      seeded = 1;
    }

  r = MRG32k3a();

  assert((r >= 0.0) && (r < 1.0));

  return (float)r;
}

static double MRG32k3a(void)
{
  long k;
  double p1, p2;

  /* Component 1 */
  p1 = A12 * s11 - A13N * s10;
  k  = (long)(p1 / M1);

  p1 -= k * M1;

  if (p1 < 0.0)
    p1 += M1;

  s10 = s11;
  s11 = s12;
  s12 = p1;

  /* Component 2 */

  p2 = A21 * s22 - A23N * s20;
  k  = (long)(p2 / M2);

  p2 -= k * M2;

  if (p2 < 0.0)
    p2 += M2;

  s20 = s21;
  s21 = s22;
  s22 = p2;

  /* Combination */

  if (p1 <= p2)
    return ((p1 - p2 + M1) * NORM);
  else
    return ((p1 - p2) * NORM);
}
#endif