Path: news.io.com!uunet!in1.uu.net!world!blanket.mitre.org!linus.mitre.org!
+ spectre!eachus
From: eachus@spectre.mitre.org (Robert I. Eachus)
Newsgroups: comp.arch.arithmetic
Subject: Re: Psuedo Random Numbers
Date: 11 Oct 1995 22:37:56 GMT
Organization: The Mitre Corp., Bedford, MA.
Lines: 77
Message-ID:
References: <44vm6r$1q5@bug.rahul.net>
+ <4593t4$2nm0@b.stat.purdue.edu>
NNTP-Posting-Host: spectre.mitre.org
In-reply-to: hrubin@b.stat.purdue.edu's message of 8 Oct 1995 13:05:24 -0500
In article <4593t4$2nm0@b.stat.purdue.edu> hrubin@b.stat.purdue.edu (Herman Rubi
n) writes:
> If this is the paper I think it is, the generator is far too slow
> to be of much use, unless one needs fantastically good
> pseudo-random numbers, and can afford the cost. Now this might
> not be the paper I am thinking of, but the one I have in mind uses
> two large primes congruent to 3 mod 4, and at each stage squares
> the seed modulo the product, and outputs the least significant
> bit. This is a lot of work for ONE bit.
Same algorithm different implementation. ;-) First of all, the
current situation is that log2 N bits not one bit are provably secure
per iteration. But the generator I use is not intended to be
cryptographically secure, just sufficiently random to pass all tests
which do not involve factoring the modulus.
The second bag of tricks and the one that makes it all practical
is to use the Chinese Remainder Theorem to allow you to do all of the
computations with machine arithmetic operations instead of
multiprecision. The biggest hit is two mods per iteration that do
require divides but a 62-bit number mod a 32-bit integer is pretty
fast on most modern hardware. (I have a version which uses
floating-point to do the mod, but writing those few instructions in
assembler is usually the best approach. Ada 95 allows you to write it
in high level code without redundant checks, but most programming
languages don't.*) Oh, the rest of the mods are implementable as test
and subtracts, since they are of the form (x+y) mod m.
If all this seems tremendously complex, the design part is. But
the generator code is very simple, and not much slower than an LCG.
Get a copy and try it out. The guarentee you get can twist your mind
into a pretzel, but don't let it worry you. Basically it says that
you can't predict the previous output value from an arbitrarily long
sequence of values. I guess you could write out a sequence and reverse
it, but for most randomness tests, and most concerns about
pseudorandom sequences, the order is unimportant. The correlation
between any two numbers adjacent or any distance apart up to 10% of
the modulus is zero.)
*Ada 95 allows you to declare modular types where the modulus is
not a power of two. So to compute X squared mod M, you can write:
type Mod_M is M; -- for M a static value less than
-- System.Max_Nonbinary_Modulus
-- now I can write:
X: Mod_M;
...
X := X*X;
-- and hopefully I get three or four instructions on most hardware:
move X to Y -- not required on all hardware but is on most
multiply 32x32--> 64
divide 64/32 --> 32 -- by a constant
discard the quotient, the remainder is the answer.
(Flame control: Some of the Power PC chips don't have the divide
in hardware, you have to emulate it in software. The SPARC
architecture now includes the divide, but not the remainder, which is
discarded. In both cases doing the mod exactly in floating point,
while tricky, is possible and reasonably fast. Of course, I have the
luxury of designing the algorithm by picking numbers such that no
precision is lost in the floating point operations. Well, the final
number generated and returned is limited by the floating point
accuracy, but it is not used explicitly to generate the next variate.
The Chinese Remainder Theorem representation kept internally is exact,
and is used to generate the next variate.
--
Robert I. Eachus
with Standard_Disclaimer;
use Standard_Disclaimer;
function Message (Text: in Clever_Ideas) return Better_Ideas is...