how your computer fakes randomness

Here is a program that asks for random numbers, and here is what it prints:

#include <stdio.h>
#include <stdlib.h>
int main(){ for(int i=0;i<5;i++) printf("%d ", rand()); }

$ ./hook   $ ./hook   $ ./hook
16807 ...  16807 ...  16807 ...

Three fresh processes printed the same five integers:

16807 282475249 1622650073 984943658 1144108930
16807 282475249 1622650073 984943658 1144108930
16807 282475249 1622650073 984943658 1144108930

The program never called srand. Nothing in its source selected 16807, yet the first process and the third process agreed down to the last bit. A process restart had apparently erased every variable in the program and restored one hidden variable to the same value.

That repeat is useful. The matrix inputs in the CPU and GPU experiment can be regenerated after a bug. A failing test can keep the input that exposed it. The same property would be disastrous if these integers became password-reset tokens. The difference is not whether the output looks irregular on a screen. It is whether another program can reconstruct the state that produces the next output.

the state survived in the source

A pseudorandom number generator needs a memory of its past. Call that memory its state. Each call transforms the current state into a new state and derives an output. No fresh uncertainty has to enter during that transformation.

The first implementation I inspected was not xoshiro. It was the function that printed 16807. The current Apple Libc repository contains a do_rand function with one stored integer named next. Its non-compatibility path computes this recurrence:

$$x_{n+1} = 16807x_n \bmod (2^{31}-1)$$

The subscript only counts calls. If the current state is $x_n$, multiply it by 16,807, divide by 2,147,483,647, and retain the remainder. That remainder becomes both the next state and the returned integer.

Apple’s source avoids overflowing a 31-bit intermediate by using this identity:

$$2^{31}-1 = 127773 \times 16807 + 2836$$

It divides the old state into a quotient and remainder, computes 16807 * lo - 2836 * hi, and adds the modulus if the result is negative. This is an implementation of the Park and Miller generator described in the source comment, not a property C promises for every rand().

The rearrangement is sometimes called Schrage’s method. Write the old state as:

$$x_n = 127773q + r$$

where this local $q$ is the quotient and $r$ is below 127,773. The source names them hi and lo. Since:

$$16807 \times 127773 = (2^{31}-1)-2836$$

multiplying the decomposed state gives:

$$16807x_n = 16807r + q((2^{31}-1)-2836)$$

Modulo $2^{31}-1$, the whole $q(2^{31}-1)$ term disappears:

$$16807x_n \equiv 16807r - 2836q \pmod{2^{31}-1}$$

That is the source expression. Its intermediate products fit in a signed long on the supported platform. The direct product of two 31-bit values may not fit in the 31-bit arithmetic the historical implementation was designed around. This is an implementation constraint shaping the code, not a different generator.

State zero is another source-level edge. A multiplicative generator would leave zero at zero forever. Apple’s function replaces a zero state with 123459876 before applying the recurrence. Consequently, srand(0) does not create the all-zero stream, and it does not behave like srand(1). That mapping is implementation behavior absent from the portable C contract.

The source initializes next to 1. The POSIX contract for rand independently says that calling rand() before srand() must produce the same sequence as srand(1). POSIX fixes the default seed. It does not fix the recurrence or the sequence across C libraries.

I put the source recurrence beside both library paths:

RAND_MAX=2147483647
implicit state:   16807 282475249 1622650073 984943658 1144108930
after srand(1):   16807 282475249 1622650073 984943658 1144108930
source model:     16807 282475249 1622650073 984943658 1144108930
all three paths match

This was measured on arm64 macOS 15.7.7 with LibSystem 1351.0.0, using Apple clang 15 at -O2. The repository’s main branch is supporting source evidence rather than a proven source tag for that installed binary. The equality test ties the relevant recurrence to this machine; it does not prove every unseen detail of the shipped library came from the current repository revision.

one output owned the rest

The recurrence does more than explain why separate processes agree. The value returned by this implementation is the complete state. Once one output is visible, the next one is:

$$(16807 \times \text{visible output}) \bmod 2147483647$$

No seed search is required. A model with one observed value predicted every later value from a separately seeded C process:

secret seed: 14542
watcher sees first 3 outputs: [244407394, 1766337894, 27048330]
watcher predicts next 5: [1482232793, 1076246751, 224385376, 263730300, 108904692]
actually were         : [1482232793, 1076246751, 224385376, 263730300, 108904692]
all correct: True

The original probe watched three values because that made the transcript easy to recognize. It only needed the third value. The first two added no information to the prediction.

An attacker does not need to decide whether those integers look random. The implementation has already exposed the one integer from which its entire future is computed.

four words hid a larger state

Replacing rand() with a larger state and a better transformation changes the statistical quality and the amount of work required to recover state. Here is the state transition for xoshiro256**, copied into a small C program from the authors’ reference implementation:

uint64_t s[4];                                    // the entire state: four 64-bit numbers
uint64_t next(){
    uint64_t r = rotl(s[1]*5,7)*9, t = s[1]<<17;  // scramble the state into an output
    s[2]^=s[0]; s[3]^=s[1]; s[1]^=s[2];           // shuffle the state for next time
    s[0]^=s[3]; s[2]^=t; s[3]=rotl(s[3],45);
    return r;
}

first 4 outputs: 11520 0 1509978240 1215971899390074240
same seed again: 11520 0 1509978240 1215971899390074240

I assigned all four state words directly to 1, 2, 3, 4. That is a full-state initialization, not a realistic one-integer seeding interface. Resetting those words resets the future.

The state transition consists of exclusive-or, shift, and rotation operations. At the level of individual bits, each next-state bit is an exclusive-or of selected current-state bits. This is a linear transformation over two possible bit values. The chosen transformation cycles through all $2^{256}-1$ nonzero states, as established in Blackman and Vigna’s scrambled linear generator paper.

The returned value is not a raw state word. It is rotl(s[1] * 5, 7) * 9, which the ** suffix names. Odd integers have multiplicative inverses modulo $2^{64}$. Multiplication by 5 and 9 can therefore be undone. Rotation can also be undone. One output reveals the value of s[1] at that step.

Four consecutive outputs reveal four successive values of s[1], or 256 observed bits. The transition between them is linear. I constructed a 256 by 256 binary matrix by placing one bit at a time into the initial state and recording which observed bits it affected. Gaussian elimination over bits recovered the complete initial state:

observed first four outputs:
  0x6666666666666c65
  0xd90633608dbae0aa
  0xdc3f17eac117d726
  0x6507a6dc8a5dc138
observation matrix rank: 256
recovered initial state:
  0x0123456789abcdef 0xfedcba9876543210
  0x0f1e2d3c4b5a6978 0x8877665544332211
state exact: True
next eight outputs predicted exactly: True

The inversion starts from the last multiplication. Arithmetic on uint64_t discards overflow, which is arithmetic modulo $2^{64}$. An integer $a$ has a multiplicative inverse in that arithmetic exactly when it is odd. Both 5 and 9 are odd, so Python can compute:

inv9 = pow(9, -1, 1 << 64)
inv5 = pow(5, -1, 1 << 64)

def recover_s1(output):
    before_times_9 = output * inv9
    before_rotate = rotate_right(before_times_9, 7)
    return before_rotate * inv5

Every operation is masked back to 64 bits. Applying recover_s1 to one output gives the exact s[1] used to create it. This does not reveal the other three state words immediately.

The next-state transformation supplies the missing equations. At the bit level, exclusive-or is addition modulo 2:

0 xor 0 = 0
0 xor 1 = 1
1 xor 0 = 1
1 xor 1 = 0

A shift or rotation only changes which old bit contributes to which new position. If the 256 unknown initial bits are collected in a vector $s$, the four recovered s[1] words can be collected in another 256-bit vector $o$. There is a binary matrix $T$ such that:

$$o = Ts \pmod 2$$

I did not derive all 65,536 matrix entries by hand. For input bit 0, initialize the state with only bit 0 set, advance it four times, and record the four observed s[1] words. That produces column 0. Repeat for every input bit:

columns = []
for bit in range(256):
    basis = [0, 0, 0, 0]
    basis[bit // 64] = 1 << (bit % 64)
    columns.append(observe_four_s1_words(basis))

Because the state transition is linear, exclusive-oring any set of these basis states exclusive-ors their observations too. The columns therefore describe every possible state, not merely the 256 test states used to construct them.

Gaussian elimination uses exclusive-or where ordinary elimination uses subtraction. It found a pivot in every one of the 256 columns. The matrix rank was 256, so the four observations determined one unique initial state for this observation layout. Had the rank been 240, sixteen state bits would have remained unconstrained and more outputs would have been necessary.

The recovery program generated twelve outputs from a state known only to the generator side, retained four as observations, solved for the original state, and regenerated all twelve. Comparing the eight later values separates reconstruction from a solver that merely fits the four values it was given.

This is not a practical claim that every observer of every generator gets state after four values. It uses the exact public xoshiro256** algorithm and full 64-bit outputs. Truncated outputs require more work. A generator with a cryptographic output function changes the state-recovery problem again.

It does establish a narrower point. A state space with 78 decimal digits and a period of $2^{256}-1$ can still be predictable after a few observations. Period length, visible statistical quality, and state-recovery resistance are separate properties.

The recovery also explains why changing only the output scrambler cannot create cryptographic security casually. A reversible scrambler hides a state word from sight but not from algebra. A non-reversible output function may leave many possible states, yet later outputs can eliminate candidates. Cryptographic generator design asks for computational resistance against that entire observation process, including partial state compromise and large output volumes. Statistical balance is only one constraint.

the low bit refused the easy explanation

The first draft blamed rand() for the most familiar linear-congruential failure: an alternating lowest bit. I had imported a warning about generators whose modulus is a power of two and attached it to a generator whose modulus is prime. Counting the bits rejected that explanation:

low bit (rand&1): 11100000110100111100
low bit flip rate over 200000 draws: 0.4985 (strict alternation would be 1.0000)
low bit fraction of 1s: 0.5026 (fair=0.5000)

Strict alternation would produce 101010... or 010101... and a flip rate of 1. The measured rate was 0.4985. The fraction of ones was 0.5026. Those two small tests found no defect in that bit.

They do not prove the bit is good. Balance asks one question: were zeros and ones equally frequent in this sample? Flip rate asks another: how often did adjacent bits differ? A stream can satisfy both while obeying a more distant linear relation. Xoshiro’s authors report that the lowest bits of the + variant fail linear-complexity tests even though the same balance and lag-one probes pass. An easy test can reject a stream. Passing an easy test says only that the stream escaped that test.

The alternating-bit result does apply to a common LCG form with modulus $2^k$ and an odd increment. Reducing the recurrence modulo 2 leaves only the lowest bit. Multiplication by an odd number preserves that bit and adding an odd increment flips it. The low bit therefore alternates. Park-Miller uses the prime modulus $2^{31}-1$ and has no additive increment. That proof never applied here.

rand() on this machine still exposes its state, still has a visible modular recurrence, and still repeats. The low-bit accusation was simply the wrong one.

a full period was still too short

A deterministic finite-state machine must eventually revisit a state. Once it does, every later output repeats too. A multiplicative LCG with prime modulus $m$ can visit at most $m-1$ nonzero states. Whether it reaches all of them depends on the multiplier. I built a version modulo 101 small enough to watch the wrap:

tiny LCG (m=101,a=2): returned to seed after 100 steps (m-1=100)
tiny LCG (m=101,a=95): period 5 (short! a bad multiplier)

With a = 2, every nonzero state appears once before the seed returns. The multiplier 2 is a primitive root modulo 101, so its powers generate the entire nonzero multiplicative group. With a = 95, the seed returns after five transitions.

For the Apple recurrence, 16,807 gives the full nonzero period:

$$2^{31}-2 = 2,147,483,646$$

The current bulk-fill measurement produced about 274 million rand() values per second. At that measured rate, exhausting one full period would take about 7.8 seconds. The loop would then reuse the exact sequence. A program does not have to store all two billion outputs for the repeat to matter. A long simulation can accidentally compare later work with an earlier copy of the same stream.

Xoshiro256** raises the period to $2^{256}-1$. Its authors also provide jump functions that advance the state by a vast fixed distance without generating every skipped value. That is useful for assigning streams to parallel workers. It is not a cryptographic defense. The four-output recovery already crossed almost the entire period in the only direction an observer cared about: from present output to future output.

the seed was smaller than the state

A seed is input to an initialization procedure. It need not contain the entire state, and expanding it does not add information. The common srand(time(NULL)) path compresses the clock to the unsigned seed accepted by srand. time() exposes whole seconds on this path, so two initializations that see the same second receive the same seed:

same-second seed 1783432282
  process A model: 1717102395 1454704379
  process B model: 1717102395 1454704379
  streams equal: yes

These are two initializations with the same modeled clock value, not two wall-clock launches. The result follows from determinism. If both calls enter the same state, every output agrees.

Knowing that a service started within a ten-second window leaves at most eleven candidate seeds if it used whole seconds. A watcher can generate the first value under all eleven and identify the matching future. Adding a process identifier or nanosecond counter increases the search space, but predictable metadata is still not entropy merely because it has more digits.

The search cost depends on what the observer knows, not on the type width. An unsigned 32-bit seed has more than four billion representations. If logs reveal a launch timestamp within one second and the code seeds from a nanosecond clock, scheduling and clock behavior may leave far fewer than one billion plausible readings. If a container starts processes with predictable identifiers, exclusive-oring the PID into the timestamp does not multiply uncertainty by the full PID range.

Several weak values can even describe the same underlying event. Wall-clock seconds, wall-clock nanoseconds, process start time, and a timestamp embedded in a request may all be correlated observations of one moment. Adding their bit widths double-counts uncertainty.

The prediction test for a time seed has a direct form:

for candidate in range(observed_second - 5, observed_second + 6):
    if first_output(candidate) == visible_output:
        print("candidate seed", candidate)

For this Park-Miller output function, one visible output usually identifies the seed directly because the multiplier is invertible modulo the prime. Searching the window is useful when an API transforms, truncates, or hides the first outputs.

A unique time seed is still not necessarily a statistically sound parallel seed. If a thousand workers initialize from successive clock readings, their seeds encode launch order and may enter related parts of a poorly initialized state space. Conversely, a fixed master seed with a sound stream-allocation scheme can be better for simulation even though it contains no fresh uncertainty. The purpose decides which property matters.

For test reproduction, deriving a seed from the clock is also self-defeating. The failing input disappears unless the selected seed is printed and retained. A better test runner accepts a seed as an argument, prints it on failure, and chooses a documented default. A separate mode can obtain a seed from the system when input diversity is desired, then immediately record that value as part of the test artifact.

The xoshiro example needs four 64-bit words. Its authors recommend SplitMix64 to initialize those words from a 64-bit seed. I ran adjacent seeds through that expansion:

SplitMix64 expansion from adjacent 64-bit seeds:
  word 0: 910a2dec89025cc1  975835de1c9756ce
  word 1: beeb8da1658eec67  bfc846100bfc1e42
  word 2: f893a2eefb32555e  987bbcbfdd7e532f
  word 3: 71c18690ee42c90b  c3f2827affe7f664

The four words look unrelated. They are still one of at most $2^{64}$ possible expansions because only $2^{64}$ seeds entered the function. SplitMix64 disperses nearby seed values so the target generator does not begin in visibly related states. It cannot turn 64 uncertain bits into 256 uncertain bits.

That distinction changes how experiments should be recorded. Reproducing an output requires the generator algorithm and version, the complete initialization procedure, the seed, and the way raw integers were mapped into the values the program consumed. “Seed 42” is incomplete if a library later changes any of the other three.

There are two different reasons to expand a seed. The first is diffusion: adjacent input seeds should not produce states whose early outputs remain visibly related. SplitMix64 is different in structure from xoshiro’s linear recurrence and spreads changed input bits across each output word. The two four-word rows above demonstrate diffusion for seeds 1 and 2, but two examples do not establish every statistical property of initialization.

The second possible reason would be adding uncertainty. Deterministic expansion cannot do that. If an initialization function maps a 64-bit seed to a 256-bit state, at most $2^{64}$ of the $2^{256}$ states are reachable through that interface. The reachable fraction is:

$$\frac{2^{64}}{2^{256}} = 2^{-192}$$

That tiny fraction is not automatically a problem for a reproducible simulation. A well-designed mapping can place those reachable states far apart and the period can remain enormous from each one. It does matter if someone argues that the 256-bit state contains 256 bits of secret uncertainty. The secret input selected only 64 bits’ worth of alternatives.

Zero needs explicit handling too. Xoshiro’s all-zero state maps to itself forever because shift, rotate, and exclusive-or of zero remain zero. Its claimed period excludes that one state. Directly zero-initialized static storage would therefore produce zeros forever. The initialization procedure must prevent all four words from becoming zero together.

Parallel seeds add another distinction. Giving worker $i$ the integer seed master + i is reproducible, but absence of stream overlap or correlation depends on the initialization function and generator family. A jump polynomial uses the known linear transition to advance by a fixed number of steps exactly, without looping over every skipped output. Xoshiro256**‘s published jump advances as if $2^{128}$ calls occurred. Beginning each worker at successive jump points divides a period of roughly $2^{256}$ into vast named regions.

No-overlap is not independence. It prevents two workers from consuming the same state positions. Their streams still arise from one algebraic recurrence. The authors analyze overlap probability for randomly selected subsequences and give the bound:

$$\Pr(\text{overlap}) \le \frac{n^2L}{P}$$

where $n$ is the number of streams, $L$ is the length consumed by each, and $P$ is the period. Units matter here. Both $L$ and $P$ count generator steps, while $n$ is a count of streams. The ratio is dimensionless, as a probability bound must be.

For one million workers, one billion values per worker, and period near $2^{256}$:

$$\frac{(10^6)^2(10^9)}{2^{256}} \approx 8.6\times10^{-57}$$

That bound concerns overlap under its stated random-start assumptions. It does not say the generated values pass a statistical battery, resist state recovery, or model the application’s randomness correctly. One formula closes one failure mode.

the child inherited the same future

State is ordinary process memory unless an implementation takes special measures. fork() creates a child whose address space initially contains a copy of the parent’s memory. I seeded rand, consumed one value, then forked. Parent and child each asked for the next five:

rand parent: 1790989824 2035175616 77048696 24794531 109854999
rand child : 1790989824 2035175616 77048696 24794531 109854999
inherited rand stream equal: yes

The result can be predicted before running it. Both processes begin after the same consumed output. Their copies of next are equal. They execute the same recurrence five times.

This matters beyond demonstrations. A worker manager may initialize a simulation once and fork eight workers. Without per-worker stream assignment or post-fork reseeding, all eight can spend their first samples on identical work. The nominal sample count grows by eight while the number of independent samples does not. Confidence intervals then shrink on paper without gaining the information they claim.

One response is to derive each worker’s stream from an explicit master seed and stable worker identifier. Another is to use a generator’s jump operation to allocate disjoint regions. A thread scheduler’s current execution order is not a stable identifier. If tasks pull from one shared stream, a timing change can send different random values to different tasks even when the global sequence is unchanged.

For a reproducible parallel experiment, I want the mapping to be a function of recorded facts:

$$\text{worker state} = \operatorname{initialize}(\text{master seed},\ \text{worker id})$$

That equation is a design contract, not proof of independent streams. The initialization or jump construction still needs analysis. It does prevent the accidental claim that “same master seed” alone determines a parallel run whose task interleaving changes.

the scheduler reassigned the same numbers

Fork duplication is visible because two processes print identical sequences. A shared mutable generator creates a subtler replay failure even when it returns one valid global sequence.

Take four consecutive generator outputs and two workers, A and B. Under schedule AABB, worker A consumes positions 0 and 1 while B consumes 2 and 3. Under ABAB, A consumes 0 and 2 while B consumes 1 and 3. The generator’s global output is identical. The values attached to each worker’s computation change.

A durable probe used SplitMix64 only as a simple deterministic sequence for this assignment experiment:

shared mutable stream, schedule AABB:
  A=[16294208416658607535, 10451216379200822465]
  B=[10905525725756348110, 2092789425003139053]
shared mutable stream, schedule ABAB:
  A=[16294208416658607535, 10905525725756348110]
  B=[10451216379200822465, 2092789425003139053]
shared worker results stable: False

Adding a mutex around next() would serialize state access and prevent a data race. It would not make the assignment reproducible. The mutex determines one total call order, but the scheduler still influences which worker wins each call.

One repair is mutable state per worker. Another is a counter mapping in which the value is a deterministic function of a master key, stable worker identifier, and per-worker draw index. The probe used:

$$\text{value} = f((\text{worker id} \ll 32)\mathbin{\vert}\text{draw index})$$

and produced the same per-worker values under both schedules:

per-worker counter mapping, schedule AABB:
  A=[14135772400868000056, 2324861979054413167]
  B=[16695506628682495282, 14160868677571091529]
per-worker counter mapping, schedule ABAB:
  A=[14135772400868000056, 2324861979054413167]
  B=[16695506628682495282, 14160868677571091529]
counter-mapped worker results stable: True

The particular SplitMix mapping is not offered as a cryptographic or parallel simulation design. It isolates the ownership property. A real counter-based generator must specify how keys and counters are partitioned, analyze collisions and correlations, and define what happens when a counter wraps.

Stable worker identifiers also require care. A transient thread number, rank assigned by launch order, or retry attempt can change between runs. A task identifier derived from logical input partitions is more durable. If work is stolen between threads, the random stream should usually travel with the logical task rather than the executing thread.

Reproduction after a checkpoint needs current state, not merely the original seed. Replaying from a seed requires consuming the exact same number of values in the exact same order before reaching the checkpoint. Adding one diagnostic sample or changing a branch shifts every later value. Saving the generator state with application state resumes from the actual position. For a parallel run, that means one state or counter position per logical stream, plus the generator and mapping versions.

different bytes did not prove unpredictability

The fork probe also filled five words through arc4random_buf:

arc4random_buf parent: 41d3c47f 6cbebe1b 204e7321 3cd32506 1faabfb0
arc4random_buf child : af8e621b 620d59bd 55c0c1ce 9b708225 dbccc448
secure buffers equal: no

That difference is a control result, not a cryptographic proof. A time-seeded rand() would usually differ too. The stronger evidence is the contract installed with the macOS 15.7 SDK. Its arc4random(3) manual describes a cryptographic pseudorandom generator, says the subsystem is reseeded from the kernel regularly and upon fork, and recommends arc4random_uniform when a bounded value is needed.

The phrase cryptographic pseudorandom generator contains the correction the earlier draft missed. Cryptographic output is still produced by deterministic state transitions between reseeds. The operating system does not need one fresh physical event for every returned word.

NIST SP 800-90A specifies deterministic random bit generators. SP 800-90B treats the entropy sources used to seed them. The split exists because collecting uncertain physical observations and expanding secret state are different jobs.

A simplified system has three paths:

entropy source observations
        |
        v
conditioning and health tests
        |
        v
instantiate or reseed secret generator state
        |
        v
deterministic generate requests ------> application bytes

The first path is allowed to be slow, irregular, and hardware-dependent. Conditioning compresses noisy observations into material suitable for initialization. Health tests try to detect an entropy source stuck in a failure mode. The generate path can then return far more bits than were collected in that moment because security relies on the secrecy and cryptographic evolution of state.

This model has a state-compromise question absent from the original article. If an attacker reads current generator state, can they reconstruct old outputs, predict future outputs, or both? Updating state through a one-way construction after generation can provide backtracking resistance for erased earlier states. Mixing fresh entropy after compromise can restore future unpredictability. The exact guarantees depend on the generator and reseed design. I did not inspect Apple’s current internal implementation deeply enough to assign it a specific state-compromise guarantee.

The local arc4random(3) contract does establish the parts used here: a cryptographic pseudorandom generator, process-wide data pool, periodic kernel reseeding, and reseeding upon fork. It also says the historical RC4 implementation was replaced by an AES-based design in OS X 10.12 and may change again. That warning is a reason to program to the interface contract rather than serialize or infer its private state.

A small ChaCha20 block implementation makes deterministic cryptographic expansion visible without claiming that it is Apple’s current internal algorithm. I used the key, nonce, counter, and expected 64-byte block from RFC 8439:

RFC 8439 block vector matched: True
same key, nonce, and counter repeated exactly: True
one key bit changed 248 of 512 output bits
deterministic output can still be computationally unpredictable when the key remains secret

The first two lines are not a contradiction. Anyone with the same key and inputs can reproduce the block. An observer without the secret key is not known to have a practical route from output back to state. That computational asymmetry is the property missing from rand() and xoshiro256**. Different output on two runs is neither necessary nor sufficient to establish it.

ChaCha20 also shows why a nonce is not a seed. A nonce may be public, but it must obey the construction’s uniqueness rules with a given key. Reusing a key and nonce restarts the same keystream. The secret uncertainty remains in the key. Combining a predictable nonce with a secret key is normal. Combining a public key with a “random-looking” nonce supplies no secret state.

The 248 changed bits in the probe are an avalanche observation from one input pair, not a security test. A deliberately weak hash can change half its output bits when one input bit changes. The evidence that matters for ChaCha20 here is conformance to the RFC vector and the construction’s cryptographic analysis, not an attractive Hamming-distance count.

the collision arrived before the sequence ended

Unpredictability and uniqueness are different requirements. A 32-bit cryptographic output can be difficult to predict before it appears and still be far too small for a fleet of session identifiers.

With $S=2^b$ equally likely values and $n$ independent draws, the first value cannot collide with an earlier one. The second avoids the first with probability $(S-1)/S$. The third must avoid two existing values, so its probability is $(S-2)/S$. Continuing gives:

$$\Pr(\text{no collision}) = \frac{S}{S} \times \frac{S-1}{S} \times \frac{S-2}{S} \times \cdots \times \frac{S-n+1}{S}$$

When $n$ is much smaller than $S$, this is approximately:

$$\Pr(\text{at least one collision}) \approx 1-\exp\left(-\frac{n(n-1)}{2S}\right)$$

Set the collision probability to one half and solve for $n$:

$$n \approx \sqrt{2S\ln 2}$$

For 32-bit values, $S=4,294,967,296$ and the halfway point is about 77,163 draws. The generator may have an immense internal period. The application threw most of that period away when it retained only 32 output bits.

An attacker’s guessing problem uses a different scale. If one valid secret token is uniformly selected from $2^b$ possibilities and each online guess tests one candidate, a single guess succeeds with probability $2^{-b}$. Roughly $2^{b-1}$ guesses are needed on average to find that one fixed token by exhaustive search. Collisions among many generated tokens appear around $2^{b/2}$ because every pair is a collision opportunity.

A biased mapping reduces both safety margins. If some tokens have more preimages, those tokens appear more often and an attacker can try high-probability values first. Counting string length without counting the distribution of generated strings overstates entropy.

NIST’s entropy-source terminology distinguishes output length from min-entropy, the negative base-two logarithm of the most likely outcome. If a 128-bit seed is selected from only one million plausible clock and process combinations, its representation is 128 bits wide but its input uncertainty is at most:

$$\log_2(1{,}000{,}000) \approx 19.93\text{ bits}$$

Hashing or running SplitMix64 over those candidates can spread them across a 128-bit or 256-bit state. An observer can still enumerate the one million inputs and reproduce the expansion.

This changes the choice of interface. Simulation wants named algorithms, recorded seeds, separate streams, and exact replay. Security wants a maintained system CSPRNG, enough output bits for both guessing and collision requirements, automatic handling of fork and reseeding, and no seed written to an experiment log. Treating both jobs as “get a random number” hides the property each one would sacrifice.

the call boundary was the expensive part

The old benchmark compared scalar xoshiro calls with scalar arc4random() calls and small buffered device reads, then described the gap as the cost of physical entropy. I changed the experiment to fill the same four-million-word buffer in five ways. Each method ran eleven times. Method order was shuffled by a separate fixed control generator on every repetition. The timed region covered the fill, and a checksum consumed values after timing.

4000000 uint32_t values (15.26 MiB), 11 repetitions, randomized order
  rand fill                median=  3.655 ns/value, MAD= 0.012
  xoshiro256** fill        median=  0.732 ns/value, MAD= 0.007
  arc4random scalar fill   median= 25.262 ns/value, MAD= 0.162
  arc4random_buf bulk      median=  0.785 ns/value, MAD= 0.011
  /dev/urandom bulk read   median=  6.047 ns/value, MAD= 0.024

These measurements describe one arm64 Mac, Apple clang 15, one 15.26 MiB output size, warm caches, and one OS release. CPU affinity and frequency were not controlled. Median absolute deviation, or MAD, is the median distance from the median. It records the typical within-run spread without pretending the eleven timings form a universal performance estimate.

Scalar arc4random() was about 34.5 times slower per stored word than the xoshiro fill. arc4random_buf was only about 7 percent slower. The bulk interface amortized the call boundary and allowed the system implementation to generate into a contiguous destination. Reading the device in bulk was also far below the earlier scalar-like number, although slower than arc4random_buf here.

The final verifier reran the complete benchmark and measured medians of 3.644, 0.733, 25.785, 0.858, and 6.033 ns per value in the same method order shown above. The bulk cryptographic path was 17 percent slower than xoshiro in that run rather than 7 percent. The scalar-to-bulk direction remained large. The exact bulk gap did not. Both runs are retained because the smaller number would otherwise make the comparison look more stable than it was.

Converting nanoseconds per 32-bit value into output bandwidth keeps the units visible:

$$\text{GiB/s} = \frac{4\ \text{bytes/value}} {\text{seconds/value}} \times \frac{1\ \text{GiB}}{2^{30}\ \text{bytes}}$$

The medians correspond to roughly 5.09 GiB/s for xoshiro, 4.75 GiB/s for arc4random_buf, 1.02 GiB/s for rand, 0.62 GiB/s for the device read, and 0.15 GiB/s for scalar arc4random. These are bytes written into the benchmark array. They are not measured DRAM bandwidth. Cache allocation, later eviction, and writeback are outside the simple division.

The source resets rand and xoshiro state before each timed fill so all repetitions perform the same state path. It opens /dev/urandom once before measurement, then reads until the entire buffer is filled, handling interrupted reads. The device descriptor is not reopened for every word. The scalar paths store every result into the same array as the bulk paths.

The checksum samples two positions after each timed fill. It prevents the entire buffer from being an obviously dead result at the program boundary, but it does not validate every byte or compare distributions. Separate functional probes establish known xoshiro and rand sequences. The cryptographic interfaces deliberately cannot be checked against a fixed expected buffer.

The measurement also conflates generation and memory stores by design. That matches the concrete task “fill this buffer.” It does not isolate the minimum latency of one generated word held in a register. Blackman and Vigna’s published microbenchmarks warn that generator timing inside an application depends on register allocation, vectorization, and how state is embedded. Replacing a generator in the actual consumer remains the more relevant experiment.

The table does not measure entropy collection cost. It measures public interfaces after the system has initialized their hidden machinery. It also does not prove arc4random_buf is interchangeable with xoshiro for simulation. One is designed for cryptographic unpredictability and automatic reseeding. The other provides explicit reproducibility, jumpable state, and stable algorithmic control. Similar throughput did not merge those contracts.

the endpoint rounded upward

A generator’s raw output is rarely the value an application wants. Simulations ask for floating-point values in $[0,1)$. Shuffles ask for an integer below a changing bound. A generator can have a uniform raw output while the conversion on top introduces bias or an illegal endpoint.

The tempting floating conversion is:

double u = (double)next64() / 0x1p64;

The intended interval excludes 1. A 64-bit unsigned integer has 64 significant bits. A double carries 53. Before the division, converting a large integer to double rounds away eleven low bits. UINT64_MAX is close enough to $2^{64}$ that round-to-nearest produces $2^{64}$ exactly:

UINT64_MAX direct / 2^64 = 1
naive reached 1.0: yes

That endpoint can turn an array index floor(u * count) into count, one past the last element. It can also break a transformation whose derivation assumes $u < 1$.

The xoshiro authors document a conversion that takes the upper 53 bits first:

double u = (next64() >> 11) * 0x1.0p-53;

Now the integer numerator ranges from 0 to $2^{53}-1$, and every one of those exactly representable values receives the same number of raw 64-bit inputs:

UINT64_MAX upper 53 bits = 0.99999999999999989
upper-53 mapping stayed below 1.0: yes

The authors explicitly warn that direct 64-bit division does not produce the desired floating distribution because of round-to-even. The conversion rule is part of the generator contract. Recording the raw generator and seed while omitting it is not enough to reproduce a floating-point experiment.

There is more than one reasonable meaning of “uniform floating-point value.” The upper-53-bit construction is uniform over the grid:

$$\left\{ 0,\frac{1}{2^{53}},\frac{2}{2^{53}},\ldots, \frac{2^{53}-1}{2^{53}} \right\}$$

It is not uniform over all representable double values in $[0,1)$. Floating-point numbers are much denser near zero than near one because their exponent changes spacing. Giving every representable value equal probability would heavily favor tiny magnitudes. Giving intervals of equal real width equal probability naturally gives more probability to widely spaced values near one.

For most simulation APIs, uniform over an evenly spaced 53-bit grid is the intended contract. A tail-sensitive application may want a construction capable of producing much smaller positive values with probabilities proportional to the intervals they represent. That choice changes inverse transforms such as -log(u).

Endpoint conventions also differ:

• $[0,1)$ includes zero and excludes one;
• $(0,1]$ excludes zero and includes one;
• $(0,1)$ excludes both;
• $[0,1]$ includes both.

Inverse exponential sampling needs a value greater than zero before applying log. An array-index mapping needs a value below one. Converting one base uniform convention into another by clamping creates point masses at the clamp value. Rejection or an explicit integer construction preserves a cleaner distribution.

These distinctions are invisible in a short histogram. They appear at endpoints and in extreme tails, exactly where a finite sample has the least power to notice them.

six did not divide the byte

For a die roll from 0 through 5, the tempting integer conversion is next() % 6. It is exact only when the raw generator has a number of possible values divisible by 6. I counted the failure first on a ten-value source mapped to three buckets:

range 0..9 mapped with %3:  bucket0=4  bucket1=3  bucket2=3   (0 is favored)

The values 0, 3, 6, and 9 map to bucket 0. The other buckets receive three values each. The modulo operation made a uniform ten-value source lopsided because $10 = 3 \times 3 + 1$.

For an eight-bit source and six faces:

$$256 = 6 \times 42 + 4$$

Four faces receive 43 preimages. Two receive 42. Before sampling, this predicts that each favored face is $43/42 - 1$, or about 2.38 percent, more likely than an unfavored face.

The toy made the shape of the bias visible; I wanted to watch it bite at scale, and watch the fix erase it. So I built a generator that draws a uniform byte, 0 to 255, and rolled a die from it with byte % 6, sixty million times. The bias should be small but not invisible: 256 isn’t a multiple of 6, 256 % 6 = 4, so the four lowest faces (0, 1, 2, 3) each get one extra byte mapping to them and faces 4 and 5 don’t:

naive % (8-bit gen, die 0..5), 60000000 rolls:
  face 0: 10238152  (2.382% off fair)
  face 1: 10235367  (2.354% off fair)
  face 2: 10241740  (2.417% off fair)
  face 3: 10245267  (2.453% off fair)
  face 4: 9994748  (-0.053% off fair)
  face 5: 9997662  (-0.023% off fair)
  chi-square=23074.1  (5 dof, >11.07 = reject fairness at p=0.05; >>that here)

Faces 0 through 3 measured between 2.35 and 2.45 percent above the fair count. The prediction arrived before the run. The chi-square statistic compares squared count errors with the counts expected under a fair six-way result. With five independent degrees of freedom, 23,074 is incompatible with ordinary sampling variation.

The values 0 through 251 form 42 complete copies of the six remainders. Values 252 through 255 form the overhang. Reject that overhang and draw again:

rejection-corrected (reject v>=252, redraw):
  face 0: 9999944  (-0.001% off fair)
  face 1: 9997763  (-0.022% off fair)
  face 2: 10002937  (0.029% off fair)
  face 3: 10006946  (0.069% off fair)
  face 4: 9994748  (-0.053% off fair)
  face 5: 9997662  (-0.023% off fair)
  chi-square=9.49  (should be a few, order dof=5)
  values rejected & redrawn: 952936 (1.56% of draws)

Chi-square fell to 9.49. The correction redrew 1.56 percent of candidate bytes, close to the predicted $4/256 = 1.5625$ percent.

For a full 64-bit source and bound 6, an overflow-safe rejection threshold can be computed as -bound % bound. Values below that threshold are discarded. In six million measured draws, the threshold was 4 and no rejection occurred. Four rejected values out of $2^{64}$ is too small a probability to expect one in that run. The algorithm remains correct even though this sample never exercised its rare branch.

The unsigned expression looks backwards until its arithmetic is written out. In 64-bit arithmetic, -bound wraps to $2^{64}-\text{bound}$. Therefore:

$$(-b)\bmod b = (2^{64}-b)\bmod b = 2^{64}\bmod b$$

Call that remainder $t$. The range from $t$ through $2^{64}-1$ contains:

$$2^{64}-t$$

values, which is divisible by $b$ by construction. Each remainder below $b$ then has exactly $(2^{64}-t)/b$ accepted preimages.

A compact implementation is:

uint64_t bounded(uint64_t bound) {
    uint64_t threshold = -bound % bound;
    for (;;) {
        uint64_t value = next64();
        if (value >= threshold)
            return value % bound;
    }
}

This function requires bound > 0. With bound == 0, the modulo operation is invalid. It also assumes next64() is uniform over every 64-bit value. If the raw source has bias, equalizing the preimage counts of % bound does not remove that upstream bias.

The expected number of candidate draws is:

$$\frac{1}{1-t/2^{64}}$$

For $b=6$, $t=4$, so the excess over one draw is about $2.17\times10^{-19}$. For a worst-case bound just above $2^{63}$, almost half the raw values may be rejected and the expected count approaches two. The same correct algorithm can therefore have a data-dependent cost that matters in a tight loop.

Modulo also has a security timing question. The number of iterations can reveal that a rejection occurred. For ordinary bounded selection this is normally acceptable because the bound is public and each raw value is secret but ephemeral. A cryptographic protocol with a secret-dependent bound or a more complex rejection condition needs its own side-channel analysis. “Unbiased” describes the output distribution, not constant-time execution.

The local macOS manual says arc4random_uniform avoids modulo bias and may iterate more than once. Other libraries may use rejection, multiply-high methods, or a combination. The necessary property is equal preimage count for every result, not one mandatory implementation.

twenty-seven paths could not share six endings

The wrong three-item shuffle swaps each position with any of the three positions. It makes three choices, each with three possibilities. That creates $3^3 = 27$ equally likely execution paths. There are $3! = 6$ possible final orderings. Twenty-seven does not divide by six.

I enumerated every path before simulating:

exact execution-path counts
  (0, 1, 2): wrong=4/27 Fisher-Yates=1/6
  (0, 2, 1): wrong=5/27 Fisher-Yates=1/6
  (1, 0, 2): wrong=5/27 Fisher-Yates=1/6
  (1, 2, 0): wrong=5/27 Fisher-Yates=1/6
  (2, 0, 1): wrong=4/27 Fisher-Yates=1/6
  (2, 1, 0): wrong=4/27 Fisher-Yates=1/6

Three orderings receive five paths. Three receive four. Their exact probabilities are therefore 18.519 percent and 14.815 percent. A six-hundred-thousand-run simulation landed between 18.441 and 18.604 percent for the favored group, and between 14.766 and 14.812 percent for the other group.

Fisher-Yates chooses among three positions for the first element, two remaining positions for the second, and one for the last. It has $3 \times 2 \times 1 = 6$ paths, one per ordering. The measured orderings ranged from 16.619 to 16.728 percent, around the exact one-sixth target.

The wrong shuffle was not repaired by a better generator. Its bias came from the structure of the transformation. Replacing rand() with cryptographic bytes would reproduce the same 4-to-5 path imbalance more securely.

the logarithm changed the shape

A uniform value assigns equal probability to equal-width intervals. Many models need something else. Suppose $X$ is a waiting time with rate 1 and the probability of waiting beyond $x$ should be $e^{-x}$. Begin with a uniform value $U$ in $(0,1]$ and set:

$$X = -\ln U$$

Then:

$$\Pr(X > x) = \Pr(-\ln U > x) = \Pr(U < e^{-x}) = e^{-x}$$

The last equality uses the definition of a uniform value: the probability of landing below a number between zero and one equals that number. The algebra predicts the complete survival curve, not merely a mean.

One million transformed values gave:

exponential inverse transform, N=1000000
  mean=0.99898, analytic=1
  variance=0.99647, analytic=1
  P(X>1)=0.36786, analytic=0.36788
  P(X>3)=0.04932, analytic=0.04979

The run used Python 3.14.6’s random.Random with seed 20260725. These are direct measurements against analytically derived targets. The last probability differs by 0.00047, which is plausible sampling error for about 49,790 expected tail observations. It is not exact agreement, and no exact agreement should be expected from a finite run.

If the raw uniform can equal zero, log(0) is negative infinity. The probe redraws zero. A different library might construct a uniform interval that excludes zero directly. That endpoint choice belongs in the transformation’s contract.

the tail ended at 8.57

The Box-Muller transform makes two standard normal values from two independent uniform values:

$$R = \sqrt{-2\ln U_1}, \qquad \theta = 2\pi U_2$$ $$Z_0 = R\cos\theta, \qquad Z_1 = R\sin\theta$$

The geometry starts from a two-dimensional Gaussian. Its direction is uniform around a circle, which supplies $\theta$. Its squared radius follows the exponential form above, which supplies $R$. The result has mean zero and variance one in each coordinate.

Four million values measured:

Box-Muller, N=4000000
  mean=-0.000805, standard deviation=1.000039
  P(|Z|>4)=5.55e-05, analytic=6.334248e-05
  P(|Z|>5)=7.5e-07, analytic=5.733031e-07
  largest |Z| observed=5.2816

Only 222 samples exceeded magnitude 4, so the relative noise in that count is much larger than the noise in the mean. Three exceeded magnitude 5. A display with six decimal places can make a three-event estimate look more decisive than it is.

The transform also inherits the finite resolution of $U_1$. For the upper-53-bit mapping, the smallest positive value is $2^{-53}$. The largest possible radius from that mapping is:

$$\sqrt{-2\ln(2^{-53})} = 8.571674\ldots$$

This is not a universal limit of Box-Muller or of double. It is the limit produced by this particular 53-bit, evenly spaced uniform construction. A method that represents probabilities closer to zero can reach farther into the tail. A method that clamps zero to an arbitrary decimal such as 1e-300 changes the tail again.

The opening conversion bug matters here from the other side. Rounding the raw maximum to exactly 1 makes $R=0$, which adds slightly too much mass at the origin. Allowing exact zero makes $R$ infinite. A distribution sampler is not independent of the bit-to-floating conversion beneath it.

two thirds of the points were discarded

For a density proportional to $x^2$ on $[0,1]$, I sampled a point $(x,y)$ uniformly in the unit square and kept it when $y \le x^2$. The probability of accepting a proposed $x$ is therefore $x^2$, exactly the relative weight the target asks for.

The area below the curve predicts the acceptance rate:

$$\int_0^1 x^2\,dx = \frac{1}{3}$$

After normalization, the target density is $3x^2$. Its predicted mean is:

$$\int_0^1 x(3x^2)\,dx = 3\int_0^1x^3\,dx = \frac{3}{4}$$

The probe accepted five hundred thousand values:

rejection sample proportional to x^2, accepted=500000
  acceptance=0.33347, analytic=1/3
  mean=0.75021, analytic=3/4

Rejection sampling trades proposed values for a transformation that can handle shapes without a convenient inverse. The proposal must cover the target everywhere. If the supposed envelope falls below the target in one region, that region is silently underrepresented.

one item survived the stream

The next transformation does not change a numeric distribution. It chooses one item uniformly from a stream whose final length is unknown. Keep the first item. When item $i$ arrives, counting from one, replace the held item with probability $1/i$.

For item $j$ to survive a stream of length $n$, it must be selected when it arrives and then survive every later replacement:

$$\Pr(j\text{ survives}) = \frac{1}{j} \times \frac{j}{j+1} \times \frac{j+1}{j+2} \times \cdots \times \frac{n-1}{n} = \frac{1}{n}$$

Every numerator cancels the previous denominator. The result is independent of $j$. A run over two hundred thousand streams of length ten produced:

one-item reservoir, 200000 streams of length 10
  rates: 0.1005 0.0985 0.1006 0.1003 0.1009
         0.0994 0.0999 0.0999 0.1003 0.0998

The ten rates wobble around 0.1. The proof establishes uniformity if each bounded choice is itself uniform. The measurement checks the implementation. Neither substitutes for the other.

Token sampling in the forward-pass investigation uses another discrete transformation. A softmax supplies one probability per token, and a cumulative bounded draw selects an index. Temperature changes those probabilities before selection. It does not create randomness, increase entropy in the generator state, or repair a biased bounded-integer routine.

one uniform value selected a token

Take three possible tokens with probabilities 0.5, 0.3, and 0.2. Their cumulative endpoints divide $[0,1)$ into three intervals:

[0.0, 0.5) -> token A
[0.5, 0.8) -> token B
[0.8, 1.0) -> token C

A uniform draw lands in each interval with probability equal to its width. The selection code can subtract probabilities until the draw crosses zero, or build cumulative sums and use a binary search. Both implement the same partition when arithmetic is exact.

The endpoint bug now has a visible consequence. A draw rounded to 1.0 lies outside every half-open interval. Code that falls through to the last token gives C a tiny extra point mass. Code that uses an array search can return an index one past the token array. The upper-53-bit conversion avoids that endpoint before sampling begins.

The probabilities usually start as model logits, arbitrary real scores $z_i$. Temperature $T>0$ changes their relative scale before normalization:

$$p_i(T) = \frac{\exp(z_i/T)} {\sum_j \exp(z_j/T)}$$

Lower $T$ enlarges score differences and concentrates probability on the largest logits. Higher $T$ compresses differences and moves the distribution toward uniform over the available tokens. The random generator underneath has not changed. The intervals built on top of its uniform draw changed width.

Numerical stability from the forward-pass experiment still applies. Subtracting the largest logit before exponentiation leaves every ratio unchanged:

$$\frac{\exp((z_i-m)/T)} {\sum_j\exp((z_j-m)/T)} = \frac{\exp(z_i/T)} {\sum_j\exp(z_j/T)}$$

because the common factor $\exp(-m/T)$ cancels. Without that shift, a large positive logit can overflow even though the final probability should be finite.

Top-k and top-p sampling change support before the cumulative intervals are built. Top-k keeps a fixed number of candidates. Top-p sorts by probability and keeps the smallest prefix whose cumulative mass reaches a threshold. The remaining probabilities are renormalized. A token outside the retained set has probability exactly zero under that policy, regardless of how long the generator runs.

This makes generation replay more demanding than storing a seed. Exact token reproduction also needs model weights, tokenizer, logits arithmetic, temperature, truncation policy, tie behavior, batch scheduling, and generator consumption order. If one request leaves a continuous batch, later requests may consume different positions from a shared GPU generator even with the same process seed.

Per-request generator state or a counter mapping can isolate request streams. That improves replay and prevents one customer’s output length from shifting another customer’s random choices. It does not guarantee the same logits across kernel versions or floating-point reduction orders. Deterministic random draws cannot compensate for nondeterministic probabilities.

There is a security boundary here too. Sampling generated text usually needs statistical quality and controlled replay, not secret unpredictability. Generating an authentication token needs a system CSPRNG and enough entropy. Reusing the model’s fast sampling generator for both jobs joins two contracts that fail in opposite directions.

the missing constant was larger than one

The quarter circle inside the unit square has area $\pi/4$. Draw $x$ and $y$ uniformly from $[0,1]$, record 1 when $x^2+y^2 \le 1$, and record 0 otherwise. If the fraction of ones is $\hat p$, then $4\hat p$ estimates $\pi$.

One C run used srand(12345) and the measured Apple rand() stream:

N=1000       pi~3.12000  err=0.02159   (1/sqrt(N)=0.03162)
N=100000     pi~3.13876  err=0.00283   (1/sqrt(N)=0.00316)
N=10000000   pi~3.14169  err=0.00010   (1/sqrt(N)=0.00032)

The individual error did not fall monotonically. It happened to be smaller at every displayed sample count, but a later run can move away from $\pi$ when one more point is added. The useful prediction concerns the distribution over repeated independent runs.

Let $I_i$ be the 0-or-1 result of dart $i$. Its probability of being 1 is $p=\pi/4$. For any 0-or-1 variable:

$$\operatorname{Var}(I_i) = p(1-p)$$

The sample fraction is:

$$\hat p = \frac{I_1 + I_2 + \cdots + I_N}{N}$$

If the darts are independent, variances add. Dividing by $N$ also divides standard deviation by $N$, so:

$$\operatorname{Var}(\hat p) = \frac{p(1-p)}{N}$$

Multiplying the estimate by four multiplies its standard deviation by four:

$$\operatorname{SD}(4\hat p) = \frac{4\sqrt{p(1-p)}}{\sqrt N} = \frac{1.642183\ldots}{\sqrt N}$$

The $1/\sqrt N$ shape came from independent variance addition. The constant came from the event being estimated. The original table compared its absolute error with $1/\sqrt N$ and silently used 1 instead of 1.642.

I repeated complete estimates with independently advanced streams and measured root-mean-square error:

derived RMS constant: 4*sqrt(p*(1-p))=1.642183
N=  1000, trials=2000, RMS=0.052143, RMS*sqrt(N)=1.64892
N= 10000, trials= 500, RMS=0.016126, RMS*sqrt(N)=1.61260
N=100000, trials= 100, RMS=0.005319, RMS*sqrt(N)=1.68209

The recovered constants fluctuate around 1.642. The largest-$N$ row has the fewest repeated trials because each trial costs more. Its recovered constant is consequently not the closest. The derivation predicts the population quantity; finite validation runs still have sampling error.

Reducing typical error by a factor of ten requires one hundred times as many independent samples:

$$\frac{C}{\sqrt{100N}} = \frac{1}{10}\frac{C}{\sqrt N}$$

This square-root rate applies to the mean of independent observations with finite variance. It is not a law that rescues every estimator. Correlated samples change the effective count. Infinite-variance distributions do not provide the same standard error. Bias does not disappear merely because variance shrinks.

eight workers still threw one set of darts

The forked rand() result breaks the independence step in the derivation. Suppose eight child processes each draw $N$ identical darts because they inherited the same state. The combined program reports $8N$ samples, but every group of eight contains copies of one observation.

If the implementation treats all $8N$ values as independent, it reports standard error proportional to:

$$\frac{1}{\sqrt{8N}}$$

The actual information came from $N$ distinct darts, so the relevant scale is:

$$\frac{1}{\sqrt N}$$

The reported interval is too narrow by $\sqrt 8$, about 2.83. More processes made the computation more expensive and the uncertainty statement less honest.

Perfect duplication is an extreme case of correlation. For a stationary sequence with lag correlations $\rho_k$, one common approximation writes the variance of the mean using an integrated autocorrelation time:

$$\operatorname{Var}(\bar X) \approx \frac{\sigma^2}{N} \left(1 + 2\sum_{k=1}^{\infty}\rho_k\right)$$

The parenthesized factor says how many nominal samples one effectively spends per independent sample. If adjacent observations have positive correlation, it exceeds one. A generator can pass a one-dimensional frequency test while an application’s state update introduces correlation between measured samples.

This is also why rerunning a deterministic simulation with the same seed does not measure uncertainty. It verifies reproducibility. To examine sampling variation, each replicate needs a documented distinct stream. To debug one surprising replicate, the seed and stream assignment need to remain recoverable.

the rare event consumed almost every draw

I estimated the probability that a standard normal exceeds 3 by direct counting. The analytic value is:

$$p = \Pr(Z>3) = \frac{1}{2}\operatorname{erfc}\left(\frac{3}{\sqrt2}\right) \approx 0.0013499$$

The measured counts were:

normal >3, N=  10000, hits=  17, estimate=0.0017000,
  relative error=0.259, 1/sqrt(Np)=0.272
normal >3, N=1000000, hits=1359, estimate=0.0013590,
  relative error=0.007, 1/sqrt(Np)=0.027

For a rare Bernoulli event, $\sqrt{p(1-p)/N}$ is the absolute standard deviation of the estimated probability. Divide by $p$ to express it relative to the small answer:

$$\frac{\sqrt{p(1-p)/N}}{p} = \sqrt{\frac{1-p}{Np}} \approx \frac{1}{\sqrt{Np}}$$

The useful count is $Np$, the expected number of hits. Ten thousand normal draws produce only about 13.5 expected hits beyond 3. Almost every random number is spent confirming that an ordinary value is ordinary.

For a target relative standard error $\epsilon$, rearranging the approximation gives:

$$N \approx \frac{1}{\epsilon^2p}$$

At $p=0.0013499$ and $\epsilon=0.1$, that is about 74,080 draws. At a six-sigma one-sided probability near $10^{-9}$, the same target would need around $10^{11}$ direct draws. A generator with a long period does not make that estimator efficient. Techniques such as importance sampling change where samples are spent, then correct the weights so the changed sampling distribution does not bias the answer.

the faster mean was not enough

Every benchmark in the timing investigation estimates behavior from a finite sample. The estimator may be a mean, median, quantile, throughput ratio, or tail probability. A smaller mean in one run set is an observation. Whether it supports a general performance claim depends on the sampling design and the variation.

Start with synthetic values whose true generating process is known:

A: mean=98.0ms std=13.5   B: mean=87.4ms std=13.6   (B looks 11% faster)

There are twenty independent samples in each group. If sample standard deviations are $s_A$ and $s_B$, the estimated standard error of the difference between independent means is:

$$\operatorname{SE}(\bar A-\bar B) = \sqrt{\frac{s_A^2}{n_A}+\frac{s_B^2}{n_B}}$$

Putting 13.5 ms and 13.6 ms into that expression gives about 4.29 ms. It is larger than either individual mean’s standard error because both estimated endpoints move.

standard error of A mean: 13.5/sqrt(20) = 3.02 ms
standard error of B mean: 13.6/sqrt(20) = 3.04 ms
standard error of difference: sqrt(3.02^2 + 3.04^2) = 4.29 ms

Using 1.96 as a multiplier assumes a normal reference distribution with known variance or a sufficiently accurate large-sample approximation. The variances here were estimated from twenty observations. Welch’s interval instead uses a Student t critical value and an approximate number of degrees of freedom:

$$\nu = \frac{\left(s_A^2/n_A+s_B^2/n_B\right)^2} {\frac{(s_A^2/n_A)^2}{n_A-1} +\frac{(s_B^2/n_B)^2}{n_B-1}}$$

For equal sample sizes and nearly equal variances, $\nu$ is close to 38, not 19. The earlier draft called a deliberately conservative 19-degree multiplier “statistically correct.” It was not the Welch calculation for these two groups.

I tested the coverage claim directly because the synthetic generator fixes the true gap at 10.6 ms. Each of one hundred thousand experiments drew two new groups of twenty values and built both intervals:

coverage over 100000 experiments, nominal 95%
  normal critical value: 94.309%
  Welch t critical value: 95.111%

The normal interval missed the true value too often. Welch’s procedure landed near its nominal long-run rate for this normally distributed synthetic case.

A 95 percent confidence interval does not mean there is a 95 percent probability that this fixed true gap lies in the particular interval already computed. In the frequentist construction, the unknown gap is fixed and the interval is random before sampling. Under the model assumptions, 95 percent of intervals produced by repeated experiments cover the fixed gap.

Excluding zero is evidence against the zero-difference model under those assumptions. It does not prove the optimization, identify the cause, establish practical importance, or guarantee the effect on another machine. A narrow interval around a benchmark corrupted by warmup or changing input shapes is a precise answer to the wrong experiment.

the order became part of the sample

Timing all A runs and then all B runs makes method identity indistinguishable from time. A temperature rise, frequency policy change, background job, cache transition, or battery event can become an apparent method difference.

I measured two loops over ten million doubles. The first uses one accumulator. Each addition depends on the previous addition. The second uses four accumulators and combines them at the end, exposing independent addition chains. Before measuring, the prediction was that four chains would reduce dependency latency and win if memory delivery did not dominate.

Each of forty pairs ran both methods. A fixed control generator randomized which method ran first inside each pair. Four unmeasured warmups preceded collection. The program was compiled with Apple clang 15 at -O3. It retained the checksum through a volatile sink.

Apple clang -O3, 40 randomized pairs, 10000000 doubles
  naive mean=5.7533 ms, sd=0.2130
  four-accumulator mean=1.5824 ms, sd=0.0725
  paired difference=4.1709 ms, 95% t CI=[4.1224, 4.2194]
  checksums remain live: 4995038.282457879

This run supports the predicted direction on the measured machine. It does not prove the dependency chain is the sole cause because no instruction-level counter or generated assembly is part of this probe. The performance mechanism remains an inference from source structure and the large timing change.

The final verifier reran those forty pairs and found a 4.1596 ms mean paired difference with interval [4.1196, 4.1995]. The two intervals overlap heavily and preserve the predicted direction. They are repeated executions, not a cross-machine validation.

Pairing matters when both measurements in a pair share environmental noise. Let $D_i=A_i-B_i$. The paired estimator is simply the mean of the $D_i$ values, with standard error:

$$\operatorname{SE}(\bar D) = \frac{s_D}{\sqrt N}$$

The common disturbance cancels when the difference is formed. To isolate that effect, I generated thirty synthetic pairs with a true difference of 4 ms. Each pair shared a large random disturbance and each method also received smaller independent noise:

synthetic randomized pairs with true A-B difference 4.0 ms
  observed difference=4.0939 ms
  unpaired Welch 95% CI=[-0.0628, 8.2505], df=58.00
  paired t 95% CI=[3.3977, 4.7900]

The unpaired calculation discarded the pairing and included zero. The paired interval retained the fact that measurements taken together moved together and became far narrower.

Pairing is not automatically valid. If running A changes the immediately following B through cache state or thermal effects, then the pair contains an order interaction. Randomizing order spreads that interaction instead of always favoring one method. A stronger design can use balanced AB and BA pairs and analyze the order effect explicitly.

the bootstrap resampled differences

The t interval uses a reference distribution for a standardized mean. A bootstrap asks a more mechanical question: if the observed paired differences approximate the population of differences, what happens when samples are drawn from those observed differences with replacement?

For the thirty synthetic pairs, each bootstrap replicate drew thirty differences with replacement and computed their mean. Twenty thousand replicate means formed an empirical distribution. Its 2.5 and 97.5 percent quantiles produced:

paired percentile bootstrap 95% CI=[3.4184, 4.7301]

That is close to the paired t interval. Agreement in one constructed case is a useful control, not evidence that both procedures always work. The basic percentile bootstrap can have poor coverage for small, biased, or highly skewed samples. Dependent time series cannot be resampled as independent points without destroying their correlation. A block bootstrap or a model of the dependence would be needed.

Bootstrap resampling also cannot repair a biased measurement design. If every B run happened after thermal throttling, resampling those B values reproduces the confounding twenty thousand times.

The random seed used by the bootstrap is operational state, just like the xoshiro seed. Keeping it makes a reported interval exactly reproducible. Repeating the bootstrap with several seeds checks whether the finite resample count itself contributes visible Monte Carlo noise.

twenty fair tests found a winner

Suppose twenty implementations are actually equal and each comparison has a 5 percent false-positive rate. Even if tests were independent, the probability that none produces a false difference would be:

$$0.95^{20} \approx 0.3585$$

The probability of at least one false difference is:

$$1 - 0.95^{20} \approx 0.6415$$

I simulated thirty thousand families. Every family compared twenty pairs of equal normal distributions with Welch intervals:

twenty equal variants, each tested at nominal 5%
  at least one false difference: 63.67%
  independent-test prediction: 64.15%

Selecting only the smallest measured time or the interval that excluded zero turns exploration into an unreported multiple-comparison procedure. Autotuners do this intentionally. Their winner needs evaluation on fresh measurements that were not used for selection.

A negative control passes identical work through both measurement labels. It should show no systematic difference. If it finds a “winner,” the measurement program, ordering, timer, or analysis can create effects without an optimization. Negative controls do not prove the main comparison is clean. They expose certain ways it is not.

the best run selected favorable noise

Taking the minimum of repeated timings can approximate an uncontended lower envelope when interruptions only add delay. It is also a selection operation. More repetitions create more chances for an unusually favorable measurement, so a best-of-100 result is not directly comparable with a best-of-10 result.

The minimum discards how often the system attains that performance. That may be appropriate for isolating a code path from scheduler interruption. It is inappropriate for request latency experienced by users, where the upper tail is part of the behavior being purchased.

Mean, median, minimum, and quantiles answer different questions:

• the mean contributes to total capacity and Little’s Law calculations when the sampling design is representative;
• the median describes a typical ordered observation but can hide a bimodal mode switch;
• the minimum approximates a lower envelope under strong assumptions about additive noise;
• the 99th percentile concerns the slow tail and needs far more than one hundred observations to estimate with useful resolution.

Reporting all of them does not rescue an uncontrolled workload. Warmup policy, cache state, inputs, compiler flags, process placement, background load, method order, failures, and discarded samples still determine what population those statistics describe.

Randomization does one specific job. It prevents a systematic nuisance trend from always aligning with one method. Repetition estimates variation under the chosen conditions. Blocking or pairing removes shared variation. Blinding can prevent an operator from stopping when the favored result appears. None of these changes a bad clock, a dead-code-eliminated loop, or a workload unrelated to production.

The seed used to randomize benchmark order belongs in the artifact. So does the realized order. If the result changes sharply under another valid order seed, order sensitivity is part of the finding rather than noise to average away.

one hundred requests could not locate the tail

A 99th-percentile latency is the value below which 99 percent of the target request population falls. With one hundred measured requests, only one observation is expected above that threshold. Depending on the quantile convention, the reported p99 becomes the largest or an interpolation between the largest few observations.

The rare-event calculation applies again. Let an exceedance beyond the true p99 be a 0-or-1 event with probability $q=0.01$. Its expected count is $Nq$. The relative standard error of the measured exceedance rate is approximately:

$$\sqrt{\frac{1-q}{Nq}}$$

At $N=100$, this is about 0.995. The tail count has relative uncertainty near 100 percent. To reduce that relative standard error to 10 percent:

$$N \approx \frac{1-q}{q(0.1)^2} = \frac{0.99}{0.01\times0.01} = 9900$$

That estimate still assumes independent representative requests and a stationary latency distribution. Autocorrelation from garbage collection, thermal throttling, queue buildup, or periodic background work reduces effective sample size.

The request generator can bias the population before the quantile is computed. A closed-loop client sends a request, waits for its response, and then sends the next. When the service slows, the client sends less load. The slowdown therefore suppresses the offered traffic that would have revealed a larger queue.

An open-loop generator schedules arrivals independently of completions. If it cannot issue an arrival on time and records latency only from the delayed actual send, it omits the waiting time caused by its own inability to keep up. This is coordinated omission. The histogram can improve as the system overloads because the missing requests never enter its clock.

Random arrival times are consequently part of the workload model, not a cosmetic source of jitter. The distribution, rate, burst structure, dropped arrivals, and timestamp origin determine which queue the measurement includes. A precise p99 from the wrong arrival process answers a different production question.

Bootstrap resampling of one hundred recorded latencies cannot invent tail events that were never observed. More resamples make the bootstrap distribution smoother while its information remains limited by the original hundred requests.

the grid looked in the wrong direction

The low-bit checks missed a known linear weakness. I then tried to see the Park-Miller lattice by putting adjacent normalized outputs into a 200 by 200 grid. Four million pairs produced:

2D pair test, 200x200 grid, N=4000000: dof=39999
   MINSTD : chi=40347  z=(chi-dof)/sqrt(2dof)=+1.23
   xoshiro: chi=39895  z=-0.37   (|z|>4 = clearly non-uniform)

For a chi-square statistic with many cells and adequate expected counts, subtracting its degrees of freedom and dividing by the approximate standard deviation $\sqrt{2\nu}$ gives a rough standard-normal scale. Neither score was unusual. The home-built test failed to distinguish the streams.

The recurrence says exactly where to look. Before the modulus:

$$x_{n+1} = Ax_n - kM$$

Here $A=16807$, $M=2^{31}-1$, and $k$ is the integer quotient discarded by the remainder operation. Divide both coordinates by $M$:

$$\frac{x_{n+1}}{M} = A\frac{x_n}{M} - k$$

If $u_n=x_n/M$, then every adjacent point satisfies:

$$u_{n+1} = 16807u_n-k$$

Each possible $k$ names one parallel line. The perpendicular spacing between neighboring lines is:

$$\frac{1}{\sqrt{16807^2+1}} \approx 0.0000595$$

A grid cell was $1/200=0.005$ wide, about 84 times that spacing. Many lines crossed each cell. Counting cells erased the structure before chi-square saw it.

I replaced the generic grid with a test aligned to the exact modular relation. For each adjacent pair, it checked whether:

$$(16807x_n-x_{n+1}) \bmod M = 0$$

One million pairs gave:

recurrence-aligned test, 1000000 adjacent pairs
  MINSTD: 1000000 exact modular matches,
          16807 observed parallel-line indices
  xoshiro values reduced modulo M: 0 matches,
          independent-uniform reference=0.000466

The result does not say xoshiro outputs have no structure. Every finite-state deterministic generator does. It says xoshiro values reduced modulo $M$ do not obey this Park-Miller recurrence in the observed sample. The test found the defect because the implementation supplied a precise alternative to uniformity.

a passed test left the question open

The recurrence-aligned test would be unfair as a general contest. It was designed after reading one generator’s recurrence. A separate test set is needed when the goal is to evaluate a new generator without tailoring every test to its source.

Statistical batteries collect tests for frequency, runs, linear complexity, matrix rank, collisions, serial relationships, Hamming weights, and other distinguishers. Blackman and Vigna report xoshiro256** results from TestU01 and additional Hamming-weight analysis on the authors’ site. Those are externally reported results. I did not install or run TestU01 or PractRand in this environment.

A test begins with a null model, usually that the relevant statistic behaves as it would for independent uniform bits. Its p-value measures how extreme the observed statistic is under that model. It is not the probability that the generator is random. A large p-value can mean the generator matches the model for that statistic, or that the test lacks power at the available sample size.

Thousands of tests also create extreme p-values by chance. A battery must account for its testing policy, expected distribution of p-values, repeated seeds, bit order, output transformations, and systematic rather than isolated failures. Searching configurations until one p-value looks alarming repeats the twenty-benchmark-winners problem.

Passing a battery does not certify a generator for cryptographic use. Xoshiro256** can pass strong statistical tests and still surrender its state through the invertible scrambler and linear transition demonstrated above. Conversely, a cryptographic generator can fail an application because the program maps its output with biased modulo or duplicates state across an unsupported process boundary.

The useful question is narrower than “is this random?”:

• can a stated observer predict future values;
• does the distribution consumed by the application match its model;
• are streams independent enough for the estimator;
• can the exact run be reproduced when reproduction is wanted;
• is the state kept secret and refreshed when secrecy is wanted;
• does the measurement design distinguish the effect from its nuisance variables.

Each question needs different evidence. One histogram cannot answer all six.

the first value had stopped being mysterious

The opening 16807 is now predictable before running the executable. The C contract supplies an implicit seed of 1. Apple’s measured recurrence multiplies that state by 16,807 and takes a modulus larger than the product, so the first result remains 16,807. The next product wraps only when it crosses the modulus. Restarting the process reconstructs the static state and repeats the calculation.

For a test matrix, that is a useful mechanism. I would record the algorithm, implementation version, complete seed-to-state procedure, stream assignment, raw-to-value mapping, and seed. A future run then has a chance to reproduce the same values even if the surrounding library changes.

For a secret, every one of those replay properties is dangerous if an observer can obtain the replay material. The system CSPRNG exists to maintain secret state, mix entropy into it, handle lifecycle boundaries such as fork, and expose bounded or bulk interfaces without asking application code to invent those mechanisms.

For a benchmark, randomness is neither decoration nor permission to average. It chooses order, constructs confidence procedures, powers negative controls, and creates the repeated samples whose independence the mathematics assumes. Its seed belongs in the evidence. Its state must not be accidentally shared. Its selected winner must face fresh data.

The five integers did not become genuinely random after any of these changes. They became accounted for. State explained the repeat. Source explained the state transition. Algebra predicted the sequence. Fork exposed the process copy. Mapping explained the bias. Variance explained how much a sample mean could move. Once every hidden choice had a name, 16807 was no longer surprising.