If you don't need very high quality randomness, and close-to-uniform distribution is good enough, you can go really fast, especially on a modern CPU with efficient SIMD integer vectors like x86 with SSE2 or AVX2.
This is like @NominalAnimal's answer since we both had the same idea, but manually vectorized for x86. (And with with worse quality random numbers, but still probably good enough for a lot of use-cases.) This runs about 15 or 30 times faster than @Nominal's code, at ~13GB/s of ASCII output on a 2.5GHz Intel Haswell CPU with AVX2. That's still less than theoretical max main memory bandwidth (dual channel DDR3-1600 is about 25.6GB/s), but I was timing writing to /dev/null so it's actually just rewriting a buffer that stays hot in cache. Skylake should run this same code significantly faster than Haswell (see the bottom of this answer).
Assuming you actually bottleneck on I/O to disk or piping this somewhere, a fast implementation means your CPU doesn't even have to clock higher than idle. It uses much less total energy to produce the result. (Battery life / heat / global warming.)
This is so fast that you probably don't want to write it to disk. Just re-generate as-needed (from the same seed if you want the same data again). Even if you want to feed it to a multi-threaded process that can use all CPUs, running this to pipe the data to it will leave it hot in L3 cache (and L2 cache on the core that wrote it), and use so very little CPU time. (But note that piping adds a lot of overhead vs. writing to /dev/null. On a Skylake i7-6700k, piping to wc -c or another program that just reads + discards its input, it's about 8x slower than writing to /dev/null, and only uses 70% of a CPU. But that's still 4.0GB/s on a 3.9GHz CPU.
Re-generating it is faster than re-reading it even from a fast PCIe-connected SSD, but IDK if it's more power efficient (the vector-integer multiplier is kept pretty busy, and it's probably pretty power-hungry, along with other AVX2 256b vector ALUs). OTOH, I don't know how much CPU time reading from disk would take away from something that was maxing out all cores processing this input. I'd guess that a context-switch to re-generate in 128k chunks might be competitive with running filesystem / pagecache code and allocating pages to read data from disk. Of course, if it's already hot in the pagecache, it's just basically memcpy. OTOH, we already write about as fast as memcpy! (which has to split main memory bandwidth between reading and writing). (Also note that writing to memory that's not already hot in cache usually triggers a read-for-ownership to maintain cache coherency, which can be avoided with non-temporal stores, or with x86's rep movsb (optimized memcpy and memset in microcode, which avoids RFO, since Andy Glew's implementation of it in P6 (Pentium Pro))).
So far this is only a proof of concept, and the newline handling is only approximately correct. It's wrong around the ends of a power-of-2 buffer. With more development time. I'm confident I could find a more efficient way to insert newlines that's also exactly correct, with overhead at least as low as this (compared to outputting only spaces). I think this is something like 10 to 20%. I'm only interested in knowing how fast we could make this run, not in actually having a polished version of it, so I'll leave that part as an exercise for the reader, with comments describing some ideas.
On a Haswell i5 at its 2.5GHz max turbo, with DDR3-1600MHz RAM, timed producing 100GiB but scaled down. (Timed on cygwin64 on Win10 with gcc5.4 -O3 -march=native, omitted -funroll-loops since I was having a hard enough time getting decent timing runs on this borrowed laptop. Should have just booted Linux on a USB).
writing to /dev/null unless otherwise specified.
- James Hollis's: (not tested)
- Nominal's fwrite version: ~2.21s
- this (SSE2): ~0.142s (unscaled times = real=14.232s, user=13.999s, sys=0.187s).
- this (AVX-128): ~0.140s
- this (AVX2): ~0.073s (unscaled: real=0m7.291s, user=0m7.125s, sys=0m0.155s).
- this (AVX2) cygwin piping to
wc -c, with 128kiB buffer size: 0.32s with CPU at 2.38GHz (max dual-core turbo). (unscaled times: real=32.466s user=11.468s sys=41.092s, including both this and wc). Only half the data was actually copied, though, because my silly program assumes that write does the full buffer, even though that's not the case and cygwin write() only does 64k per call into a pipe.
So with SSE2 this is about 15 times faster than @Nominal Animal's scalar code. With AVX2, it's about 30 times faster. I didn't try a version of Nominal's code which just uses write() instead of fwrite(), but presumably for large buffers stdio mostly stays out of the way. If it is copying the data, that would account for a lot of slowdown.
Times to produce 1GB of data on a Core2Duo E6600 (Merom 2.4GHz, 32kiB private L1, 4MiB shared L2 caches), DDR2-533MHz in 64-bit Linux 4.2 (Ubuntu 15.10). Still using a 128kiB buffer size for write(), haven't explored that dimension.
writing to /dev/null unless otherwise specified.
- (SSE2) this with newline handling and 4 vectors of digits from each vector of random bytes: 0.183s (timed doing 100GiB in 18.3s, but similar results for 1GiB runs). 1.85 instructions per cycle.
- (SSE2) this, piping to
wc -c: 0.593s (unscaled: real=59.266s
user=20.148s sys=1m6.548s, including wc's CPU time). Same number of write() system calls as with cygwin, but actually piping all the data because Linux handles all 128k of a write() to a pipe.
- NominalAnimal's
fwrite() version (gcc5.2 -O3 -march=native), run with ./decdig 100 $((1024*1024*1024/200)) > /dev/null: 3.19s +/- 0.1%, with 1.40 instruction per cycle. -funroll-loops made maybe a tiny difference. clang-3.8 -O3 -march=native: 3.42s +/- 0.1%
- Nominal-
fwrite piping to wc -c: real=3.980s user=3.176s sys=2.080s
- James Hollis's line-at-a-time version (
clang++-3.8 -O3 -march=native): 22.885s +/- 0.07%, with 0.84 instructions per cycle. (g++5.2 was slightly slower: 22.98s). Writing only one line at a time probably hurt significantly.
- Stéphane Chazelas's
tr < /dev/urandom | ...: real=41.430s user=26.832s
sys=40.120s. tr was getting all of a CPU core to itself most of the time, spending nearly all its time in the kernel driver generating random bytes and copying them to a pipe. The other core on this dual core machine was running the rest of the pipeline.
time LC_ALL=C head -c512M </dev/urandom >/dev/null: i.e. just reading that much randomness with no piping: real=35.018s user=0.036s sys=34.940s.- Lưu Vĩnh Phúc's perl program (perl v5.20.2 from Ubuntu15.10):
LANG=en_CA.UTF-8: real=4m32.634s user=4m3.288s sys=0m29.364.
LC_ALL=C LANG=C: real=4m18.637s user=3m50.324s sys=0m29.356s. Still very slow.
- (SSE2) this with no newline handling, and either 3 or 4 vectors of digits from each vector of random bytes (almost exactly the same speed: the
dig3 = v%10 step is about break-even on this HW): 0.166s (1.82 instructions per cycle). This is basically the lower limit for what we can come close to with perfectly efficient newline handling.
- (SSE2) Old version of this with no newline handling, but only getting one digit per uint16_t element using
v%10, 0.222 seconds +/- 0.4%, 2.12 instructions per cycle. (Compiled with gcc5.2, -march=native -O3 -funroll-loops. Unroll loops does happen to help for this code on this hardware. Don't use it blindly, especially for large programs).
- (SSE2) Old version of this, writing to a file (on a RAID10f2 of 3 fast magnetic hard drives, not very optimized for writes): ~4 seconds. Could go faster by tweaking kernel I/O buffer settings to allow a lot more dirty data before write() blocks. "System" time is still ~1.0 seconds, much higher than "user" time. On this old system with slow DDR2-533 RAM, it takes ~4x longer for the kernel to memcpy the data into the pagecache and run XFS functions than it does for my loop to keep rewriting it in-place in a buffer that stays hot in cache.
How it's done
A fast PRNG is obviously essential. xorshift128+ can be vectorized, so you have two or four 64-bit generators in parallel, in elements of a SIMD vector. Each step produces a full vector of random bytes. (256b AVX2 implementation here with Intel intrinsics). I picked it over Nominal's choice of xorshift*, because 64-bit vector integer multiplication is only possible in SSE2/AVX2 with extended-precision techniques.
Given a vector of random bytes, we can chop up each 16-bit element into multiple decimal digits. We produce multiple vectors of 16-bit elements that are each one ASCII digit + ASCII space. We store that directly into our output buffer.
My original version just used x / 6554 to get one random digit from every uint16_t element of a vector. It's always between 0 and 9, inclusive. It's biased away from 9, because (2^16 -1 ) / 6554 is only 9.99923. (6554 = ceil((2^16-1)/10), which ensures that the quotient is always < 10.)
x/6554 can be computed with one multiply by a "magic" constant (the fixed-point reciprocal) and a right shift of the high-half result. This is the best case for division by a constant; some divisors take more operations, and signed division takes extra work. x % 10 has a similar bias and isn't as cheap to compute. (gcc's asm output is equivalent to x - 10*(x/10), i.e. an extra multiply and subtract on top of the division using a modular multiplicative inverse.) Also, the lowest bit of xorshift128+ is not as high quality, so dividing to take entropy from high bits is better (for quality as well as speed) than modulo to take entropy from low bits.
However, we can use more of the entropy in each uint16_t by looking at the low decimal digits, like @Nominal's digit() function. For maximum performance, I decided to take the low 3 decimal digits and x/6554, to save one PMULLW and PSUBW (and probably some MOVDQA) vs. the higher quality option of taking the 4 low decimal digits. x/6554 is slightly affected by the low 3 decimal digits, so there is some correlation between digits from the same element (8 or 16 digits separation in the ASCII output, depending on vector width).
I think gcc is dividing by 100 and by 1000, rather than a longer chain that successively divides by 10, so it's probably not significantly shortening the length of the non-loop-carried dependency chain that produces 4 results from each PRNG output. port0 (vector multiply and shift) is the bottleneck because of the modular multiplicative inverses, and the shifts in xorshift+, so it's definitely useful to save a vector-multiply.
xorshift+ is so fast that even using only ~3.3 bits of randomness from every 16 (i.e. 20% efficiency) is not a lot slower than chopping it up into multiple decimal digits. We only approximate the uniform distribution, because this answer is focused on speed as long as the quality isn't too bad.
Any kind of conditional behaviour that keeps a variable number of elements would take much more work. (But could maybe still be done somewhat efficiently using SIMD left-packing techniques. However, that gets less efficient for small element sizes; giant shuffle-mask lookup tables are not viable, and there's no AVX2 lane-crossing shuffle with smaller than 32-bit elements. A 128b PSHUFB version might still be able to generate a mask on the fly with BMI2 PEXT/PDEP, like you can for AVX2 with larger elements, but it's tricky because a 64-bit integer only holds 8 bytes. The godbolt link on that answer has some code that might work for higher element counts.)
If latency of the RNG is a bottleneck, we could go even faster by running two vectors of generators in parallel, alternating which one we use. The compiler can still easily keep everything in registers in an unrolled loop, and that lets the two dependency chains run in parallel.
In the current version, chopping up the output of the PRNG, we actually bottleneck on port 0 throughput, not PRNG latency, so there's no need for that.
The code: AVX2 version
Full version with more comments on the Godbolt compiler explorer.
Not very tidy, sorry I have to get to sleep and want to get this posted.
To get the SSE2 version, s/_mm256/_mm, s/256/128/, s/v16u/v8u/, and change vector_size(32) to 16. Also change the newline increment from 4*16 to 4*8. (Like I said, code is messy, and not well set up for compiling two versions. Didn't originally plan on making an AVX2 version, but then I really wanted to test on a Haswell CPU I had access to.)
#include <immintrin.h>
#include <unistd.h>
#include <stdint.h>
#include <stdio.h>
//#include <string.h>
// This would work equally fast 128b or 256b at a time (AVX2):
// https://stackoverflow.com/questions/24001930/avx-sse-version-of-xorshift128
struct rngstate256 {
__m256i state0;
__m256i state1;
};
static inline __m256i xorshift128plus_avx2(struct rngstate256 *sp)
{
__m256i s1 = sp->state0;
const __m256i s0 = sp->state1;
sp->state0 = s0;
s1 = _mm256_xor_si256(s1, _mm256_slli_epi64(s1, 23));
__m256i state1new = _mm256_xor_si256(_mm256_xor_si256(_mm256_xor_si256(s1, s0),
_mm256_srli_epi64(s1, 18)),
_mm256_srli_epi64(s0, 5));
sp->state1 = state1new;
return _mm256_add_epi64(state1new, s0);
}
// GNU C native vectors let us get the compiler to do stuff like %10 each element
typedef unsigned short v16u __attribute__((vector_size(32)));
__m256i* vec_store_digit_and_space(__m256i vec, __m256i *restrict p)
{
v16u v = (v16u)vec;
v16u ten = (v16u)_mm256_set1_epi16(10);
v16u divisor = (v16u)_mm256_set1_epi16(6554); // ceil((2^16-1) / 10.0)
v16u div6554 = v / divisor; // Basically the entropy from the upper two decimal digits: 0..65.
// Probably some correlation with the modulo-based values, especially dig3, but we do this instead of
// dig4 for more ILP and fewer instructions total.
v16u dig1 = v % ten;
v /= ten;
v16u dig2 = v % ten;
v /= ten;
v16u dig3 = v % ten;
// dig4 would overlap much of the randomness that div6554 gets
const v16u ascii_digitspace = (v16u)_mm256_set1_epi16( (' '<<8) | '0');
v16u *vecbuf = (v16u*)p;
vecbuf[0] = div6554 | ascii_digitspace;
vecbuf[1] = dig1 | ascii_digitspace;
vecbuf[2] = dig2 | ascii_digitspace;
vecbuf[3] = dig3 | ascii_digitspace;
return p + 4; // always a constant number of full vectors
}
void random_decimal_fill_buffer(char *restrict buf, size_t len, struct rngstate256 *restrict rngstate)
{
buf = __builtin_assume_aligned(buf, 32);
// copy to a local so clang can keep state in register, even in the non-inline version
// restrict works for gcc, but apparently clang still thinks that *buf might alias *rngstate
struct rngstate256 rng_local = *rngstate;
__m256i *restrict p = (__m256i*restrict)buf;
__m256i *restrict endbuf = (__m256i*)(buf+len);
static unsigned newline_pos = 0;
do {
__m256i rvec = xorshift128plus_avx2(&rng_local);
p = vec_store_digit_and_space(rvec, p); // stores multiple ASCII vectors from the entropy in rvec
#if 1
// this is buggy at the end or start of a power-of-2 buffer:
// usually there's a too-short line, sometimes a too-long line
const unsigned ncols = 100;
newline_pos += 4*16;
if (newline_pos >= ncols) {
newline_pos -= ncols;
char *cur_pos = (char*)p;
*(cur_pos - newline_pos*2 - 1) = '\n';
}
#endif
// Turning every 100th space into a newline.
// 1) With an overlapping 1B store to a location selected by a counter. A down-counter would be more efficient
// 2) Or by using a different constant for ascii_digitspace to put a newline in one element
// lcm(200, 16) is 400 bytes, so unrolling the loop enough to produce two full lines makes a pattern of full vectors repeat
// lcm(200, 32) is 800 bytes
// a power-of-2 buffer size doesn't hold a whole number of lines :/
// I'm pretty sure this can be solved with low overhead, like maybe 10% at worst.
} while(p <= endbuf-3);
*rngstate = rng_local;
}
#define BUFFER_SIZE (128 * 1024)
const static size_t bufsz = BUFFER_SIZE;
__attribute__((aligned(64))) static char static_buf[BUFFER_SIZE];
int main(int argc, char *argv[])
{
// TODO: choose a seed properly. (Doesn't affect the speed)
struct rngstate256 xorshift_state = {
_mm256_set_epi64x(123, 456, 0x123, 0x456),
_mm256_set_epi64x(789, 101112, 0x789, 0x101112)
};
for (int i=0; i < 1024ULL*1024*1024 / bufsz * 100; i++) {
random_decimal_fill_buffer(static_buf, bufsz, &xorshift_state);
size_t written = write(1, static_buf, bufsz);
(void)written;
//fprintf(stderr, "wrote %#lx of %#lx\n", written, bufsz);
}
}
Compile with gcc, clang, or ICC (or hopefully any other compiler that understands the GNU C dialect of C99, and Intel's intrinsics). GNU C vector extensions are highly convenient to get the compiler to generate the magic numbers for division/modulo using modular multiplicative inverses, and occasional __attribute__s are useful.
This could be written portably, but it would take more code.
Performance notes:
The overlapping-store to insert newlines has significant overhead to decide where to place it (branch mispredictions, and frontend bottlenecks on Core2), but the store itself has no impact on performance. Commenting out just that store instruction in the compiler's asm (leaving all the branching the same) left the performance on Core2 completely unchanged, with repeated runs giving the same time to +/- less than 1%. So I conclude that the store buffer / cache handle it just fine.
Still, using some kind of rotating window of ascii_digitspace with one element having a newline might be even faster, if we unroll enough that any counters/branching go away.
Writing to /dev/null is basically a no-op, so the buffer probably stays hot in L2 cache (256kiB per core on Haswell). The perfect speedup from 128b vectors to 256b vectors is expected: there are no extra instructions, and everything (including the stores) happens with twice the width. The newline-insertion branch is taken twice as often, though. I unfortunately didn't time on my Haswell cygwin setup with that part #ifdefed out.
2.5GHz * 32B / 13.7GB/s = 5.84 cycles per AVX2-store on Haswell. That's pretty good, but could be faster. Maybe there's some overhead in the cygwin system calls than I thought. I didn't try commenting those out in the compiler's asm output (which would ensure that nothing optimized away.)
L1 cache can sustain one 32B store per clock, and L2 is not much lower bandwidth (higher latency, though).
When I looked at IACA a few versions ago (without the branching for newlines, but only getting one ASCII vector per RNG vector), it was predicting something like one 32B vector store per 4 or 5 clocks.
I was hoping to get more of a speedup from extracting more data from each RNG result, based on looking at the asm myself, considering Agner Fog's guides and other optimization resources which I've added links for in the SO x86 tag wiki.)
Likely it would be significantly faster on Skylake, where vector integer multiply and shift can run on twice as many ports (p0 / p1) compared to Haswell (p0 only). xorshift and the digit extraction both use a lot of shifts and multiplies. (Update: Skylake runs it at 3.02 IPC, giving us 3.77 cycles per 32-byte AVX2 store, timed at 0.030s per 1GB iteration, writing to /dev/null on Linux 4.15 on i7-6700k at 3.9GHz.
It doesn't require 64-bit mode to work well. The SSE2 version is just as fast when compiled with -m32, because it doesn't need very many vector registers, and all the 64-bit math is done in vectors, not general-purpose registers.
It's actually slightly faster in 32-bit mode on Core2, because compare/branch macro-fusion only works in 32-bit mode, so there are fewer uops for the out-of-order core (18.3s (1.85 Instructions Per Clock) vs. 16.9s (2.0 IPC)). The smaller code-size from having no REX prefixes also helps Core2's decoders.
Also, some reg-reg vector moves are replaced with loads, since not all the constants fix in vector regs anymore. Since load throughput from L1 cache isn't a bottleneck, this actually helps. (e.g. multiplying by a constant vector of set1(10): movdqa xmm0, xmm10 / pmullw xmm0, xmm1 turns into movdqa xmm0, [constant] / pmullw xmm0, xmm1.) Since reg-reg MOVDQA requires an ALU port, it competes with the real work being done, but a MOVDQA load only competes for front-end decode bandwidth. (Having a 4-byte address inside many instructions cancels out a lot of the gain from saving REX prefixes.
I wouldn't be surprised if saving ALU MOVDQA uops is where the real gains are coming from, since the frontend should be keeping up with the average of 2.0 IPC pretty well.
All these differences disappear on Haswell, where the whole thing should run from the decoded-uop cache, if not the loopback buffer. ALU+branch macro-fusion works in both modes since Nehalem.
yes 4 | tr '\n' ' ' | fold -w 200 | head -c1G– Matthew Crumley Nov 17 '16 at 21:00