相同除数时的快速 AVX512 模
Posted
技术标签:
【中文标题】相同除数时的快速 AVX512 模【英文标题】：Fast AVX512 modulo when same divisor 【发布时间】：2018-05-31 00:32:09 【问题描述】：我试图找到潜在阶乘素数（n!+-1 形式的数）的除数，因为我最近购买了 Skylake-X 工作站，我认为使用 AVX512 指令可以加快速度。
算法简单，主要步骤是对同一个除数重复取模。主要是循环大范围的 n 值。这是用 c 编写的幼稚方法（P 是素数表）：
uint64_t factorial_naive(uint64_t const nmin, uint64_t const nmax, const uint64_t *restrict P)
uint64_t n, i, residue;
for (i = 0; i < APP_BUFLEN; i++)
residue = 2;
for (n=3; n <= nmax; n++)
residue *= n;
residue %= P[i];
// Lets check if we found factor
if (nmin <= n)
if( residue == 1)
report_factor(n, -1, P[i]);
if(residue == P[i]- 1)
report_factor(n, 1, P[i]);
return EXIT_SUCCESS;
这里的想法是检查大范围的 n，例如1,000,000 -> 10,000,000 针对同一组除数。所以我们将对同一个除数取模数百万次。使用 DIV 非常慢，因此根据计算范围有几种可能的方法。在我的例子中，n 很可能小于 10^7，潜在除数 p 小于 10,000 G (
我发现了一些用于 power pc 的旧代码，其中 FMA 用于在使用双精度时获得高达 106 位的精确乘积（我猜）。所以我将这种方法转换为 AVX 512 汇编器（英特尔内部）。这是 FMA 方法的一个简单版本，它基于 Dekker (1971) 的工作，Dekker 产品和 TwoProduct 的 FMA 版本在尝试查找/谷歌搜索背后的基本原理时是有用的词。这个论坛也讨论过这种方法（例如 here）。
int64_t factorial_FMA(uint64_t const nmin, uint64_t const nmax, const uint64_t *restrict P)
uint64_t n, i;
double prime_double, prime_double_reciprocal, quotient, residue;
double nr, n_double, prime_times_quotient_high, prime_times_quotient_low;
for (i = 0; i < APP_BUFLEN; i++)
residue = 2.0;
prime_double = (double)P[i];
prime_double_reciprocal = 1.0 / prime_double;
n_double = 3.0;
for (n=3; n <= nmax; n++)
nr = n_double * residue;
quotient = fma(nr, prime_double_reciprocal, rounding_constant);
quotient -= rounding_constant;
prime_times_quotient_high= prime_double * quotient;
prime_times_quotient_low = fma(prime_double, quotient, -prime_times_quotient_high);
residue = fma(residue, n, -prime_times_quotient_high) - prime_times_quotient_low;
if (residue < 0.0) residue += prime_double;
n_double += 1.0;
// Lets check if we found factor
if (nmin <= n)
if( residue == 1.0)
report_factor(n, -1, P[i]);
if(residue == prime_double - 1.0)
report_factor(n, 1, P[i]);
return EXIT_SUCCESS;
这里我使用了魔法常数
static const double rounding_constant = 6755399441055744.0;
那是 2^51 + 2^52 双打的幻数。
我将其转换为 AVX512（每个循环 32 个潜在除数）并使用 IACA 分析结果。它告诉吞吐量瓶颈：由于分配资源不可用，后端和后端分配停滞不前。 我对汇编程序不是很有经验，所以我的问题是我可以做些什么来加快速度并解决这个后端瓶颈？
AVX512 代码在这里，也可以从 github 找到
uint64_t factorial_AVX512_unrolled_four(uint64_t const nmin, uint64_t const nmax, const uint64_t *restrict P)
// we are trying to find a factor for a factorial numbers : n! +-1
//nmin is minimum n we want to report and nmax is maximum. P is table of primes
// we process 32 primes in one loop.
// naive version of the algorithm is int he function factorial_naive
// and simple version of the FMA based approach in the function factorial_simpleFMA
const double one_table[8] __attribute__ ((aligned(64))) =1.0, 1.0, 1.0,1.0,1.0,1.0,1.0,1.0;
uint64_t n;
__m512d zero, rounding_const, one, n_double;
__m512i prime1, prime2, prime3, prime4;
__m512d residue1, residue2, residue3, residue4;
__m512d prime_double_reciprocal1, prime_double_reciprocal2, prime_double_reciprocal3, prime_double_reciprocal4;
__m512d quotient1, quotient2, quotient3, quotient4;
__m512d prime_times_quotient_high1, prime_times_quotient_high2, prime_times_quotient_high3, prime_times_quotient_high4;
__m512d prime_times_quotient_low1, prime_times_quotient_low2, prime_times_quotient_low3, prime_times_quotient_low4;
__m512d nr1, nr2, nr3, nr4;
__m512d prime_double1, prime_double2, prime_double3, prime_double4;
__m512d prime_minus_one1, prime_minus_one2, prime_minus_one3, prime_minus_one4;
__mmask8 negative_reminder_mask1, negative_reminder_mask2, negative_reminder_mask3, negative_reminder_mask4;
__mmask8 found_factor_mask11, found_factor_mask12, found_factor_mask13, found_factor_mask14;
__mmask8 found_factor_mask21, found_factor_mask22, found_factor_mask23, found_factor_mask24;
// load data and initialize cariables for loop
rounding_const = _mm512_set1_pd(rounding_constant);
one = _mm512_load_pd(one_table);
zero = _mm512_setzero_pd ();
// load primes used to sieve
prime1 = _mm512_load_epi64((__m512i *) &P[0]);
prime2 = _mm512_load_epi64((__m512i *) &P[8]);
prime3 = _mm512_load_epi64((__m512i *) &P[16]);
prime4 = _mm512_load_epi64((__m512i *) &P[24]);
// convert primes to double
prime_double1 = _mm512_cvtepi64_pd (prime1); // vcvtqq2pd
prime_double2 = _mm512_cvtepi64_pd (prime2); // vcvtqq2pd
prime_double3 = _mm512_cvtepi64_pd (prime3); // vcvtqq2pd
prime_double4 = _mm512_cvtepi64_pd (prime4); // vcvtqq2pd
// calculates 1.0/ prime
prime_double_reciprocal1 = _mm512_div_pd(one, prime_double1);
prime_double_reciprocal2 = _mm512_div_pd(one, prime_double2);
prime_double_reciprocal3 = _mm512_div_pd(one, prime_double3);
prime_double_reciprocal4 = _mm512_div_pd(one, prime_double4);
// for comparison if we have found factors for n!+1
prime_minus_one1 = _mm512_sub_pd(prime_double1, one);
prime_minus_one2 = _mm512_sub_pd(prime_double2, one);
prime_minus_one3 = _mm512_sub_pd(prime_double3, one);
prime_minus_one4 = _mm512_sub_pd(prime_double4, one);
// residue init
residue1 = _mm512_set1_pd(2.0);
residue2 = _mm512_set1_pd(2.0);
residue3 = _mm512_set1_pd(2.0);
residue4 = _mm512_set1_pd(2.0);
// double counter init
n_double = _mm512_set1_pd(3.0);
// main loop starts here. typical value for nmax can be 5,000,000 -> 10,000,000
for (n=3; n<=nmax; n++) // main loop
// timings for instructions:
// _mm512_load_epi64 = vmovdqa64 : L 1, T 0.5
// _mm512_load_pd = vmovapd : L 1, T 0.5
// _mm512_set1_pd
// _mm512_div_pd = vdivpd : L 23, T 16
// _mm512_cvtepi64_pd = vcvtqq2pd : L 4, T 0,5
// _mm512_mul_pd = vmulpd : L 4, T 0.5
// _mm512_fmadd_pd = vfmadd132pd, vfmadd213pd, vfmadd231pd : L 4, T 0.5
// _mm512_fmsub_pd = vfmsub132pd, vfmsub213pd, vfmsub231pd : L 4, T 0.5
// _mm512_sub_pd = vsubpd : L 4, T 0.5
// _mm512_cmplt_pd_mask = vcmppd : L ?, Y 1
// _mm512_mask_add_pd = vaddpd : L 4, T 0.5
// _mm512_cmpeq_pd_mask = vcmppd L ?, Y 1
// _mm512_kor = korw L 1, T 1
// nr = residue * n
nr1 = _mm512_mul_pd (residue1, n_double);
nr2 = _mm512_mul_pd (residue2, n_double);
nr3 = _mm512_mul_pd (residue3, n_double);
nr4 = _mm512_mul_pd (residue4, n_double);
// quotient = nr * 1.0/ prime_double + rounding_constant
quotient1 = _mm512_fmadd_pd(nr1, prime_double_reciprocal1, rounding_const);
quotient2 = _mm512_fmadd_pd(nr2, prime_double_reciprocal2, rounding_const);
quotient3 = _mm512_fmadd_pd(nr3, prime_double_reciprocal3, rounding_const);
quotient4 = _mm512_fmadd_pd(nr4, prime_double_reciprocal4, rounding_const);
// quotient -= rounding_constant, now quotient is rounded to integer
// countient should be at maximum nmax (10,000,000)
quotient1 = _mm512_sub_pd(quotient1, rounding_const);
quotient2 = _mm512_sub_pd(quotient2, rounding_const);
quotient3 = _mm512_sub_pd(quotient3, rounding_const);
quotient4 = _mm512_sub_pd(quotient4, rounding_const);
// now we calculate high and low for prime * quotient using decker product (FMA).
// quotient is calculated using approximation but this is accurate for given quotient
prime_times_quotient_high1 = _mm512_mul_pd(quotient1, prime_double1);
prime_times_quotient_high2 = _mm512_mul_pd(quotient2, prime_double2);
prime_times_quotient_high3 = _mm512_mul_pd(quotient3, prime_double3);
prime_times_quotient_high4 = _mm512_mul_pd(quotient4, prime_double4);
prime_times_quotient_low1 = _mm512_fmsub_pd(quotient1, prime_double1, prime_times_quotient_high1);
prime_times_quotient_low2 = _mm512_fmsub_pd(quotient2, prime_double2, prime_times_quotient_high2);
prime_times_quotient_low3 = _mm512_fmsub_pd(quotient3, prime_double3, prime_times_quotient_high3);
prime_times_quotient_low4 = _mm512_fmsub_pd(quotient4, prime_double4, prime_times_quotient_high4);
// now we calculate new reminder using decker product and using original values
// we subtract above calculated prime * quotient (quotient is aproximation)
residue1 = _mm512_fmsub_pd(residue1, n_double, prime_times_quotient_high1);
residue2 = _mm512_fmsub_pd(residue2, n_double, prime_times_quotient_high2);
residue3 = _mm512_fmsub_pd(residue3, n_double, prime_times_quotient_high3);
residue4 = _mm512_fmsub_pd(residue4, n_double, prime_times_quotient_high4);
residue1 = _mm512_sub_pd(residue1, prime_times_quotient_low1);
residue2 = _mm512_sub_pd(residue2, prime_times_quotient_low2);
residue3 = _mm512_sub_pd(residue3, prime_times_quotient_low3);
residue4 = _mm512_sub_pd(residue4, prime_times_quotient_low4);
// lets check if reminder < 0
negative_reminder_mask1 = _mm512_cmplt_pd_mask(residue1,zero);
negative_reminder_mask2 = _mm512_cmplt_pd_mask(residue2,zero);
negative_reminder_mask3 = _mm512_cmplt_pd_mask(residue3,zero);
negative_reminder_mask4 = _mm512_cmplt_pd_mask(residue4,zero);
// we and prime back to reminder using mask if it was < 0
residue1 = _mm512_mask_add_pd(residue1, negative_reminder_mask1, residue1, prime_double1);
residue2 = _mm512_mask_add_pd(residue2, negative_reminder_mask2, residue2, prime_double2);
residue3 = _mm512_mask_add_pd(residue3, negative_reminder_mask3, residue3, prime_double3);
residue4 = _mm512_mask_add_pd(residue4, negative_reminder_mask4, residue4, prime_double4);
n_double = _mm512_add_pd(n_double,one);
// if we are below nmin then we continue next iteration
if (n < nmin) continue;
// Lets check if we found any factors, residue 1 == n!-1
found_factor_mask11 = _mm512_cmpeq_pd_mask(one, residue1);
found_factor_mask12 = _mm512_cmpeq_pd_mask(one, residue2);
found_factor_mask13 = _mm512_cmpeq_pd_mask(one, residue3);
found_factor_mask14 = _mm512_cmpeq_pd_mask(one, residue4);
// residue prime -1 == n!+1
found_factor_mask21 = _mm512_cmpeq_pd_mask(prime_minus_one1, residue1);
found_factor_mask22 = _mm512_cmpeq_pd_mask(prime_minus_one2, residue2);
found_factor_mask23 = _mm512_cmpeq_pd_mask(prime_minus_one3, residue3);
found_factor_mask24 = _mm512_cmpeq_pd_mask(prime_minus_one4, residue4);
if (found_factor_mask12 | found_factor_mask11 | found_factor_mask13 | found_factor_mask14 |
found_factor_mask21 | found_factor_mask22 | found_factor_mask23|found_factor_mask24)
// we find factor very rarely
double *residual_list1 = (double *) &residue1;
double *residual_list2 = (double *) &residue2;
double *residual_list3 = (double *) &residue3;
double *residual_list4 = (double *) &residue4;
double *prime_list1 = (double *) &prime_double1;
double *prime_list2 = (double *) &prime_double2;
double *prime_list3 = (double *) &prime_double3;
double *prime_list4 = (double *) &prime_double4;
for (int i=0; i <8; i++)
if( residual_list1[i] == 1.0)
report_factor((uint64_t) n, -1, (uint64_t) prime_list1[i]);
if( residual_list2[i] == 1.0)
report_factor((uint64_t) n, -1, (uint64_t) prime_list2[i]);
if( residual_list3[i] == 1.0)
report_factor((uint64_t) n, -1, (uint64_t) prime_list3[i]);
if( residual_list4[i] == 1.0)
report_factor((uint64_t) n, -1, (uint64_t) prime_list4[i]);
if(residual_list1[i] == (prime_list1[i] - 1.0))
report_factor((uint64_t) n, 1, (uint64_t) prime_list1[i]);
if(residual_list2[i] == (prime_list2[i] - 1.0))
report_factor((uint64_t) n, 1, (uint64_t) prime_list2[i]);
if(residual_list3[i] == (prime_list3[i] - 1.0))
report_factor((uint64_t) n, 1, (uint64_t) prime_list3[i]);
if(residual_list4[i] == (prime_list4[i] - 1.0))
report_factor((uint64_t) n, 1, (uint64_t) prime_list4[i]);
return EXIT_SUCCESS;
【问题讨论】：
赞成一个详细且问得很好的问题。欢迎使用 Stack Overflow！ 出于好奇，这个if(residue == prime_double - 1.0) 是否可靠地工作（==）？仅通过阅读源代码对我来说并不明显，这些值将仅保持整数并在双尾数限制内，因此不会丢失低位数字。但它可能是，取决于fma 的实现......对我来说仍然足够脆弱，值得额外的源评论，为什么它应该工作。
@Nuutti：FMA 吞吐量的后端瓶颈很好，这意味着您正在使机器的 FMA 吞吐量饱和，而不是延迟或前端的瓶颈。 （我认为这就是“分配资源”的意思，但发布 IACA 摘要输出。）总会有某种瓶颈。就正确应用蛮力而言，FMA 吞吐量（端口 0 / 端口 5 饱和）是您想要达到的瓶颈。运行得更快需要重新组合您的操作以执行更多的 FMA 和更少的 add / mul，或者以其他方式节省操作，但这可能无法获得准确的结果。
IACA_trace_analysis：github.com/NudeSurfer/Factoring/blob/master/…IACA 分析：github.com/NudeSurfer/Factoring/blob/master/IACA_analysis.txt
另外，您不需要那么快地进行分支。假设某个特定因素成功的概率非常低，您可以将所有掩码组合在一起并每千次检查一次吗？迭代？然后，如果它显示成功，您可以重新运行该块以准确找出它是哪个因素。
【参考方案1】：
正如一些评论者所建议的那样：“后端”瓶颈是您对该代码的期望。这表明你吃得很好，这正是你想要的。
看报告，这部分应该有机会：
// Lets check if we found any factors, residue 1 == n!-1
found_factor_mask11 = _mm512_cmpeq_pd_mask(one, residue1);
found_factor_mask12 = _mm512_cmpeq_pd_mask(one, residue2);
found_factor_mask13 = _mm512_cmpeq_pd_mask(one, residue3);
found_factor_mask14 = _mm512_cmpeq_pd_mask(one, residue4);
// residue prime -1 == n!+1
found_factor_mask21 = _mm512_cmpeq_pd_mask(prime_minus_one1, residue1);
found_factor_mask22 = _mm512_cmpeq_pd_mask(prime_minus_one2, residue2);
found_factor_mask23 = _mm512_cmpeq_pd_mask(prime_minus_one3, residue3);
found_factor_mask24 = _mm512_cmpeq_pd_mask(prime_minus_one4, residue4);
if (found_factor_mask12 | found_factor_mask11 | found_factor_mask13 | found_factor_mask14 |
found_factor_mask21 | found_factor_mask22 | found_factor_mask23|found_factor_mask24)
来自 IACA 分析：
| 1 | 1.0 | | | | | | | | kmovw r11d, k0
| 1 | 1.0 | | | | | | | | kmovw eax, k1
| 1 | 1.0 | | | | | | | | kmovw ecx, k2
| 1 | 1.0 | | | | | | | | kmovw esi, k3
| 1 | 1.0 | | | | | | | | kmovw edi, k4
| 1 | 1.0 | | | | | | | | kmovw r8d, k5
| 1 | 1.0 | | | | | | | | kmovw r9d, k6
| 1 | 1.0 | | | | | | | | kmovw r10d, k7
| 1 | | 1.0 | | | | | | | or r11d, eax
| 1 | | | | | | | 1.0 | | or r11d, ecx
| 1 | | 1.0 | | | | | | | or r11d, esi
| 1 | | | | | | | 1.0 | | or r11d, edi
| 1 | | 1.0 | | | | | | | or r11d, r8d
| 1 | | | | | | | 1.0 | | or r11d, r9d
| 1* | | | | | | | | | or r11d, r10d
处理器将生成的比较掩码 (k0-k7) 移至常规寄存器以进行“或”运算。您应该能够消除这些动作，并在 6ops 与 8 中进行“或”汇总。
注意：found_factor_mask 类型定义为__mmask8，它们应该是__mask16（512 位fector 中的16x 双浮点数）。这可能会让编译器进行一些优化。如果没有，请按照评论者的说明进行组装。
还有相关的：有多少部分的迭代会触发这个 or-mask 子句？正如另一位评论者所观察到的，您应该能够通过累积“或”操作来展开它。在每次展开迭代结束时（或 N 次迭代后）检查累积的“或”值，如果为“真”，则返回并重新计算值以找出触发它的 n 值。
（并且，您可以在“滚动”中进行二分搜索以找到匹配的 n 值——这可能会有所收获）。
接下来，你应该可以摆脱这个循环中的检查了：
// if we are below nmin then we continue next iteration, we
if (n < nmin) continue;
这里显示：
| 1* | | | | | | | | | cmp r14, 0x3e8
| 0*F | | | | | | | | | jb 0x229
这可能不是一个巨大的收益，因为预测器将（可能）得到这个（大部分）正确，但您应该通过为两个“阶段”设置两个不同的循环来获得一些收益：
n=3 到 n=nmin-1 n=nmin 及以上即使你获得一个周期，也就是 3%。而且由于这通常与上面的大“或”运算有关，因此可能会发现更多的聪明之处。
【讨论】：
删除分支并将循环分成两个阶段可能根本无济于事，如果代码真的是 后端 绑定，即使它被采用并且可能创建一些前端-结束气泡。cmp/jcc 在端口 6 上运行，该端口没有任何向量 ALU。但值得一试，较低的 uop 吞吐量将使其对超线程更友好，而代价是稍大的 uop-cache 占用空间。以上是关于相同除数时的快速 AVX512 模的主要内容，如果未能解决你的问题，请参考以下文章