在探索性活动的背景下，我开始研究32位平台上整数和定点算术构建块。我的主要目标是ARM32(具体来说是
在过去的三十年中，我已经多次实现了这个和其他算术构建块，但总是在各种架构的汇编语言中实现。然而，此时此刻，我希望这些可以用纯ISO-C编程，而不需要使用intrinsic。理想情况下，单个版本的C代码将在各种架构上生成良好的机器代码。我知道操纵HLL代码以控制机器代码的方法通常很脆弱，但希望处理器架构和工具链已经成熟到使这种方法可行。
一些在汇编语言实现中使用的方法不适合移植到 C 语言。在下面的示例代码中，我选择了六种变体，这些变体似乎适合 HLL 实现。除了生成部分积之外，这些方法有两个基本方法：（1）使用 64 位算术求和部分积，让编译器处理 32 位半之间的进位传播。在这种情况下，存在多种选择来对部分积进行求和。（2）使用 32 位算术进行求和，直接模拟进位标志。在这种情况下，我们可以选择在加法后生成进位 (
在编译器资源管理器上，我专注于gcc 12.2和clang 15.0这两个平台的工具链。我使用了
特别值得注意的是三种情况。虽然我以前曾与编译器工程师合作过，并对基本的代码转换、阶段排序问题等有一定的了解，但我所知道的唯一可能适用于此代码的技术是成语识别，而我不认为这本身就能解释这些观察结果。第一个情况是变体3，其中clang 15.0生成了非常紧凑的代码，仅包含10条指令，我认为这是无法改进的：
相比之下，gcc 生成了两倍的指令并需要两倍的寄存器。我假设它不知道如何在这里使用乘加指令
生成的代码将模拟进位转换为真实进位传播。clang 15.0 使用了九个指令和五个寄存器，如果不花更多时间进行分析，我无法真正理解它试图做什么。第三个观察结果是关于在变体5与变体6产生的机器代码中看到的差异，特别是使用clang时。这些使用相同的基本算法，其中一个变体在加法之前计算模拟进位，另一个在加法之后计算。最终我发现，一种变体，即变体4，在两个工具链和两个架构上都很有效。然而，在我继续处理其他构建块并面临类似的挑战之前，我想询问：
armv7)，并且关注RISC-V32，我预计它会在嵌入式领域中变得越来越重要。我选择检查的第一个示例构建块是无符号64x64->128位整数乘法。该网站上关于该构建块的其他问题没有提供有关32位平台的详细覆盖。在过去的三十年中，我已经多次实现了这个和其他算术构建块，但总是在各种架构的汇编语言中实现。然而，此时此刻，我希望这些可以用纯ISO-C编程，而不需要使用intrinsic。理想情况下，单个版本的C代码将在各种架构上生成良好的机器代码。我知道操纵HLL代码以控制机器代码的方法通常很脆弱，但希望处理器架构和工具链已经成熟到使这种方法可行。
一些在汇编语言实现中使用的方法不适合移植到 C 语言。在下面的示例代码中，我选择了六种变体，这些变体似乎适合 HLL 实现。除了生成部分积之外，这些方法有两个基本方法：（1）使用 64 位算术求和部分积，让编译器处理 32 位半之间的进位传播。在这种情况下，存在多种选择来对部分积进行求和。（2）使用 32 位算术进行求和，直接模拟进位标志。在这种情况下，我们可以选择在加法后生成进位 (
a = a + b; carry = a < b;) 或在加法之前生成进位 (carry = ~a < b; a = a + b;)。下面的变体 1 到 3 属于前一类，变体 5 和 6 属于后一类。在编译器资源管理器上，我专注于gcc 12.2和clang 15.0这两个平台的工具链。我使用了
-O3进行编译。总体发现是，平均而言，clang生成的代码比gcc更高效，并且在clang中，变体之间（指令数和寄存器使用量）的差异更加明显。虽然这在RISC-V作为较新的架构的情况下可以理解，但在armv7的情况下，这让我感到惊讶，因为它已经存在了十多年。特别值得注意的是三种情况。虽然我以前曾与编译器工程师合作过，并对基本的代码转换、阶段排序问题等有一定的了解，但我所知道的唯一可能适用于此代码的技术是成语识别，而我不认为这本身就能解释这些观察结果。第一个情况是变体3，其中clang 15.0生成了非常紧凑的代码，仅包含10条指令，我认为这是无法改进的：
umul64wide:
push {r4, lr}
umull r12, r4, r2, r0
umull lr, r0, r3, r0
umaal lr, r4, r2, r1
umaal r0, r4, r3, r1
ldr r1, [sp, #8]
strd r0, r4, [r1]
mov r0, r12
mov r1, lr
pop {r4, pc}
相比之下，gcc 生成了两倍的指令并需要两倍的寄存器。我假设它不知道如何在这里使用乘加指令
umaal，但这是全部原因吗？在变体6中出现了相反的情况，但不太戏剧化，其中gcc 12.2生成了这个包含18条指令的序列，并且寄存器使用率较低：umul64wide:
mov ip, r0
push {r4, r5, lr}
mov lr, r1
umull r0, r1, r0, r2
ldr r4, [sp, #12]
umull r5, ip, r3, ip
adds r1, r1, r5
umull r2, r5, lr, r2
adc ip, ip, #0
umull lr, r3, lr, r3
adds r1, r1, r2
adc r2, ip, #0
adds r2, r2, r5
adc r3, r3, #0
adds r2, r2, lr
adc r3, r3, #0
strd r2, r3, [r4]
pop {r4, r5, pc}
生成的代码将模拟进位转换为真实进位传播。clang 15.0 使用了九个指令和五个寄存器，如果不花更多时间进行分析，我无法真正理解它试图做什么。第三个观察结果是关于在变体5与变体6产生的机器代码中看到的差异，特别是使用clang时。这些使用相同的基本算法，其中一个变体在加法之前计算模拟进位，另一个在加法之后计算。最终我发现，一种变体，即变体4，在两个工具链和两个架构上都很有效。然而，在我继续处理其他构建块并面临类似的挑战之前，我想询问：
(1) 在下面的代码中，是否有我没有考虑到的编码习惯或算法，可能会导致更优秀的结果？(2) 是否有特定的优化开关，例如假设的-ffrobnicate（请参见here），它没有包含在-O3中，可以帮助编译器更一致地生成有效的代码，以处理这些位操作场景？对于观察到的代码生成显著差异的编译器机制的解释，以及如何影响或解决它们，也可能是有帮助的。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#define VARIANT (3)
#define USE_X64_ASM_REF (0)
/* Multiply two unsigned 64-bit operands a and b. Returns least significant 64
bits of product as return value, most significant 64 bits of product via h.
*/
uint64_t umul64wide (uint64_t a, uint64_t b, uint64_t *h)
{
uint32_t a_lo = (uint32_t)a;
uint32_t a_hi = a >> 32;
uint32_t b_lo = (uint32_t)b;
uint32_t b_hi = b >> 32;
uint64_t p0 = (uint64_t)a_lo * b_lo;
uint64_t p1 = (uint64_t)a_lo * b_hi;
uint64_t p2 = (uint64_t)a_hi * b_lo;
uint64_t p3 = (uint64_t)a_hi * b_hi;
#if VARIANT == 1
uint32_t c = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);
*h = p3 + (p1 >> 32) + (p2 >> 32) + c;
return p0 + ((p1 + p2) << 32);
#elif VARIANT == 2
uint64_t s = (p0 >> 32) + p1;
uint64_t t = (uint32_t)s + p2;
*h = (s >> 32) + (t >> 32) + p3;
return (uint32_t)p0 + (t << 32);
#elif VARIANT == 3
*h = (p1 >> 32) + (((p0 >> 32) + (uint32_t)p1 + p2) >> 32) + p3;
return p0 + ((p1 + p2) << 32);
#elif VARIANT == 4
uint64_t t = (p0 >> 32) + p1 + (uint32_t)p2;
*h = (p2 >> 32) + (t >> 32) + p3;
return (uint32_t)p0 + (t << 32);
#elif VARIANT == 5
uint32_t r0, r1, r2, r3, r4, r5, r6;
r0 = (uint32_t)p0;
r1 = p0 >> 32;
r5 = (uint32_t)p1;
r2 = p1 >> 32;
r1 = r1 + r5;
r6 = r1 < r5;
r2 = r2 + r6;
r6 = (uint32_t)p2;
r5 = p2 >> 32;
r1 = r1 + r6;
r6 = r1 < r6;
r2 = r2 + r6;
r4 = (uint32_t)p3;
r3 = p3 >> 32;
r2 = r2 + r5;
r6 = r2 < r5;
r3 = r3 + r6;
r2 = r2 + r4;
r6 = r2 < r4;
r3 = r3 + r6;
*h = ((uint64_t)r3 << 32) | r2;
return ((uint64_t)r1 << 32) | r0;
#elif VARIANT == 6
uint32_t r0, r1, r2, r3, r4, r5, r6;
r0 = (uint32_t)p0;
r1 = p0 >> 32;
r5 = (uint32_t)p1;
r2 = p1 >> 32;
r4 = ~r1;
r4 = r4 < r5;
r1 = r1 + r5;
r2 = r2 + r4;
r6 = (uint32_t)p2;
r5 = p2 >> 32;
r4 = ~r1;
r4 = r4 < r6;
r1 = r1 + r6;
r2 = r2 + r4;
r4 = (uint32_t)p3;
r3 = p3 >> 32;
r6 = ~r2;
r6 = r6 < r5;
r2 = r2 + r5;
r3 = r3 + r6;
r6 = ~r2;
r6 = r6 < r4;
r2 = r2 + r4;
r3 = r3 + r6;
*h = ((uint64_t)r3 << 32) | r2;
return ((uint64_t)r1 << 32) | r0;
#else
#error unsupported VARIANT
#endif
}
#if defined(__SIZEOF_INT128__)
uint64_t umul64wide_ref (uint64_t a, uint64_t b, uint64_t *h)
{
unsigned __int128 prod = ((unsigned __int128)a) * b;
*h = (uint64_t)(prod >> 32);
return (uint64_t)prod;
}
#elif defined(_MSC_VER) && defined(_WIN64)
#include <intrin.h>
uint64_t umul64wide_ref (uint64_t a, uint64_t b, uint64_t *h)
{
*h = __umulh (a, b);
return a * b;
}
#elif USE_X64_ASM_REF
uint64_t umul64wide_ref (uint64_t a, uint64_t b, uint64_t *h)
{
uint64_t res_l, res_h;
__asm__ (
"movq %2, %%rax;\n\t" // rax = a
"mulq %3;\n\t" // rdx:rax = a * b
"movq %%rdx, %0;\n\t" // res_h = rdx
"movq %%rax, %1;\n\t" // res_l = rax
: "=rm" (res_h), "=rm"(res_l)
: "rm"(a), "rm"(b)
: "%rax", "%rdx");
*h = res_h;
return res_l;
}
#else // generic (and slow) reference implementation
#define ADDCcc(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), cy=t0<cy, t0=t0+t1, t1=t0<t1, cy=cy+t1, t0=t0)
#define ADDcc(a,b,cy,t0,t1) \
(t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)
#define ADDC(a,b,cy,t0,t1) \
(t0=(b)+cy, t1=(a), t0+t1)
uint64_t umul64wide_ref (uint64_t a, uint64_t b, uint64_t *h)
{
uint32_t cy, t0, t1;
uint32_t a_lo = (uint32_t)a;
uint32_t a_hi = a >> 32;
uint32_t b_lo = (uint32_t)b;
uint32_t b_hi = b >> 32;
uint64_t p0 = (uint64_t)a_lo * b_lo;
uint64_t p1 = (uint64_t)a_lo * b_hi;
uint64_t p2 = (uint64_t)a_hi * b_lo;
uint64_t p3 = (uint64_t)a_hi * b_hi;
uint32_t p0_lo = (uint32_t)p0;
uint32_t p0_hi = p0 >> 32;
uint32_t p1_lo = (uint32_t)p1;
uint32_t p1_hi = p1 >> 32;
uint32_t p2_lo = (uint32_t)p2;
uint32_t p2_hi = p2 >> 32;
uint32_t p3_lo = (uint32_t)p3;
uint32_t p3_hi = p3 >> 32;
uint32_t r0 = p0_lo;
uint32_t r1 = ADDcc (p0_hi, p1_lo, cy, t0, t1);
uint32_t r2 = ADDCcc (p1_hi, p2_hi, cy, t0, t1);
uint32_t r3 = ADDC (p3_hi, 0, cy, t0, t1);
r1 = ADDcc (r1, p2_lo, cy, t0, t1);
r2 = ADDCcc (r2, p3_lo, cy, t0, t1);
r3 = ADDC (r3, 0, cy, t0, t1);
*h = ((uint64_t)r3 << 32) + r2;
return ((uint64_t)r1 << 32) + r0;
}
#endif
/*
https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
From: geo <gmars...@gmail.com>
Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
Subject: 64-bit KISS RNGs
Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)
This 64-bit KISS RNG has three components, each nearly
good enough to serve alone. The components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64 (kiss64_t = (kiss64_x << 58) + kiss64_c, \
kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64 (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
kiss64_y ^= (kiss64_y << 43))
#define CNG64 (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
int main (void)
{
uint64_t a, b, res_hi, res_lo, ref_hi, ref_lo, count = 0;
printf ("Smoke test of umul64wide variant %d\n", VARIANT);
do {
a = KISS64;
b = KISS64;
ref_lo = umul64wide_ref (a, b, &ref_hi);
res_lo = umul64wide (a, b, &res_hi);
if ((res_lo ^ ref_lo) | (res_hi ^ ref_hi)) {
printf ("!!!! error: a=%016llx b=%016llx res=%016llx_%016llx ref=%016llx_%016llx\n",
a, b, res_hi, res_lo, ref_hi, ref_lo);
return EXIT_FAILURE;
}
if (!(count & 0xfffffff)) printf ("\r%llu", count);
count++;
} while (count);
return EXIT_SUCCESS;
}
__UINT_FAST64_TYPE__，可用于32 x 32 => 64，让编译器完成工作 - 就像uint64_t一样 - 但这可能会导致函数调用。您是否只想使用纯粹的习惯用法C，只能依赖于uint32_t，没有特殊的32 x 32 => 64或类似于ARM上的umulh的内联汇编，并将其用作64 x 64的构建块？ - Brett Halestdint.h可用。我在问题中提到了“没有内置函数”。同样地，没有内联汇编。编译器中的供应商扩展不受欢迎（可移植性）。40年前，我手工组装了Z80指令并输入了结果十六进制代码。现在，无论是国际象棋还是围棋，没有人能够与计算机相比，我希望我的编译器能够将C代码转换为精选和排列良好的机器指令。我最后一次使用汇编语言编写这个特定的构建块是大约五年前为一个Maxwell级GPU。 - njuffa_ExtInt(128)，这应该可以在32位机器上工作（不像unsigned __int128）。ISO C23可能会标准化unsigned _BitInt(128) foo。（https://en.wikipedia.org/wiki/C2x）。但目前只有clang支持，即使使用`-std=gnu2x`也不支持GCC trunk：https://godbolt.org/z/r3cejeKsr。（一个这么宽的值在32位ISA上通过引用返回，但除此之外，clang似乎对ARM32编译得相当高效，但x86 32看起来很混乱，过早地设置CF并且必须跨mul保存它。） - Peter Cordes