WV.28-大数阶乘算法8-入门篇之三汇编的威力

入门篇之三汇编的威力

在上一篇文章《大数阶乘之计算从入门到精通-入门篇之二》中，我们给出两个计算大数阶乘的程序，其中第2个程序由于用到64位整数的除法，速度反而更慢。在本文中，我们采用在C语言中嵌入汇编代码的方法，改进瓶颈部分，使计算速度提高到原先3倍多。

我们首先看一下计算大数阶乘的核心代码(见下)，可以看到，在每次循环中，需要计算1次64位的乘法，1次64位的加法，2次64位的除法(求余视为除法)。我们知道，除法指令是慢于乘法指令的，对于64位的除法更是如此。在vc中，两个64位数的除法是调aulldiv函数来实现的，两个64位数的求余运算是调用aullrem来实现的。通过分析这两个函数的源代码(汇编代码)得知，在做除法运算时，首先检查除数，如果它小于2的32次方,则需两次除法以得到一个64位商，如果除数大于等于2³²，运算将更为复杂。同样在做求余运算时，如果除数小于2³²，为了得到一个32位的余数，也需2次除法。在我们这个例子中，除数为1000000000，小于2³²，因此，这段代码需4次除法。

for (carry=0,p=pTail;p>=pHead;p--)
{
     prod=(UINT64)(*p) * (UINT64)i +carry;
     *p=(DWORD)(prod % TEN9);
         carry=prod / TEN9;
}

?

我们注意到，在这段代码中，在进行除法运算时，其商总是小于2的32次方,因此，这个64位数的除法和求余运算可以用一条除法指令来实现。下面我们用C函数中嵌入汇编代码的方法，执行一次除法指令同时得到商和余数。函数原形为：DWORD Div_TEN9_2 (UINT64 x,DWORD *pRemainder );这个函数返加x 除以1000000000的商，除数存入pRemainder。下面是函数Div_TEN9_2和计算阶乘的代码，用C中嵌入式汇编实现。关于C代码中嵌入汇编的用法，详见MSDN。

?

inline DWORD Div_TEN9_2(UINT64 x,DWORD *pRemainder )
{
         _asm
         {
                  mov eax,dword ptr [x]   // low DWORD
                  mov edx,dword ptr [x+4] // high DWORD
                 
                  mov ebx,TEN9
                  div ebx
                  mov ebx,pRemainder
                  mov [ebx],edx                 // remainder-> *remainder
                                                     // eax,return value
         }
}
void calcFac2(DWORD n)
{
     DWORD *buff,*pHead,*pTail,*p;
     DWORD t,i,len,carry;
     UINT64 prod;
   
     if (n==0)
     { printf("%d!=1",n); return; }
 //---计算并分配所需的存储空间
     t=GetTickCount();
        
     len=calcResultLen(n,TEN9);
     buff=(DWORD*)malloc( sizeof(DWORD)*len);
     if (buff==NULL)
             return ;
 //以下代码计算n!
     pHead=pTail=buff+len-1;
     for (*pTail=1,i=2;i<=n;i++)
     {
                  for (carry=0,p=pTail;p>=pHead;p--)
                  {
                      prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry;
                      carry=Div_TEN9_2(prod,p );
                  }
                 
                  while (carry>0)
                  {
                          pHead--;
                          *pHead=(DWORD)(carry % TEN9);
                          carry /= TEN9;
                  }
       }
  t=GetTickCount()-t;
 
     //显示计算结果 
         printf("It take %d ms/n",t);
     printf("%d!=%d",n,*pHead);
     for (p=pHead+1;p<=pTail;p++)
             printf("%09d",*p);
     printf("/n");
 free(buff);                  //释放分配的内存
}

注意，本文题名虽为汇编的威力，用汇编语言并不总是那么有效，一般情况下，将关键代码用汇编改写后，性能提升的幅度并不像上面的例子那么明显，可能至多提高30%，本程序是特例。

?

最后的改进：如果我们分析一下这几个计算阶乘的函数，就会发现，计算阶乘其实际上是一个二重循环，内循环部分计算出(i-1)! *?i,?外循环则依次计算2!,3!，4!，直到n!,假如们己计算出来r=(i-1)!,可否先算出 prod=?i*(i+1)* …m,?使得i*(i+1)* …刚好小于2^32,而i*(i+1)* …m*(m+1)则 >=2^32,再计算r * prod,如此一来，可减少外循环的次数，从而提高速度。理论和测试结果都表明，当计算30000以下的阶乘时，速度可提高1倍以上。下面给出代码。

void calcFac3(DWORD n)
{
       DWORD *buff,*p;
       DWORD t,carry;
       UINT64 prod;
   
        if (n==0)
        { printf("%d!=1",n); return; }
 //---计算并分配所需的存储空间
        t=GetTickCount();
  
        len=calcResultLen(n,TEN9);
        buff=(DWORD*)malloc( sizeof(DWORD)*len);
        if (buff==NULL)
            return ;
   //以下代码计算n!
                pHead=pTail=buff+len-1;
        for (*pTail=1,i=2;i+15<n;)
        {
              UINT64 t=i++;
               while (t<4294967296I64)
               {
                      t *= (UINT64)i;                 i++;
               }
               i--;
               t/=i;
          
               for (carry=0,p=pTail;p>=pHead;p--)
               {
                       prod=(UINT64)(*p) * t +(UINT64)carry;
                       carry=Div_TEN9_2(prod,p );
               }
          
               while (carry>0)
               {
                      pHead--;
                      *pHead=(DWORD)(carry % TEN9);
                       carry /= TEN9;
               }
        }
  
        for (;i<=n;i++)
        {
                    for (carry=0,p=pTail;p>=pHead;p--)
                    {
                             prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry;
                             carry=Div_TEN9_2(prod,p );
                    }
          
                    while (carry>0)
                    {
                             pHead--;
                             *pHead=(DWORD)(carry % TEN9);
                             carry /= TEN9;
                    }
        }
  
      printf("It take %d ms/n",GetTickCount()-t);
   //显示计算结果
     printf("%d!=%d",*pHead);
     for (p=pHead+1;p<=pTail;p++)
           printf("%09d",*p);
     printf("/n");
 free(buff);//释放分配的内存
}