naras@stat.uucp (B. Narasimhan) (02/24/89)
As I have seen quite a clamour for good random number generators I am enclosing
the code for the UNIVERSAL RANDOM NUMBER GENERATOR coded in Microsoft C. It is
not the best implementation, but I hope it is clear enough for the user to
improve upon. For a reference see:
"Toward a Universal Random Number Generator"
by
G. Marsaglia, Arif Zaman and Wai Wan Tsang
to appear in Letters in Statistics and Probability sometime this year.
In my opinion, this is the best RNG existing today and has passed ALL
stringent tests put to it besides being completely portable from ETA10
down to CP/M machines.
B. Narasimhan e-mail: naras@stat.fsu.edu on Internet.
/*******************************************************
* This program tests the Universal Random Number *
* generator. *
* *
* SOFTWARE PROVIDED BY *
* Dr. G. Marsaglia and B. Narasimhan *
* Dept. of Statistics, Florida State University, *
* Tallahassee, FL 32306. *
* Ph. (904)-644-3218. *
*******************************************************/
/*----------------------------------------------------*/
/* Application Notes: */
/* This program has been compiled and tested with */
/* Microsoft C 5.0. Hopefully, it is portable. Note */
/* the code can be made more efficient in easy-to-see */
/* spots. We have avoided these obvious changes in */
/* interest of readability. */
/*----------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define TAB_SIZE 97
#define NO_OF_BITS 24
#define MOD_1 179
#define MOD_2 169
#define CONST_1 53
#define CONST_2 362436
#define CONST_3 7654321
#define CONST_4 16777213
#define CONST_5 96
#define CONST_6 32
#define CONST_7 16777216.0
#define MASK 0xffffff
unsigned long UNI_TABL [TAB_SIZE];
signed long C_VAL;
short I_IND, J_IND;
SET_UNI(a, b, c, d)
int a, b, c, d;
{
int i, j;
ldiv_t div_result;
unsigned long sum;
for (i = 0; i < TAB_SIZE ; i++) {
sum = 0;
for (j = NO_OF_BITS; j > 0; j--) {
div_result = ldiv(a*b,MOD_1);
div_result = ldiv(div_result.rem * c,MOD_1);
a = b;
b = c;
c = div_result.rem;
div_result = ldiv(CONST_1 * d + 1, MOD_2);
d = div_result.rem;
sum *= 2;
div_result = ldiv(d*c,64);
if (div_result.rem >= 32)
sum += 1;
}
UNI_TABL[i] = sum;
}
C_VAL = CONST_2;
I_IND = CONST_5;
J_IND = CONST_6;
}
float UNI(void)
{
unsigned long uni_val;
uni_val = UNI_TABL[I_IND] - UNI_TABL[J_IND];
uni_val &= MASK;
UNI_TABL[I_IND] = uni_val;
if (I_IND == 0)
I_IND = CONST_5;
else
I_IND--;
if (J_IND == 0)
J_IND = CONST_5;
else
J_IND--;
C_VAL -= CONST_3;
if (C_VAL < 0)
C_VAL += CONST_4;
uni_val -= C_VAL;
uni_val &= MASK;
return (float) uni_val / CONST_7;
}
main()
{
int i;
float temp;
unsigned long stemp;
SET_UNI(12,34,56,78);
for (i = 20000; i > 0; i--)
temp = UNI(); /* Generate 20000 unis. */
for (i = 6; i > 0; i--) {
stemp = UNI() * CONST_7;
printf("\t%lu\n",stemp); /* Print next 6 of them as integers. */
}
}dhesi@bsu-cs.UUCP (Rahul Dhesi) (02/24/89)
In article <7354@pyr.gatech.EDU> naras@stat.fsu.edu (B. Narasimhan) presents a universal random number generator that includes the definitions: >#define CONST_2 362436 >#define CONST_3 7654321 >#define CONST_4 16777213 Better make these long, with an L suffix. Else some compilers I've encountered will assume they are 16-bit ints, and truncate them accordingly. -- Rahul Dhesi UUCP: <backbones>!{iuvax,pur-ee}!bsu-cs!dhesi ARPA: bsu-cs!dhesi@iuvax.cs.indiana.edu
gwyn@smoke.BRL.MIL (Doug Gwyn ) (02/25/89)
In article <5873@bsu-cs.UUCP> dhesi@bsu-cs.UUCP (Rahul Dhesi) writes:
-In article <7354@pyr.gatech.EDU> naras@stat.fsu.edu (B. Narasimhan) presents
-a universal random number generator that includes the definitions:
->#define CONST_2 362436
->#define CONST_3 7654321
->#define CONST_4 16777213
-Better make these long, with an L suffix. Else some compilers I've
-encountered will assume they are 16-bit ints, and truncate them
-accordingly.
Since that is in direct contradiction to K&R1's specs for integer
constants, one could make a strong case that such a compiler is not
worth using. Could you tell us whose compilers are thus broken,
so we could steer clear of them?mark@jhereg.Jhereg.MN.ORG (Mark H. Colburn) (02/26/89)
In article <9704@smoke.BRL.MIL> gwyn@brl.arpa (Doug Gwyn (VLD/VMB) <gwyn>) writes: +In article <5873@bsu-cs.UUCP> dhesi@bsu-cs.UUCP (Rahul Dhesi) writes: +-In article <7354@pyr.gatech.EDU> naras@stat.fsu.edu (B. Narasimhan) presents +-a universal random number generator that includes the definitions: +->#define CONST_2 362436 +->#define CONST_3 7654321 +->#define CONST_4 16777213 +-Better make these long, with an L suffix. Else some compilers I've +-encountered will assume they are 16-bit ints, and truncate them +-accordingly. + +Since that is in direct contradiction to K&R1's specs for integer +constants, one could make a strong case that such a compiler is not +worth using. Could you tell us whose compilers are thus broken, +so we could steer clear of them? The XENIX compiler, which is a descendant of a Microsoft compiler (I think), does this. It has been a problem in PAX, and some other code, where I had to explicitly cast the constant 1024 to to a long, so that the compiler wouldn't puke. Yuk! -- Mark H. Colburn "Look into a child's eye; Minnetech Consulting, Inc. there's no hate and there's no lie; mark@jhereg.mn.org there's no black and there's no white."
chip@ateng.ateng.com (Chip Salzenberg) (03/02/89)
According to mark@jhereg.Jhereg.MN.ORG (Mark H. Colburn): >The XENIX compiler, which is a descendant of a Microsoft compiler (I think), >[does not automatically consider large constants to be long]. It has been >a problem in PAX, and some other code, where I had to explicitly cast >the constant 1024 to to a long, so that the compiler wouldn't puke. Although it's true that the Xenix compiler does not promote the constant 1024 to long, that's not the issue. The issue is whether a *constant* too large for an int is automatically considered a long. The Microsoft compiler, whatever its foibles, does handle long constants correctly. It also gives a warning when it does so, which is a nice feature in my estimation. For the curious: The bug in the pax source was the calculation of "1024 * 1024", which should have been coded as "1024L * 1024L". -- Chip Salzenberg <chip@ateng.com> or <uunet!ateng!chip> A T Engineering Me? Speak for my company? Surely you jest! "It's no good. They're tapping the lines."