Topic — Hardware random numbers

Needs: Lesson 2 — Values, control flow and functions, Lesson 3 — How a PE runs a loop: vector instructions and cycles.

Each PE has four pseudo-random number generators in hardware. Tungsten's rand() reads the active one. This topic seeds the generators, chooses which one rand() reads and how it steps, changes that choice while the program runs, and maps the raw bits to uniform values in (0, 1). Every value the PE produces is checked against a host model of the generator, prng_ref.py.

What you will be able to do

• Draw 16-bit and 32-bit pseudo-random values with rand(), and say how many generator steps a loop of draws uses.
• Seed the four generators, pick the one rand() reads, and change its recurrence at run time.
• Space a store to the generator's registers from the next draw, so that the draw sees it.
• Turn 16-bit draws into uniform values in (0, 1) that are exact in sp, and predict a kernel's sequence on the host.

The program

One PE, no colors. The kernel hwrand.w imports the module hwprng.w and runs these steps in order:

seed[8] (from the image) ──> prng_seed ──> rnd16[N]  N draws into xp
                                       ──> rnd32[N]  N draws into u32, continuing the sequence
                             prng_select_poly(prng, poly2)
                                       ──> rnd16b[N] N draws into xp, under the new recurrence
rnd16[N] ──> u[N] = rnd16 · 2⁻¹⁶ + (½ + 2⁻¹⁷)

The three draw loops read one generator, so they consume one continuous sequence. check.py rebuilds that sequence from the seed with prng_ref.py and requires every value to match exactly.

The code

The generator's registers: hwprng.w

xp param prng;   // the generator rand() reads, 0..3
xp param poly;   // the recurrence it starts with, 0..3

xp prng_cfg[n=1] address(0x7e27) = {
    (prng << 8) | (0xe4 & ~(3 << (2 * prng))) | (poly << (2 * prng))
};
xp prng_state[n=8] address(0x7dd0);

address() places an array at a fixed word address, as the timestamp counter of lesson 3 is. These two arrays are the generator's registers.

prng_cfg is the configuration register. Its bits are:

Bits	Meaning
9:8	The generator rand() reads in the program's main code
7:6	The recurrence of generator 3 (default 3)
5:4	The recurrence of generator 2 (default 2)
3:2	The recurrence of generator 1 (default 1)
1:0	The recurrence of generator 0 (default 0)

The initializer sets bits 9:8 to prng, and gives generator prng the recurrence poly. The other generators keep their defaults, which are the 0xe4 in the expression. prng and poly are parameters, so the whole expression is a constant, and the value is in the image when the program starts. The parentheses matter: &, | and ^ share one precedence level in Tungsten (see the operator precedence).

prng_state holds the four generators' state, two words each.

Seeding: prng_seed

function prng_seed(xp s[8]) {
    ∀i in [0, 8) prng_state[i] = s[i];
    // two instructions before the next draw
    __builtin_nop();
    __builtin_nop();
}

Each generator has 23 bits of state. s[2k] holds the low 16 bits of generator k, and the low 7 bits of s[2k+1] hold the other 7. The function writes all four generators, whichever one the kernel reads.

The two __builtin_nop() calls are each one instruction that does nothing. They are there because of a hardware rule, explained under What the hardware does.

Changing the recurrence: prng_select_poly

function prng_select_poly(xp id, xp sel) {
    xp shift = 2 * (id & 3);
    xp mask = 3 << shift;
    prng_cfg[0] = (prng_cfg[0] & ~mask) | ((sel & 3) << shift);
    // three instructions before the next draw
    __builtin_nop();
    __builtin_nop();
    __builtin_nop();
}

It replaces the two bits of generator id in the configuration register with sel. The generator's state stays as it is, so the sequence continues from where it was, under the new rule.

Drawing: hwrand.w

main {
    let i in [0, N);

    prng_seed(seed);

    ∀i rnd16[i] = rand();
    ∀i rnd32[i] = rand();

    prng_select_poly(prng, poly2);

    ∀i rnd16b[i] = rand();

    ∀i u[i] = rnd16[i];
    ∀i u[i] = u[i] * c + f;
}

rand() takes its width from its destination: 16 bits into an xp, 32 bits into a u32. Each draw loop is one instruction. The listing (test-tutorial/hwrand_artefacts/hwrand1.yeti) shows the first two:

	imov16 [T27] = lfsr;  //@../hwrand.w:24
	imov32 [T28] = lfsr;  //@../hwrand.w:25

