如何加速棘手的随机数生成

Question

我有一些代码执行许多日志tan和cos双打操作.我需要这个尽可能快.目前我使用的代码如

#include <stdio.h>
#include <stdlib.h>
#include "mtwist.h"
#include <math.h>


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\n",x); 

   return EXIT_SUCCESS;
}

我正在使用gcc.

有两种明显的方法可以加快速度.首先是选择更快的RNG.第二是加快先验功能.
要做到这一点,我想知道

1. 如何在x86上的程序集中实现tan和cos？我的CPU是AMD FX-8350,如果它有所作为.(答案fcos为cos和fptan的tan.)
2. 如何使用查找表来加速计算？我只需要32位的精度.例如,你可以使用一个大小为2 ^ 16的表来加速tan和cos操作吗？

在英特尔优化手册说

如果不需要使用80位的扩展精度来评估超越函数,则应用程序应考虑使用基于软件的替代方法,例如使用插值技术的基于查找表的算法.通过选择所需的数值精度和查找表的大小,并利用SSE和SSE2指令的并行性,可以通过这些技术提高超越性能.

根据这个非常有用的表,fcos有延迟154并fptan具有延迟166-231.

---

您可以使用编译我的代码

gcc -O3 -Wall random.c mtwist-1.5/mtwist.c -lm -o random

我的C代码使用的梅森倍捻机RNG的C代码来自 这里.您应该能够运行我的代码来测试它.如果你不能,请告诉我.

---

更新 @rhashimoto加速了我的代码从20秒到6秒!

RNG似乎应该可以加速.然而,在我的测试中,http://www.math.sci.hiroshima-u.ac.jp/~%20m-mat/MT/SFMT/index.html#dSFMT花费了相同的时间(有没有人看到任何不同的东西) ).如果有人能找到更快的RNG(通过所有死硬测试),我将非常感激.

请显示您建议的任何改进的实际时间,因为这有助于确定哪些有效或无效.

Answer 1

你可以重写

tan(w1)*(M_PI_2-w1)+log(w2*cos(w1)/(M_PI_2-w1))

如

tan(w1)*(M_PI_2-w1) + log(cos(w1)/(M_PI_2-w1)) + log(w2).

根据w1这里的不同,您可以使用minimax多项式来处理这些东西.组成64个左右,每个组成1/64的范围,你可能只需要3级或4级.

你计算w2为

w2 = -log(u2);

对于均匀u2在(0,1).所以你真的是在计算log(log(1/u2)).我打赌你可以使用类似的技巧来获得log(log(1/x))分块的多项式近似值(0,1).(该功能的作用吓人近0和1,所以你可能需要做一些花哨的存在吧.)

Answer 2

我喜欢@ tmyklebu的建议,即为整体计算创建一个minimax近似值.有一些很好的工具可以帮助解决这个问题,包括Remez函数逼近工具包

为了提高速度,你可以比MT做得更好; 例如,参见Dr. Dobbs文章:快速,高质量,并行随机数发生器

另请参阅ACML - AMD核心数学库以利用SSE和SSE2.

Answer 3

你可以试试log(x)我用 SSE2 内在函数写的这个替换：

#include <assert.h>
#include <immintrin.h>

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]
}

该代码实现了德州仪器简介中描述的算法，即反复平方参数并连接指数位。它不适用于非规范化输入。它提供至少 30 位的精度。

这比log()我的一台机器运行得快，而在另一台机器上运行得更慢，所以你的里程可能会有所不同；我并不认为这一定是最好的方法。但是，此代码实际上使用 128 位 SSE2 字的两半并行计算两个对数（尽管函数原样仅返回一个结果），因此它可以改编为整个 SIMD 计算的构建块功能（我认为这log是困难的部分，因为cos表现得非常好）。此外，您的处理器支持 AVX，它可以将 4 个双精度元素打包到一个 256 位字中，并且将此代码扩展到 AVX 应该很简单。

如果您选择不走全SIMD，你仍然可以通过流水线同时使用对数位-即计算log(w2*cos(w1)/(M_PI_2-w1))与当前迭代log(u2)的下一个迭代。

