Arata: Fast Quality Pseudo-Random Number Generator Based on ARX

• Fast - Much faster than the common PCG and slightly faster than the xoshiro256 series.
• Quality - Passed terabytes of PractRand easily, as well as BigCrush and SmokeRand tests.
• Minimum Period - 2^64 (or 2^32 for the 32-bit variant), avoids excessively short periods.
• Domain Separation - For RNG instances with different domain IDs, the output sequences are guaranteed not to overlap. Suitable for large-scale parallel applications.
• ARX-Based - No multiplication and offers better portability in terms of efficiency, e.g., via SIMD implementation or porting to a GPU.
• Non-Crypto - Not intended for cryptographic purposes.

Provide C reference implementation and Rust crate

Variants

Capacity refers to how many bytes an RNG instance can produce before statistical flaws become detectable. This estimate is derived by extrapolating from the performance of variants with a smaller bit-width.

• We primarily recommend using the default Arata, as it delivers the best throughput and possesses ample capacity to handle any non-crypto applications.

• If 64-bit arithmetic is inefficient, switch to Arata32 or ArataQuad32 depending on the required capacity; the former should suffice for most applications. However, be mindful of the period issue: they only have 32-bit counters, which means that in the worst-case, they would cycle after only outputting 2^34 bytes (16 GiB), though this is very unlikely to happen.

  Solutions

  • Directly adding the carry from w + φ to domain_id, this creates a 64-bit counter; however, this approach makes the domain separation unusable.

  • A more balanced approach is to use only the lower 16 bits of the domain_id as counter extension, as shown in the code:

    let (w_, carry) = w.overflowing_add(0x9e3779b9);
     w = w_;
    domain_id = 0x0000FFFF & domain_id.wrapping_add(carry as u32)
              | 0xFFFF0000 & domain_id;
    

    This creates a 48-bit counter and a total of 65,536 non-overlapping streams. Each stream guarantees an output of at least 2^50 bytes (1,024 TiB) without cycling, truly sufficient for the vast majority of applications.

    However, neither the C reference implementation nor this Rust crate provides such a solution yet. You may need to write the implementation yourself.

• As for ArataQuad, that's for the paranoid who crave endless randomness.

Benchmarks

Single call

NAME                LATENCY         NOTE
Arata               0.877 ns/u64    without domain separation
Arata.domain        0.962 ns/u64    with domain separation
ArataQuad           1.039 ns/u64    without domain separation
ArataQuad.domain    1.054 ns/u64    with domain separation
RomuTrio            1.104 ns/u64
Pcg64               1.675 ns/u64
Pcg32               1.840 ns/u64
Xoshiro256++        1.072 ns/u64
Xoshiro256**        1.133 ns/u64
Xoroshiro128++      1.054 ns/u64
Xoroshiro128**      1.055 ns/u64
ChaCha8             1.865 ns/u64    SSE2 and AVX are used

Fill bytes

NAME                THROUGHPUT      NOTE
Arata               12.338 GiB/s    without domain separation
Arata.domain        10.449 GiB/s    with domain separation
ArataQuad           11.263 GiB/s    without domain separation
ArataQuad.domain     9.691 GiB/s    with domain separation
RomuTrio            12.345 GiB/s
Pcg64                5.822 GiB/s
Pcg32                4.261 GiB/s
Xoshiro256++         9.710 GiB/s
Xoshiro256**        10.372 GiB/s
Xoroshiro128++       7.046 GiB/s
Xoroshiro128**       7.054 GiB/s
ChaCha8              6.439 GiB/s    SSE2 and AVX are used

Tested on AMD Ryzen 7 5700G desktop with criterion.rs.

Statistical tests

PractRand

Note: an older version of the rotation constants was used during the testing of ArataQuad16.

Based on the performance of variants with smaller bit-widths, and taking into account the masking effect of the Weyl counter on defects in narrow-width variants, we conservatively set $k = 0.6875$ and estimate the capacity using the following formula: $\text{Est. Capacity} = 2^{k \cdot \text{bits}}$ where $\text{bits}$ is the footprint bits.

This estimation is feasible because the number of state bits in LCG and MCG generators has been observed to have a linear relationship with the exponent of the failure point. See M.E. O'Neill, Does It Beat the Minimal Standard?

TestU01

SmokeRand

Note: an older version of the rotation constants was used during the coll64dec test.

The design

Fig. 1: The next-state function of Arata.

Fig. 2: The next-state function of ArataQuad.

Both output a + x; this is not shown in the figures.

• w is the Weyl counter, used to provide the minimum period for the generator.
• a and b alternately employ modular addition and XOR operations, providing critical nonlinear confusion.
• A linear diffusion step is applied to x to propagate the confusion results throughout the entire state as rapidly as possible.
• id is the optional domain ID; setting it to a non-zero value enables domain separation. For RNG instances with different domain IDs, the output sequences are guaranteed not to overlap, regardless of their initial states.
• $\varphi$ is the fractional part of the golden ratio, 0x9e3779b97f4a7c15.

As can be seen, the Quad variant simply inserts an extra state-word into the standard variant. Honestly, I feel that the current linear diffusion step does not "make efficient use" of this additional state-word. Even so, the extra capacity this provides to ArataQuad32 should enable it to handle all non-crypto applications.

I guess taking Arata's output once every 32 calls to next (discarding the first 31 entirely) should yield a keystream with true cryptographic strength :p

Disclaimer: Do not actually use this method to encrypt data. A proper encryption algorithm requires rigorous theoretical analysis, not just statistical test; furthermore, the implementation of encryption and decryption must be correct—for instance, allowing the decryption of arbitrary data (not just valid ciphertext) is a massive security risk in itself. Do not attempt to design and implement your own encryption unless you truly know what you are doing. Be responsible for your own information security and that of others.

The mix function

A mix function is required to eliminate inter-stream correlation between RNGs with similar initial states. We define it explicitly to ensure portability. Refer to the C reference implementation and the test for details.

Rationale

Selection of incremental constant

Prefixes of the Weyl sequence generated by $\frac 1 \varphi$ are known to be "well spread out", and there is no reason not to choose them. See Donald E. Knuth, Sorting and Searching (second edition), Volume 3 of The Art of Computer Programming (Addison-Wesley, Reading, Massachusetts, 1998), Exercise 6.4-9.

Selection of rotation constants

The initial rotation constants were selected empirically: (7, 12), (13, 24), and (25, 48) for 16-, 32-, and 64-bit variants, respectively. However, we later noticed that these sets of rotation constants did not perform optimally during the linear diffusion step. We evaluated their performance by tracking the trajectories of their Hamming weights:

$$ F(x) = x \oplus (x \lll {\tt r1}) \oplus (x \lll {\tt r2}) $$ $$ x_{n+1} = F(x_n), \ x_0 = 1 $$ $$ w_n = \text{HammingWeight}(x_n) $$

In an ideal diffusion process, the influence of a single bit should rapidly propagate to multiple bits, rather than collapsing into a state where it no longer affects other bits after only a few iterations. Observing the trajectories of the old constants reveals that they collapse back to a Hamming weight of 1 after exactly four iterations; at this point, the operation is equivalent to a single rotation of the input with no diffusion whatsoever.

We subsequently wrote a program to search for theoretically superior rotation constants based on the following criteria:

• Maximizing the number of iterations required for the Hamming weight to collapse back to 1 for the first time;
• Maximizing the total sum of Hamming weights over the complete trajectory (a total of ${\tt W} - 1$ iterations);
• Minimizing Hamming weight fluctuations during the first $\frac {\tt W} 2$ iterations of the trajectory.

And we applied the following filtering conditions:

• $\frac 1 8 {\tt W} < {\tt r1} < \frac 1 2 {\tt W}$
• ${\tt r1} < {\tt r2} < \frac 7 8 {\tt W}$
• ${\tt r1} - {\tt r2} \equiv 1 \pmod 2$
• $\frac 3 8 {\tt W} < {\tt r2} - {\tt r1} < \frac 1 2 {\tt W}$

where r1 and r2 are the rotation constants to be searched, and W is the bit-width.

Multiple sets of constants meeting the criteria were found, and one set was selected based on personal preference. These are the final rotation constants selected and used in the released version: (6, 13), (10, 25), and (17, 42) for 16-, 32-, and 64-bit variants, respectively.

8-bit, r1=2, r2=5:
1 ->  3 ->  3 ->  5 ->  1 ->  3 ->  3 ->  5 ->  1

--- old ---
16-bit, r1=7, r2=12:
1 ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1
  ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1

32-bit, r1=13, r2=24:
1 ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1
  ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1
  ... 2 same lines.

64-bit, r1=25, r2=48:
1 ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1
  ->  3 ->  3 ->  9 ->  1 ->  3 ->  3 ->  9 ->  1
  ... 6 same lines.

--- new ---
16-bit, r1=6, r2=13:
1 ->  3 ->  3 ->  7 ->  3 ->  9 ->  5 ->  9 ->  1
  ->  3 ->  3 ->  7 ->  3 ->  9 ->  5 ->  9 ->  1

32-bit, r1=10, r2=25:
1 ->  3 ->  3 ->  9 ->  3 ->  9 ->  9 -> 19 ->  3
  ->  9 ->  9 -> 19 ->  5 -> 15 ->  7 -> 17 ->  1
  ->  3 ->  3 ->  9 ->  3 ->  9 ->  9 -> 19 ->  3
  ->  9 ->  9 -> 19 ->  5 -> 15 ->  7 -> 17 ->  1

64-bit, r1=17, r2=42:
1 ->  3 ->  3 ->  9 ->  3 ->  9 ->  9 -> 25 ->  3
  ->  9 ->  9 -> 19 ->  9 -> 27 -> 19 -> 35 ->  3
  ->  9 ->  9 -> 27 ->  9 -> 27 -> 19 -> 31 ->  5
  -> 15 -> 15 -> 37 ->  7 -> 21 -> 17 -> 37 ->  1
  ->  3 ->  3 ->  9 ->  3 ->  9 ->  9 -> 25 ->  3
  ->  9 ->  9 -> 19 ->  9 -> 27 -> 19 -> 35 ->  3
  ->  9 ->  9 -> 27 ->  9 -> 27 -> 19 -> 31 ->  5
  -> 15 -> 15 -> 37 ->  7 -> 21 -> 17 -> 37 ->  1

Is there further evidence that the new constants are superior to the old ones? In the scaled-down test for Arata16:

• With the old constants, the failure occurred at DC6-9x1Bytes-1, p ≈ 1e-21;
• With the new constants, the failure occurred at [Low4/16]FPF/16:all, p ≈ 1e-16.

The former represents a fundamental structural issue, whereas the latter relates to the distribution of the lower bits, and the reported p-value is less extreme. This indicates that the new rotation constants facilitate more effective diffusion. This holds even though the latter reports more "suspicious" items—because PractRand includes multiple non-orthogonal variants of the FPF test, all of which are reported together.

License

The algorithm itself and its C implementation are released into the public domain under the CC0 1.0 license.

This Rust crate is dual-licensed under the MIT and Apache 2.0 licenses, at your option.