有用吧？ 写了个比sqrtf()快一倍的sqrt()函数

最近在写界面时需要根据一幅图片中指定颜色作为透明色生成带alpha透明度的32位图片, 比较特殊的要求是透明色不是纯色，需要判断颜色相似度在一定范围内都作为透明色并根据相似度设成不同的alpha， 原来打算将RGB值转换成HSL然后判断H和L计算相似度和透明度再转回RGB, 尝试之后发觉效果并不理想而且速度慢了点. 后来把主意打到了3D距离公式上，即计算sqrt(R*R+G*G+B*B)再和阈值比较来计算实际透明度。试下来效果不错。接下来就是优化了， 显然优化首选当然是计算中的sqrt()。 在网上逛了一天， 又冥思苦想了几小时，终于写出了自己的float sqrt(float)函数，测试下来速度比CRT的sqrtf()快了整一倍，平均错误不超过0.0005, 最大错误0.0015左右。 用到程序中对图片的处理速度也提高了近75%。

个人还是比较满意这个结果的。 现在公布出来一起分享一下。

先说说原理， 实际上和著名的Quake3的invsqrt()差不多， 都是通过直接对float浮点格式进行位移操作获取近似值，然后用逼近公式迭代一次以获得满足一定精度的最终结果。

比较头疼的是迭代公式中原本存在除法, 而除法在CPU/FPU中都是相对比较慢而且并行很差的. 最后不得已我在函数中只能用很粗暴的查表法来解决这个问题了. 我实际应用中程序里保存了处理过的1.0~2.0之间9bit精度的512个倒数. 也试过只存8bit精度的倒数,不过相应地错误率也扩大了一倍--虽然不影响我的这个应用. 既然1k内存和2k内存的差别对L1 cache命中率几乎没有影响(实际使用得出的结论), 那么当然精度更高点好了.

下面的测试程序比较了3种快速sqrt()方法的差别, 分别列出计算0-0x1000000之间所有整数的sqrt后, 相对CRT的sqrtf()的速度比, 平均错误率, 最大错误率, 最大错误发生的原始数与具体偏差的值. 

其中: 
float qsqrt(float)是我写的sqrt()函数
float qsqrt2(float)是改造的Quake3的invsqrt()
float qsqrt3(float)是只位移, 不逼近, 直接获得近似结果

release版本下测试结果:

C/C++ code

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

使用 8bit 倒数表时:

 

           SPEEDUP     AVG(%)     MAX(%)        @VAL

 

via_lut():    2.00     0.0798     0.2655    4538365( 2130.3438 +5.655273)

via_inv():    1.74     0.1001     0.1752   15643569( 3955.1953 -6.930664)

via_bin():    3.13     1.5444     3.5164    4505600( 2122.6399 -74.639893)

 

使用 9bit 倒数表时：

 

           SPEEDUP     AVG(%)     MAX(%)        @VAL

 

via_lut():    2.00     0.0480     0.1648    4521982( 2126.4951 +3.504150)

via_inv():    1.74     0.1001     0.1752   15643569( 3955.1953 -6.930664)

via_bin():    3.13     1.5444     3.5164    4505600( 2122.6399 -74.639893)

生成倒数表的代码:

C/C++ code

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

void fillrecip9(DWORD* arr)

{

    union {

        int i;

        float f;

    } u;

    for( int i = 0; i < 512; ++i )

    {

        u.i = (i<<14)+0x3f800000;

        u.f = 1.0f/u.f;

        arr[i] = u.i-0x3f800000;

    }

 

}

 

void fillrecip8(DWORD* arr)

{

    union {

        int i;

        float f;

    } u;

    for( int i = 0; i < 256; ++i )

    {

        u.i = (i<<15)+0x3f800000;

        u.f = 1.0f/u.f;

        arr[i] = u.i-0x3f800000;

    }

 

}

具体测试代码：

C/C++ code

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

#include <windows.h>

#include "stdafx.h"

#include <conio.h>

#include <math.h>

 

//static DWORD recip8[256] = {

//    0x00000000, 0xffff00ff, 0xfffe03f8, 0xfffd08e5, 0xfffc0fc1, 0xfffb1885, 0xfffa232d, 0xfff92fb2,

//    0xfff83e10, 0xfff74e40, 0xfff6603e, 0xfff57404, 0xfff4898d, 0xfff3a0d5, 0xfff2b9d6, 0xfff1d48c,

//    0xfff0f0f1, 0xfff00f01, 0xffef2eb7, 0xffee500f, 0xffed7304, 0xffec9791, 0xffebbdb3, 0xffeae564,

//    0xffea0ea1, 0xffe93965, 0xffe865ac, 0xffe79373, 0xffe6c2b4, 0xffe5f36d, 0xffe52598, 0xffe45933,

//    0xffe38e39, 0xffe2c4a7, 0xffe1fc78, 0xffe135aa, 0xffe07038, 0xffdfac1f, 0xffdee95c, 0xffde27eb,