即使这个函数log单独进行基准测试比较慢，它仍然值得用你的实际函数测试。这段代码根本不会对数据缓存施加压力，因此它可能对其他有压力的代码更友好（例如cos，使用查找表的 ）。此外，微指令调度可以通过周围的代码改进（或不改进），这取决于其他代码是否使用 SSE。

我的其他建议（从评论中重复）是：

• 尝试-march=native -mtune=native让 gcc 优化您的 CPU。
• 避免在同一个参数上同时调用tan和cos- 使用sincos或触发标识。
• 考虑使用 GPU（例如 OpenCL）。

似乎计算sin而不是计算更好cos- 原因是您可以将它用于tan_w1 = sin_w1/sqrt(1.0 - sin_w1*sin_w1). cos我最初建议使用，在计算 时会丢失正确的符号tan。正如其他回答者所说，通过在 [-pi/2, pi/2] 上使用极小极大多项式，您似乎可以获得很好的加速。此链接上的 7 项函数（确保您得到minimaxsin，而不是taylorsin）似乎工作得很好。

所以这是用所有 SSE2 超越近似值重写的程序：

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <immintrin.h>
#include "mtwist.h"

#if defined(__AVX__)
#define VECLEN 4
#elif defined(__SSE2__)
#define VECLEN 2
#else
#error // No SIMD available.
#endif

#if VECLEN == 4
#define VBROADCAST(K) { K, K, K, K };
typedef double vdouble __attribute__((vector_size(32)));
typedef long vlong __attribute__((vector_size(32)));
#elif VECLEN == 2
#define VBROADCAST(K) { K, K };
typedef double vdouble __attribute__((vector_size(16)));
typedef long vlong __attribute__((vector_size(16)));
#endif

static const vdouble FALLBACK_THRESHOLD = VBROADCAST(1.0 - 0.001);

vdouble sse_sin(vdouble x) {
   static const vdouble a0 = VBROADCAST(1.0);
   static const vdouble a1 = VBROADCAST(-1.666666666640169148537065260055e-1);
   static const vdouble a2 = VBROADCAST( 8.333333316490113523036717102793e-3);
   static const vdouble a3 = VBROADCAST(-1.984126600659171392655484413285e-4);
   static const vdouble a4 = VBROADCAST( 2.755690114917374804474016589137e-6);
   static const vdouble a5 = VBROADCAST(-2.502845227292692953118686710787e-8);
   static const vdouble a6 = VBROADCAST( 1.538730635926417598443354215485e-10);

   vdouble xx = x*x;
   return x*(a0 + xx*(a1 + xx*(a2 + xx*(a3 + xx*(a4 + xx*(a5 + xx*a6))))));
}

static inline vlong shiftRight(vlong x, int bits) {
#if VECLEN == 4
   __m128i lo = (__m128i)_mm256_extractf128_si256((__m256i)x, 0);
   __m128i hi = (__m128i)_mm256_extractf128_si256((__m256i)x, 1);
   return (vlong)
      _mm256_insertf128_si256(
         _mm256_castsi128_si256(_mm_srli_epi64(lo, bits)),
         _mm_srli_epi64(hi, bits),
         1);
#elif VECLEN == 2
   return (vlong)_mm_srli_epi64((__m128i)x, bits);
#endif
}

static inline vlong shiftLeft(vlong x, int bits) {
#if VECLEN == 4
   __m128i lo = (__m128i)_mm256_extractf128_si256((__m256i)x, 0);
   __m128i hi = (__m128i)_mm256_extractf128_si256((__m256i)x, 1);
   return (vlong)
      _mm256_insertf128_si256(
         _mm256_castsi128_si256(_mm_slli_epi64(lo, bits)),
         _mm_slli_epi64(hi, bits),
         1);
#elif VECLEN == 2
   return (vlong)_mm_slli_epi64((__m128i)x, bits);
#endif
}

static const vlong EXPONENT_MASK = VBROADCAST(0x7ff0000000000000L);
static const vlong EXPONENT_BIAS = VBROADCAST(0x00000000000003ffL);
static const vlong EXPONENT_ZERO = VBROADCAST(0x3ff0000000000000L);
static const vdouble FIXED_SCALE = VBROADCAST(9.31322574615478515625e-10); // 2^-30
static const vdouble LOG2ERECIP = VBROADCAST(0.6931471805599453094172);
static const int EXPONENT_SHIFT = 52;

// Extract IEEE754 double exponent as integer.
static inline vlong extractExponent(vdouble x) {
   return shiftRight((vlong)x & EXPONENT_MASK, EXPONENT_SHIFT) - EXPONENT_BIAS;
}

// Set IEEE754 double exponent to zero.
static inline vdouble clearExponent(vdouble x) {
   return (vdouble)(((vlong)x & ~EXPONENT_MASK) | EXPONENT_ZERO);
}

// Compute log(x) using SSE2 intrinsics to >= 30 bit precision, except
// denorms.
vdouble sse_log(vdouble base) {
   vlong iLog = extractExponent(base);
   vlong fLog = VBROADCAST(0);

   base = clearExponent(base);
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   fLog = shiftLeft(fLog, 10);
   fLog |= extractExponent(base);

   base = clearExponent(base);
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   fLog = shiftLeft(fLog, 10);
   fLog |= extractExponent(base);

   base = clearExponent(base);
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   base = base*base;
   fLog = shiftLeft(fLog, 10);
   fLog |= extractExponent(base);

   // No _mm_cvtepi64_pd() exists so use 32-bit conversion to double.
#if VECLEN == 4
   __m128i iLogLo = _mm256_extractf128_si256((__m256i)iLog, 0);
   __m128i iLogHi = _mm256_extractf128_si256((__m256i)iLog, 1);
   iLogLo = _mm_srli_si128(_mm_shuffle_epi32(iLogLo, 0x80), 8);
   iLogHi = _mm_slli_si128(_mm_shuffle_epi32(iLogHi, 0x08), 8);

   __m128i fLogLo = _mm256_extractf128_si256((__m256i)fLog, 0);
   __m128i fLogHi = _mm256_extractf128_si256((__m256i)fLog, 1);
   fLogLo = _mm_srli_si128(_mm_shuffle_epi32(fLogLo, 0x80), 8);
   fLogHi = _mm_slli_si128(_mm_shuffle_epi32(fLogHi, 0x08), 8);

   vdouble result = _mm256_cvtepi32_pd(iLogHi | iLogLo) +
      FIXED_SCALE*_mm256_cvtepi32_pd(fLogHi | fLogLo);
#elif VECLEN == 2
   iLog = (vlong)_mm_shuffle_epi32((__m128i)iLog, 0x8);
   fLog = (vlong)_mm_shuffle_epi32((__m128i)fLog, 0x8);

   vdouble result = _mm_cvtepi32_pd((__m128i)iLog) +
      FIXED_SCALE*_mm_cvtepi32_pd((__m128i)fLog);
#endif

   // Convert from base 2 logarithm and extract result.
   return LOG2ERECIP*result;
}

// Original computation.
double fallback(double u1, double u2) {
   double w1 = M_PI*(u1-1/2.0);
   double w2 = -log(u2);
   return tan(w1)*(M_PI_2-w1)+log(w2*cos(w1)/(M_PI_2-w1));
}

int main() {
   static const vdouble ZERO = VBROADCAST(0.0)
   static const vdouble ONE = VBROADCAST(1.0);
   static const vdouble ONE_HALF = VBROADCAST(0.5);
   static const vdouble PI = VBROADCAST(M_PI);
   static const vdouble PI_2 = VBROADCAST(M_PI_2);

   int i,j;
   vdouble x = ZERO;
   for(i = 0; i < 100000000; i += VECLEN) {
      vdouble u1, u2;
      for (j = 0; j < VECLEN; ++j) {
         ((double *)&u1)[j] = mt_drand();
         ((double *)&u2)[j] = mt_drand();
      }

      vdouble w1 = PI*(u1 - ONE_HALF);
      vdouble w2 = -sse_log(u2);

      vdouble sin_w1 = sse_sin(w1);
      vdouble sin2_w1 = sin_w1*sin_w1;

#if VECLEN == 4
      int nearOne = _mm256_movemask_pd(sin2_w1 >= FALLBACK_THRESHOLD);
#elif VECLEN == 2
      int nearOne = _mm_movemask_pd(sin2_w1 >= FALLBACK_THRESHOLD);
#endif
      if (!nearOne) {
#if VECLEN == 4
         vdouble cos_w1 = _mm256_sqrt_pd(ONE - sin2_w1);
#elif VECLEN == 2
         vdouble cos_w1 = _mm_sqrt_pd(ONE - sin2_w1);
#endif
         vdouble tan_w1 = sin_w1/cos_w1;

         x += tan_w1*(PI_2 - w1) + sse_log(w2*cos_w1/(PI_2 - w1));
      }
      else {
         vdouble result;
         for (j = 0; j < VECLEN; ++j)
            ((double *)&result)[j] = fallback(((double *)&u1)[j], ((double *)&u2)[j]);
         x += result;
      }
   }

   double sum = 0.0;
   for (i = 0; i < VECLEN; ++i)
      sum += ((double *)&x)[i];

   printf("%lf\n", sum);
   return 0;
}

我遇到了一个烦人的问题 - sin±pi/2 附近的近似误差会使值稍微超出 [-1, 1] 并且 (1) 导致 的计算tan无效并且 (2) 在 log 参数时导致过大的影响接近 0。为了避免这种情况，代码测试是否sin(w1)^2“接近”到 1，如果是，则返回到原始的全双精度路径。“关闭”的定义FALLBACK_THRESHOLD在程序的顶部——我随意将它设置为 0.999，这似乎仍然返回 OP 原始程序范围内的值，但对性能几乎没有影响。

我已经编辑了代码以使用针对 SIMD的特定于gcc 的语法扩展。如果您的编译器没有这些扩展，那么您可以返回编辑历史记录。如果在编译器中启用，代码现在使用 AVX 一次处理 4 个双精度（而不是 SSE2 的 2 个双精度）。

在我的机器上没有调用mt_seed()以获得可重复结果的结果是：

Version   Time         Result
original  14.653 secs  -1917488837.945067
SSE        7.380 secs  -1917488837.396841
AVX        6.271 secs  -1917488837.422882

由于先验近似，SSE/AVX 结果与原始结果不同是有道理的。我认为您应该能够进行调整FALLBACK_THRESHOLD以权衡精度和速度。我不确定为什么 SSE 和 AVX 结果彼此略有不同。

Answer 4

首先，稍微改造一下。这是原来的总和：

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));
}

这个总和在数学上是等价的：

for(i = 0; i < 100000000; ++i) {
    u1 = M_PI - mt_drand()* M_PI;
    u2 = mt_drand();
    x += u1 / tan (u1) + log (sin (u1) / u1) + log (- log (u2));
}

并且既然应该等价于用1.0-mt_rand()替换mt_drand()，我们可以让u1 = mt_drand() * M_PI。

for(i = 0; i < 100000000; ++i) {
    u1 = mt_drand()* M_PI;
    u2 = mt_drand();
    x += u1 / tan (u1) + log (sin (u1) / u1) + log (- log (u2));
}

因此，现在可以很好地将其分为可以单独处理的随机变量的两个函数；x += f (u1) + g (u2)。这两个函数在长距离上都非常平滑。对于 u1 > 0.03，f 非常平滑，对于较小的值，1 / f 非常平滑。除了接近 0 或 1 的值外，g 非常平滑。所以我们可以对区间 [0 .. 0.01]、[0.01 .. 0.02] 等使用 100 种不同的近似值。除了选择正确的近似值是耗时的。

要解决这个问题：区间 [0 .. 1] 中的线性随机函数将在区间 [0 .. 0.01] 中具有一定数量的值，在 [0.01 .. 0.02] 中具有其他数量的值，依此类推。我认为您可以通过假设正态分布来计算 100,000,000 中有多少随机数落入区间 [0 .. 0.01] 。然后计算剩余的有多少落入 [0.01 .. 0.02] 等等。如果您计算出有 999,123 个数字落入 [0.00, 0.01]，那么您将在区间中生成该数量的随机数，并对区间中的所有数字使用相同的近似值。

要在区间 [0.33 .. 0.34] 中找到 f (x) 的近似值，仅作为示例，您可以对 [-1 .. 1] 中的 x 近似 f (0.335 + x / 200)。通过采用 n 次插值多项式，在 Chebysev 节点 xk = cos (pi * (2k - 1) / 2n) 进行插值，您将获得相当好的结果。

顺便说一句，旧的 x87 三角函数和对数运算的性能很慢。绝对远不及评估低次多项式。并且间隔足够小，您不需要很高的多项式次数。