在linux内核中使用整数运算模拟浮点运算计算三角函数

 作者声明：
本文系作者独创作品，作者保留所有相关权利。本文首次发表于百度空间——猎人的家园，欢迎网友引用及批评指正。如需转载，请注明出处，并连同本声明一并转载。禁止将本文全部或部分内容用于商业目的，特别鄙视某些源码下载站利用他人作品牟利的行为。

最近在搞一个设备在linux下的驱动，在驱动中需要进行浮点运算和正弦sin、余弦cos和反余弦acos的运算，用来进行设备输出数据和系统消息数据的变换。本来这些运算是在应用程序中进行的，然后再将运算结果反馈给系统。为了提高处理效率，决定将这些变换运算放在内核中进行，以减小系统开销。

看了一些资料，得到的结论是linux内核中不能使用协处理器进行浮点运算（不知正确否？），如果一定要进行浮点运算，需要在建立内核时选上math－emu，使用软件模拟浮点运算。据说这样作的代价有两个：
1.用户在安装驱动时需要重建内核；
2.可能会影响到其他应用程序，使得这些应用程序在进行浮点运算时也使用math－emu，会严重将低效率。

由于没有太多的精力去研究math－emu对应用程序到底有没有影响，于是决定不用math－emu。好在我的驱动里面用到的浮点运算只是加减乘除，还有上面提到的那三个三角函数，那我们就用软件方法在我的驱动里模拟三角函数吧，还不会对其他程序造成影响。

说干就干，我们来实现它。我们的数据变换公式如下：

int d1,d2,BL  //这三个数是已知的
double g,h,I  //这三个数是中间结果
int u,v       //这两个数是最终结果

g=d1*0.5667
h=d2*0.5667
I=acos((h*h+BL*BL-g*g)/(2*h*BL))
u=h*sin(2.3562-I)
v=h*cos(2.3562-I)+989.9

在应用程序中，照上面这么写就可以了。但是在内核中，是不可以的，因为用到了浮点四则运算和三角函数。我们必须用整数运算来进行上述浮点运算。

进行整数运算，首先需要考虑的是处理器的字长，它决定了CPU可以直接进行的整数运算的操作数和结果的大小。我们通常的PC机CPU都是32位的，可以用来表示0～2^32-1范围的无符号整数，也可以用来表示－2^31~2^31-1范围的有符号整数，正好对应C语言的32位整数（int或者long）。如果在运算过程中操作数、中间结果或者最终结果超出了上述的范围，将会产生溢出，从而导致结果不正确。因此，我们要使用来代替浮点数的整数小一些，使得操作数、中间结果和最终结果都不超出上述范围。

要考虑的第二个问题是运算精度。我们知道，要使运算精度高，就要增加有效数字的位数，而对于整数来说，增加有效数字就意味着数值变大，太大了就会产生溢出。这与第一个问题是矛盾的。

现在就要研究在不溢出的情况下，我们能做到的最大精度是多少。

先来说说我们的三个已知数d1、d2和BL，d1和d2的最大值不会超过30000，BL大约在2000左右，都不超过2^15＝32768。而我们要进行的是有符号运算，也就是说我们要保证运算的结果绝对值小于2^31。对于g和h，由于比d1和d2小，计算h*h和g*g的时候也不会溢出，BL*BL肯定也不会溢出。但是h*h+BL*BL-g*g会不会溢出呢？由于h和g都小于2^15，所以h*h+g*g<2^29+2^29=2*2^29=2^30，不会溢出。但是如果给d1和d2再增加一位有效数字呢？结论是会溢出，网友可自行推算。

这样我们就得出了结论：在不发生溢出的情况下，g、h和acos括号里面的东西可以有5位有效数字，不能再多了。

基于同样的理由，sin、cos、acos这三个函数也只能保留5位有效数字，后面还会讲到这三个函数我们也只能做到5位有效数字。我们知道sin、cos的值域是[-1,1],acos的值域是[0,3.1415926...]，这意味着我们可以将这三个函数的自变量和函数值都扩大到10000倍，然后用整数来表示，比如15078表示1.5078。

而对于g=d1*0.5667这样的运算，我们可以用

g=d1*5667/10000

代替。

为了避免在运算过程中产生溢出，需要注意的是，在多个数进行乘除混合运算时，要乘除交替进行，绝对不能先连乘再连除，或者先连除在连乘。安排乘除运算的原则是，只要预计到不溢出就可以乘，否则就算除法。好在乘除法的次序可以交换，这个小学就学过的内容可以保证结果的正确性。同样的，对于加减法也是同样的安排原则。

有了这些基础，我们就可以用整数运算代替浮点运算了，只不过还有一个问题需要注意，即运算过程中产生的舍入误差，其实只有舍，没有入。这个误差会减少有效数字的位数，严重降低精度。为了尽量避免或者减小这种误差，对四则运算的次序需要作仔细的安排。安排运算次序的原则是：
1.尽量避免两个接近的数相减；
2.尽量避免用很大的数去除很小的数。

