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

最近在搞一个设备在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位整数的四则运算,具体算法网友可自行研究。
Logo

更多推荐