lfsr is the active generator, used as a source operand. The values are raw bits. Read as an xp, they cover -32768 to 32767.

seed comes from the image. Its initializer is a fallback: the Makefile writes the eight words that gen.py derives from SEED over it with miset, so a new seed needs no recompile.

Uniform values in (0, 1)

sp c = 1.0 / 65536.0;
sp f = 0.5 + 0.5 / 65536.0;

In main, ∀i u[i] = rnd16[i]; converts each xp draw to sp, and ∀i u[i] = u[i] * c + f; maps [-32768, 32767] to [2⁻¹⁷, 1 − 2⁻¹⁷]. Every step is exact in sp: c is a power of two, so the product is exact, and the sum needs at most 17 significant bits, well under the 24 that sp has. The result therefore never rounds to 0.0 or 1.0, which matters before a logarithm or a division by 1 - u. Because it is exact, check.py compares the bit patterns.

What the hardware does

Four generators, one active. One generator is active at a time, and rand() reads it. The program's main code starts with the generator that bits 9:8 of the configuration register name. A task that a wavelet starts (lesson 7) begins instead with the generator that its input queue's configuration names.

One draw. A draw takes the low 20 bits of the 23-bit state and stretches them to 64 bits with a second, fixed 20-bit shift register, x²⁰ + x¹⁹ + x¹⁶ + x¹⁴ + 1. The 20 state bits become bits 63:44, and the 44 generated bits fill bits 43:0. The state then advances 128 steps of the selected recurrence. A step shifts the state left by one bit and feeds the XOR of two tap bits into bit 0.

Lanes per draw. One draw feeds four 16-bit destinations or two 32-bit ones, from the low bits up. A 16-bit draw loop therefore fills at most four elements per cycle, and a 32-bit one two. A loop of length L into xp uses ⌈L/4⌉ draws, and into u32 ⌈L/2⌉; a partial last group discards the rest of its draw. A single assignment such as a[0] = rand(); uses a whole draw for one value. All 64 bits of a draw come from 20 state bits, so the four 16-bit values of one draw are not independent of each other: together they take at most 2²⁰ of the 2⁶⁴ possible values.

The four recurrences. Each is a primitive trinomial, so it visits every non-zero state of its width before it repeats:

poly	Recurrence	Width	Period
0	x²³ + x¹⁸ + 1	23	8,388,607
1	x²² + x²¹ + 1	22	4,194,303
2	x²¹ + x¹⁹ + 1	21	2,097,151
3	x²⁰ + x¹⁷ + 1	20	1,048,575

A draw advances 128 steps, and 128 shares no factor with these odd periods, so the sequence of draws has the same period as the state. These are short linear generators: consecutive outputs are linearly related, so they fail statistical tests that modern generators pass. They suit dithering, jitter and randomized rounding.

Width is part of the state. A recurrence narrower than 23 bits clears the state bits at and above its width on every step. While a generator keeps one recurrence, this cannot be seen. After a switch to a wider recurrence, the new taps read bits that the old one had cleared, so the sequence differs from a fresh start. prng_ref.py models this, and rows g2, g3 and odd switch to a wider recurrence.

Register writes take effect late. A store to the generator's state or configuration registers lands a few cycles after the store instruction, and the next draw does not wait for it:

• after a store to prng_state, leave at least two instructions before the next draw;
• after a store to prng_cfg, leave at least three.

These are the smallest counts that work in the simulator. The hardware documentation gives no count for these registers.

With fewer, the draw goes wrong in a way no message reports. With one instruction after the seed, the first draw can read a partly written state, and the seed's last word then overwrites the step that draw made, so every later value comes one draw late. prng_seed writes generator 3 last, so a kernel that reads generator 3 is the one hit. With two instructions after a recurrence change, the whole next draw loop still steps with the old recurrence. The compiler does not add this spacing: write the __builtin_nop() calls yourself, and check in the listing that nothing else moved between the store and the draw.

The starting state. When a program uses rand(), the compiler writes a fixed state, 0x211234, into all four generators in the image. Without prng_seed, every run, and every PE that runs the same kernel, therefore draws the same numbers.

Stochastic rounding. The same generators can also supply random rounding bits to the floating-point instructions. That mode is set in the floating-point control word. The compiler never sets that word, so the programs in this tutorial run with its default, round to nearest even.

The test

