经过若干次试验修改,研究出下面这个快速/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)”。