Генератор случайных чисел по алгоритму Парка-Миллера.
По заверениям дает белый шум.
Код:
#include <stdio.h>
/* Minimal portable random generator by Park and Miller */
/* Lewis-Goodman-Miller constants */
#define IA 16807
#define IM 2147483647
#define AM (1./IM)
/* Scharge constants */
#define IQ 12773
#define IR 2836
#define NTAB 32
#define NWUP 8
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
/* Special mask to be explained below */
#define MASK 123456789
static long dummy;
/* initial seed, for all the generators here */
void Seed(long dum) {
dummy = dum;
}
/* returns random uniformly distributed between 0 and 1 */
float unirand0(void) {
long k;
float ans;
dummy ^= MASK; /* avoid dummy==0 */
k = dummy / IQ;
if ((dummy=IA*(dummy-k*IQ)-IR*k)<0)
dummy += IM;
ans = AM * dummy;
dummy ^= MASK; /* restore unmasked dummy */
return(ans);
}
float unirand1(void) {
int j;
long k;
static long iy=0,iv[NTAB];
float temp;
/* initialize */
if (dummy <= 0 || !iy) {
/* avoid negative or zero seed */
if (dummy < 0) dummy =- dummy; else
if (dummy == 0) dummy = 1;
/* after NWUP warmups, initialize shuffle table */
for(j = NTAB + NWUP - 1; j >= 0; j--) {
k = dummy / IQ;
if ((dummy = IA * (dummy - k * IQ) - IR * k) < 0)
dummy += IM;
if (j < NTAB)
iv[j] = dummy;
}
/* first specimen from the table */
iy=iv[0];
}
/* regular work: generate new number */
k = dummy / IQ;
if ((dummy = IA * (dummy - k * IQ) - IR * k) < 0)
dummy += IM;
/* shuffle output */
iy = iv[j = iy / NDIV];
iv[j] = dummy;
/* return */
if ((temp = AM * iy) > RNMX)
return(RNMX);
else
return(temp);
}
void main(){
int i;
Seed(6723);
for (i = 0; i < 100; i++)
printf("%f\n", unirand1());
}
>>bcc32 parkmiller.c
З.Ы. Код не совсем нормализован, некоторые места можно сократить.