Row	N	prng	poly → poly2	What it covers
tutorial	16	1	1 → 3	generator 1; a switch to a narrower recurrence
g0	16	0	0 → 0	generator 0; the switch changes nothing
g2	16	2	2 → 0	generator 2; a switch to a wider recurrence
g3	16	3	3 → 1	generator 3, the last seed words; a wider recurrence
odd	7	0	2 → 1	partial lane groups at both widths
big	256	0	0 → 2	a long sequence; enough values for the uniformity check

check.py prints each array's comparison with the model. With 64 or more values, it also bins the values of u into 8 equal bins and fails if the counts are grossly uneven. That check only catches a broken mapping; the exact comparison already fixes every value.

check.py builds the host model from the row's seed and build parameters:

    cfg = prng_ref.cfg_for_build(args.prng, args.poly)
    seed = np.fromfile(args.seed_file, dtype=np.uint16)
    machine = prng_ref.Prng(seed, cfg)

and then draws from it in the kernel's order:

    ref16 = machine.lanes(n, 16)
    ref32 = machine.lanes(n, 32)
    machine.select_poly(args.prng, args.poly2)

machine.lanes(n, 16) gives the values that a draw loop of length n stores into an xp, as unsigned 16-bit integers, and machine.lanes(n, 32) those it stores into a u32. Each call continues the sequence where the one before stopped, as the kernel's loops do. select_poly makes the change that prng_select_poly makes in the kernel, and prng_ref.uniform_from_int16 maps 16-bit values to (0, 1) as hwrand.w does. A host program that imports prng_ref predicts a kernel's sequence the same way.

Run it

In your copy of the tutorial (see Lesson 1):

cd tutorial/ext-numerics/random-numbers
make TEST=tutorial          # one row
make TEST='*'               # every row
make clean
make TEST=tutorial SEED=42  # the same row with another seed

A passing row prints one line ending in PASS tutorial [ N=16 PRNG=1 POLY=1 POLY2=3 ARCH=sdr ]. test-tutorial/make-log then contains:

comparison:
  rnd16 : ok
...
  u     : ok
...
Ok!

A failing row says FAIL, and the log shows the first differing values, for example rnd16 : MISMATCH (12 of 16), followed by [ 0] device 0x... model 0x... lines.

After you edit a source file, run make clean before the next test: make reuses an existing test-* directory.

Pitfalls

• Too little spacing. The spacing rule is not checked by the compiler or by the hardware. A seed or a recurrence change that is one instruction too close gives plausible-looking numbers that are wrong.
• A zero seed is stuck. Zero is a fixed point of every recurrence, so a generator seeded with 0 emits the same value forever. A seed whose only set bits are at or above the recurrence's width is also stuck, because the first step clears them: {0, 0x40} sets only bit 22, which recurrences 1 to 3 clear. Give every generator a non-zero bit among state bits 19:0; gen.py draws its seeds from those 20 bits.
• rand() into an sp reads the 32 random bits as a float. The value can have any exponent and can be a NaN. Draw an integer and convert it, as u does.
• The same numbers on every PE. Without a seed per PE, all PEs that run one kernel draw the same sequence. Write a different seed into each PE, as lesson 8 writes per-PE data with miset order="x n".
• Every draw advances the generator. The three loops in hwrand.w share one sequence. Inserting a draw anywhere shifts every later value.

Exercises

1. Remove one __builtin_nop() from prng_seed, then run make clean and make TEST='*'. Only g3 fails, because it reads generator 3, whose words prng_seed writes last. In its log, rnd16 differs from element 4 on, and the device's values 8 to 11 are the model's values 4 to 7: after the first draw, every value comes one draw late. Restore the line.
2. Estimate π. Draw 2·N 16-bit values, map them to points (u, v) in (0, 1)², and count the points with u² + v² < 1 in an xp counter. Use if, not the value of the comparison, to count: a comparison stored as a value is -1 when true (lesson 2). Four times the count divided by N approaches π. Take u from the first N values and v from the last N: values next to each other share a draw, and are not independent. Add the expected count to check.py: after the calls that replay the kernel's earlier draws, machine.lanes(2 * n, 16) gives the 2·N 16-bit values of one vector read of length 2·N, which takes ⌈2·N/4⌉ draws.

Next

This is the last topic of the group. The next group, Single-PE Programming Techniques, starts with Topic — Vector-step loops, predication and memory iterators, which runs a loop body of several statements, conditions included, as hardware vector instructions.

References

• Unary functions (rand)
• Hardware register access (address())
• Imports
• Other built-ins (__builtin_nop)
• Architecture notes