How to speed up tricky random number generation

Question

I have some code which performs many log, `tan` and `cos` operations on doubles. I need this to be as fast as possibly. Currently I use code such as ``` #include #include #include "mtwist.h" #include int main(void) { int i; double x; mt_seed(); double u1; double u2; double w1; double w2; x = 0; for(i = 0; i < 100000000; ++i) { u1 = mt_drand(); u2 = mt_drand(); w1 = M_PI*(u1-1/2.0); w2 = -log(u2); x += tan(w1)*(M_PI_2-w1)+log(w2*cos(w1)/(M_PI_2-w1)); } printf("%f ",x); return EXIT_SUCCESS; } ``` I am using gcc. There are two obvious ways to speed this up. The first is to choose a faster RNG. The second is to speed up the transcendental functions. To do this I would like to know - How are tan and cos implemented in assembly on x86? My CPU is the AMD FX-8350 if that makes a difference. (Answer `fcos` for `cos` and `fptan` for `tan`.) - How can you use a lookup table to speed up the calculations? I only need 32 bits of accuracy. Can you make use a table of s…

Accepted Answer

You can try this `log(x)` replacement I wrote using SSE2 intrinsics: ``` #include #include static __m128i EXPONENT_MASK; static __m128i EXPONENT_BIAS; static __m128i EXPONENT_ZERO; static __m128d FIXED_SCALE; static __m128d LOG2ERECIP; static const int EXPONENT_SHIFT = 52; // Required to initialize constants. void sselog_init() { EXPONENT_MASK = _mm_set1_epi64x(0x7ff0000000000000UL); EXPONENT_BIAS = _mm_set1_epi64x(0x00000000000003ffUL); EXPONENT_ZERO = _mm_set1_epi64x(0x3ff0000000000000UL); FIXED_SCALE = _mm_set1_pd(9.31322574615478515625e-10); // 2^-30 LOG2ERECIP = _mm_set1_pd(0.693147180559945309417232121459); // 1/log2(e) } // Extract IEEE754 double exponent as integer. static inline __m128i extractExponent(__m128d x) { return _mm_sub_epi64( _mm_srli_epi64( _mm_and_si128(_mm_castpd_si128(x), EXPONENT_MASK), EXPONENT_SHIFT), EXPONENT_BIAS); } // Set IEEE754 double exponent to zero. static inline __m128d clearExponent(__m128d x) { return _mm_castsi128_pd( _mm_or_si128( _mm_andnot_si128( EXPONENT_MASK, _mm_castpd_si128(x)), EXPONENT_ZERO)); } // Compute log(x) using SSE2 intrinsics to >= 30 bit precision, except denorms. double sselog(double x) { assert(x >= 2.22507385850720138309023271733e-308); // no denormalized // Two independent logarithms could be computed by initializing // base with two different values, either with independent // arguments to _mm_set_pd() or from contiguous memory with // _mm_load_pd(). No other changes should be needed other than to // extract both results at the end of the function (or just return // the packed __m128d). __m128d base = _mm_set_pd(x, x); __m128i iLog = extractExponent(base); __m128i fLog = _mm_setzero_si128(); base = clearExponent(base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); // fLog = _mm_slli_epi64(fLog, 10); // Not needed first time through. fLog = _mm_or_si128(extractExponent(base), fLog); base = clearExponent(base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); fLog = _mm_slli_epi64(fLog, 10); fLog = _mm_or_si128(extractExponent(base), fLog); base = clearExponent(base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); base = _mm_mul_pd(base, base); fLog = _mm_slli_epi64(fLog, 10); fLog = _mm_or_si128(extractExponent(base), fLog); // No _mm_cvtepi64_pd() exists so use _mm_cvtepi32_pd() conversion. iLog = _mm_shuffle_epi32(iLog, 0x8); fLog = _mm_shuffle_epi32(fLog, 0x8); __m128d result = _mm_mul_pd(_mm_cvtepi32_pd(fLog), FIXED_SCALE); result = _mm_add_pd(result, _mm_cvtepi32_pd(iLog)); // Convert from base 2 logarithm and extract result. result = _mm_mul_pd(result, LOG2ERECIP); return ((double *)&result)[0]; // other output in ((double *)&result)[1] } ``` The code implements an algorithm described in this Texas Instruments brief, of repeatedly squaring the argument and concatenating exponent bits. It will not work with denormalized inputs. It provides at least 30 bits of precision. This runs faster than `log()` on one of my machines and slower on the other, so your mileage may vary; I don't claim this is necessarily the best approach. This code, however, actually computes two logarithms in parallel using both halves of a 128-bit SSE2 word (although the function as-is only returns one result), so it could be adapted into one building block of a SIMD calculation of your entire function (and I think `log` is the difficult part as `cos` is pretty well-behaved). Furthermore, your processor supports AVX which can pack 4 double-precision elements into a 256-bit word, and extending this code to AVX should be straightforward. If you choose not to go full SIMD, you could still use both logarithm slots by pipelining - i.e. compute `log(w2*cos(w1)/(M_PI_2-w1))` for the current iteration with `log(u2)` for the next iteration. Even if this function benchmarks slower against `log` in isolation, it may still be worth testing with your actual function. This code does not stress th…

How to speed up tricky random number generation

Solution

Problem

Related problems