Skip to content

[FEA]: Multiply-with-carry engines for <cuda/random> #12287

Open

Description

@vigna

Is this a duplicate?

Area

libcu++

Is your feature request related to a problem? Please describe.

I would like to generate random numbers on GPUs with a generator that has a small state
(registers are scarce), high throughput, good statistical quality, and cheap independent streams.
The options available today in CUDA trade these against each other:

  • XORWOW, the default generator of cuRAND (also available in cuRANDDx), fails PractRand
    at its first checkpoint, after 8 MB of output: its lowest bit is a linear function of the
    state, plus the alternating lowest bit of the Weyl sequence.
  • Philox4x32-10 is counter-based, which is ideal when random access to the sequence matters,
    but it needs 20 32x32->64-bit multiplications for four outputs.
  • MRG32k3a has a long, well-studied history, but it is eight times slower.

XORWOW (step of curand(), Marsaglia's initial state; PractRand 0.96, expanded test set, extra
folding) fails more and more tests as the output grows, from 10 failed results at 8 MB to 29 at
1 GB:

First failure Test
8 MB BRank (binary rank)
8 MB BDayS2 (birthday spacings)
64 MB DC7
256 MB BCFN
1 GB NS3

Describe the solution you'd like

I propose to add to the non-standard header <cuda/random> two engines from my paper
“Scrambled Multiply-with-Carry Pseudorandom Number Generators”.

  • cuda::mwc64x32: 32-bit outputs, 64 bits of state, period about 2^63;
  • cuda::mwc128x32: 32-bit outputs, 128 bits of state, period about 2^127.

On GPUs they are the fastest generators I measured: on an H100, more than 50% faster than
cuRAND's XORWOW and three and a half times as fast as Philox4x32-10; they pass PractRand to 32 TB
of output, and they can jump ahead by any number of steps in logarithmic time, which is what
discard() needs to give each thread its own stream.

// x: digits (one for mwc64x32, three for mwc128x32), c: carry
uint32_t mwc64x32::operator()() {
    const uint32_t result = x ^ c;            // the scrambler X
    const uint64_t t = 0xfff145f8ULL * x + c;
    x = uint32_t(t);
    c = uint32_t(t >> 32);
    return result;
}

uint32_t mwc128x32::operator()() {
    const uint32_t result = z ^ c;
    const uint64_t t = 0xffcccc66ULL * x + c;
    x = y;
    y = z;
    z = uint32_t(t);
    c = uint32_t(t >> 32);
    return result;
}

The interface is that of the engines of <cuda/std/random>, such as cuda::std::philox4x32,
following the RandomNumberEngine requirements:
result_type (uint32_t), min(), max(), default_seed, constructors and seed(),
operator(), discard(), comparison and stream operators.

I can provide known-answer tests from the reference implementation and a pull request with
the two engines and their tests. I am also submitting a request for enhancement to cuRAND, which
is not open source; an engine in <cuda/random> would be useful anyway, as it composes with the
distributions of libcu++ and Thrust.

Performance

Throughput of the whole GPU in billions of 32-bit outputs per second (higher is better): each
thread keeps its state in registers and accumulates 2^24 outputs; one wave of threads sized by
the occupancy of each kernel; best over five runs for each block size in {128, 256, 512, 1024};
SM clock locked at its maximum (on the H100 SXM, it stays at 1830 MHz).
64-bit outputs count as two 32-bit outputs. Sources of the
benchmark: https://github.com/vigna/prng-speed (directory gpu/).

Generator V100 (CUDA 11.4) H100 (CUDA 13.2)
MWC64X/32 1929 4397
MWC128X/32 1828 4193
SFC32 1399 3376
XORWOW (cuRAND) 1391 2783
MWC128X (b = 2^64) 1171 2618
xoshiro128++ 1164 2278
PCG32 (pcg32_random_r()) 879 1719
Philox4x32-10 (cuRAND, curand4()) 613 1242
MRG32k3a (cuRAND) 233 563

Quality

PractRand 0.96, expanded test set, extra folding:

Generator Result
MWC64X/32 no anomalies at 32 TB
MWC128X/32 no anomalies at 32 TB
XORWOW fails at 8 MB (first checkpoint)

Links

Describe alternatives you've considered

No response

Additional context

No response

Activity

  1. github-actions commented on Oct 27, 2026

    @github-actions
    Contributor

    Hi @vigna!

    Thanks for submitting this issue - the CCCL team has been notified and we'll get back to you as soon as we can!
    In the mean time, feel free to add any relevant information to this issue.

  2. fbusato commented on Oct 27, 2026

    @fbusato
    Contributor

    @vigna, thanks for suggesting your algorithm in CCCL!!
    both speed and statistical tests look amazing. My understanding is that the paper is not peer-reviewed yet. We need to evaluate the claims before going ahead with the implementation

  3. vigna commented on Oct 27, 2026

    @vigna
    Author

    Yeah, that'll take a while. Note that the paper is more about a theoretical analysis: Marsaglia proposed MWCs in the '90s. The paper chooses good multipliers. The scrambler has appeared before, too (albeit floating around the net, not in a paper).

    You can run PractRand on your own. I have also logs.

    For some reason MWCs were underestimated, except in certain circles (e.g., people in videogames use the mother-of-all MWC from Marsaglia). They were used to generate immense periods nobody needs, whereas their strength is to produce very high quality random bits with just a multiplication and a sum, which makes them perfect for GPUs and modern 64-bit processors with fast multiplications.

  4. fbusato commented on Nov 6, 2026

    @fbusato
    Contributor

    We completed our analysis.
    The algorithm has excellent performance, small state, and use few registers. We could improve it further with some low-level optimizations. Our plan is to expose the engine under cuda::experimental as the paper has not yet been peer-reviewed.

    We have also some concerns about stream correlations. On the other hand, it is fine for several applications that prioritize speed (including LLM's applications) and that are not in cryptographic, simulation, or scientific computing contexts.

  5. vigna commented on Nov 7, 2026

    @vigna
    Author

    Note that MWC generators have been around for more than 30 years (GMWC for about 20), and their mathematics is well understood. What the paper provides is some more mathematical analysis to justify the choice of multipliers in a situation in which some scores are defective, and the X scrambler (which, in the case of MWC, degenerates to Thomas's scrambler). The scrambler can only be tested empirically, which we do at large scale using the coll-birth crate and PractRand—and that part is easily reproducible.

    Stream correlation happens if you use states whose associated integers (carry · bᵏ + digits) are related by a small multiplier—e.g., c = 0 and last digit 1, 2, 3, …, but also c = 1 and last digit 1, 2, 3, …. I was suggesting to fill the digits with SplitMix and set c = 1, but in the case of MWC64X/32 the state is so small that this strategy restricts the pool of initial states to 2^32, and in such a small pool two random states related by a small multiplier are not unlikely. I didn't realize it—you made me think about it.

    So now the suggestion is to fill the whole state, carry included, with SplitMix, mapping the carry to 1 ≤ c ≤ a − 2. In this way, MWC64X/32 has 2^63 possible starting points instead of 2^32. I applied the same change to all generators, even if for the others the advantage is marginal. Note also that with MWC64X/32 many randomly seeded streams can overlap (the period is about 2^63), so for many streams (e.g., one per GPU thread) the right choice is jumping, which is cheap.

    Thank you for the feedback!

  6. moved this from Todo to In Progress in CCCLon Nov 11, 2026
  7. moved this from In Progress to In Review in CCCLon Nov 11, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

Labels

needs triageIssues that require the team's attention

Type

No type

Projects

  • Status
    In Review

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions