libc: Reimplement the *rand48 family of functions

Rather than implementing the recurrence using 3 16-bit integers, as was
done in _dorand48() before this patch, provide an equivalent
implementation using 64-bit integers.

For drand48() and erand48(), replace the use of ldexp() with
bit-twiddling assuming IEEE 754 double-precision float layout.

This implementation is significantly faster and requires less code,
while producing identical outputs on supported platforms.

While here, add a STANDARDS section to rand48.3.

Obtained from:	https://github.com/apple-oss-distributions/libc
MFC after:	3 weeks
Sponsored by:	Klara, Inc.
Differential Revision:	https://reviews.freebsd.org/D52429

(cherry picked from commit 2ba20004ef7649db7654520e8376927c4410d9c3)
This commit is contained in:
Mark Johnston 2025-09-04 19:34:17 +00:00
parent 79fbc4ccb6
commit 36224ed427
12 changed files with 95 additions and 90 deletions

View file

@ -13,34 +13,6 @@
#include "rand48.h"
unsigned short _rand48_seed[3] = {
RAND48_SEED_0,
RAND48_SEED_1,
RAND48_SEED_2
};
unsigned short _rand48_mult[3] = {
RAND48_MULT_0,
RAND48_MULT_1,
RAND48_MULT_2
};
unsigned short _rand48_add = RAND48_ADD;
void
_dorand48(unsigned short xseed[3])
{
unsigned long accu;
unsigned short temp[2];
accu = (unsigned long) _rand48_mult[0] * (unsigned long) xseed[0] +
(unsigned long) _rand48_add;
temp[0] = (unsigned short) accu; /* lower 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += (unsigned long) _rand48_mult[0] * (unsigned long) xseed[1] +
(unsigned long) _rand48_mult[1] * (unsigned long) xseed[0];
temp[1] = (unsigned short) accu; /* middle 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0];
xseed[0] = temp[0];
xseed[1] = temp[1];
xseed[2] = (unsigned short) accu;
}
uint48 _rand48_seed = RAND48_SEED;
uint48 _rand48_mult = RAND48_MULT;
uint48 _rand48_add = RAND48_ADD;

View file

@ -13,10 +13,10 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
double
drand48(void)
{
return erand48(_rand48_seed);
ERAND48_BEGIN;
_DORAND48(_rand48_seed);
ERAND48_END(_rand48_seed);
}

View file

@ -16,8 +16,9 @@
double
erand48(unsigned short xseed[3])
{
_dorand48(xseed);
return ldexp((double) xseed[0], -48) +
ldexp((double) xseed[1], -32) +
ldexp((double) xseed[2], -16);
uint48 tmp;
ERAND48_BEGIN;
DORAND48(tmp, xseed);
ERAND48_END(tmp);
}

View file

@ -11,14 +11,13 @@
* to anyone/anything when using this software.
*/
#include <stdint.h>
#include "rand48.h"
long
jrand48(unsigned short xseed[3])
{
uint48 tmp;
_dorand48(xseed);
return ((int32_t)(((uint32_t)xseed[2] << 16) | (uint32_t)xseed[1]));
DORAND48(tmp, xseed);
return ((int)((tmp >> 16) & 0xffffffff));
}

View file

@ -13,18 +13,10 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
extern unsigned short _rand48_mult[3];
extern unsigned short _rand48_add;
void
lcong48(unsigned short p[7])
{
_rand48_seed[0] = p[0];
_rand48_seed[1] = p[1];
_rand48_seed[2] = p[2];
_rand48_mult[0] = p[3];
_rand48_mult[1] = p[4];
_rand48_mult[2] = p[5];
LOADRAND48(_rand48_seed, &p[0]);
LOADRAND48(_rand48_mult, &p[3]);
_rand48_add = p[6];
}

View file

@ -13,11 +13,9 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
long
lrand48(void)
{
_dorand48(_rand48_seed);
return ((long) _rand48_seed[2] << 15) + ((long) _rand48_seed[1] >> 1);
_DORAND48(_rand48_seed);
return (_rand48_seed >> 17) & 0x7fffffff;
}

View file

@ -15,13 +15,9 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
long
mrand48(void)
{
_dorand48(_rand48_seed);
return ((int32_t)(((uint32_t)_rand48_seed[2] << 16) |
(uint32_t)_rand48_seed[1]));
_DORAND48(_rand48_seed);
return ((int)((_rand48_seed >> 16) & 0xffffffff));
}

View file

@ -16,6 +16,8 @@
long
nrand48(unsigned short xseed[3])
{
_dorand48(xseed);
return ((long) xseed[2] << 15) + ((long) xseed[1] >> 1);
uint48 tmp;
DORAND48(tmp, xseed);
return ((tmp >> 17) & 0x7fffffff);
}

View file

@ -9,7 +9,7 @@
.\" of any kind. I shall in no event be liable for anything that happens
.\" to anyone/anything when using this software.
.\"
.Dd September 4, 2012
.Dd September 11, 2025
.Dt RAND48 3
.Os
.Sh NAME
@ -183,5 +183,8 @@ generator calls.
.Xr arc4random 3 ,
.Xr rand 3 ,
.Xr random 3
.Sh STANDARDS
The functions described in this page are expected to conform to
.St -p1003.1-2008 .
.Sh AUTHORS
.An Martin Birgmeier

View file

@ -14,10 +14,11 @@
#ifndef _RAND48_H_
#define _RAND48_H_
#include <sys/types.h>
#include <math.h>
#include <stdlib.h>
void _dorand48(unsigned short[3]);
#include "fpmath.h"
#define RAND48_SEED_0 (0x330e)
#define RAND48_SEED_1 (0xabcd)
@ -27,4 +28,62 @@ void _dorand48(unsigned short[3]);
#define RAND48_MULT_2 (0x0005)
#define RAND48_ADD (0x000b)
typedef uint64_t uint48;
extern uint48 _rand48_seed;
extern uint48 _rand48_mult;
extern uint48 _rand48_add;
#define TOUINT48(x, y, z) \
((uint48)(x) + (((uint48)(y)) << 16) + (((uint48)(z)) << 32))
#define RAND48_SEED TOUINT48(RAND48_SEED_0, RAND48_SEED_1, RAND48_SEED_2)
#define RAND48_MULT TOUINT48(RAND48_MULT_0, RAND48_MULT_1, RAND48_MULT_2)
#define LOADRAND48(l, x) do { \
(l) = TOUINT48((x)[0], (x)[1], (x)[2]); \
} while (0)
#define STORERAND48(l, x) do { \
(x)[0] = (unsigned short)(l); \
(x)[1] = (unsigned short)((l) >> 16); \
(x)[2] = (unsigned short)((l) >> 32); \
} while (0)
#define _DORAND48(l) do { \
(l) = (l) * _rand48_mult + _rand48_add; \
} while (0)
#define DORAND48(l, x) do { \
LOADRAND48(l, x); \
_DORAND48(l); \
STORERAND48(l, x); \
} while (0)
#define ERAND48_BEGIN \
union { \
union IEEEd2bits ieee; \
uint64_t u64; \
} u; \
int s
/*
* Optimization for speed: assume doubles are IEEE 754 and use bit fiddling
* rather than converting to double. Specifically, clamp the result to 48 bits
* and convert to a double in [0.0, 1.0) via division by 2^48. Normalize by
* shifting the most significant bit into the implicit one position and
* adjusting the exponent accordingly. The store to the exponent field
* overwrites the implicit one.
*/
#define ERAND48_END(x) do { \
u.u64 = ((x) & 0xffffffffffffULL); \
if (u.u64 == 0) \
return (0.0); \
u.u64 <<= 5; \
for (s = 0; !(u.u64 & (1LL << 52)); s++, u.u64 <<= 1) \
; \
u.ieee.bits.exp = 1022 - s; \
return (u.ieee.d); \
} while (0)
#endif /* _RAND48_H_ */

