Defects in ga68 random number generators
Nelson H. F. Beebe
beebe@math.utah.edu
Sat May 23 16:05:46 GMT 2026
The argument-less procedures for accessing floating-point random
numbers in Algol 68 are defined extremely sparsely in the Revised
Report in Algol68_revised_report-AB-600dpi.pdf, with a little more
detail in the Very Informal Introduction to Algol 68 in
Lindsey_van_der_Meulen-IItA68-Revised.pdf, which says
>> ...
>> there is a proc real random, which returns the next pseudo-random real
>> value from a uniformly distributed sequence on the interval [0,1)
>> (i.e. 0 <= random < 1).
>> ...
However, there are no further details on the quality, or period, of
the generator, which essentially leaves it up to the compiler and
library implementation.
We know far more about random number generation today than was known
in the 1950s and 1960s when programming languages, including the Algol
family, were first defined. I devoted the 57 pages of Chapter 7 of
The Mathematical Function Handbook to random number generation,
including a summary of their desirable properties, and generation of
numbers corresponding to various distributions, of which uniform is by
far the most common. That chapter also discusses reproducible
sequences needed for simulation, and true random sequences needed for
cryptography and some other applications.
I also actively maintain a large bibliography of publications in this
area at
https://www.math.utah.edu/pub/tex/bib/prng.bib
I was therefore interested in learning how ga68 supplies random
numbers: it recognizes shortrandom, random, longrandom,
longlongrandom, longlonglongrandom, and possibly versions with more
"long" prefixes. I wrote test code, and dug into compiler sources to
find in the file gcc-17-20260517/libga68/ga68-standenv.c this block:
float
_libga68_random (void)
{
float res = (float) rand () / (float) (RAND_MAX);
return res;
}
double
_libga68_longrandom (void)
{
double res = (double) rand () / (float) (RAND_MAX);
return res;
}
long double
_libga68_longlongrandom (void)
{
long double res = (long double) rand () / (float) (RAND_MAX);
return res;
}
The ISO C99 Standard says:
>> ...
>> 7.20.2 Pseudo-random sequence generation functions
>>
>> 7.20.2.1 The rand function
>>
>> Synopsis
>>
>> #include <stdlib.h>
>> int rand(void);
>>
>> Description
>>
>> The rand function computes a sequence of pseudo-random integers in
>> the range 0 to RAND_MAX.
>>
>> The implementation shall behave as if no library function calls the
>> rand function.
>>
>> Returns
>>
>> The rand function returns a pseudo-random integer.
>>
>> Environmental limits
>>
>> The value of the RAND_MAX macro shall be at least 32767.
>> ...
Notice that nothing is said about the algorithm(s) used for random
number generation, and testing quickly shows that the sequences
generated differ across platforms. Neither does it say anything about
the initial seed. That makes it impossible to use rand() for
reproducible simulations, and consequently, rand() should generally be
avoided in any nontrivial use of random numbers.
Most systems will have RAND_MAX = 2**n - 1, where n is the number of
bits in a C int, ignoring the sign bit: historically, this has been
2**15 - 1 = 0x7fff, 2**31 - x = 0x7fff_ffff, and on 36-bit systems,
2**35 - 1 = +0x7_ffff_ffff. I have in the past encountered
Unix-family systems that have 32-bit int, but where RAND_MAX = 0x7fff,
with a dreadfully small period of only 32_768 numbers.
Floating-point on modern systems is almost universally based on IEEE
754 arithmetic, though sometimes with subsetted features. Thus,
(float)RAND_MAX, with RAND_MAX = 0x7fff_ffff converts that value with
rounding to a 24-bit number = 0x8000_0000 = 2**31, which exceeds
RAND_MAX. That produces an exact scaling of the numerator by 0x1p-31.
With 64-bit, 80-bit, and 128-bit IEEE 754 results, the largest value
of the quotient is strictly less than 1.0. However, with a 32-bit
float, the result can be exactly 1.0, because with M = RAND_MAX =
0x7fff_ffff, any value of the form (float)(M - k)/M, for k = 0 to 64,
rounds to exactly 1.0, violating the Algol 68 Revised Report
specification.
A rewrite like this restores conformance:
float
_libga68_random (void)
{
float res;
do
res = (float) rand () / (float) (RAND_MAX);
while (res == 1.0F);
return res;
}
Even on the same platform, however, the ga68 random family will not
necessarily use the same sequence on repeated runs. It does on most,
but OpenBSD chooses a unique initial seed on each run.
To guarantee identical cross-platform behavior, we need a common
generator across all systems, and one way to do so without supplying
private random-number code is to use the POSIX rand48 family member
lrand48(), which returns a value in the interval [0, 2**31 - 1].
However, in addition, the runtime startup code before the user program
is entered MUST set a common initial seed, such as with a call to
srand48((long int)DEFAULT_SEED).
On OpenBSD, one has to go further, because they wrongly chose to
change the behavior of the POSIX rand48 family. The code that I use
to defend against that is this wrapper:
#if defined(__OpenBSD__) && !defined(__Bitrig__) && !defined(__MirBSD__)
/* OpenBSD replaces POSIX-standard *rand48() deterministic functions
with nondeterministic ones of the same name. We need to guarantee
the same random sequence on every system, so revert to the OpenBSD
function that obeys POSIX requirements. */
#undef srand48
#define srand48(x) srand48_deterministic(x)
#undef seed48
#define seed48(x) seed48_deterministic(x)
#endif /* defined(__OpenBSD__) && !defined(__Bitrig__) && !defined(__MirBSD__) */
The safe solution for the ga68 source code is probably to include a
private copy of the open-source rand48 code, which is natively
available on modern systems, but as noted, its good intents are
thwarted by at least OpenBSD.
Suitable macros in the ga68 source tree should therefore redefine all
of the public names in the rand48 code with private ones, such as
a68_lrand(), to avoid conflict with system library versions of the
rand48 family.
ga68 also needs to supply, and document, a standard seed-setting
function for its random family; at present, it does not.
There remains however another serious problem in _libga68_longrandom
and _libga68_longlongrandom: because they use only a single call to
rand(), they get at best 31 random bits, instead of the 53, 64, or 113
that their types can hold.
They should therefore be reimplemented something like this in IEEE 754
arithmetic:
double
_libga68_longrandom (void)
{
// double res = (double) rand () / (float) (RAND_MAX);
double res;
int r1, r2;
assert(RAND_MAX >= 0x7fffffff);
(r1 = rand(), r2 = rand());
res = ((double)r1 +
0x1p-31 * (double)r2) /
(double)(float) (RAND_MAX);
return res;
}
long double
_libga68_longlongrandom (void)
{
// long double res = (long double) rand () / (float) (RAND_MAX);
long double res;
assert(RAND_MAX >= 0x7fffffff);
int r1, r2, r3, r4;
(r1 = rand(), r2 = rand(), r3 = rand(), r4 = rand());
res = ((long double)r1 +
0x1p-31L * (long double)r2 +
0x1p-62L * (long double)r3 +
0x1p-93L * (long double)r4) /
(long double)(float) (RAND_MAX);
return res;
}
Notice the critical use of the C comma operator to enforce strict
order evaluation of (a, b), and (a, b, c, d), to prevent optimizers
from rearranging the order of calls to rand().
Also notice that the longlongrandom function might correspond to
either 80-bit or 128-bit arithmetic, and on deficient platforms, like
Apple macOS , some of the BSD family, and Linux on HPPA (PA-RISC),
just to 64-bit arithmetic. We therefore must use FOUR random numbers
to guarantee sufficient random bits, and to generate the same
sequences (modulo rounding errors in the evaluation of res) across
platforms.
Of course, RAND_MAX and rand() in those examples need to be changed to
the recommended rand48 family values.
----------------------------------------
P.S. We can discuss later for ga68g whether the linear congruential
generator family rand48 should be avoided entirely, in favor of
portable code with vastly larger periods than the 2**47 - 1 of
lrand48(); there are many such in the random-number literature.
See
https://www.math.utah.edu/~beebe/random
for some recent excellent generators. Random number streams can be
thoroughly tested with the excellent open-source dieharder suite that
is in the package systems of many O/Ses. The jones-1.0.0/TESTING file
at that site supplies further commentary on dieharder.
Alternatively, the Bays--Durham shuffle buffer technique (entry
Bays:1976:IPR in prng.bib) can be used to substantially increase the
period of ANY generator, and improve randomness. In GNU gawk, we use
that in the file support/random.c with a 512-entry buffer to increase
the period of the original fairly decent generator by a factor of
about 10**569, at essentially zero run-time cost, and 4096 bytes of
additional buffer memory.
-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: +1 801 581 5254 -
- University of Utah -
- Department of Mathematics, 110 LCB Internet e-mail: beebe@math.utah.edu -
- 155 S 1400 E RM 233 beebe@acm.org beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA URL: https://www.math.utah.edu/~beebe -
-------------------------------------------------------------------------------
More information about the Algol68
mailing list