Skywind Inside

经过若干次试验修改，研究出下面这个快速/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的时间代价属于可以接受的范围。

下面是测试程序（点击more展开）：

#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)”。