View file

@ -13,24 +13,14 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
extern unsigned short _rand48_mult[3];
extern unsigned short _rand48_add;
unsigned short *
seed48(unsigned short xseed[3])
{
static unsigned short sseed[3];
sseed[0] = _rand48_seed[0];
sseed[1] = _rand48_seed[1];
sseed[2] = _rand48_seed[2];
_rand48_seed[0] = xseed[0];
_rand48_seed[1] = xseed[1];
_rand48_seed[2] = xseed[2];
_rand48_mult[0] = RAND48_MULT_0;
_rand48_mult[1] = RAND48_MULT_1;
_rand48_mult[2] = RAND48_MULT_2;
STORERAND48(_rand48_seed, sseed);
LOADRAND48(_rand48_seed, xseed);
_rand48_mult = RAND48_MULT;
_rand48_add = RAND48_ADD;
return sseed;
return (sseed);
}

View file

@ -13,18 +13,11 @@
#include "rand48.h"
extern unsigned short _rand48_seed[3];
extern unsigned short _rand48_mult[3];
extern unsigned short _rand48_add;
void
srand48(long seed)
{
_rand48_seed[0] = RAND48_SEED_0;
_rand48_seed[1] = (unsigned short) seed;
_rand48_seed[2] = (unsigned short) (seed >> 16);
_rand48_mult[0] = RAND48_MULT_0;
_rand48_mult[1] = RAND48_MULT_1;
_rand48_mult[2] = RAND48_MULT_2;
_rand48_seed = TOUINT48(RAND48_SEED_0, (unsigned short)seed,
(unsigned short)(seed >> 16));
_rand48_mult = RAND48_MULT;
_rand48_add = RAND48_ADD;
}