www.pudn.com > opengpssim.zip > boxmuller.c

/* boxmuller.c            
 
  Implements the Polar form of the Box-Muller Transformation 
 
  (c) Copyright 1994, Everett F. Carter Jr. 
      Permission is granted by the author to use 
      this software for any application provided this 
      copyright notice is preserved. 
 
*/ 
 
#include  
#include  
 
#define BOXMULLR_H 
 
#include "boxmuller.h" 
 
#undef BOXMULLR_H 
 
 
inline static float ranf( void) 
{ 
//  return (float)( rand()) / 32767.0; 
  return (float)( rand()) / RAND_MAX; 
} 
 
/*  
 *  normal random variate generator  
 *  mean m, standard deviation s  
 */ 
float randn( float m, float s)   
{                                    
  float x1, x2, w, y1; 
  static float y2; 
  static int use_last   = 0; 
  static float rand_max = (float)( RAND_MAX); 
 
  if ( use_last)            /* use value from previous call */ 
  { 
    y1 = y2; 
    use_last = 0; 
  } 
  else 
  { 
    do  
    { 
//      x1 = 2.0 * ranf() - 1.0; 
//      x2 = 2.0 * ranf() - 1.0; 
      x1 = 2.0 * (float)( rand()) / rand_max - 1.0; 
      x2 = 2.0 * (float)( rand()) / rand_max - 1.0; 
      w = x1 * x1 + x2 * x2; 
    } while ( w >= 1.0); 
 
    w  = sqrt( (-2.0 * log( w) ) / w ); 
    y1 = x1 * w; 
    y2 = x2 * w; 
    use_last = 1; 
  } 
 
  return ( m + y1 * s); 
} 
 
 
/*  
 *  normal random variate generator  
 *  mean m, standard deviation s  
 */ 
double randndbl( double m, double s)   
{                                    
  double x1, x2, w, y1; 
  static double y2; 
  static int use_last   = 0; 
  static double rand_max = (float)( RAND_MAX); 
 
  if ( use_last)            /* use value from previous call */ 
  { 
    y1 = y2; 
    use_last = 0; 
  } 
  else 
  { 
    do  
    { 
//      x1 = 2.0 * ranf() - 1.0; 
//      x2 = 2.0 * ranf() - 1.0; 
      x1 = 2.0 * (double)( rand()) / rand_max - 1.0; 
      x2 = 2.0 * (double)( rand()) / rand_max - 1.0; 
      w = x1 * x1 + x2 * x2; 
    } while ( w >= 1.0 ); 
 
    w = sqrt( (-2.0 * log( w ) ) / w ); 
    y1 = x1 * w; 
    y2 = x2 * w; 
    use_last = 1; 
  } 
 
  return( m + y1 * s ); 
} 
 
/* ============================== End of File ============================== */