Terje Mathisen 发明了一种非常快的 itoa(),不需要查找表。如果您对其工作原理的解释不感兴趣,请跳至性能或实现。
15 多年前,Terje Mathisen 提出了以 10 为基数的并行 itoa()。其想法是采用 32 位值并将其分解为两个 5 位数字的块。 (在谷歌上快速搜索“Terje Mathisen itoa”给出了这篇文章:http://computer-programming-forum.com/46-asm/7aa4b50bce8dd985.htm http://computer-programming-forum.com/46-asm/7aa4b50bce8dd985.htm)
我们这样开始:
void itoa(char *buf, uint32_t val)
{
lo = val % 100000;
hi = val / 100000;
itoa_half(&buf[0], hi);
itoa_half(&buf[5], lo);
}
现在我们只需要一个可以将域 [0, 99999] 中的任何整数转换为字符串的算法。一个天真的方法可能是:
// 0 <= val <= 99999
void itoa_half(char *buf, uint32_t val)
{
// Move all but the first digit to the right of the decimal point.
float tmp = val / 10000.0;
for(size_t i = 0; i < 5; i++)
{
// Extract the next digit.
int digit = (int) tmp;
// Convert to a character.
buf[i] = '0' + (char) digit;
// Remove the lead digit and shift left 1 decimal place.
tmp = (tmp - digit) * 10.0;
}
}
我们将使用 4.28 定点数学,而不是使用浮点,因为在我们的例子中它要快得多。也就是说,我们将二进制小数点固定在第 28 位位置,这样 1.0 就表示为 2^28。要转换为定点,我们只需乘以 2^28 即可。我们可以通过使用 0xf0000000 掩码轻松地向下舍入到最接近的整数,并且可以通过使用 0x0fffffff 掩码来提取小数部分。
(注:Terje 的算法在定点格式的选择上略有不同。)
所以现在我们有:
typedef uint32_t fix4_28;
// 0 <= val <= 99999
void itoa_half(char *buf, uint32_t val)
{
// Convert `val` to fixed-point and divide by 10000 in a single step.
// N.B. we would overflow a uint32_t if not for the parentheses.
fix4_28 tmp = val * ((1 << 28) / 10000);
for(size_t i = 0; i < 5; i++)
{
int digit = (int)(tmp >> 28);
buf[i] = '0' + (char) digit;
tmp = (tmp & 0x0fffffff) * 10;
}
}
此代码的唯一问题是 2^28 / 10000 = 26843.5456,它被截断为 26843。这会导致某些值不准确。例如,itoa_half(buf, 83492) 生成字符串“83490”。如果我们在转换为 4.28 定点时应用一个小修正,那么该算法适用于域 [0, 99999] 中的所有数字:
// 0 <= val <= 99999
void itoa_half(char *buf, uint32_t val)
{
fix4_28 const f1_10000 = (1 << 28) / 10000;
// 2^28 / 10000 is 26843.5456, but 26843.75 is sufficiently close.
fix4_28 tmp = val * ((f1_10000 + 1) - (val / 4);
for(size_t i = 0; i < 5; i++)
{
int digit = (int)(tmp >> 28);
buf[i] = '0' + (char) digit;
tmp = (tmp & 0x0fffffff) * 10;
}
}
Terje 将 itoa_half 部分交错用于低半部和高半部:
void itoa(char *buf, uint32_t val)
{
fix4_28 const f1_10000 = (1 << 28) / 10000;
fix4_28 tmplo, tmphi;
lo = val % 100000;
hi = val / 100000;
tmplo = lo * (f1_10000 + 1) - (lo / 4);
tmphi = hi * (f1_10000 + 1) - (hi / 4);
for(size_t i = 0; i < 5; i++)
{
buf[i + 0] = '0' + (char)(tmphi >> 28);
buf[i + 5] = '0' + (char)(tmplo >> 28);
tmphi = (tmphi & 0x0fffffff) * 10;
tmplo = (tmplo & 0x0fffffff) * 10;
}
}
如果循环完全展开,还有一个额外的技巧可以使代码稍微快一些。乘以 10 被实现为 LEA+SHL 或 LEA+ADD 序列。我们可以通过乘以 5 来节省 1 条指令,这仅需要一个 LEA。这与每次通过循环将 tmphi 和 tmplo 右移 1 个位置具有相同的效果,但我们可以通过调整移位计数和掩码来进行补偿,如下所示:
uint32_t mask = 0x0fffffff;
uint32_t shift = 28;
for(size_t i = 0; i < 5; i++)
{
buf[i + 0] = '0' + (char)(tmphi >> shift);
buf[i + 5] = '0' + (char)(tmplo >> shift);
tmphi = (tmphi & mask) * 5;
tmplo = (tmplo & mask) * 5;
mask >>= 1;
shift--;
}
这仅在循环完全展开时才有用,因为您可以预先计算每次迭代的移位和掩码值。
最后,该例程产生补零结果。如果 val == 0,您可以通过返回指向非 0 的第一个字符或最后一个字符的指针来消除填充:
char *itoa_unpadded(char *buf, uint32_t val)
{
char *p;
itoa(buf, val);
p = buf;
// Note: will break on GCC, but you can work around it by using memcpy() to dereference p.
if (*((uint64_t *) p) == 0x3030303030303030)
p += 8;
if (*((uint32_t *) p) == 0x30303030)
p += 4;
if (*((uint16_t *) p) == 0x3030)
p += 2;
if (*((uint8_t *) p) == 0x30)
p += 1;
return min(p, &buf[15]);
}
还有一个适用于 64 位(即 AMD64)代码的额外技巧。额外、更宽的寄存器可以有效地在寄存器中累加每个 5 位数字组;计算出最后一位数字后,您可以将它们与 SHRD 粉碎在一起,或与 0x3030303030303030 一起粉碎,然后存储到内存中。这使我的性能提高了约 12.3%。
矢量化
我们可以在 SSE 单元上按原样执行上述算法,但性能几乎没有任何提升。然而,如果我们将值分割成更小的块,我们就可以利用 SSE4.1 32 位乘法指令。我尝试了三种不同的分割:
- 2组5位数字
- 3组4位数字
- 4组3位数字
最快的变体是 4 组 3 位数字。结果见下文。
表现
除了 vitaut 和 Inge Henriksen 建议的算法之外,我还测试了 Terje 算法的许多变体。我通过对输入的详尽测试来验证每个算法的输出与 itoa() 匹配。
我的数据取自运行 Windows 7 64 位的 Westmere E5640。我以实时优先级进行基准测试并锁定到核心 0。我将每个算法执行 4 次,以将所有内容强制放入缓存中。我使用 RDTSCP 对 2^24 个调用进行计时,以消除任何动态时钟速度变化的影响。
我对 5 种不同的输入模式进行了计时:
- itoa(0 .. 9) -- 接近最佳情况的性能
- itoa(1000 .. 1999) -- 更长的输出,没有分支错误预测
- itoa(100000000 .. 999999999) -- 最长的输出,没有分支错误预测
- itoa(256 个随机值) -- 不同的输出长度
- itoa(65536 随机值) -- 不同的输出长度and冲击 L1/L2 缓存
数据:
ALG TINY MEDIUM LARGE RND256 RND64K NOTES
NULL 7 clk 7 clk 7 clk 7 clk 7 clk Benchmark overhead baseline
TERJE_C 63 clk 62 clk 63 clk 57 clk 56 clk Best C implementation of Terje's algorithm
TERJE_ASM 48 clk 48 clk 50 clk 45 clk 44 clk Naive, hand-written AMD64 version of Terje's algorithm
TERJE_SSE 41 clk 42 clk 41 clk 34 clk 35 clk SSE intrinsic version of Terje's algorithm with 1/3/3/3 digit grouping
INGE_0 12 clk 31 clk 71 clk 72 clk 72 clk Inge's first algorithm
INGE_1 20 clk 23 clk 45 clk 69 clk 96 clk Inge's second algorithm
INGE_2 18 clk 19 clk 32 clk 29 clk 36 clk Improved version of Inge's second algorithm
VITAUT_0 9 clk 16 clk 32 clk 35 clk 35 clk vitaut's algorithm
VITAUT_1 11 clk 15 clk 33 clk 31 clk 30 clk Improved version of vitaut's algorithm
LIBC 46 clk 128 clk 329 clk 339 clk 340 clk MSVCRT12 implementation
我的编译器(VS 2013 Update 4)生成了令人惊讶的糟糕代码; Terje 算法的汇编版本只是一个简单的翻译,而且速度足足快了 21%。我还对 SSE 实现的性能感到惊讶,我预计它会更慢。最令人惊讶的是 INGE_2、VITAUT_0 和 VITAUT_1 的速度有多快。值得称赞的是,vitaut 提出了一种便携式解决方案,甚至在装配级别上也超越了我的最佳努力。
注意:INGE_1 是 Inge Henriksen 的第二个算法的修改版本,因为原始算法有一个错误。
INGE_2 基于 Inge Henriksen 给出的第二种算法。它将字符串本身存储在 char[][5] 数组中,而不是在 char*[] 数组中存储指向预先计算的字符串的指针。另一个重大改进是它在输出缓冲区中存储字符的方式。它存储比所需更多的字符,并使用指针算术返回指向第一个非零字符的指针。结果明显更快——即使与 Terje 算法的 SSE 优化版本相比也具有竞争力。应该指出的是,微基准测试有点偏向于这种算法,因为在实际应用中,600K 数据集会不断地破坏缓存。
VITAUT_1 基于 vitaut 算法,有两个小变化。第一个变化是它在主循环中复制字符对,从而减少存储指令的数量。与INGE_2类似,VITAUT_1复制两个最终字符并使用指针算术返回指向字符串的指针。
执行
这里我给出了 3 个最有趣的算法的代码。
TERJE_ASM:
; char *itoa_terje_asm(char *buf<rcx>, uint32_t val<edx>)
;
; *** NOTE ***
; buf *must* be 8-byte aligned or this code will break!
itoa_terje_asm:
MOV EAX, 0xA7C5AC47
ADD RDX, 1
IMUL RAX, RDX
SHR RAX, 48 ; EAX = val / 100000
IMUL R11D, EAX, 100000
ADD EAX, 1
SUB EDX, R11D ; EDX = (val % 100000) + 1
IMUL RAX, 214748 ; RAX = (val / 100000) * 2^31 / 10000
IMUL RDX, 214748 ; RDX = (val % 100000) * 2^31 / 10000
; Extract buf[0] & buf[5]
MOV R8, RAX
MOV R9, RDX
LEA EAX, [RAX+RAX] ; RAX = (RAX * 2) & 0xFFFFFFFF
LEA EDX, [RDX+RDX] ; RDX = (RDX * 2) & 0xFFFFFFFF
LEA RAX, [RAX+RAX*4] ; RAX *= 5
LEA RDX, [RDX+RDX*4] ; RDX *= 5
SHR R8, 31 ; R8 = buf[0]
SHR R9, 31 ; R9 = buf[5]
; Extract buf[1] & buf[6]
MOV R10, RAX
MOV R11, RDX
LEA EAX, [RAX+RAX] ; RAX = (RAX * 2) & 0xFFFFFFFF
LEA EDX, [RDX+RDX] ; RDX = (RDX * 2) & 0xFFFFFFFF
LEA RAX, [RAX+RAX*4] ; RAX *= 5
LEA RDX, [RDX+RDX*4] ; RDX *= 5
SHR R10, 31 - 8
SHR R11, 31 - 8
AND R10D, 0x0000FF00 ; R10 = buf[1] << 8
AND R11D, 0x0000FF00 ; R11 = buf[6] << 8
OR R10D, R8D ; R10 = buf[0] | (buf[1] << 8)
OR R11D, R9D ; R11 = buf[5] | (buf[6] << 8)
; Extract buf[2] & buf[7]
MOV R8, RAX
MOV R9, RDX
LEA EAX, [RAX+RAX] ; RAX = (RAX * 2) & 0xFFFFFFFF
LEA EDX, [RDX+RDX] ; RDX = (RDX * 2) & 0xFFFFFFFF
LEA RAX, [RAX+RAX*4] ; RAX *= 5
LEA RDX, [RDX+RDX*4] ; RDX *= 5
SHR R8, 31 - 16
SHR R9, 31 - 16
AND R8D, 0x00FF0000 ; R8 = buf[2] << 16
AND R9D, 0x00FF0000 ; R9 = buf[7] << 16
OR R8D, R10D ; R8 = buf[0] | (buf[1] << 8) | (buf[2] << 16)
OR R9D, R11D ; R9 = buf[5] | (buf[6] << 8) | (buf[7] << 16)
; Extract buf[3], buf[4], buf[8], & buf[9]
MOV R10, RAX
MOV R11, RDX
LEA EAX, [RAX+RAX] ; RAX = (RAX * 2) & 0xFFFFFFFF
LEA EDX, [RDX+RDX] ; RDX = (RDX * 2) & 0xFFFFFFFF
LEA RAX, [RAX+RAX*4] ; RAX *= 5
LEA RDX, [RDX+RDX*4] ; RDX *= 5
SHR R10, 31 - 24
SHR R11, 31 - 24
AND R10D, 0xFF000000 ; R10 = buf[3] << 24
AND R11D, 0xFF000000 ; R11 = buf[7] << 24
AND RAX, 0x80000000 ; RAX = buf[4] << 31
AND RDX, 0x80000000 ; RDX = buf[9] << 31
OR R10D, R8D ; R10 = buf[0] | (buf[1] << 8) | (buf[2] << 16) | (buf[3] << 24)
OR R11D, R9D ; R11 = buf[5] | (buf[6] << 8) | (buf[7] << 16) | (buf[8] << 24)
LEA RAX, [R10+RAX*2] ; RAX = buf[0] | (buf[1] << 8) | (buf[2] << 16) | (buf[3] << 24) | (buf[4] << 32)
LEA RDX, [R11+RDX*2] ; RDX = buf[5] | (buf[6] << 8) | (buf[7] << 16) | (buf[8] << 24) | (buf[9] << 32)
; Compact the character strings
SHL RAX, 24 ; RAX = (buf[0] << 24) | (buf[1] << 32) | (buf[2] << 40) | (buf[3] << 48) | (buf[4] << 56)
MOV R8, 0x3030303030303030
SHRD RAX, RDX, 24 ; RAX = buf[0] | (buf[1] << 8) | (buf[2] << 16) | (buf[3] << 24) | (buf[4] << 32) | (buf[5] << 40) | (buf[6] << 48) | (buf[7] << 56)
SHR RDX, 24 ; RDX = buf[8] | (buf[9] << 8)
; Store 12 characters. The last 2 will be null bytes.
OR R8, RAX
LEA R9, [RDX+0x3030]
MOV [RCX], R8
MOV [RCX+8], R9D
; Convert RCX into a bit pointer.
SHL RCX, 3
; Scan the first 8 bytes for a non-zero character.
OR EDX, 0x00000100
TEST RAX, RAX
LEA R10, [RCX+64]
CMOVZ RAX, RDX
CMOVZ RCX, R10
; Scan the next 4 bytes for a non-zero character.
TEST EAX, EAX
LEA R10, [RCX+32]
CMOVZ RCX, R10
SHR RAX, CL ; N.B. RAX >>= (RCX % 64); this works because buf is 8-byte aligned.
; Scan the next 2 bytes for a non-zero character.
TEST AX, AX
LEA R10, [RCX+16]
CMOVZ RCX, R10
SHR EAX, CL ; N.B. RAX >>= (RCX % 32)
; Convert back to byte pointer. N.B. this works because the AMD64 virtual address space is 48-bit.
SAR RCX, 3
; Scan the last byte for a non-zero character.
TEST AL, AL
MOV RAX, RCX
LEA R10, [RCX+1]
CMOVZ RAX, R10
RETN
INGE_2:
uint8_t len100K[100000];
char str100K[100000][5];
void itoa_inge_2_init()
{
memset(str100K, '0', sizeof(str100K));
for(uint32_t i = 0; i < 100000; i++)
{
char buf[6];
itoa(i, buf, 10);
len100K[i] = strlen(buf);
memcpy(&str100K[i][5 - len100K[i]], buf, len100K[i]);
}
}
char *itoa_inge_2(char *buf, uint32_t val)
{
char *p = &buf[10];
uint32_t prevlen;
*p = '\0';
do
{
uint32_t const old = val;
uint32_t mod;
val /= 100000;
mod = old - (val * 100000);
prevlen = len100K[mod];
p -= 5;
memcpy(p, str100K[mod], 5);
}
while(val != 0);
return &p[5 - prevlen];
}
VITAUT_1:
static uint16_t const str100p[100] = {
0x3030, 0x3130, 0x3230, 0x3330, 0x3430, 0x3530, 0x3630, 0x3730, 0x3830, 0x3930,
0x3031, 0x3131, 0x3231, 0x3331, 0x3431, 0x3531, 0x3631, 0x3731, 0x3831, 0x3931,
0x3032, 0x3132, 0x3232, 0x3332, 0x3432, 0x3532, 0x3632, 0x3732, 0x3832, 0x3932,
0x3033, 0x3133, 0x3233, 0x3333, 0x3433, 0x3533, 0x3633, 0x3733, 0x3833, 0x3933,
0x3034, 0x3134, 0x3234, 0x3334, 0x3434, 0x3534, 0x3634, 0x3734, 0x3834, 0x3934,
0x3035, 0x3135, 0x3235, 0x3335, 0x3435, 0x3535, 0x3635, 0x3735, 0x3835, 0x3935,
0x3036, 0x3136, 0x3236, 0x3336, 0x3436, 0x3536, 0x3636, 0x3736, 0x3836, 0x3936,
0x3037, 0x3137, 0x3237, 0x3337, 0x3437, 0x3537, 0x3637, 0x3737, 0x3837, 0x3937,
0x3038, 0x3138, 0x3238, 0x3338, 0x3438, 0x3538, 0x3638, 0x3738, 0x3838, 0x3938,
0x3039, 0x3139, 0x3239, 0x3339, 0x3439, 0x3539, 0x3639, 0x3739, 0x3839, 0x3939, };
char *itoa_vitaut_1(char *buf, uint32_t val)
{
char *p = &buf[10];
*p = '\0';
while(val >= 100)
{
uint32_t const old = val;
p -= 2;
val /= 100;
memcpy(p, &str100p[old - (val * 100)], sizeof(uint16_t));
}
p -= 2;
memcpy(p, &str100p[val], sizeof(uint16_t));
return &p[val < 10];
}