根据上面提到的原则，我们可以很仔细的安排四则运算，以保证精度，而需要的数学基础也就是加法交换律、乘法交换律和分配律。

用整数运算模拟三角函数就需要点高等数学的基础了，先来复习一下。

sin、cos、arcsin、arccos可以展开为幂级数，表达式如下：

sin(x)=x-x^3/3!+x^5/5!-x^7/7!+...

cos(x)=1-x^2/2!+x^4/4!-x^6/6!+...

arcsin(x)=x+(1/2)*x^3/3+((1*3)/(2*4))*x^5/5+((1*3*5)/(2*4*6))*x^7/7+...

arccos(x)=pi/2-arccos(x)

对于arcsin和arccos要求-1<=x<=1。

这就是我们模拟三角函数的数学基础，即把三角函数变成四则运算，其实数学家们计算数学用表的时候就是这么算的。

现在我们来看看，达到要求精度的运算量，以sin为例。sin函数的幂级数拉格朗日余项

Rn(x)<=x^(n+1)/(n+1)!

　　我们要求5位有效数字，也就是小数点后4为小数，允许万分之一的误差。在x不大于pi/2（约等于1.6）的情况下，n+1=10就可以了。也就是说我们只需要计算到x^9/9!这项就可以了。这说明sin的幂级数收敛还是很快的，cos也是一样。对于x大于pi/2的情况可以通过三角函数性质变换为x小于pi/2的情形。

再看看arcsin，这个函数的幂级数分子是n,而不是n的阶乘(不要看那个2*4*6...增长得很快，但那个分子1*3*5...增长也很快)，因此这个收敛速度比较慢，我没有仔细计算，只是测试了一下，要一直算到n等于289才能达到要求的精度。加快收敛的算法我还没仔细研究，以后有时间再说吧。

公式已经有了，现在我们来研究怎么扩大10000倍。

设y=sin(x),其中x在-pi/2和pi/2之间，因为sin是奇函数，我们只需要考虑x大于0的情形，即0<=x<=1.5708,此时0<=y<=1.

我们令BY＝10000y，BX＝10000x，则有

x=BX/10000
y=BY/10000

代如sin的幂级数表达式，得到

BY=10000*sin(x)=10000*sin(BX/10000)
=10000*((BX/10000)-(BX/10000)^3/3!+(BX/10000)^5/5!-(BX/10000)^7/7!+...)
=BX/10000-BX^3/(10000^3*3!)+BX^5/(10000^5*5!)-BX^7/(10000^7*7!)+...

这样我们就得到了自变量和函数值都扩大了10000倍的sin函数的近似公式，只要我们按照前面提到的原则来合理安排运算次序，就能得到精确到万分之一的正弦函数值，而这是通过整数运算得到的。通过这个近似公式我们看到，想精确到万分之一，需要将自变量括达到10000倍，如果将再精确，就要扩大10万倍或者更大倍数，这会在运算过程中产生乘法溢出，进而出错。这就是前面提到的只能做到万分之一精度的原因。

类似的，还可以得到cos的近似公式：
10000-BX^2/(10000^2*2!)+BX^4/(10000^4*4!)-BX^6/(10000^6*6!)+...

arcsin的近似公式：
BX＋1/2*BX^3/10000^2/3+1/2*3/4*BX^5/10000^4/5+1/2*3/4*5/6*BX^7/10000^6/7...

注意观察上面公式的每一项，可以发现可以利用前一项来计算当前项的数值，这样我们只要保留前一项的结果，就可以大大减少乘除法的次数。

下面给出具体实现函数名分别为i_sin,i_cos,i_asin,i_cos。在具体实现中，进行了乘积位数的判断，目的是判断是否将会溢出，以便调整动态乘除法的次序。

int Bits_of_Num(int n)    //这是计算二进制数位数的函数
{
int i=0;
while(n){
i++;
n>>=1;
}
return i;
}



int i_sin(int ix)
{
int i;
int j;
int f;
int sum;
int mkj,mk10000,mkn;

if(ix<0){
ix=-ix;
mkn=1;
}
else mkn=0;

j=1;
f=ix;
sum=f;
for(i=1;i<6;i++){
mkj=2;
mk10000=2;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
f/=(++j);
mkj--;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
f/=(++j);
mkj--;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
f/=10000;
mk10000--;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
f/=10000;
mk10000--;
}
}
}
}
f*=ix;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
if(mkj){
f/=(++j);
mkj--;
}
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
if(mkj){
f/=(++j);
mkj--;
}
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
if(mk10000){
f/=10000;
mk10000--;
}
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
if(mk10000){
f/=10000;
mk10000--;
}
}
}
}
}
f*=ix;
if(mkj){
f/=(++j);
mkj--;
}
if(mkj){
f/=(++j);
mkj--;
}
if(mk10000){
f/=10000;
mk10000--;
}
if(mk10000){
f/=10000;
mk10000--;
}
if(i%2) sum-=f;
else sum+=f;
}
if(mkn) return -sum;
return sum;
}