//    0xffdd67c9, 0xffdca8f1, 0xffdbeb62, 0xffdb2f17, 0xffda740e, 0xffd9ba42, 0xffd901b2, 0xffd84a5a,

//    0xffd79436, 0xffd6df44, 0xffd62b81, 0xffd578e9, 0xffd4c77b, 0xffd41733, 0xffd3680d, 0xffd2ba08,

//    0xffd20d21, 0xffd16154, 0xffd0b6a0, 0xffd00d01, 0xffcf6475, 0xffcebcf9, 0xffce168a, 0xffcd7127,

//    0xffcccccd, 0xffcc2978, 0xffcb8728, 0xffcae5d8, 0xffca4588, 0xffc9a634, 0xffc907da, 0xffc86a79,

//    0xffc7ce0c, 0xffc73294, 0xffc6980c, 0xffc5fe74, 0xffc565c8, 0xffc4ce08, 0xffc43730, 0xffc3a13e,

//    0xffc30c31, 0xffc27806, 0xffc1e4bc, 0xffc15250, 0xffc0c0c1, 0xffc0300c, 0xffbfa030, 0xffbf112b,

//    0xffbe82fa, 0xffbdf59d, 0xffbd6910, 0xffbcdd53, 0xffbc5264, 0xffbbc841, 0xffbb3ee7, 0xffbab656,

//    0xffba2e8c, 0xffb9a786, 0xffb92144, 0xffb89bc3, 0xffb81703, 0xffb79301, 0xffb70fbb, 0xffb68d31,

//    0xffb60b61, 0xffb58a48, 0xffb509e7, 0xffb48a3a, 0xffb40b41, 0xffb38cfa, 0xffb30f63, 0xffb2927c,

//    0xffb21643, 0xffb19ab6, 0xffb11fd4, 0xffb0a59b, 0xffb02c0b, 0xffafb322, 0xffaf3ade, 0xffaec33e,

//    0xffae4c41, 0xffadd5e6, 0xffad602b, 0xffaceb10, 0xffac7692, 0xffac02b0, 0xffab8f6a, 0xffab1cbe,

//    0xffaaaaab, 0xffaa392f, 0xffa9c84a, 0xffa957fb, 0xffa8e83f, 0xffa87917, 0xffa80a81, 0xffa79c7b,

//    0xffa72f05, 0xffa6c21e, 0xffa655c4, 0xffa5e9f7, 0xffa57eb5, 0xffa513fd, 0xffa4a9cf, 0xffa44029,

//    0xffa3d70a, 0xffa36e72, 0xffa3065e, 0xffa29ecf, 0xffa237c3, 0xffa1d13a, 0xffa16b31, 0xffa105a9,

//    0xffa0a0a1, 0xffa03c17, 0xff9fd80a, 0xff9f747a, 0xff9f1166, 0xff9eaecd, 0xff9e4cad, 0xff9deb07,

//    0xff9d89d9, 0xff9d2922, 0xff9cc8e1, 0xff9c6917, 0xff9c09c1, 0xff9baadf, 0xff9b4c70, 0xff9aee73,

//    0xff9a90e8, 0xff9a33cd, 0xff99d723, 0xff997ae7, 0xff991f1a, 0xff98c3bb, 0xff9868c8, 0xff980e41,

//    0xff97b426, 0xff975a75, 0xff97012e, 0xff96a850, 0xff964fda, 0xff95f7cc, 0xff95a025, 0xff9548e5,

//    0xff94f209, 0xff949b93, 0xff944581, 0xff93efd2, 0xff939a86, 0xff93459c, 0xff92f114, 0xff929cec,

//    0xff924925, 0xff91f5bd, 0xff91a2b4, 0xff915009, 0xff90fdbc, 0xff90abcc, 0xff905a38, 0xff900901,

//    0xff8fb824, 0xff8f67a2, 0xff8f177a, 0xff8ec7ab, 0xff8e7835, 0xff8e2918, 0xff8dda52, 0xff8d8be3,

//    0xff8d3dcb, 0xff8cf009, 0xff8ca29c, 0xff8c5584, 0xff8c08c1, 0xff8bbc51, 0xff8b7034, 0xff8b246b,

//    0xff8ad8f3, 0xff8a8dcd, 0xff8a42f8, 0xff89f874, 0xff89ae41, 0xff89645c, 0xff891ac7, 0xff88d181,

//    0xff888889, 0xff883fde, 0xff87f781, 0xff87af70, 0xff8767ab, 0xff872033, 0xff86d905, 0xff869223,

//    0xff864b8a, 0xff86053c, 0xff85bf37, 0xff85797c, 0xff853408, 0xff84eedd, 0xff84a9fa, 0xff84655e,

//    0xff842108, 0xff83dcf9, 0xff839930, 0xff8355ad, 0xff83126f, 0xff82cf75, 0xff828cc0, 0xff824a4e,

//    0xff820821, 0xff81c636, 0xff81848e, 0xff814328, 0xff810204, 0xff80c122, 0xff808081, 0xff804020

//};

 

static DWORD recip9[512] = {

 

... 

 

(我KAO!)

 

 提示信息

 

帖子内容过长