• Please review our updated Terms and Rules here

Pseudo-random Real Numbers

CP/M User

CP/M User

Veteran Member
Joined
May 2, 2003
Messages
2,986
Location
Back of Burke (Guday!), Australia
I was just wondering if anyone knows how a Pseudo-random Real Number formula is constructed. I have one though it appears to be Integer based which uses the linear congruential method which looks like this:

Rn=(multiplier*Rn-1+increment) mod modulus

Though I thought it would be handy to have one which returns Random values between 0.0 & 1.0

Thanks.
 
If all you want is randomness, just stack enough random integers to add up to the size of your floating-point type. If you want it to be positive, set the sign bit manually. If it's IEEE, just make sure that the combination field (the exponent) isn't 30 or 31 (infinity or NaN.) If you specifically want a value from 0-1, only use random values for the coefficient, and set the exponent such that the highest bit of the coefficient is the ones place.

Granted, this is not so quick-'n-easy as your multiplication approach, but it saves you a multiplication and division, and on the kind of systems we use 'round here that can be pretty significant.
 
...among the most clever that I've seen involved using noise from a diode to generate truly random numbers.

I recall an article written by someone from LANL or LLL some years ago severely critical of most builtin random number generators. I'll see if I can locate it--no promises, however.

I also have the old IBM document describing the algorithms used in their FORTRAN libraries from the 1960s. I can look that one up if you're interested.
 
Thanks Folks, I've gone through all your comments and came up with this routine:

Code: 
   PROGRAM realrnd;
   
   VAR result : real;
       num1   : integer;
       num2   : integer;
   
   BEGIN
     num1:=random(0);
     num2:=random(0);
     IF num1<=num2 THEN
         result:=num1/num2;
     IF num1>num2 THEN
         result:=num2/num1;
     write(result:9:5);
   END.

Unfortunately the Pascal I'm using has it's own way of doing things. Integers for example range from -32768 to 32767 (which is the same as Turbo Pascal Integers), and Reals range from 5.9E-39 (smallest) to 3.4E38 (largest).
What I can do in it is generate a random seed with random(0), so I've done that twice, and then I work out which number is larger and divide that with the first and put that into result.

After doing that though I noticed my write statement was making it look like my values were incorrect, so I have applied a condition, I think I used 9 because I was counting 9 digits before the Exp sign came up and I've specified 5 to return a number with 5 decimal places, which seems to have fixed it.


But I'd be interested to see your FORTRAN algorithms just to compare how they went about it. :D
 
I was going for the S/360 FORTRAN book, when I spied on the shelf below, the old "bible" for such stuff; "Software Manual for the Elementary Functions", by William J. Cody Jr. and William Waite (Prentice-Hall, 1980). It's part of the PH "Computational Mathematics" series and was an excellent source book when I was coding math packs (you'll find my handiwork in a couple of Sorcim products).

At any rate, here's what they give for a PRN generator (in FORTRAN):

Code: 
      REAL FUNCTION RAN(K)
C       
C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 256 BY PIKE AND
C     HILL (MODIFIED BY HANSSON); COMMUNICATIONS OF THE ACM,
C     VOL 8, NO. 10, OCTOBER, 1965.
C
C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTER WITH
C     FIXED-POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS BEST IF
C     THE FLOATING-POINT SIGNIFICAND HAS AT MOST 29 BITS.
C
      INTEGER IY, J, K
      DATA IY/100001/
C
      J = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      RAN = FLOAT(IY) / 2796203.0E0
      RETURN
      END

I do have ACM CALGO 256 if you're interested.
 
Chuck(G) said:
At any rate, here's what they give for a PRN generator (in FORTRAN):

Code: 
      REAL FUNCTION RAN(K)
C       
C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 256 BY PIKE AND
C     HILL (MODIFIED BY HANSSON); COMMUNICATIONS OF THE ACM,
C     VOL 8, NO. 10, OCTOBER, 1965.
C
C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTER WITH
C     FIXED-POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS BEST IF
C     THE FLOATING-POINT SIGNIFICAND HAS AT MOST 29 BITS.
C
      INTEGER IY, J, K
      DATA IY/100001/
