快速正弦余弦运算

文章来源：http://blog.csdn.net/KyosukeNo1/archive/2006/05/28/759445.aspx

作者：裕作

近年来的新游戏，越来越多的采用了三维的画面――虽然比起二维游戏来，三维游戏只是多了一个深度Z坐标，但它给游戏所加入的深度、广度则是以前的游戏所远远无法比拟的。在三维空间的运算中，计算量最大的是矩阵乘法，然而，正弦余弦的计算频率则是紧跟其后的。虽然几乎任何编程语言都会自带有计算正弦余弦的函数库、类库，但在实际的使用中，它们的运算效率却往往无法达到使用者的需求――在这个时候，快速的正弦余弦查找表，则成为了程序员的救星。查找表本身并非什么神秘的东西，任何上过初中数学的朋友，都曾会拥有过那本不算太大，但却印满了密密麻麻数字的数学查找表。然而，计算机的内存是有限的，而且在精度需求不高、但运算速度却有极高需求的情况下，程序员们就需要一种高效、省内存的查找表及其算法。以下的算法适用于C＋＋和Java的快速查找表：

 

在机械运算中，加法减法比乘法除法要快，而数字的位运算，又要比加减法做得迅速，所以怎么最大限度的应用位运算来完成查找，则是本次算法的关键。首先，这个算法的出现基于以下两个理由：

 

1. 众所周知，一个圆的圆心角是360度，或者表示成2*PI，即2*3.1415926。然而，所谓的角度只是人为确定的因素而已，而且无论360还是2*PI，都并非2的次方，因此这两种方式的计算都并不迅速。

 

2. 其次，考虑到一个浮点型的数值在C＋＋和Java里都是占用32位字节，而双精度则是占用64位字节，然而，由于所有的正弦余弦结果都不可能超过1（换而言之，我们不需要整数部分），因而，32甚至64位的小数对于我们珍贵的内存空间而言，确实有点奢侈了。我们可以考虑使用16位，甚至8位的数据来表示所有的小数，在绝大多数场合，这都是已经足够的了。

 

综合上述的因素，我们可以试着把圆心角的角度设置成某个2的次方的数字，例如4096（它的16进制表达式是0x1000，在下面的部分，我们会继续讨论这个角度的取值范围问题），而把值域的小数位设成16位（即以0xC001到0x3FFF表示从-1到1之间的范围，最高位为符号位――最高位为1时是负数，最高位为0时是正数）。

 

因此，我们可以得到一个有4096个数字的查找表：表的定义域从0到4096，即对应从0到360度的圆周；表的值域从-16383到16383之间，对应从-1到1的Sin三角函数值域范围。当需要在这查找表里查找某数值时，可以使用以下这个宏：

#define QSINE( x ) (SINE_TABLE[x & 0xfff])

 

以上公式中，用0xfff与x相与的，是为了保证x的取值范围不会大于4095。

 

现在来观察正弦余弦的函数波形，由于余弦与正弦函数波形刚好相差了1/4个周期，因此只要简单的把当前的余弦值加上1/4周期（在我们的表里，是圆周的1/4，即1024）。用以下的宏，就可以查表计算余弦值：

#define QCOS( x ) (SINE_TABLE[(x + 1024) & 0xfff])

 

在实际使用时，可以代入自己所需要的角度数目相应的数值（自己预先用计算器计算好就行，可以省下将来程序里运算的时间。例如需要5度，则代入的数据是0x1000*5/360。注意，最好先乘以5再除以360，这样可以保证运算的精度）。那么，这个宏的运算结果又怎么去运用呢？大家不必急着把它转化回浮点数，事实上，运算的时候用整型的速度会更快的――因此在计算时，先把需要和该Sin值相乘的数值（一般来说，这个数值也采用整型以保证速度）和Sin查找表的结果相乘，然后把结果右移14位，结果就出来了。至于为什么是右移14位呢？其实很简单：当我们使用16位的Short类型作为Sin函数的值域时，有效的位数只有15位，最高位是用作正负号的。此外，由于剩下的取值范围需要分成正负2等份，所以取值范围需要再除以2。以下的例子通过查找表计算一个已知斜边是12345，顶角是22度（即在我们的定义域中为4096*22/360 = 250）的三角形的对边长度：

12345 * QSINE(250) >> 14 = 12345 * 0x17F1 >> 14 = 4618

 

而同样以三角库函数计算的结果是：

12345 * sin( 22 ) = 12345 * 0.3746 = 4624.437

 

由以上两公式所得结果的误差在0.139%，误差并不严重。

 

现在我们回过头来，讨论以下查找表的定义域和值域的问题。在上文中，我们把查找表的定义域规定在0到4095之间，而值域则在-16383到16383之间，即：平均一个输入值对应8个单位的输出值。由于正弦曲线上的X与Y轴的增减趋势并非正比例，因此有某些区间内会出现几个输入值对应同一输出值的情况：仔细观察刚才产生的查找表，在90度和270度处，都有近16个输入对应一个输出的情况――不过值得庆幸的是，除了这两处外，其他的输出值的都比较正常，而且16个输入所对应的角度也只是1.4度，所以对整体精度影响并不大。假如在某些对精度要求更高的场合，我们可以扩大定义域范围到原来的N倍，而值域范围扩大到32位（必须取16位，32位和64位等值，因为这是int和long类型的位数），则新的查找表大小为原来的2N倍，所以程序员必须小心控制查找表的总大小。

 

在某些编程环境中（例如J2ME JSR184的3D规范中），矩阵运算需要使用浮点数的输入，在这种情况下，使用上述的查找表反而会导致运算量的增大。因此，更为简单有效的方法是，按照所需要的精度把浮点数值直接写入到查找表中：例如设定圆周为360度，则记录360个角度各自对应的正弦值，利用递加1/4周期的特性，可以在同一张表中查出余弦值；另外，也可以为了节省空间，只记录第一象限的值，在运算过程中推出余下3个象限的值――在至于采取哪种方法比较简单有效，则需要程序员在实际的应用中自己衡量了。

我做的事情：QSINE(250)为什么等于0x17F1？

这个表大致是根据这样的原理建立的：

自变量：250/1024 = 22/90

函数值：sin(22)/sin(90) = x/16383

x = 6137.0718 = 0x17F9 跟博文中计算的相当接近了。