【原创】快速除以255的方法
经过若干次试验修改,研究出下面这个快速/255的宏,可以在 X属于[0,65536]的范围内误差为零:
#define div_255_fast(x) (((x) + (((x) + 257) >> 8)) >> 8)
传统来说,人们习惯于将 /255改为 >> 8,但是这样误差挺大的,比如先乘以255再除以255,连续做十次,如果用>>8来代替除法,那么十次之后,误差为10. 另外一种常见的近似法是((x) + 255) >> 8,这种累积误差也挺厉害的。 因此>>8代替/255结果是比较粗糙的。而这个宏的开销比起>>8来说成本大12%。 经过测试65536000次计算中,使用/255的时间是325ms, 使用div_255_fast的时间是70ms,使用>>8的时间是62ms,div_255_fast的时间代价属于可以接受的范围。 下面是测试程序
#include <stdio.h>
#include <windows.h>
#include <mmsystem.h>
#define div_255_accurate(x) ((x) / 255)
#define div_255_approximate(x) ((x) >> 8)
// x >= 0 && x <= 65536
#define div_255_fast(x) (((x) + (((x) + 257) >> 8)) >> 8)
// 验证快速除以255的正确性
void test1(void)
{
int error = 0;
int i;
for (i = 0; i <= 65536; i++) {
if (i / 255 != div_255_fast(i)) error++;
}
printf("div_255_fast(x) ERROR = %d !!\n\n", error);
}
int A[0x10000], B[0x10000], C[0x10000], D[0x10000];
#define N 1000
// 评测三种方法的速度
void test2(void)
{
DWORD t1, t2, t3;
int i, k;
printf("benchmark:\n");
for (i = 0; i < 0x10000; i++) A[i] = rand();
// test: B = A / 255
Sleep(100);
t1 = timeGetTime();
for (k = 0; k < N; k++) {
for (i = 0; i < 0x10000; i += 4) {
B[i + 0] = div_255_accurate(A[i + 0]); // U
B[i + 1] = div_255_accurate(A[i + 1]); // V
B[i + 2] = div_255_accurate(A[i + 2]); // U
B[i + 3] = div_255_accurate(A[i + 3]); // V
}
}
t1 = timeGetTime() - t1;
printf("div_255_accurate: %dms\n", (int)t1);
// test: B = (A + ((A + 1) >> 8)) >> 8
Sleep(100);
t2 = timeGetTime();
for (k = 0; k < N; k++) {
for (i = 0; i < 0x10000; i += 4) {
B[i + 0] = div_255_fast(A[i + 0]); // U
B[i + 1] = div_255_fast(A[i + 1]); // V
B[i + 2] = div_255_fast(A[i + 2]); // U
B[i + 3] = div_255_fast(A[i + 3]); // V
}
}
t2 = timeGetTime() - t2;
printf("div_255_fast: %dms\n", (int)t2);
// test: B = A >> 8
Sleep(100);
t3 = timeGetTime();
for (k = 0; k < N; k++) {
for (i = 0; i < 0x10000; i += 4) {
B[i + 0] = div_255_approximate(A[i + 0]); // U
B[i + 1] = div_255_approximate(A[i + 1]); // V
B[i + 2] = div_255_approximate(A[i + 2]); // U
B[i + 3] = div_255_approximate(A[i + 3]); // V
}
}
t3 = timeGetTime() - t3;
printf("div_255_approximate: %dms\n", (int)t3);
}
int main(void)
{
test1();
test2();
return 0;
}
/*
OUTPUT:
div_255_fast(x) ERROR = 0 !!
benchmark:
div_255_accurate: 325ms
div_255_fast: 70ms
div_255_approximate: 62ms
*/
PS: 后来我在某网站看到另外一个快速做法:(((x) * 257 + 257) >> 16),也是在0-65536里无误差的,虽然 *257
可以化解成移位和加法,但是感觉它的实现不太适合SIMD,因为他乘以257把原来16位的数字扩大到了24位,占用太多位置了,而我提供的这个方法,确可以很好的用一个分段的长整数一次性做n个/255,并且每段诗句只占用16位,适合SIMD。
最后,帖一下 SSE2 的版本:
// (x + ((x + 257) >> 8)) >> 8
static inline __m128i _mm_fast_div_255_epu16(__m128i x)
{
return _mm_srli_epi16(_mm_adds_epu16(x,
_mm_srli_epi16(_mm_adds_epu16(x, _mask_8x0101), 8)), 8);
}
其中全局常量 _mask_8x0101就是 257这个常量,初始化为:“_mm_set1_epi16(0x0101)”。