C
      J = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      RAN = FLOAT(IY) / 2796203.0E0
      RETURN
      END
Click to expand...

Hmmm, sadly that little bit of code has be stumped, I've gone through my FORTRAN book, though it's FORTRAN 77, which I'm presuming is a bit different from the earlier FORTRANs. I don't think that even Pascal has something like FLOAT which converts an INTEGER into a REAL, even though there are other functions which can take an INTEGER and make the result an REAL. I'm a bit confused what happens with J & K, assuming that when you apply the FUNCTION RAN K holds a Seed? Which then goes into J, though the trail seems to finish there unless I'm missing something?
 
I double-checked Cody and Waite and I got it right. Since functions in FORTRAN require at least one argument, K would be a dummy. I wonder if J prevents some compiler diagnostic for an unused variable. No comments in the book about it. So...

I went to my copy of CALGO Volume II and here's both Pike and Hill's Algorithm 266 and Hansson's modification. There's also a certification of the algorithm by one Walter Sullins of Indiana State U, Terre Haute, but it's just a table of frequencies.

Here's Pike and Hill in all of its glorious ALGOL-60:
Code: 
real procedure random(a, b, y);
  real a,b; integer y;
comment Generates a random number in the interval (a,b);
begin
  y :=3125 × y; y : = y - (y ÷67108864) × 67108864;
  random := y/67108864.0 × (b - a) + a
end random

And here's Hansson's modification:

Code: 
real procedure random(a, b, y);
  real a,b; integer y;
begin
  y :=125 × y; y : = y - (y ÷ 2796203) × 2796203;
  random := y / 2796203 × (b - a) + a
end random

Hansson adds that any starting value for y between 1 and 2796202 will serve as a seed.

ACM CALGO is a fascinating history in mathematics of computing. FWIW, Algorithm 1 is from 1960 and concerns quadrature integration.
 
There are a few random number generation routines in the SWAG under miscellaneous (should be under 'math'). Don't know how good they are, but when looking for ideas on different ways of doing thing's it's sometimes a decent starting point.

http://kd5col.info/swag/MISC/index.html
 
Chuck(G) said:
I double-checked Cody and Waite and I got it right. Since functions in FORTRAN require at least one argument, K would be a dummy. I wonder if J prevents some compiler diagnostic for an unused variable. No comments in the book about it. So...
Hi, By the looks of it, it looks like a Variable FORTRAN wants, though not necessarily needed for that routine and I vaguely recall how fussy those early languages like FORTRAN and COBOL were where when it comes to where they were merely spaced out. :D

I went to my copy of CALGO Volume II and here's both Pike and Hill's Algorithm 266 and Hansson's modification. There's also a certification of the algorithm by one Walter Sullins of Indiana State U, Terre Haute, but it's just a table of frequencies.

Here's Pike and Hill in all of its glorious ALGOL-60:
Code: 
real procedure random(a, b, y);
  real a,b; integer y;
comment Generates a random number in the interval (a,b);
begin
  y :=3125 × y; y : = y - (y ÷67108864) × 67108864;
  random := y/67108864.0 × (b - a) + a
end random

And here's Hansson's modification:

Code: 
real procedure random(a, b, y);
  real a,b; integer y;
begin
  y :=125 × y; y : = y - (y ÷ 2796203) × 2796203;
  random := y / 2796203 × (b - a) + a
end random

Hansson adds that any starting value for y between 1 and 2796202 will serve as a seed.

ACM CALGO is a fascinating history in mathematics of computing. FWIW, Algorithm 1 is from 1960 and concerns quadrature integration.
Click to expand...

I can follow the Algol examples a bit better, though I'm still stumped. :( The problem I'm having must be something to do with the value of the "y" variable. If I'd assume that the "y" value was all the same then "y : = y - (y ÷ 2796203) × 2796203;" appears to cancel itself. :( Though with "y :=125 × y;" does that mean the 2nd "y" is the same as the first "y".

