如何实现 SWAR 无符号小于？
How to implement SWAR unsigned less-than?
我正在尝试使用 uint64_t 就好像它是 uint8_t 的 8 个车道；我的目标是实现逐车道小于。给定 x 和 y，如果 x 中相应车道的值小于该车道的值，则此操作应在车道中产生 0xFF 的结果在 y 中，否则 0x00。逐条车道小于或等于也可以。
根据我所见，我猜我需要一个车道差或零操作（定义为 doz(x, y) = if (x < y) then 0 else (x - y)），然后用它来构造一个选择掩码。然而，我看到的所有车道减法方法都是有符号的，我不确定我将如何使用它们来完成这种任务。
有没有一种方法可以做到这一点，使用差或零或其他方式？
事实证明，以 DOZ 为基础是错误的做法。这些都是没有意义的，不要使用它。
However, all the lane-wise subtraction approaches I've seen are signed
这很奇怪，因为减法既不是有符号的也不是无符号的，只有一种减法并且可以两种方式解释。至少，这就是它在 2 的补码世界中的工作方式。
作为参考，SWAR 减法看起来像这样：（来源：SIMD and SWAR Techniques）
SWAR sub z = x - y
z = ((x | H) - (y &~H)) ^ ((x ^~y) & H)
DOZ 可以以此为基础。 full DOZ 有点矫枉过正，如果它是一个有意义的原语。但是 SWAR DOZ 会通过计算差异来工作，然后在 x < y 的情况下将其归零，这是我们一直想要的条件。所以我们只计算它而不是整个 DOZ。那个条件是根据这个：什么时候有借位高位？
- 如果x的高位为0，y的高位为1。
- 如果x和y的高位相同，且高位差为1。等价地：如果x和y的高位相同，并且有次高位的借位。
SWAR sub 的第一部分，((x | H) - (y &~H))，计算（除其他外）第二高位的借位。 SWAR 差异的高位是第二高位借位的倒数（来自 H 的位要么被借位“吃掉”，要么不被借位）。
将它们放在一起，SWAR unsigned-less-than 可以像这样工作：
tmp = ((~x ^ y) & ~((x | H) - (y &~H)) | (~x & y)) & H
less_than_mask = (tmp << 1) - (tmp >> 7)
零件：
(~x ^ y) = “位相同”的掩码，用于“高位相同”~((x | H) - (y &~H)) = 元素低位差，用于“借用次高位”(~x & y) = “x为零，y为一”的掩码，用于“x的高位为零，y的高位为一”& H接近尾部，用于只抓取高位借位对应的位(tmp << 1) - (tmp >> 7) 将上一步抓取的位展开到车道掩码中。备选方案：(tmp >> 7) * 255。这是 SWAR 逻辑明确取决于通道大小的唯一步骤，并且每个通道都需要相同，即使对于 SWAR 子您可以混合通道大小。
可以通过应用 De Morgan 规则在表达式级别删除一个操作：
tmp = (~(x ^ y | (x | H) - (y & ~H)) | ~x & y) & H
但是 ~x 无论如何都需要计算，因此在汇编级别可能无济于事，具体取决于它的编译方式。
也许可以进行一些简化。
这是一种独立于体系结构的方法。我确定它可以使用改进，但它似乎工作正常。使用 x86 gcc/clang，它编译为 20/19 条指令。
想法是首先解决两个字节是否小于128 时的问题，将每个字节的第7 位设置为该结果。然后修补其他情况。最后把bit 7往下抹。
#include <stdio.h>
#include <stdint.h>
uint64_t bwlt(uint64_t a, uint64_t b) {
uint64_t lo7 = ~0ull / 255 * 127, // low 7 bits set in each byte
alo7 = a & lo7, // mask low 7 bits in a
blo7 = b & lo7, // mask low 7 bits in b
r = (lo7 - alo7 + blo7) & ~lo7, // set 8th bits with a < b
diff = (a ^ b) & ~lo7; // 8th bits that differ
r &= ~(a & diff); // unset if a[i]_7=1,b[i]_7=0
r |= b & diff; // set if a[i]_7=0,b[i]_7=1
return (r << 1) - (r >> 7);
}
int main(void) {
uint64_t a = 0x11E1634052A6B7CB;
uint64_t b = 0x1EAEF1E85F26734E;
printf("r=%016llx\n", bwlt(a, b));
return 0;
}
一个测试用例：
$ gcc foo.c -o foo
$ ./foo
r=ff00ffffff000000
我很高兴弄清楚如何使用 64bit unsigned int 并仅使用 logical operators 和算术 + 和 -.[=29 创建 SWAR x LT (Less Than) y 函数=]
我在网上看了一些资料(https://www.chessprogramming.org/SIMD_and_SWAR_Techniques)，从那里我知道这个函数可以从减法开始完成(x - y)。
查看最高位的含义：x、y和(x - y)在使用unsigned int时，我创造了以下真相table 其中：
当 LT 条件发生时，R（结果）为 1。
D是减法(x-y)的最高位，
X为待测X值的最高位，
Y为待测Y值的最高位
D X Y | R
0 0 0 | 0
0 0 1 | 1
0 1 0 | 0
0 1 1 | 0
1 0 0 | 1
1 0 1 | 1
1 1 0 | 0
1 1 1 | 1
将卡诺图 (https://getcalc.com/karnaugh-map/3variable-kmap-solver.htm) 应用于上面的 table 我们得到以下公式：
(~ X & Y) | (D & ~ X) | (D & Y)
宏 SWARLTU(x, y) 的起源（参见下面的文件 swar.h）。
由于不满意，我观察了编译器如何生成宏 SWARLTU 的汇编代码，然后按照该代码编写了宏 SWARLTU2(x, y)（参见文件 swar.h以下）。最后一个宏应该在逻辑上进行优化。
此代码的限制是 LT 结果的值为 0x80 而不是问题中要求的 0xFF。
程序可以通过三种不同的方式启动：
没有参数，在本例中它将对随机数执行10次测试。
只有一个参数，该参数表示随机测试的次数
有两个参数，两个0xnnnnn形式的数字，在这种情况下只会显示输入值的控件。
这里是代码：
文件 swar.h（此文件还包含其他 SWAR 宏，例如：SHL、SHR）
#ifndef SWAR_H
#define SWAR_H
/*
https://www.chessprogramming.org/SIMD_and_SWAR_Techniques
SWAR add z = x + y
z = ((x &~H) + (y &~H)) ^ ((x ^ y) & H)
SWAR sub z = x - y
z = ((x | H) - (y &~H)) ^ ((x ^~y) & H)
SWAR average z = (x+y)/2 based on x + y = (x^y) + 2*(x&y)
z = (x & y) + (((x ^ y) & ~L) >> 1)
*/
// 0 1 2 3 4 5 6 7
#define SWARH 0x8080808080808080LL
#define SWARL 0x0101010101010101LL
#define SWARADD(x,y) \
((( (x) &~SWARH) + ( (y) &~SWARH)) ^ (( (x) ^ (y) ) & SWARH))
#define SWARSUB(x,y) \
((( (x) | SWARH) - ( (y) &~SWARH)) ^ (( (x) ^~(y) ) & SWARH))
#define SWARAVE(x,y) \
(( (x) & (y) ) + ((( (x) ^ (y)) & ~SWARL) >> 1))
#define SWARLTI(x,y) \
( SWARSUB(x,y) & SWARH )
#define SWARSHL(x) \
(((x)&(~SWARH))<<1)
#define SWARSHR(x) \
(((x)&(~SWARL))>>1)
/*** Computing unsigned less than
Truth table considering the HIGH bit setting of
Differece, X Value, Y Value
D X Y | R
0 0 0 | 0
0 0 1 | 1
0 1 0 | 0
0 1 1 | 0
1 0 0 | 1
1 0 1 | 1
1 1 0 | 0
1 1 1 | 1
***/
#define _SWARDH (SWARSUB(x,y) & SWARH)
#define _SWARXH ((x)&SWARH)
#define _SWARYH ((y)&SWARH)
#define SWARLTU(x,y) \
((~_SWARXH & _SWARYH) | (_SWARDH & ~_SWARXH) | (_SWARDH & _SWARYH))
// Elaborated from the generated ASM of the previous.
#define SWARLTU2(X,Y) \
((((~(X & SWARH)) & ((((~(X ^ Y)) & SWARH) ^ ((X | SWARH) - Y)) | Y)) | \
((((~(X ^ Y)) & SWARH) ^ ((X | SWARH) - Y)) & Y)) & SWARH)
#endif // SWAR_H
文件main.c
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>
#include "swar.h"
char * verifyltu(char * rs,uint64_t x, uint64_t y, uint64_t v);
void printvalues(uint64_t x,uint64_t y,uint64_t r,uint64_t r1);
uint64_t rand64();
int main(int argc, char *argv[])
{
int rnd=1;
size_t i,n=10;
uint64_t x=0,y=0,r,r1;
srand(time(NULL));
if (argc>1) {
if (argc==2) {
n=strtoul(argv[1],NULL,0);
} else {
x=strtoull(argv[1],NULL,0);
y=strtoull(argv[2],NULL,0);
rnd=0;
}
}
if (rnd) {
for(i=0;i<n;i++) {
x=rand64();
y=rand64();
r=SWARLTU(x,y);
r1=SWARLTU2(x,y);
printvalues(x,y,r,r1);
}
} else {
r=SWARLTU(x,y);
r1=SWARLTU2(x,y);
printvalues(x,y,r,r1);
}
return 0;
}
char * verifyltu(char * rs,uint64_t x, uint64_t y, uint64_t v)
{
size_t i;
uint8_t *xs, *ys, *vs;
xs=(uint8_t *)&x; ys=(uint8_t *)&y;
vs=(uint8_t *)&v;
for(i=0;i<sizeof(uint64_t);i++) {
if ( ( xs[i]<ys[i] && vs[i]&0x80) ||
( !(xs[i]<ys[i]) && !(vs[i]&0x80) ) )
{
rs[i*2]='*';rs[i*2+1]=' ';
} else {
rs[i*2]='-';rs[i*2+1]=' ';
}
}
rs[i*2]=0;
return rs;
}
void printvalues(uint64_t x,uint64_t y,uint64_t r,uint64_t r1)
{
char rs[17],rs1[17];
printf(
"X %016" PRIX64 " <\n"
"Y %016" PRIX64 "\n"
" ----------------\n"
"LTU %016" PRIX64 "\n"
"*=Ok %s\n"
"LTU2 %016" PRIX64 "\n"
"*=Ok %s\n\n",
x,y,
r,verifyltu(rs,x,y,r),
r1,verifyltu(rs1,x,y,r1)
);
}
uint64_t rand64()
{
uint64_t x;
x=rand(); x=(x<<32)+rand();
return x;
}
哪种方法最快取决于目标平台的处理器架构中可用的指令类型，例如移位加法、加法、三输入加法、三输入逻辑指令。它还取决于人们是否需要吞吐量或延迟优化版本，以及处理器架构的超标量。
以下 ISO C 99 代码提供了两种选择。一个使用直接取自文献 (LTU_VARIANT = 1) 的无符号字节比较，另一个 (LTU_VARIANT = 0) 我在减半加法的基础上设计了自己（即两个整数的总和除以二, 向下舍入）。这是基于以下事实：对于 [0,255] 中的补码整数 a、b，a < u b ⇔ ~a + b >= 256.
然而，这将需要 9 位求和，因此我们可以使用 a < u b ⇔ ((~a + b) >> 1) >= 128 代替，其中可以在 8 以内计算平均值通过众所周知的位旋转技术位。我所知道的唯一提供 SIMD 减半加法作为硬件指令 vhadd 的处理器架构是 Arm NEON.
我已经包含了一个用于功能测试的测试框架，但需要进行基准测试以确定哪个版本在给定平台上的性能更好。
此答案可能与其他答案部分重叠；剥金橘皮的方法只有这么多。
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#define LTU_VARIANT (1) // 0 or 1
#define UINT64_H8 (0x8080808080808080U) // byte-wise sign bits (MSBs)
uint64_t sign_to_mask8 (uint64_t a)
{
a = a & UINT64_H8; // isolate sign bits
a = a + a - (a >> 7); // extend them to full byte to create mask
return a;
}
uint64_t vhaddu8 (uint64_t a, uint64_t b)
{
/* Peter L. Montgomery's observation (newsgroup comp.arch, 2000/02/11,
https://groups.google.com/d/msg/comp.arch/gXFuGZtZKag/_5yrz2zDbe4J):
(A+B)/2 = (A AND B) + (A XOR B)/2.
*/
return (a & b) + (((a ^ b) >> 1) & ~UINT64_H8);
}
uint64_t ltu8_core (uint64_t a, uint64_t b)
{
/* Sebastiano Vigna, "Broadword implementation of rank/select queries."
In: International Workshop on Experimental and Efficient Algorithms,
pp. 154-168, Springer Berlin Heidelberg, 2008.
*/
return (((a | UINT64_H8) - (b & ~UINT64_H8)) | (a ^ b)) ^ (a | ~b);
}
uint64_t vcmpltu8 (uint64_t a, uint64_t b)
{
#if LTU_VARIANT==1
return sign_to_mask8 (ltu8_core (a, b));
#else // LTU_VARIANT
return sign_to_mask8 (vhaddu8 (~a, b));
#endif // LTU_VARIANT
}
uint64_t ref_func (uint64_t a, uint64_t b)
{
uint8_t a0 = (uint8_t)((a >> 0) & 0xff);
uint8_t a1 = (uint8_t)((a >> 8) & 0xff);
uint8_t a2 = (uint8_t)((a >> 16) & 0xff);
uint8_t a3 = (uint8_t)((a >> 24) & 0xff);
uint8_t a4 = (uint8_t)((a >> 32) & 0xff);
uint8_t a5 = (uint8_t)((a >> 40) & 0xff);
uint8_t a6 = (uint8_t)((a >> 48) & 0xff);
uint8_t a7 = (uint8_t)((a >> 56) & 0xff);
uint8_t b0 = (uint8_t)((b >> 0) & 0xff);
uint8_t b1 = (uint8_t)((b >> 8) & 0xff);
uint8_t b2 = (uint8_t)((b >> 16) & 0xff);
uint8_t b3 = (uint8_t)((b >> 24) & 0xff);
uint8_t b4 = (uint8_t)((b >> 32) & 0xff);
uint8_t b5 = (uint8_t)((b >> 40) & 0xff);
uint8_t b6 = (uint8_t)((b >> 48) & 0xff);
uint8_t b7 = (uint8_t)((b >> 56) & 0xff);
uint8_t r0 = (a0 < b0) ? 0xff : 0x00;
uint8_t r1 = (a1 < b1) ? 0xff : 0x00;
uint8_t r2 = (a2 < b2) ? 0xff : 0x00;
uint8_t r3 = (a3 < b3) ? 0xff : 0x00;
uint8_t r4 = (a4 < b4) ? 0xff : 0x00;
uint8_t r5 = (a5 < b5) ? 0xff : 0x00;
uint8_t r6 = (a6 < b6) ? 0xff : 0x00;
uint8_t r7 = (a7 < b7) ? 0xff : 0x00;
return ( ((uint64_t)r0 << 0) +
((uint64_t)r1 << 8) +
((uint64_t)r2 << 16) +
((uint64_t)r3 << 24) +
((uint64_t)r4 << 32) +
((uint64_t)r5 << 40) +
((uint64_t)r6 << 48) +
((uint64_t)r7 << 56) );
}
/*
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, ref, n = 0;
printf ("Testing vcmpltu8: byte-wise unsigned comparison with mask result\n");
printf ("using LTU variant %d\n", LTU_VARIANT);
do {
a = KISS64;
b = KISS64;
res = vcmpltu8 (a, b);
ref = ref_func (a, b);
if (res != ref) {
printf ("\nerr @ a=%016" PRIx64 " b=%016" PRIx64 " : res=%016" PRIx64 " ref=%016" PRIx64 "\n",
a, b, res, ref);
return EXIT_FAILURE;
}
n++;
if (!(n & 0xffffff)) printf ("\r%016" PRIx64, n);
} while (a);
printf ("\ntest passed\n");
return EXIT_SUCCESS;
}
我想到了
uint64_t magic(uint64_t a, uint64_t b) {
auto H = 0x8080808080808080ull;
auto c = (a|H) - (b&(~H));
auto d = a^b;
auto e = ((a & d) | (c & (~d))) & H;
return e ^ H;
}
逻辑与哈罗德的逻辑几乎相同；不同之处在于将最高位解释为
c = 1aaaaaaa -> the carry to the H bit is 0 only IFF a<b
0bbbbbbb
a = 80, b = 00 different sign -> select b, i.e. ~a
a = 00, b = 80 different sign -> select b, i.e. ~a
a = 80, b = 80 same sign -> select ~c
a = 00, b = 00 same sign -> select ~c
如果可以使用反转掩码（即 b >= a），那么最后一个 ^H 也可以省略。
对于 arm64 和 x64 使用 clang / godbolt 结果的指令数将是有和（没有）sign_to_mask。
instructions arm64 x64
-------------+----------+---------
vhaddu8 | 5 (8) | 8 (13)
ltu8_core | 7 (10)| 11 (15)
magic | 9 (10)| 12 (15)
harold | 9 (11)| 14 (17)
bwlt | 10 (12)| 15 (18)
SirJoBlack | 11 (13)| 16 (19)
我正在尝试使用 uint64_t 就好像它是 uint8_t 的 8 个车道；我的目标是实现逐车道小于。给定 x 和 y，如果 x 中相应车道的值小于该车道的值，则此操作应在车道中产生 0xFF 的结果在 y 中，否则 0x00。逐条车道小于或等于也可以。
根据我所见，我猜我需要一个车道差或零操作（定义为 doz(x, y) = if (x < y) then 0 else (x - y)），然后用它来构造一个选择掩码。然而，我看到的所有车道减法方法都是有符号的，我不确定我将如何使用它们来完成这种任务。
有没有一种方法可以做到这一点，使用差或零或其他方式？
事实证明，以 DOZ 为基础是错误的做法。这些都是没有意义的，不要使用它。
However, all the lane-wise subtraction approaches I've seen are signed
这很奇怪，因为减法既不是有符号的也不是无符号的，只有一种减法并且可以两种方式解释。至少，这就是它在 2 的补码世界中的工作方式。
作为参考，SWAR 减法看起来像这样：（来源：SIMD and SWAR Techniques）
SWAR sub z = x - y
z = ((x | H) - (y &~H)) ^ ((x ^~y) & H)
DOZ 可以以此为基础。 full DOZ 有点矫枉过正，如果它是一个有意义的原语。但是 SWAR DOZ 会通过计算差异来工作，然后在 x < y 的情况下将其归零，这是我们一直想要的条件。所以我们只计算它而不是整个 DOZ。那个条件是根据这个：什么时候有借位高位？
- 如果x的高位为0，y的高位为1。
- 如果x和y的高位相同，且高位差为1。等价地：如果x和y的高位相同，并且有次高位的借位。
SWAR sub 的第一部分，((x | H) - (y &~H))，计算（除其他外）第二高位的借位。 SWAR 差异的高位是第二高位借位的倒数（来自 H 的位要么被借位“吃掉”，要么不被借位）。
将它们放在一起，SWAR unsigned-less-than 可以像这样工作：
tmp = ((~x ^ y) & ~((x | H) - (y &~H)) | (~x & y)) & H
less_than_mask = (tmp << 1) - (tmp >> 7)
零件：
(~x ^ y)= “位相同”的掩码，用于“高位相同”~((x | H) - (y &~H))= 元素低位差，用于“借用次高位”(~x & y)= “x为零，y为一”的掩码，用于“x的高位为零，y的高位为一”& H接近尾部，用于只抓取高位借位对应的位(tmp << 1) - (tmp >> 7)将上一步抓取的位展开到车道掩码中。备选方案：(tmp >> 7) * 255。这是 SWAR 逻辑明确取决于通道大小的唯一步骤，并且每个通道都需要相同，即使对于 SWAR 子您可以混合通道大小。
可以通过应用 De Morgan 规则在表达式级别删除一个操作：
tmp = (~(x ^ y | (x | H) - (y & ~H)) | ~x & y) & H
但是 ~x 无论如何都需要计算，因此在汇编级别可能无济于事，具体取决于它的编译方式。
也许可以进行一些简化。
这是一种独立于体系结构的方法。我确定它可以使用改进，但它似乎工作正常。使用 x86 gcc/clang，它编译为 20/19 条指令。
想法是首先解决两个字节是否小于128 时的问题，将每个字节的第7 位设置为该结果。然后修补其他情况。最后把bit 7往下抹。
#include <stdio.h>
#include <stdint.h>
uint64_t bwlt(uint64_t a, uint64_t b) {
uint64_t lo7 = ~0ull / 255 * 127, // low 7 bits set in each byte
alo7 = a & lo7, // mask low 7 bits in a
blo7 = b & lo7, // mask low 7 bits in b
r = (lo7 - alo7 + blo7) & ~lo7, // set 8th bits with a < b
diff = (a ^ b) & ~lo7; // 8th bits that differ
r &= ~(a & diff); // unset if a[i]_7=1,b[i]_7=0
r |= b & diff; // set if a[i]_7=0,b[i]_7=1
return (r << 1) - (r >> 7);
}
int main(void) {
uint64_t a = 0x11E1634052A6B7CB;
uint64_t b = 0x1EAEF1E85F26734E;
printf("r=%016llx\n", bwlt(a, b));
return 0;
}
一个测试用例：
$ gcc foo.c -o foo
$ ./foo
r=ff00ffffff000000
我很高兴弄清楚如何使用 64bit unsigned int 并仅使用 logical operators 和算术 + 和 -.[=29 创建 SWAR x LT (Less Than) y 函数=]
我在网上看了一些资料(https://www.chessprogramming.org/SIMD_and_SWAR_Techniques)，从那里我知道这个函数可以从减法开始完成(x - y)。
查看最高位的含义：x、y和(x - y)在使用unsigned int时，我创造了以下真相table 其中：
当 LT 条件发生时，R（结果）为 1。
D是减法(x-y)的最高位，
X为待测X值的最高位，
Y为待测Y值的最高位
D X Y | R
0 0 0 | 0
0 0 1 | 1
0 1 0 | 0
0 1 1 | 0
1 0 0 | 1
1 0 1 | 1
1 1 0 | 0
1 1 1 | 1
将卡诺图 (https://getcalc.com/karnaugh-map/3variable-kmap-solver.htm) 应用于上面的 table 我们得到以下公式：
(~ X & Y) | (D & ~ X) | (D & Y)
宏 SWARLTU(x, y) 的起源（参见下面的文件 swar.h）。
由于不满意，我观察了编译器如何生成宏 SWARLTU 的汇编代码，然后按照该代码编写了宏 SWARLTU2(x, y)（参见文件 swar.h以下）。最后一个宏应该在逻辑上进行优化。
此代码的限制是 LT 结果的值为 0x80 而不是问题中要求的 0xFF。
程序可以通过三种不同的方式启动：
没有参数，在本例中它将对随机数执行10次测试。
只有一个参数，该参数表示随机测试的次数
有两个参数，两个0xnnnnn形式的数字，在这种情况下只会显示输入值的控件。
这里是代码：
文件 swar.h（此文件还包含其他 SWAR 宏，例如：SHL、SHR）
#ifndef SWAR_H
#define SWAR_H
/*
https://www.chessprogramming.org/SIMD_and_SWAR_Techniques
SWAR add z = x + y
z = ((x &~H) + (y &~H)) ^ ((x ^ y) & H)
SWAR sub z = x - y
z = ((x | H) - (y &~H)) ^ ((x ^~y) & H)
SWAR average z = (x+y)/2 based on x + y = (x^y) + 2*(x&y)
z = (x & y) + (((x ^ y) & ~L) >> 1)
*/
// 0 1 2 3 4 5 6 7
#define SWARH 0x8080808080808080LL
#define SWARL 0x0101010101010101LL
#define SWARADD(x,y) \
((( (x) &~SWARH) + ( (y) &~SWARH)) ^ (( (x) ^ (y) ) & SWARH))
#define SWARSUB(x,y) \
((( (x) | SWARH) - ( (y) &~SWARH)) ^ (( (x) ^~(y) ) & SWARH))
#define SWARAVE(x,y) \
(( (x) & (y) ) + ((( (x) ^ (y)) & ~SWARL) >> 1))
#define SWARLTI(x,y) \
( SWARSUB(x,y) & SWARH )
#define SWARSHL(x) \
(((x)&(~SWARH))<<1)
#define SWARSHR(x) \
(((x)&(~SWARL))>>1)
/*** Computing unsigned less than
Truth table considering the HIGH bit setting of
Differece, X Value, Y Value
D X Y | R
0 0 0 | 0
0 0 1 | 1
0 1 0 | 0
0 1 1 | 0
1 0 0 | 1
1 0 1 | 1
1 1 0 | 0
1 1 1 | 1
***/
#define _SWARDH (SWARSUB(x,y) & SWARH)
#define _SWARXH ((x)&SWARH)
#define _SWARYH ((y)&SWARH)
#define SWARLTU(x,y) \
((~_SWARXH & _SWARYH) | (_SWARDH & ~_SWARXH) | (_SWARDH & _SWARYH))
// Elaborated from the generated ASM of the previous.
#define SWARLTU2(X,Y) \
((((~(X & SWARH)) & ((((~(X ^ Y)) & SWARH) ^ ((X | SWARH) - Y)) | Y)) | \
((((~(X ^ Y)) & SWARH) ^ ((X | SWARH) - Y)) & Y)) & SWARH)
#endif // SWAR_H
文件main.c
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>
#include "swar.h"
char * verifyltu(char * rs,uint64_t x, uint64_t y, uint64_t v);
void printvalues(uint64_t x,uint64_t y,uint64_t r,uint64_t r1);
uint64_t rand64();
int main(int argc, char *argv[])
{
int rnd=1;
size_t i,n=10;
uint64_t x=0,y=0,r,r1;
srand(time(NULL));
if (argc>1) {
if (argc==2) {
n=strtoul(argv[1],NULL,0);
} else {
x=strtoull(argv[1],NULL,0);
y=strtoull(argv[2],NULL,0);
rnd=0;
}
}
if (rnd) {
for(i=0;i<n;i++) {
x=rand64();
y=rand64();
r=SWARLTU(x,y);
r1=SWARLTU2(x,y);
printvalues(x,y,r,r1);
}
} else {
r=SWARLTU(x,y);
r1=SWARLTU2(x,y);
printvalues(x,y,r,r1);
}
return 0;
}
char * verifyltu(char * rs,uint64_t x, uint64_t y, uint64_t v)
{
size_t i;
uint8_t *xs, *ys, *vs;
xs=(uint8_t *)&x; ys=(uint8_t *)&y;
vs=(uint8_t *)&v;
for(i=0;i<sizeof(uint64_t);i++) {
if ( ( xs[i]<ys[i] && vs[i]&0x80) ||
( !(xs[i]<ys[i]) && !(vs[i]&0x80) ) )
{
rs[i*2]='*';rs[i*2+1]=' ';
} else {
rs[i*2]='-';rs[i*2+1]=' ';
}
}
rs[i*2]=0;
return rs;
}
void printvalues(uint64_t x,uint64_t y,uint64_t r,uint64_t r1)
{
char rs[17],rs1[17];
printf(
"X %016" PRIX64 " <\n"
"Y %016" PRIX64 "\n"
" ----------------\n"
"LTU %016" PRIX64 "\n"
"*=Ok %s\n"
"LTU2 %016" PRIX64 "\n"
"*=Ok %s\n\n",
x,y,
r,verifyltu(rs,x,y,r),
r1,verifyltu(rs1,x,y,r1)
);
}
uint64_t rand64()
{
uint64_t x;
x=rand(); x=(x<<32)+rand();
return x;
}
哪种方法最快取决于目标平台的处理器架构中可用的指令类型，例如移位加法、加法、三输入加法、三输入逻辑指令。它还取决于人们是否需要吞吐量或延迟优化版本，以及处理器架构的超标量。
以下 ISO C 99 代码提供了两种选择。一个使用直接取自文献 (LTU_VARIANT = 1) 的无符号字节比较，另一个 (LTU_VARIANT = 0) 我在减半加法的基础上设计了自己（即两个整数的总和除以二, 向下舍入）。这是基于以下事实：对于 [0,255] 中的补码整数 a、b，a < u b ⇔ ~a + b >= 256.
然而，这将需要 9 位求和，因此我们可以使用 a < u b ⇔ ((~a + b) >> 1) >= 128 代替，其中可以在 8 以内计算平均值通过众所周知的位旋转技术位。我所知道的唯一提供 SIMD 减半加法作为硬件指令 vhadd 的处理器架构是 Arm NEON.
我已经包含了一个用于功能测试的测试框架，但需要进行基准测试以确定哪个版本在给定平台上的性能更好。
此答案可能与其他答案部分重叠；剥金橘皮的方法只有这么多。
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#define LTU_VARIANT (1) // 0 or 1
#define UINT64_H8 (0x8080808080808080U) // byte-wise sign bits (MSBs)
uint64_t sign_to_mask8 (uint64_t a)
{
a = a & UINT64_H8; // isolate sign bits
a = a + a - (a >> 7); // extend them to full byte to create mask
return a;
}
uint64_t vhaddu8 (uint64_t a, uint64_t b)
{
/* Peter L. Montgomery's observation (newsgroup comp.arch, 2000/02/11,
https://groups.google.com/d/msg/comp.arch/gXFuGZtZKag/_5yrz2zDbe4J):
(A+B)/2 = (A AND B) + (A XOR B)/2.
*/
return (a & b) + (((a ^ b) >> 1) & ~UINT64_H8);
}
uint64_t ltu8_core (uint64_t a, uint64_t b)
{
/* Sebastiano Vigna, "Broadword implementation of rank/select queries."
In: International Workshop on Experimental and Efficient Algorithms,
pp. 154-168, Springer Berlin Heidelberg, 2008.
*/
return (((a | UINT64_H8) - (b & ~UINT64_H8)) | (a ^ b)) ^ (a | ~b);
}
uint64_t vcmpltu8 (uint64_t a, uint64_t b)
{
#if LTU_VARIANT==1
return sign_to_mask8 (ltu8_core (a, b));
#else // LTU_VARIANT
return sign_to_mask8 (vhaddu8 (~a, b));
#endif // LTU_VARIANT
}
uint64_t ref_func (uint64_t a, uint64_t b)
{
uint8_t a0 = (uint8_t)((a >> 0) & 0xff);
uint8_t a1 = (uint8_t)((a >> 8) & 0xff);
uint8_t a2 = (uint8_t)((a >> 16) & 0xff);
uint8_t a3 = (uint8_t)((a >> 24) & 0xff);
uint8_t a4 = (uint8_t)((a >> 32) & 0xff);
uint8_t a5 = (uint8_t)((a >> 40) & 0xff);
uint8_t a6 = (uint8_t)((a >> 48) & 0xff);
uint8_t a7 = (uint8_t)((a >> 56) & 0xff);
uint8_t b0 = (uint8_t)((b >> 0) & 0xff);
uint8_t b1 = (uint8_t)((b >> 8) & 0xff);
uint8_t b2 = (uint8_t)((b >> 16) & 0xff);
uint8_t b3 = (uint8_t)((b >> 24) & 0xff);
uint8_t b4 = (uint8_t)((b >> 32) & 0xff);
uint8_t b5 = (uint8_t)((b >> 40) & 0xff);
uint8_t b6 = (uint8_t)((b >> 48) & 0xff);
uint8_t b7 = (uint8_t)((b >> 56) & 0xff);
uint8_t r0 = (a0 < b0) ? 0xff : 0x00;
uint8_t r1 = (a1 < b1) ? 0xff : 0x00;
uint8_t r2 = (a2 < b2) ? 0xff : 0x00;
uint8_t r3 = (a3 < b3) ? 0xff : 0x00;
uint8_t r4 = (a4 < b4) ? 0xff : 0x00;
uint8_t r5 = (a5 < b5) ? 0xff : 0x00;
uint8_t r6 = (a6 < b6) ? 0xff : 0x00;
uint8_t r7 = (a7 < b7) ? 0xff : 0x00;
return ( ((uint64_t)r0 << 0) +
((uint64_t)r1 << 8) +
((uint64_t)r2 << 16) +
((uint64_t)r3 << 24) +
((uint64_t)r4 << 32) +
((uint64_t)r5 << 40) +
((uint64_t)r6 << 48) +
((uint64_t)r7 << 56) );
}
/*
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, ref, n = 0;
printf ("Testing vcmpltu8: byte-wise unsigned comparison with mask result\n");
printf ("using LTU variant %d\n", LTU_VARIANT);
do {
a = KISS64;
b = KISS64;
res = vcmpltu8 (a, b);
ref = ref_func (a, b);
if (res != ref) {
printf ("\nerr @ a=%016" PRIx64 " b=%016" PRIx64 " : res=%016" PRIx64 " ref=%016" PRIx64 "\n",
a, b, res, ref);
return EXIT_FAILURE;
}
n++;
if (!(n & 0xffffff)) printf ("\r%016" PRIx64, n);
} while (a);
printf ("\ntest passed\n");
return EXIT_SUCCESS;
}
我想到了
uint64_t magic(uint64_t a, uint64_t b) {
auto H = 0x8080808080808080ull;
auto c = (a|H) - (b&(~H));
auto d = a^b;
auto e = ((a & d) | (c & (~d))) & H;
return e ^ H;
}
逻辑与哈罗德的逻辑几乎相同；不同之处在于将最高位解释为
c = 1aaaaaaa -> the carry to the H bit is 0 only IFF a<b
0bbbbbbb
a = 80, b = 00 different sign -> select b, i.e. ~a
a = 00, b = 80 different sign -> select b, i.e. ~a
a = 80, b = 80 same sign -> select ~c
a = 00, b = 00 same sign -> select ~c
如果可以使用反转掩码（即 b >= a），那么最后一个 ^H 也可以省略。
对于 arm64 和 x64 使用 clang / godbolt 结果的指令数将是有和（没有）sign_to_mask。
instructions arm64 x64
-------------+----------+---------
vhaddu8 | 5 (8) | 8 (13)
ltu8_core | 7 (10)| 11 (15)
magic | 9 (10)| 12 (15)
harold | 9 (11)| 14 (17)
bwlt | 10 (12)| 15 (18)
SirJoBlack | 11 (13)| 16 (19)