int i_cos(int ix)
{
int i;
int j;
int f;
int sum;
int mkj,mk10000;

if(ix<0) ix=-ix;
j=0;
f=10000;
sum=f;
for(i=1;i<6;i++){
mkj=2;
mk10000=2;
if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
f/=(++j);
mkj--;
　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 f/=(++j);
　　 　　 　　 　　 mkj--;
　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 　　 }
　　 　　 　　 　　 }
　　 　　 　　 }
　　 　　 }
　　 　　 f*=ix;
　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 if(mkj){
　　 　　 　　 　　 f/=(++j);
　　 　　 　　 　　 mkj--;
　　 　　 　　 }
　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 if(mkj){
　　 　　 　　 　　 　　 f/=(++j);
　　 　　 　　 　　 　　 mkj--;
　　 　　 　　 　　 }
　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 if(mk10000){
　　 　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 　　 }
　　 　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 　　 if(mk10000){
　　 　　 　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 　　 　　 }
　　 　　 　　 　　 　　 }
　　 　　 　　 　　 }
　　 　　 　　 }
　　 　　 }
　　 　　 f*=ix;
　　 　　 if(mkj){
　　 　　 　　 f/=(++j);
　　 　　 　　 mkj--;
　　 　　 }
　　 　　 if(mkj){
　　 　　 　　 f/=(++j);
　　 　　 　　 mkj--;
　　 　　 }
　　 　　 if(mk10000){
　　 　　 　　 f/=10000;
　　 　　 　　 mk10000--;
　　 　　 }
　　 　　 if(mk10000){
　　 　　 　　 f/=10000;
　　 　　 　　 mk10000--;
　　 　　 }
　　 　　 if(i%2) sum-=f;
　　 　　 else sum+=f;
　　 }
　　 return sum;
}

int i_asin(int ix)
{
　　 int i;
　　 int j;
　　 int k;
　　 int f;
　　 int sum;
　　 int mkk,mk10000,mkn;

　　 if(ix<0){
　　 　　 ix=-ix;
　　 　　 mkn=1;
　　 }
　　 else mkn=0;

　　 j=1;
　　 k=2;
　　 f=ix;
　　 sum=f;
　　 for(i=1;i<144;i++){
　　 　　 mkk=1;
　　 　　 mk10000=2;
　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 f/=k;
　　 　　 　　 k+=2;
　　 　　 　　 mkk--;
　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 }
　　 　　 　　 }
　　 　　 }
　　 　　 f*=ix;
　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 if(mkk){
　　 　　 　　 　　 f/=k;
　　 　　 　　 　　 k+=2;
　　 　　 　　 　　 mkk--;
　　 　　 　　 }
　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 if(mk10000){
　　 　　 　　 　 　　 　 f/=10000;
　　 　　 　　 　　 　　 mk10000--;
　　 　　 　　 　　 }
　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(ix)>31){
　　 　　 　　 　　 　　 if(mk10000){
　　 　　 　　 　 　　 　 　　 f/=10000;
　　 　　 　 　　 　 　　 　　 mk10000--;
　　 　　 　 　　 　 　　 }
　　 　　 　　 　　 }
　　 　 　　 　 }
　 　　 　 }
　 　　 　 f*=ix;
　 　　 　 if(Bits_of_Num(f)+Bits_of_Num(j)>31){
　　 　　 　　 if(mkk){
　　 　　 　　 　　 f/=k;
　　 　　 　　 　　 k+=2;
　　 　　 　　 　　 mkk--;
　　 　 　　 　 }
　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(j)>31){
　　 　　 　 　　 　 if(mk10000){
　　 　　 　 　　 　 　　 f/=10000;
　　 　　 　 　　 　 　　 mk10000--;
　　 　　 　 　　 　 }
　　 　　 　　 　　 if(Bits_of_Num(f)+Bits_of_Num(j)>31){
　　 　　 　　 　 　　 　 if(mk10000){
　　 　　 　　 　　 　　 　　 f/=10000;
　　 　　 　　 　　 　 　　 　 mk10000--;
　　 　　 　 　　 　 　　 }
　　 　　 　　 　　 }
　　 　　 　　 }
　　 　　 }
　　 　　 f*=j;
　　 　　 j+=2;
　　 　　 if(mkk){
　　 　 　　 　 f/=k;
　　 　　 　　 k+=2;
　　 　　 　　 mkk--;
　　 　　 }
　　 　　 if(mk10000){
　　 　　 　　 f/=10000;
　　 　　 　　 mk10000--;
　　 　　 }
　　 　　 if(mk10000){
　　 　　 　　 f/=10000;
　 　　 　　 　 mk10000--;
　　 　　 }
　　 　sum+=f/j;
　　 }
　　 if(mkn) return -sum;
return sum;
}

int i_acos(int ix){
return 15708-i_asin(ix);
}

前面提到的算法，是在避免处理器产生溢出的约束条件下设计的。在这个条件下，如果CPU字长是64位，采用64位整数进行运算，则可以得到更高的精度。

如果采用两个32位整数表示一个64位整数的办法，也可以进行64位整数的四则运算，具体算法网友可自行研究。