Be good if there is some Algol-60 online, I found one for "Algol-68" though I've been reading about how notoriously complicated that is, so I don't really want to try that one and my Computer hasn't got Algol-60. :(
Click to expand...
Click to expand...
 
deathshadow said:
There are a few random number generation routines in the SWAG under miscellaneous (should be under 'math'). Don't know how good they are, but when looking for ideas on different ways of doing thing's it's sometimes a decent starting point.

http://kd5col.info/swag/MISC/index.html
Click to expand...

Heh, I used to some fun with those SWAG routines, converting the Inline Assembly into Inline M/C, porting it to TP3 and running it in CP/M-86. :D

The Mandelbrot routine on my Website was once one of those SWAG routines, though much of it has been transformed to work on my Amstrad, the Setup program (which takes forever to setup the files with the Mandelbrot Data) still has the algorithms to the Mandelbrot. :)
 
CP/M User said:
I can follow the Algol examples a bit better, though I'm still stumped. :( The problem I'm having must be something to do with the value of the "y" variable. If I'd assume that the "y" value was all the same then "y : = y - (y ÷ 2796203) × 2796203;" appears to cancel itself. :( Though with "y :=125 × y;" does that mean the 2nd "y" is the same as the first "y".
Click to expand...

Note that there's a big difference between the / and ÷ operators. "/" returns a floating point value (i.e. with fractional part) and "÷" performs a division and returns an integer with no fraction. If you'd like, I can code it in BASIC if that helps.
 
Chuck(G) said:
Note that there's a big difference between the / and ÷ operators. "/" returns a floating point value (i.e. with fractional part) and "÷" performs a division and returns an integer with no fraction. If you'd like, I can code it in BASIC if that helps.
Click to expand...

Oh okay, so if I did an "y:= y - int(y div 2796203) * 2796203" I should get a value. I've been trying with decimal numbers on my Calculator, but all I get is 0, but I'll try on my computer this weekend.
 
The laugh being that formula is just trying to pull the modulo / remainder.

in Pascal:
y:= y mod 2796203;

in C:
y = y % 2796203;

Think about it, it's a integer divide throwing away the remainder, times the same value as the divisor.. you subtract that from Y, you get the remainder!
 
The Algol-60 standard did not define a modulus operator.

A little CS history here--as far as I know, the definition of the Algol-60 language was one of the first, if not the first, official specification for a language in BNF (note that the "n" is Peter Naur, the author of the report, although he initially called it Backus "normal" form. The "B" being John Backus, he of FORTRAN among other things).

That really was a big thing, because you now had a regular way to describe a programming language--and, ultimately, describe the parser for it. The divide in programming languages back then was generally US = FORTRAN, Europe = Algol.

Here's the Algol-60 report
 
... Fortran didn't get it until 1990 right? Checking... yeah, was added to Fortran-90.

You're right on the divide across the pond too; I don't even remember hearing about anyone even using Algol for any real programming by the time I started programming in the late '70's... from what I understood it was strictly historical at that point so far as use here in the colonies.
 
deathshadow said:
The laugh being that formula is just trying to pull the modulo / remainder.

in Pascal:
y:= y mod 2796203;

in C:
y = y % 2796203;

Think about it, it's a integer divide throwing away the remainder, times the same value as the divisor.. you subtract that from Y, you get the remainder!
Click to expand...

I tried simulating this in BASIC, though it looks like a number like 2796203 is too big for the "mod" function to support, probably being Interpreted BASIC and only supporting numbers of a limited range. The Pascal I use would only allow that number if it were a Real, I don't have the luxury of Long Integers in the Pascal I'm using (even TP on the 8bits didn't allow this).

So I'm probably just better off doing what I had by obtaining 2 Seed numbers, working out the larger and Dividing the smaller number by the larger? The routine I was going to test this on is a 3D Star Field, and the Random Number field is testing distance and using a dull colour to make it look like a star is further away.
 
You must log in or register to reply here.
Back
Top