均勻分布和高斯分布

/* random.c */

/* Created on Apr 5, 1999 Shanghai Jiao Tong University
    by PENG Zhaowang
*****************************************************/
/* 
This file contains routines for generating random numbers. The routines are from
  Chapter 7 of the book by Press, Teukolsky, et al. 

 William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery.
Numerical Recipes in C, The Art of Scientific Computing, Second Edition. 
CAMBRIDGE UNIVERSITY PRESS, 1992, Cambridge. 

  中譯本: 傅祖蕓, 趙梅娜, 丁巖 等譯. C語言數值算法程序大全(第二版). 
  北京: 電子工業出版社, 1995. 

Small modification made by PENG Zhaowang on 5 Apr 1999.
*/

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

#include <math.h>
#include <time.h>
#include <stdlib.h> /* Needed for Borland compilers. */
#include "randnum.h"

static long random_number_seed;
static int flag_seeded=0;

double rand1(void)  /* Generating Uniform random number between 0.0 and 1.0*/
{
        if(!flag_seeded)
        {
                random_number_seed = (long)time(NULL);
                random_number_seed = 1664525L*random_number_seed + 1013904223L;
                if(random_number_seed>0)
                        random_number_seed = -random_number_seed;
                flag_seeded=1;
        }
        return ran1(&random_number_seed);
}

double gassdev(void)
{
        if(!flag_seeded)
        {
                random_number_seed = (long)time(NULL);
                random_number_seed = 1664525L*random_number_seed + 1013904223L;
                if(random_number_seed>0)
                        random_number_seed = -random_number_seed;
                flag_seeded=1;
        }
        return gasdev(&random_number_seed);
}

/* "Minimal" random number generator of Park and Miller with Bays-Durham shuffle and added
safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint
values). Call with idum a negative integer to initialize; thereafter, do not alter idum between
successive deviates in a sequence. RNMX should approximate the largest floating value that is
less than 1. */
double ran1(long *idum)
{
        int j;
        long k;
        static long iy=0;
        static long iv[NTAB];
        double temp;
        if (*idum <= 0 || !iy)
        {                                                                       /*Initialize.*/
                if (-(*idum) < 1)
                        *idum=1;                                /*Be sure to prevent idum = 0.*/
                else
                        *idum = -(*idum);
                for (j=NTAB+7;j>=0;j--)
                {                                                                /*Load the shuffle table (after 8 warm-ups).*/
                        k=(*idum)/IQ;
                        *idum=IA*(*idum-k*IQ)-IR*k;
                        if (*idum < 0)
                                *idum += IM;
                        if (j < NTAB)
                        iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ; /* Start here when not initializing. */
        *idum=IA*(*idum-k*IQ)-IR*k; /* Compute idum=(IA*idum) % IM without overflows */
        if (*idum < 0)
                *idum += IM;            /* by Schrage's method. */
        j=iy/NDIV; /* Will be in the range 0..NTAB-1. */
        iy=iv[j];  /* Output previously stored value and rell the shuffle table. */
        iv[j] = *idum;
        if ((temp=AM*iy) > RNMX)
                return RNMX; /* Because users don't expect endpoint values. */
        else
                return temp;
}



/* Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum)
as the source of uniform deviates. */
double gasdev(long *idum)
{
        static int iset=0;
        static double gset;
        double fac,rsq,v1,v2;
        if (*idum < 0)
                iset=0;                         /* Reinitialize */
        if (iset == 0)
        {
                do
                {
                        v1=2.0*ran1(idum)-1.0; /* pick two uniform numbers in the square extending */
                        v2=2.0*ran1(idum)-1.0;          /*  from -1 to +1 in each direction, */
                        rsq=v1*v1+v2*v2;  /* see if they are in the unit circle */
                }while (rsq >= 1.0 || rsq == 0.0); /* and if they are not, try again. */
                fac=sqrt(-2.0*log(rsq)/rsq); 
        /* Now make the Box-Muller transformation to get two normal deviates. Return one and
                save the other for next time. */
                gset=v1*fac;
                iset=1; /* Set flag. */
                return (v2*fac); 
                /* NOTE: Without the brackets, this function cannot work well with its release edition 
                        generated by Microsoft VC5.0, so it is strongly suggested to keep them here.
                Why? I don't know! Maybe Mr. Bill Gates knows. 
                by PENG Zhaowang on Apr 6 1999 */
        }
        else
        {                                               /* We have an extra deviate handy, */
                iset=0;                         /* so unset the flag, */
                return gset;            /* and return it. */
        }
}

--
* random.h */

/* Created on Apr 5, 1999 Shanghai Jiao Tong University
    by Peng, Zhaowang
*****************************************************/
/*      The Head File for Generating random number
*****************************************************/

#ifndef __RANDNUM_H99_4_5_22080409
#define __RANDNUM_H99_4_5_22080409

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */


/* 
This file contains declaration of the routines for generating random numbers. 
The routines are from Chapter 7 of the book by Press, Teukolsky, et al. 

 William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery.
Numerical Recipes in C, The Art of Scientific Computing, Second Edition. 
CAMBRIDGE UNIVERSITY PRESS, 1992, Cambridge. 

  中譯本: 傅祖蕓, 趙梅娜, 丁巖 等譯. C語言數值算法程序大全(第二版). 
  北京: 電子工業出版社, 1995. 

*/

/* 為簡單起見, 建議使用rand1()
Automatic seeded with the ANSI C time() function. You are suggested to use rand1(void) 
 or gassdev(void) for your application since they are automatically seeded with the ANSI C
 time() function. However, you may call ran1 or gasdev after seeding idum with something.
 */

double rand1(void); /* Uniform deviation within [0, 1] */
double gassdev(void); /* Gaussian deviation with zero mean 
                        and unit standard deviation  */

/* Not seed random number generator. */
double ran1(long *idum); /* Uniform deviation within [0, 1] */
double gasdev(long *idum); /* Gaussian deviation with zero mean 
                        and unit standard deviation */

#ifdef __cplusplus
}
#endif /* __cplusplus */

#endif /* __RANDNUM_H99_4_5_22080409 */