NO IMAGE

原文地址:http://blog.csdn.net/yxnk/article/details/1665052

我的感言:首先,有一個概念上的認知,即根據階乘定義而來的常規演算法,如果是long int型只能正確計算到12左右的階乘,如果用double型只能正確計算170左右的階乘,當然這些只是大概,需要結合實際平臺進行驗證。

下面是原文:

大數階乘的計算是一個有趣的話題,從中學生到大學教授,許多人都投入到這個問題的探索和研究之中,並發表了他們自己的研究成果。如果你用階乘作關鍵字在google上搜尋,會找到許多此類文章,另外,如果你使用google學術搜尋,也能找到一些計算大數階乘的學術論文。但這些文章和論文的深度有限,並沒有給出一個高速的演算法和程式。

 

我和許多對大數階乘感興趣的人一樣,很早就開始編制大數階乘的程式。從2000年開始寫第一個大數階乘程式算起,到現在大約己有6-7年的時光,期間我寫了多個版本的階乘計算器,在階乘計算器的演算法探討和程式的編寫和優化上,我花費了很大的時間和精力,品嚐了這一過程中的種種甘苦,我曾因一個新演算法的實現而帶來速度的提升而興奮,也曾因費了九牛二虎之力但速度反而不及舊的版本而懊惱,也有過因解一個bug而通宵敖夜的情形。我寫的大數階乘的一些程式碼片段散見於網際網路絡,而演算法和構想則常常縈繞在我的腦海。自以為,我對大數階乘計算器的演算法探索在深度和廣度上均居於先進水平。我常常想,應該寫一個關於大數階乘計算的系列文章,一為整理自己的勞動成果,更重要的是可以讓同行分享我的知識和經驗,也算為IT界做一點兒貢獻吧。

 我的第一個大數階乘計算器始於2000年,那年夏天,我買了一臺電腦,開始在家專心學習VC,同時寫了我的第一個VC程式,一個仿製windows介面的計算器。該計算器的特點是高精度和高速度,它可以將四則運算的結果精確到6萬位以內,將三角、對數和指數函式的結果精確到300位以內,也可以計算開方和階乘等。當時,我碰巧看到一個叫做實用計器的軟體。值得稱頌的是,該計算器的作者姜邊竟是一個高中生,他的這個計算器功能強大,獲得了各方很高的評價。該計算的功能之一是可以計算9000以內階乘的精確值,且速度很快。在佩服之餘,也激起我寫一個更好更快的程式的決心,經過幾次改進,終於使我的計算器在做階乘的精確計算時
(以9000!為例),可比實用計算器快10倍。而當精確到30多位有效數字時,可比windows自帶的計算器快上7500倍。其後的2001年1月,我在csdn上看到一個貼子,題目是“有誰可以用四行程式碼求出1000000的階乘”,我回復了這個貼子,給出一個相對簡潔的程式碼,這是我在網上公佈的第一個大數階乘的程式。這一時期,可以看作是我寫階乘計算器的第一個時期。

 

  我寫階乘計算器的第二個時期始於2003年9月,在那時我寫了一組專門計算階乘的程式,按執行速度來分,分為三個級別的版本,初級版、中級版和高階版。初級版本的演算法許多人都能想到,中級版則採用大數乘以大數的硬乘法,高階版本在計算大數乘法時引入分治法。期間在csdn社群發了兩個貼子,“擂臺賽:計算n!(階乘)的精確值,速度最快者2000分送上”和“擂臺賽:計算n!(階乘)的精確值,速度最快者2000分送上(續)”。其高階算的版本完成於2003年11月。此後,郭先強於2004年5月10日也發表了系列貼子,“擂臺:超大整數高精度快速演算法”、“擂臺:超大整數高精度快速演算法(續)”和“擂臺:超大整數高精度快速演算法(續2)”,該貼重點展示了大數階乘計算器的速度。這個貼子一經發表即引起了熱列的討論,除了我和郭先強先生外,郭雄輝也寫了同樣功能的程式,那時,大家都在持續改進自己的程式,看誰的程式更快。初期,郭先強的稍稍領先,中途郭子將apfloat的程式碼應用到階乘計算器中,使得他的程式勝出,後期(2004年8月後)在我將程式作了進一步的改進後,其速度又稍勝於他們。在這個貼子中,arya提到一個開放原始碼的程式,它的大數乘法採用FNT+CRT(快速數論變換+中國剩餘定理)。郭雄輝率先採用apflot來計算大數階乘,後來郭先強和我也參於到apfloat的學習和改進過程中。在這點上,郭先強做得非常好,他在apfloat的基礎上,成功地優化和改時演算法,並應用到大數階乘計算器上,同時他也將FNT演算法應用到他的<超大整數高精度快速演算法庫>中,並在2004.10.18正式推出V3.0.2.1版。此後,我在2004年9月9日也完成一個採用FNT演算法的版本,但卻不及郭先強的階乘計算器快。那時,郭先強的程式是我們所知的運算速度最快的階乘計算器,其速度超過久負盛名的數學軟體Mathematica和Maple,以及通用高精度演算法庫GMP。

我寫階乘計算器的第三個時間約開始於2006年,在2005年8月收到北大劉楚雄老師的一封e-mail,他提到了他寫的一個程式在計算階乘時比我們的更快。這使我非常吃驚,在詢問後得知,其核心部分使用的是ooura寫的FFT函式。ooura的FFT程式碼完全公開,是世界上執行的最快的FFT程式之一,從這點上,再次看到了我們和世界先進水平的差距。佩服之餘,我決定深入學習FFT演算法,看看能否寫出和ooura速度相當或者更快的程式,同時一個更大計劃開始形成,即寫一組包括更多演算法的階乘計算器,包括使用FFT演算法的終極版和使用無窮級數的stirling公式來計算部分精度的極速版,除此之外,我將重寫和優化以前的版本,力爭使速度更快,程式碼更優。這一計劃的進展並不快,曾一度停止。

  目前,csdn上blog數量正在迅速地增加,我也萌生了寫blog的計劃,藉此機會,對大數階乘之計算作一個整理,用文字和程式碼詳述我的各個版本的演算法和實現,同時也可能分析一些我在網上看到的別人寫的程式,當然在這一過程中,我會繼續編寫未完成的版本或改寫以前己經實現的版本,爭取使我公開的第一份程式碼都是精品,這一過程可能是漫長的,但是我會盡力做下去。

 菜鳥篇

程式1,一個最直接的計算階乘的程式

#include "stdio.h"

#include "stdlib.h"

 

int main(int argc, char* argv[])

{

         long i,n,p;

         printf("n=?");

         scanf("%d",&n);

         p=1;

         for (i=1;i<=n;i++)

                 p*=i;

         printf("%d!=%d/n",n,p);

        return 0;

}

 

程式2,稍微複雜了一些,使用了遞迴,一個c++初學者寫的程式

#include   <iostream.h>  
  long   int   fac(int   n);  
  void   main()  
  {  
          int   n;  
          cout<<"input   a   positive   integer:";  
          cin>>n;  
          long   fa=fac(n);  
          cout<<n<<"!   ="<<fa<<endl;  
  }  
  long   int   fac(int   n)  
  {  
          long   int   p;  
          if(n==0)   p=1;  
          else  
              p=n*fac(n-1);  
          return   p;  
  }  
  

程式點評,這兩個程式在計算12以內的數是正確,但當n>12,程式的計算結果就完全錯誤了,單從演算法上講,程式並沒有錯,可是這個程式到底錯在什麼地方呢?看來程式作者並沒有意識到,一個long型整數能夠表示的範圍是很有限的。當n>=13時,計算結果溢位,在C語言,整數相乘時發生溢位時不會產生任何異常,也不會給出任何警告。既然整數的範圍有限,那麼能否用範圍更大的資料型別來做運算呢?這個主意是不錯,那麼到底選擇那種資料型別呢?有人想到了double型別,將程式1中long型換成double型別,結果如下:

#include "stdio.h"

#include "stdlib.h"

 

int main(int argc, char* argv[])

{

   double i,n,p;

   printf("n=?");

   scanf("%lf",&n);

   p=1.0;

   for (i=1;i<=n;i++)

           p*=i;

   printf("%lf!=%.16g/n",n,p);

  return 0;

}

 

執行這個程式,將運算結果並和windows計算器對比後發現,當於在170以內時,結果在誤差範圍內是正確。但當N>=171,結果就不能正確顯示了。這是為什麼呢?和程式1類似,資料發生了溢位,即運算結果超出的資料型別能夠表示的範圍。看來C語言提供的資料型別不能滿足計算大數階乘的需要,為此只有兩個辦法。1.找一個能表示和處理大數的運算的類庫。2.自己實現大數的儲存和運算問題。方法1不在本文的討論的範圍內。本系列的後續文章將圍繞方法2來展開。

大數的表示

 

1.大數,這裡提到的大數指有效數字非常多的數,它可能包含少則幾十、幾百位十進位制數,多則幾百萬或者更多位十進位制數。有效數字這麼多的數只具有數學意義,在現實生活中,並不需要這麼高的精度,比如銀河系的直徑有10萬光年,如果用原子核的直徑來度量,31位十進位制數就可使得誤差不超過一個原子核。

 

2.大數的表示:

 2.1定點數和浮點數

 我們知道,在計算機中,數是存貯在記憶體(RAM)中的。在記憶體中儲存一個數有兩類格式,定點數和浮點數。定點數可以精確地表示一個整數,但數的範圍相對較小,如一個32位元的無符號整數可表示0-4294967295之間的數,可精確到9-10位數字(這裡的數字指10進位制數字,如無特別指出,數字一律指10進位制數字),而一個8位元組的無符號整數則能精確到19位數字。浮點數能表示更大的範圍,但精度較低。當表示的整數很大的,則可能存在誤差。一個8位元組的雙精度浮點數可表示2.22*10^-308到
1.79*10^308之間的數,可精確到15-16位數字.

2.2日常生活中的數的表示:

 對於這裡提到的大數,上文提到的兩種表示法都不能滿足需求。為此,必需設計一種表示法來儲存大數。我們以日常生活中的十進位制數為例,看看是如何表示的。如一個數N被寫成"12345",則這個數可以用一個陣列a來表示,a[0]=1,a[1]=2,a[2]=3,a[3]=4,a[4]=5,這時數N=a[4]*10^0+a[3]*10^1+a[2]*10^2+a[1]*10^3+a[0]*10^4,(10^4表示10的4次方,下同),10^i可以叫做權,在日常生活中,a[0]被稱作萬位,也說是說它的權是10000,類似的,a[1]被稱作千位,它的權是1000。

 

   2.3 大數在計算機語言表示:

  在日常生活中,我們使用的阿拉伯數字只有0-9共10個,按照書寫習慣,一個字元表示1位數字。計算機中,我們常用的最小資料儲存單位是位元組,C語言稱之為char,多個位元組可表示一個更大的儲存單位。習慣上,兩個相鄰位元組組合起來稱作一個短整數,在32位的C語言編譯器中稱之為short,彙編語語言一般記作word,4個相鄰的位元組組合起來稱為一個長整數,在32位的C語言編譯器中稱之為long,組合語言一般記作DWORD。在計算機中,按照權的不同,數的表示可分為兩種,2進位制和10進位制,嚴格說來,應該是2^k進位制和10^K進位制,前者具佔用空間少,運算速度快的優點。後者則具有容易顯示的優點。我們試舉例說明:

   例1:若一個大數用一個長為len的short型陣列A來表示,並採用權從大到小的順序依次存放,數N表示為A[0]
* 65536^(len-1)+A[1] * 65536^(len-2)+…A[len-1] * 65536^0,這時65536稱為基,其進位制2的16次方。

   例2:若一個大數用一個長為len的short型陣列A來表示並採用權從大到小的順序依次存放,數N=A[0]
* 10000^(len-1)+A[1] * 10000^(len-2)+…A[len-1] * 10000^0,這裡10000稱為基,其進製為10000,即:10^4,陣列的每個元素可表示4位數字。一般地,這時陣列的每一個元素為小於10000的數。類似的,可以用long型陣列,基為2^32=4294967296來表示一個大數;當然可以用long型組,基為1000000000來表示,這種表示法,陣列的每個元素可表示9位數字。當然,也可以用char型陣列,基為10。最後一種表示法,在新手寫的計算大數階乘程式最為常見,但計算速度卻是最慢的。使用更大的基,可以充分發揮CPU的計算能力,計算量將更少,計算速度更快,佔用的儲存空間也更少。

 2.4大尾序和小尾序,我們在書寫一個數時,總是先寫權較大的數字,後寫權較小的數字,但計算機中的數並不總是按這個的順序存放。小尾(Little Endian)就是低位位元組排放在記憶體的低端,高位位元組排放在記憶體的高階。例如對於一個4位元組的整數0x12345678,將在記憶體中按照如下順序排放,
Intel處理器大多數使用小尾(Little Endian)位元組序。

    Address[0]: 0x78

    Address[1]: 0x56

Address[2]: 0x34

    Address[3]:0x12

大尾(Big Endian)就是高位位元組排放在記憶體的低端,低位位元組排放在記憶體的高階。例如對於一個4位元組的整數0x12345678,將在記憶體中按照如下順序排放,
Motorola處理器大多數使用大尾(Big Endian)位元組序。

Address[0]: 0x12

    Address[1]: 0x34

Address[2]: 0x56

    Address[3]:0x78

  類似的,一個大數的各個元素的排列方式既可以採用低位在前的方式,也可以採用高位在前的方式,說不上那個更好,各有利弊吧。我習慣使用高位在前的方式。   

2.5不完全精度的大數表示:

 儘管以上的表示法可準確的表示一個整數,但有時可能只要求計算結果只精確到有限的幾位。如用 windows自帶的計算器計算1000的階乘時,只能得到大約32位的數字,換名話說,windows計算器的精度為32位。1000的階乘是一個整數,但我們只要它的前幾位有效數字,象windows計算器這樣,只能表示部分有效數字的表示法叫不完全精度,不完全精度不但佔用空間省,更重要的是,在只要求計算結果為有限精度的情況下,可大大減少計算量。大數的不完全精度的表示法除了需要用陣列儲存有數數字外,還需要一個數來表示第一個有效數字的權,10的階乘約等於4.023872600770937e+2567,則第一個有效數字的權是10^2567,這時我們把2567叫做階碼。在這個例子中,我們可以用一個長為16的char型陣列和一個數來表示,前者表示各位有效數字,陣列的各個元素依次為:4,0,2,3,8,7,2,6,0,0,7,7,0,9,3,7,後者表示階碼,值為2567。

 

  

2.6大數的鏈式儲存法

 如果我們搜尋大數階乘的原始碼,就會發現,有許多程式採用連結串列儲存大數。儘管這種儲存方式能夠表示大數,也不需要事先知道一個特定的數有多少位有效數字,可以在運算過程中自動擴充套件連結串列長度。但是,如果基於運算速度和記憶體的考慮,強烈不建議採用這種儲存方式,因為:

1.     這種儲存方式的記憶體利用率很低。基於大數乘法的計算和顯示,一般需要定義雙連結串列,假如我們用1個char表示1位十進位制數,則可以這樣定義連結串列的節點:

          struct _node

          {

struct _node* pre;

struct _node* next;

char n; 

};

當編譯器採用預設設定,在通常的32位編譯器,這個結構體將佔用12位元組。但這並不等於說,分配具有1000個節點的連結串列需要1000*12位元組。不要忘記,作業系統或者庫函式在從記憶體池中分配和釋放記憶體時,也需要維護一個連結串列。實驗表明,在VC編譯的程式,一個節點總的記憶體佔用量為
sizeof(struct _node) 向上取16的倍數再加8位元組。也就是說,採用這種方式表示n位十進位制數需要
n*24位元組,而採用1個char型陣列僅需要n位元組。

2採用連結串列方式表示大數的執行速度很慢.

2.1如果一個大數需要n個節點,需要呼叫n次malloc(C)或new(C++)函式,採用動態陣列則不要用呼叫這麼多次malloc.

2.2 存取陣列表示的大數比連結串列表示的大數具有更高的cache命中率。陣列的各個元素的地址是連續的,而連結串列的各個節點在記憶體中的地址是不連續的,而且具有更大的資料量。因此前者的cache的命中率高於後者,從而導致執行速度高於後者。

2.3對陣列的順序訪問也比連結串列快,如p1表示陣列當前元素的地址,則計算陣列的下一個地址時一般用p1++,而對連結串列來說則可能是p2=p2->next,毫無疑問,前者的執行速度更快。

 

近似計算之一

 

在<階乘之計算從入門到精通-菜鳥篇>中提到,使用double型數來計算階乘,當n>170,計算結果就超過double數的最大範圍而發生了溢位,故當n>170時,就不能用這個方法來計算階乘了,果真如此嗎?No,只要肯動腦筋,辦法總是有的。

通過windows計算器,我們知道,171!=1.2410180702176678234248405241031e+309,雖然這個數不能直接用double型的數來表示,但我們可以用別的方法來表示。通過觀察這個數,我們發現,這個數的表示法為科學計演算法,它用兩部分組成,一是尾數部分1.2410180702176678234248405241031,另一個指數部分309。不妨我們用兩個數來表示這個超大的數,用double型的數來表示尾數部分,用一個long型的數來表示指數部分。這會涉及兩個問題:其一是輸出,這好說,在輸出時將這兩個部分合起來就可以了。另一個就是計算部分了,這是難點所在(其實也不難)。下面我們分析一下,用什麼方法可以保證不會溢位呢?

我們考慮170!,這個數約等於7.257415e+306,可以用double型來表示,但當這個數乘以171就溢位了。我們看看這個等式:

7.257415e+306

=7.257415e+306 * 10^0 (注1)(如用兩個數來表示,則尾數部分7.257415e+306,指數部分0)

=(7.257415e+306 / 10^300 )* (10^0*10^300)

=(7.257415e6)*(10
^ 300) (如用兩個數來表示,則尾數部分7.257415e+6,指數部分300)

 

依照類似的方法,在計算過程中,當尾數很大時,我們可以重新調整尾數和指數,縮小尾數,同時相應地增大指數,使其表示的數的大小不變。這樣由於尾數很小,再乘以一個數就不會溢位了,下面給出完整的程式碼。

 

程式3.

#include "stdafx.h"

#include "math.h"

#define MAX_N 10000000.00          //能夠計算的最大的n值,如果你想計算更大的數對數,可將其改為更大的值

#define MAX_MANTISSA   (1e308/MAX_N)    //最大尾數 

struct bigNum

{

    double n1;     //表示尾數部分

     int n2;   //表示指數部分

  };

 

void calcFac(struct bigNum *p,int n)

{

     int i;

     doubleMAX_POW10_LOG=(floor(log10(1e308/MAX_N))); //最大尾數的常用對數的整數部分,

     double MAX_POW10=    (pow(10.00, MAX_POW10_LOG)); // 10 ^ MAX_POW10_LOG

        p->n1=1;

     p->n2=0;

     for (i=1;i<=n;i++)

     {

          if (p->n1>=MAX_MANTISSA)

          {

               p->n1 /= MAX_POW10;

               p->n2 += MAX_POW10_LOG;

          }

          p->n1 *=(double)i;

     }

}

 

void printfResult(struct bigNum *p,char buff[])

{

     while (p->n1 >=10.00 )

     {p->n1/=10.00; p->n2++;}

     sprintf(buff,"%.14fe%d",p->n1,p->n2);

}

 

int main(int argc, char* argv[])

{

     struct bigNum r;

     char buff[32];

     int n;

    

     printf("n=?");

     scanf("%d",&n);

     calcFac(&r,n);           //計算n的階乘

     printfResult(&r,buff);   //將結果轉化一個字串

        printf("%d!=%s/n",n,buff);

     return 0;

}

 

以上程式碼中的數的表示方式為:數的值等於尾數乘以 10 ^指數部分,當然我們也可以表示為:尾數乘以 2 ^指數部分,這們將會帶來這樣的好處,在調整尾數部分和指數部分時,不用除法,可以依據浮點數的格式直讀取階碼和修改階碼(上文提到的指數部分的標準稱呼),同時也可在一定程式上減少誤差。為了更好的理解下面的程式碼,有必要了解一下浮點數的格式。浮點數主要分為32bit單精度和64bit雙精度兩種。本方只討論64bit雙精度(double型)浮點數的格式,一個double型浮點數包括8個位元組(64bit),我們把最低位記作bit0,最高位記作bit63,則一個浮點數各個部分定義為:第一部分尾數:bit0至bit51,共計52bit,第二部分階碼:bit52-bit62,共計11bit,第三部分符號位:bit63,0:表示正數,1表示負數。如一個數為0.xxxx
* 2^ exp,則exp表示指數部分,範圍為-1023到1024,實際儲存時採用移碼的表示法,即將exp的值加上0x3ff,使其變為一個0到2047範圍內的一個值。函式GetExpBase2中各語句含義如下:1.“(*pWord
& 0x7fff)”,得到一個bit48-bit63這個16bit數,最高位清0。2.“>>4”,將其右移4位以清除最低位的4bit尾數,變成一個11bit的數(最高位5位為零)。3(rank-0x3ff)”,減去0x3ff還原成真實的指數部分。以下為完整的程式碼。

 

程式4:

#include "stdafx.h"

#include "math.h"

#define MAX_N 10000000.00      //能夠計算的最大的n值,如果你想計算更大的數對數,可將其改為更大的值

#define MAX_MANTISSA   (1e308/MAX_N) //最大尾數

typedef unsigned short WORD;

 

struct bigNum

{

double n1;     //表示尾數部分

int n2;   //表示階碼部分

};

short GetExpBase2(double a) //獲得 a的階碼
  {

// 按照IEEE 754浮點數格式,取得階碼,僅僅適用於Intel系列
cpu

        WORD *pWord=(WORD *)(&a)+3;

    WORD rank = ( (*pWord & 0x7fff) >>4 );

    return (short)(rank-0x3ff);

}

 double GetMantissa(double a) //獲得 a的尾數

{

// 按照IEEE 754浮點數格式,取得尾數,僅僅適用於Intel系列
cpu

       WORD *pWord=(WORD *)(&a)+3;

*pWord &= 0x800f; //清除階碼

*pWord |= 0x3ff0; //重置階碼

return a;

}

 

void calcFac(struct bigNum *p,int n)

{

int i;

p->n1=1;

p->n2=0;

for (i=1;i<=n;i++)

{

if (p->n1>=MAX_MANTISSA)//繼續相乘可能溢位,調整之

      {

           p->n2 += GetExpBase2(p->n1);

           p->n1 = GetMantissa(p->n1);

      }

      p->n1 *=(double)i;

}

}

 

void printfResult(struct bigNum *p,char buff[])

{

double logx=log10(p->n1)+ p->n2 * log10(2);//求計算結果的常用對數

int logxN=(int)(floor(logx)); //logx的整數部分

sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);//轉化為科學計演算法形式的字串

}

 

int main(int argc, char* argv[])

{

struct bigNum r;

char buff[32];

int n;

printf("n=?");

scanf("%d",&n);

calcFac(&r,n);           //計算n的階乘

printfResult(&r,buff);   //將結果轉化一個字串

printf("%d!=%s/n",n,buff);

return 0;

}

 

程式優化的威力

程式4是一個很好的程式,它可以計算到1千萬(當n更大時,p->n2可能溢位)的階乘,但是從執行速度上講,它仍有潛力可挖,在採用了兩種優化技術後,(見程式5),速度竟提高5倍,甚至超出筆者的預計。

第一種優化技術,將頻繁呼叫的函式定義成inline函式,inline是C++關鍵字,如果使用純C的編譯器,可採用MACRO來代替。

第二種優化技術,將迴圈體內的程式碼展開,由一個迴圈步只做一次乘法,改為一個迴圈步做32次乘法。

實際速度:計算10000000!,程式4需要0.11667秒,程式5只需要0.02145秒。測試環境為迅馳1.7G,256M
RAM。關於程式優化方面的內容,將在後續文章專門講述。下面是被優化後的部分程式碼,其餘程式碼不變。

 

程式5的部分程式碼:

 inline short GetExpBase2(double a) //獲得 a的階碼

{

// 按照IEEE 754浮點數格式,取得階碼,僅僅適用於Intel系列
cpu

    WORD *pWord=(WORD *)(&a)+3;

    WORD rank = ( (*pWord & 0x7fff) >>4 );

    return (short)(rank-0x3ff);

}

inline double GetMantissa(double a) //獲得 a的尾數 

{

// 按照IEEE 754浮點數格式,取得尾數,僅僅適用於Intel系列
cpu

    WORD *pWord=(WORD *)(&a)+3;

   *pWord &= 0x800f; //清除階碼

    *pWord |= 0x3ff0; //重置階碼

 

    return a;

 

}

 

void calcFac(struct bigNum *p,int n)

{

   int i,t;

   double f,max_mantissa;

   p->n1=1;p->n2=0;t=n-32;

   for (i=1;i<=t;i+=32)

   {

        p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1);

        f=(double)i;

        p->n1 *=(double)(f+0.0);      p->n1 *=(double)(f+1.0);

        p->n1 *=(double)(f+2.0);      p->n1 *=(double)(f+3.0);

        p->n1 *=(double)(f+4.0);      p->n1 *=(double)(f+5.0);

        p->n1 *=(double)(f+6.0);      p->n1 *=(double)(f+7.0);

        p->n1 *=(double)(f+8.0);      p->n1 *=(double)(f+9.0);

        p->n1 *=(double)(f+10.0);          p->n1 *=(double)(f+11.0);

        p->n1 *=(double)(f+12.0);          p->n1 *=(double)(f+13.0);

        p->n1 *=(double)(f+14.0);          p->n1 *=(double)(f+15.0);

        p->n1 *=(double)(f+16.0);          p->n1 *=(double)(f+17.0);

        p->n1 *=(double)(f+18.0);          p->n1 *=(double)(f+19.0);

        p->n1 *=(double)(f+20.0);          p->n1 *=(double)(f+21.0);

        p->n1 *=(double)(f+22.0);          p->n1 *=(double)(f+23.0);

        p->n1 *=(double)(f+24.0);          p->n1 *=(double)(f+25.0);

        p->n1 *=(double)(f+26.0);          p->n1 *=(double)(f+27.0);

        p->n1 *=(double)(f+28.0);          p->n1 *=(double)(f+29.0);

        p->n1 *=(double)(f+30.0);          p->n1 *=(double)(f+31.0);

   }

  

   for (;i<=n;i++)

   {

        p->n2 += GetExpBase2(p->n1);

        p->n1 = GetMantissa(p->n1);

        p->n1 *=(double)(i);

   }

}

 

注1:10^0,表示10的0次方

近似計算之二

 

在《階乘之計算從入門到精通-近似計算之一》中,我們採用兩個數來表示中間結果,使得計算的範圍擴大到1千萬,並可0.02秒內完成10000000!的計算。在保證接近16位有效數字的前提下,有沒有更快的演算法呢。很幸運,有一個叫做Stirling的公式,它可以快速計算出一個較大的數的階乘,而且數越大,精度越高。有http://mathworld.wolfram.com查詢Stirling’s
Series可找到2個利用斯特林公式求階乘的公式,一個是普通形式,一個是對數形式。前一個公式包含一個無窮級數和的形式,包含級數前5項的公式為n!=sqrt(2*PI*n)*(n/e)n*(1+
1/12/n+ 1/288/n2 –139/51840/n3-571/2488320/n4+…),這裡的PI為圓周率,e為自然對數的底。後一個公式也為一個無窮級數和的形式,一般稱這個級數為斯特林(Stirling)級數,包括級數前3項的n!的對數公式為:ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n
-1/360/n3 + 1/1260/n5)-…,下面我們分別用這兩個公式計算n!.

完整的程式碼如下:

#include "stdafx.h"

#include "math.h"

#define PI       3.1415926535897932384626433832795

#define E 2.7182818284590452353602874713527

struct bigNum

{

       double n;//科學計數法表示的尾數部分

       int    e;//科學計數法表示的指數部分,若一個bigNum為x,這x=n
* 10^e

};

const double a1[]=

{     1.00, 1.0/12.0, 1.0/288, -139/51840,-571.0/2488320.0 };

const double a2[]=

{     1.0/12.0, -1.0/360, 1.0/1260.0 };

void printfResult(struct bigNum *p,char buff[])

{

       while (p->n >=10.00 )

       {p->n/=10.00; p->e++;}

       sprintf(buff,"%.14fe%d",p->n,p->e);

}

// n!=sqrt(2*PI*n)*(n/e)^n*(1+ 1/12/n+ 1/288/n^2 -139/51840/n^3-571/2488320/n^4+…)

void calcFac1(struct bigNum *p,int n)

{

       double logx,

              s,           //級數的和

              item;  //級數的每一項

       int i;

       // x^n= pow(10,n * log10(x));

       logx= n* log10((double)n/E);

       p->e= (int)(logx);   p->n= pow(10.0, logx-p->e);

       p->n *= sqrt( 2* PI* (double)n);

      

      //求(1+ 1/12/n+ 1/288/n^2 -139/51840/n^3-571/2488320/n^4+…)

       for (item=1.0,s=0.0,i=0;i<sizeof(a1)/sizeof(double);i++)

       {

              s+= item * a1[i];

              item /= (double)n; //item= 1/(n^i)

       }

       p->n *=s;

}

//ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n^3 + 1/1260/n^5 -…)

void calcFac2(struct bigNum *p,int n)

{

       double logR;

       double s,//級數的和

       item;      //級數的每一項

       int i;

      

       logR=0.5*log(2.0*PI)+((double)n+0.5)*log(n)-(double)n;

      

       //s= (1/12/n -1/360/n^3 + 1/1260/n^5)

       for (item=1/(double)n,s=0.0,i=0;i<sizeof(a2)/sizeof(double);i++)

       {

              s+= item * a2[i];

              item /= (double)(n)* (double)n; //item= 1/(n^(2i+1))

       }

       logR+=s;

      

       //根據換底公式,log10(n)=ln(n)/ln(10)

       p->e = (int)( logR / log(10));

       p->n = pow(10.00, logR/log(10) – p->e);

}

int main(int argc, char* argv[])

{

       struct bigNum r;

       char buff[32];

       int n;

       printf("n=?");

       scanf("%d",&n);

      

       calcFac1(&r,n);           //用第一種方法,計算n的階乘

       printfResult(&r,buff);   //將結果轉化一個字串

       printf("%d!=%s by method 1/n",n,buff);

       calcFac2(&r,n);           //用第二種方法,計算n的階乘

       printfResult(&r,buff);   //將結果轉化一個字串

       printf("%d!=%s by method 2/n",n,buff);

       return 0;

}

程式執行時間的測量

 

演算法的好壞有好多評價指標,其中一個重要的指標是時間複雜度。如果兩個程式完成一個同樣的任務,即功能相同,處理的資料相同,那麼執行時間較短者為優。作業系統和庫函式一般都提供了對時間測量的函式,這麼函式一般都會返回一個代表當前時間的數值,通過在執行某個程式或某段程式碼之前呼叫一次時間函式來得到一個數值,在程式或者程式碼段執行之後再呼叫一次時間函式來得到另一個數值,將後者減去前者即為程式的執行時間。

 在windwos平臺(指windwow95及以後的版本,下同),常用的用於測量時間的函式或方法有三種:1.API函式GetTickCount或C函式clock,2.API函式QueryPerformanceCounter, 3:彙編指令RDSTC

1.API函式GetTickCount:

函式原形:DWORD GetTickCount(VOID);

該函式取回從電腦開機至現在的毫秒數,即每個時間單位為1毫秒。他的解析度比較低,常用在測量用時較長程式,如果你的程式用時為100毫秒以上,可以使用這個函式.另一個和GetTickCount類似的函式是clock,該函式的回的時間的單位的是CLOCKS_PER_SEC,在windows95/2000作業系統,該值是1000,也就是說,在windows平臺,這兩個函式的功能幾乎完全相同。

2.API函式QueryPerformanceCounter:

函式原形:BOOL QueryPerformanceCounter( LARGE_INTEGER *lpPerformanceCount);該函式取回當前的高分辨值performance計數器,用一個64bit數來表示。如果你的硬體不支援高分辨performance計數器,系統可能返加零。不像GetTickCount,每個計數器單位表示一個固定的時間1毫秒,為了知道程式確切的執行時間,你需要呼叫函式QueryPerformanceFrequency來得到每秒performance計數器的次數,即頻率。

3.彙編指令RDTSC:

RDTSC 指令讀取CPU內部的“時間戳(TimeStamp)",它是一個64位無符號數計數器,這條指令執行完畢後,儲存在EDX:EAX寄存中。該指令從intel奔騰CPU開始引入,一些老式的CPU不支援該指令,奔騰後期的CPU包括AMD的CPU均支援這條指令。和QueryPerformanceCounter類似,要想知道程式的確實的執行時間,必須知道CPU的頻率,即平常所說的CPU的主頻。不幸的是沒有現成的函式可以得到CPU的頻率。一種辦法可行的辦法延時一段指定時間,時間的測量可以用QueryPerformanceCounter來做,在這段時間的開始和結束呼叫RDTSC,得到其時鐘數的差值,然後除以這段時間的的秒數就可以了。

 

下面的程式碼給出使用3個函式封裝和測試程式碼,用RDTSC指令來計時的程式碼參考了一個Ticktest的原始碼,作者不詳。

getTime1,使用GetTickCount返回一個表示當前時間的值,單位秒。

getTime2,和getTime1類似,精度更高。

getTime3,返回一個64bit的一個計數器,欲轉換為秒,需除以CPU頻率。示例程式碼見函式test3.

#include "stdafx.h"

#include "windows.h"

#include "tchar.h"

double getTime1()

{

       DWORD t=GetTickCount();

       return (double)t/1000.00;

}

double getTime2()//使用高精度計時器

{     

       static LARGE_INTEGER s_freq;

       LARGE_INTEGER performanceCount;

       double t;

       if (s_freq.QuadPart==0)

       {

              if ( !QueryPerformanceFrequency( &s_freq))

                     return 0;

       }

      

       QueryPerformanceCounter( &performanceCount );

       t=(double)performanceCount.QuadPart / (double)s_freq.QuadPart;

       return t;

}

void test1()

{

       double t1,t2;

       t1=getTime1();

       Sleep(1000);

       t2=getTime1()-t1;

       printf("It take %.8f second/n",t2);

}

void test2()

{

       double t1,t2;

       t1=getTime2();

       Sleep(1000);

       t2=getTime2()-t1;

       printf("It take %.8f second/n",t2);

}

inline BOOL isNTOS() //檢測是否執行在NT作業系統

{

       typedef BOOL (WINAPI *lpfnGetVersionEx) (LPOSVERSIONINFO);

       static int bIsNT=-1;

       if (bIsNT!=1)

              return (BOOL)bIsNT;

      // Get Kernel handle

       HMODULE hKernel32 = GetModuleHandle(_T("KERNEL32.DLL"));

       if (hKernel32 == NULL)

              return FALSE;

 #ifdef _UNICODE

    lpfnGetVersionEx lpGetVersionEx = (lpfnGetVersionEx) GetProcAddress(hKernel32, _T("GetVersionExW"));

 #else

    lpfnGetVersionEx lpGetVersionEx = (lpfnGetVersionEx) GetProcAddress(hKernel32, _T("GetVersionExA"));

 #endif

 

       if (lpGetVersionEx)

       {

              OSVERSIONINFO osvi;

              memset(&osvi, 0, sizeof(OSVERSIONINFO));

              osvi.dwOSVersionInfoSize = sizeof(OSVERSIONINFO);

              if (!lpGetVersionEx(&osvi))

                     bIsNT=FALSE;

              else

                     bIsNT=(osvi.dwPlatformId == VER_PLATFORM_WIN32_NT);

       }

       else

       {

             //Since GetVersionEx is not available we known that

              //we are running on NT 3.1 as all modern versions of NT and

              //any version of Windows 95/98 support GetVersionEx

              bIsNT=TRUE;

       }

       return bIsNT;

}

inline static BOOL checkRDSTC()//檢測CPU是否支援RDSTC指令

{

       static int bHasRDSTC= -1;

       SYSTEM_INFO sys_info;

      

       if ( bHasRDSTC !=-1 )

              return (BOOL)bHasRDSTC;

       GetSystemInfo(&sys_info);

       if (sys_info.dwProcessorType==PROCESSOR_INTEL_PENTIUM)

       {

              try

              {

                     _asm

                     {

                     _emit 0x0f                           ; rdtsc    

                     _emit 0x31

                     }

              }

              catch (…)      // Check to see if the opcode is defined.

              {

                     bHasRDSTC=FALSE; return FALSE;

              }

             // Check to see if the instruction ticks accesses something that changes.

              volatile ULARGE_INTEGER ts1,ts2;

              _asm

              {

                     xor eax,eax

                     _emit 0x0f                                  ; cpuid

                    _emit 0xa2

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                    mov ts1.HighPart,edx

                     mov ts1.LowPart,eax

                     xor eax,eax

                     _emit 0x0f                                  ; cpuid

                    _emit 0xa2

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                    mov ts2.HighPart,edx

                     mov ts2.LowPart,eax

              }

             // If we return true then there’s a very good chance it’s a real RDTSC instruction!

              if (ts2.HighPart>ts1.HighPart)

                  bHasRDSTC=TRUE;

              else if (ts2.HighPart==ts1.HighPart && ts2.LowPart>ts1.LowPart)

                     bHasRDSTC=TRUE;

              else

              {

                     printf("RDTSC instruction NOT present./n");

                     bHasRDSTC=FALSE;

              }

       }

       else

              bHasRDSTC=FALSE;

       return bHasRDSTC;

}

//***********************************************

void getTime3( LARGE_INTEGER *pTime) //返加當前CPU的內部計數器

{

       if (checkRDSTC())

       {

              volatile ULARGE_INTEGER ts;

      //on NT don’t bother disabling interrupts as doing

       //so will generate a priviledge instruction exception

              if (!isNTOS())

                     _asm cli

          //—————-     

        _asm

              {

                     xor eax,eax

           //————-save rigister

                     push ebx

                     push ecx

                    

                     _emit 0x0f                                  ; cpuid – serialise the processor

                    _emit 0xa2

                    

                    //————

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                    

                     mov ts.HighPart,edx

                     mov ts.LowPart,eax

                    

                    pop ecx

                     pop ebx

              }

       //—————–

              if (!isNTOS())

                     _asm      sti

             //———      

        pTime->QuadPart=ts.QuadPart;

       }

       else

           pTime->QuadPart=0;

}

// maxDetermainTime:最大測定時間,單位毫秒,在首次呼叫該函式時,

// 將花費maxDetermineTime的時間來測定CPU頻率,以後的呼叫將直接返加靜態變數的值

double GetCPUFrequency(DWORD maxDetermineTime )

{

       static double CPU_freq;

       LARGE_INTEGER period,t1,t2;

       register LARGE_INTEGER goal,current;

      

       if (CPU_freq>1000)     //this value have been initilization

              return CPU_freq;

       if (!QueryPerformanceFrequency(&period) || !checkRDSTC())

       {

              CPU_freq=-1.00;

              return CPU_freq;

       }

       QueryPerformanceCounter(&goal);

       goal.QuadPart += period.QuadPart * maxDetermineTime/1000;

       getTime3( &t1);  //開始計時

       do    //延時maxDetermineTime毫秒

       {

              QueryPerformanceCounter(&current);

       } while(current.QuadPart<goal.QuadPart);

       getTime3(&t2);     //結束計時

      

       CPU_freq=double((t2.QuadPart-t1.QuadPart)*1000/maxDetermineTime);

      

       char buff[100];

       sprintf(buff,"Estimated the processor clock frequency =%gHz/n",CPU_freq);

    ::MessageBox(NULL,buff,"",MB_OK);

       return CPU_freq;

}

void test3()

{

       LARGE_INTEGER t,t1,t2;

       double f1,f2;

      

       GetCPUFrequency(100);//花費0.1秒時間計算CPU頻率

       f1=getTime2();

       getTime3(&t1);

       Sleep(1000);

       getTime3(&t2);

       f2=getTime2();

       t.QuadPart=t2.QuadPart-t1.QuadPart;

       printf("It take %.8f second by getTime3/n",(double)t.QuadPart/GetCPUFrequency(100));

       printf("It take %.8f second by getTime2/n",f2-f1);

}

int main(int argc, char* argv[])

{

       test1();

       test2();

       test3();

       return 0;

}

入門篇之一

 

在《大數階乘之計算從入門到精通-大數的表示》中,我們學習瞭如何表示和儲存一個大數。在這篇文章中,我們將討論如何對大數做乘法運算,並給出一個可以求出一個大整數階乘的所有有效數字的程式。

大整數的儲存和表示已經在上一篇文章做了詳細的介紹。其中最簡單的表示法是:大數用一個字元型陣列來表示,陣列的每一個元素表示一位十進位制數字,高位在前,低位在後。那麼,用這種表示法,如何做乘法運算呢?其實這個問題並不困難,可以使用模擬手算的方法來實現。回憶一下我們在小學時,是如何做一位數乘以多位數的乘法運算的。例如:2234*8。

  

 

我們將被乘數表示為一個陣列A[], a[1],a[2],a[3],a[4]分別為2,2,3,4,a[0]置為0。

 

Step1: 從被乘數的個位a[4]起,取出一個數字4.

Step2: 與乘數8相乘,其積是兩位數32,其個位數2作為結果的個位,存入a[4],十位數3存入進位c。

Step3: 取被乘數的上一位數字a[3]與乘數相乘,並加上上一步計算過程的進位C,得到27,將這個數的個位7作為結果的倒數第二位,存入a[3],十位數2存入進位c。

Step4:重複Step3,取a[i](i依次為4,3,2,1)與乘數相乘並加上c,其個位仍存入a[i], 十位數字存入c,直到i等於1為止。

Step5:將最後一步的進位c作為積的最高位a[0]。

 

這一過程其實和小學學過的多位數乘以一位數的珠算乘法一模一樣,學過珠算乘法的朋友還有印象嗎?

 

在計算大數階乘的過程中,乘數一般不是一位數字,那麼如何計算呢?我們可以稍作變通,將上次的進位加上本次的積得到數P, 將P除以10的餘數做為結果的本位,將P除以10的商作為進位。當被乘數的所有數字都和乘數相乘完畢後,將進位C放在積的最前面即可。下面給出C語言程式碼。

一個m位數乘以n位數,其結果為m+n-1,或者m+n位,所以需首先定義一個至少m+n個元素的陣列,並置前n位為0。

 

計算一個m位的被乘數乘以一個n位的整數k,積仍儲存於陣列a

 

 

void mul(unsigned char a[],unsigned long k,int m,int n)

{

     int i;

    unsigned long p;

    unsigned long c=0;

 

    for ( i=m+n-1; i>=n;i–)

    {

           p= a[i] * k +c;

           a[i]=(unsigned char)( p % 10);

           c= p / 10;

    }

  

    while (c>0)

    {

           a[i]=(unsigned char)( c % 10);

           i–;

           c /=10;

    }

}

 

int main(int argc, char* argv[])

{

    int i;

    unsigned char a[]={0,0,0,2,3,4,5};

    mul(a,678,4,3);

    i=0;

    while ( a[i]==0)

           i++;

    for (;i<4+3;i++)

        printf("%c",a[i]+’0’); //由於數a[i](0<=a[i] <=9)對應的可列印字任符為’0’到’9’,所以顯示為i+’0’

               return 0;

}

 

從上面的例子可知,在做乘法之前,必須為陣列保留足夠的空間。具體到計算n!的階乘時,必須準備一個能容納的n!的所有位數的陣列或者記憶體塊。即陣列採有靜態分配或者動態分配。前者程式碼簡潔,但只適應於n小於一個固定的值,後者靈活性強,只要有足夠的記憶體,可計算任意n的階乘,我們這裡討論後一種情況,如何分配一塊大小合適的記憶體。

n!有多少位數呢?我們給出一個近似的上限值:n! <(n+1)/2的n次方,下面是推導過程。

Caes 1: n是奇數,則中間的那個數mid= (n+1)/2,
除了這個數外,我們可以將1到n之間的數分成n/2組,每組的兩個數為
mid-i和mid+i (i=1到mid-1),如1,2,3,4,5,6,7可以分為數4,和3對數,它們是(3,5),(2,6)和(1,7),容易知道,每對數的積都於小mid*mid,故n!小於(n+1)/2的n的次方。

 

Case 2: n是個偶數,則中間的兩個數(n-1)/2和(n+1)/2,我們將(n+1)/2記做mid,則其它的幾對數是(mid-2,mid+1),(mid-3)(mid+2)等等,容易看出,n!小於mid的n次方。

由以上兩種情況可知,對於任意大於1的正整數n, n!<(n+1)/2的n次方。

如果想求出n!更準確的上限,可以使用司特林公式,參見該系列文章《階乘之計算從入門到精通-近似計算之二》。

 

到此,我們已經解決大數階乘之計算的主要難題,到了該寫出一個完整程式的時候了,下面給出一個完整的程式碼。

 

#include "stdio.h"

#include "stdlib.h"

#include "memory.h"

#include "math.h"

#include "malloc.h"

void calcFac(unsigned long n)

{

   unsigned long i,j,head,tail;

int blkLen=(int)(n*log10((n+1)/2)); //計算n!有數數字的個數

   blkLen+=4;  //保險起見,多加4位

  if (n<=1)

  {       printf("%d!=0/n",n);     return;}

 

  char *arr=(char *)malloc(blkLen);       

   if (arr==NULL)

   {       printf("alloc memory fail/n"); return ;}

 

   memset(arr,0,sizeof(char)*blkLen);

   head=tail=blkLen-1;

   arr[tail]=1;

  

   for (i=2;i<=n;i++)

   {

           unsigned long c=0;

           for (j=tail;j>=head;j–)

           {

                    unsigned long prod=arr[j] * i +c;

                    arr[j]=(char)( prod % 10);

                    c= prod / 10;

           }

           while (c>0)

           {

                    head–;   

                    arr[head]=(char)(c % 10);

                    c/=10;

           }

   }

   printf("%d!=",n);

   for (i=head;i<=tail;i++)

           printf("%c",arr[i]+’0′);

   printf("/n");

 

   free(arr);

}

 

void testCalcFac()

{

    int n;

    while (1)

    {

            printf("n=?");

            scanf("%ld",&n);

            if (n==0)       

                   break;

            calcFac(n);

    }

}

 

int main(int argc, char* argv[])

{

    testCalcFac();

    return 0;

}

用Stirling逼近近似計算階乘的探討與應用

江蘇省贛榆高階中學仲晨

[email protected]

【關鍵詞】:Stirling逼近,階乘,極限論,微積分,數學實驗,計算機演算法

 

“階乘”(factorial)在資訊學競賽中具有重要角色,更廣泛的說,“階乘”在數學領域也是佔有重要地位。在許多人剛剛學習計算機語言的時候,大多會被要求寫一個算階乘的程式,而在學習高精度演算法的時候,也會寫一個計算較大數字階乘的程式。不過,在實際的運用之中,可能遇到更大數字的階乘計算和不同要求的階乘結果,例如:TOJ(同濟大學ACM網路題庫,http://acm.tongji.edu.cn/problem.php)的1016題——“求N!左邊第二位的數字”,這就需要一定的精度思考了。

可是我們通常對於較大數字階乘的要求是求結果位數或前幾位數字,這怎麼辦呢?

劉汝佳《演算法藝術與資訊學競賽》一書中,(Page241)介紹了Stirling公式:

其中的符號是指“同階”或“相當”,即兩者隨n增加的大致速度相同,在n較大時,兩者極其相近。

這是一個極限的概念(現行教材高二下學期數學內容),屬於微分學內容,準確寫法為:

但遺憾的是在《演算法藝術與資訊學競賽》書中只提供了這個算式,並無他物!

本人近日看到一本數學科普讀物——《好玩的數學——不可思議的e(陳任政著,科學出版社),其中5.12節簡介了Stirling逼近近似算階乘,本人感到好奇,於是對這種演算法的具體步驟進行了分析,並研究了它的精確度,故為本文。在2005年8月7日完工之日,筆者上網搜尋了一下,找到了一些關於Stirling逼近的文章,偶然地在臺灣亞洲聚合公司蔡永裕《談Stirling公式的改良》(刊自臺灣《數學傳播》20卷4期,民國85年12月)一文中找到同感,蔡先生的做法於筆者方向不同,作出來的結果比筆者的演算法精確一個數量級左右,慚愧,於是,筆者又再次研究,尋找更好演算法,寫於本文後部。

 

在1730年,棣莫弗(棣,音Dì)(法國數學家,Abraham
De Moiver,1667~1754)發表的《分析雜論》中首先對n!地一個無窮級數展開式給出了近似公式:

但是,現在我們叫這個式子為“Stirling逼近”,中文叫做“斯特林逼近”,這是為什麼呢?

因為棣莫弗的朋友蘇格蘭數學家斯特林(James Stirling,1696~1770)在同年的《微分法或無窮級數的簡述》中也給出了等價的級數。

事實上,棣莫弗首先得到的式子是,
但是,他沒有把C求出來。而斯特林則利用棣莫弗的發現做了探討,求出了。

這些式子的來源是一個無窮級數展開式:

其中B2=1/6,B4=-1/30,B6=1/42
… B2k是雅格布·伯努力數。(具體內容請參見後文介紹)

 

這裡介紹一下,還沒上高中的同學還沒有學到,“乘方”的逆運算有兩種:開方和對數。

對於一個冪:,其中a成為底數,n成為指數,b成為冪。已知a和n求b,就是乘方運算;已知b和n,求a,就是開方運算;而已知a和b求n,就是對數運算,寫做:,這裡n就稱為以a為底b的對數(logarithm)。

當底數為10的時候,叫做常用底數,簡寫做e的時候,叫做自然對數,簡寫做;當底數為

至於e的含義:e是重要性僅次於π的數,是極限論的主要內容,具體的說,即:

意思是當n趨向於正的無限大的時候,趨向於e。

e是無理數,即無限不迴圈小數,也是超越數,即不能滿足某個整數係數代數方程的數(不能滿足某個整數係數代數方程的數叫做代數數)。

目前e只算到了幾千位。

e=2.718281828459045235360287471352662497757247093…

特別說明的是,在Pascal語言中,exp(n)函式就是e的n次方。

 

另外,有個著名的公式被成為“整個數學中最卓越的公式”:

其中的i為虛數的單位,。來自算術的0、1,來自代數的i,來自幾何的π,來自分析學的e,奇妙的組成了一個公式!

這是尤拉(瑞士數學家,Leonhard Euler,1707~1783)發現的!所以稱作“尤拉公式”

不過,真正的尤拉公式是:

那個“最卓越的公式”只是尤拉公式的一個推倒公式。

 

言歸正傳,由公式

兩邊同時取e的冪,得
即:

再經過近似等處理(這些處理比較麻煩,而且牽扯到微積分內容),我們得到了Stirling公式:

注:在本文後部有Stirling公式的推倒過程。

當然,我們可以得到它得更具體形式:

其中的θ是不定的,在(0,1)區間之內。

 

講到這裡,我們有了公式,可以開始計算了。

但是,我們計算什麼呢?難道我們要計算N!的值嗎?

雖然這個式子比階乘原始計算式簡單,但是實際計算的時候仍然得到的將是上百上千位的數字,而且這個式子畢竟是近似公式,無法真正得到確切得數!

難道我們辛苦的到的式子無用了嗎?

答案當然是否定的!我們可以求N!的位數!

 

求位數的方法是取常用對數(以10為底的對數)。那,這是為什麼呢?看看下錶:

n

n位數

 

1

0.000000

1

10

1.000000

2

100

2.000000

3

1000

3.000000

4

10000

4.000000

5

100000

5.000000

6

1000000

6.000000

7

10000000

7.000000

8

100000000

8.000000

9

好了!我們知道了n的位數等於,對嗎?

不對!不一定是整數,比如,所以,我們得到的公式是:

n的位數=

其中是取整符號,指不大於n的最大整數,在Pascal語言中為trunc(n)函式。

 

如果我們用Stirling逼近算出n!再算它的位數就太不值得了,這裡我們就要運用對數的一些特性來簡化公式。

我們兩邊取常用對數:

進而化簡為:

現在我們可以來求值了!

以求1000!的位數為例,

所以1000!的位數為2568位!

我們可以看到,這裡的計算複雜度很低,都是O(n)級別的!對於其計算機程式設計運用來說,計算過程中數字不大不會溢位,取對數後的精確度也可以保證。

這樣我們就解決了計算階乘位數問題!而且可以絕對放心,這個公式在位數方面的準確率是99.9999999999999%以上。(而且這個僅有的可能性也只是從正態性來推測的)

 

用Pascal語言來書寫:

 Trunc(0.5*Ln(2*n*3.1415926)
/
Ln(10)
+
n*Ln(n/Exp(1))
/
Ln(10))
+
1

建議使用分步計算,這樣能避免計算過程中的溢位現象。

 

這樣,筆者使用批處理計算了幾個值:

n

n!位數

1

1

10

7

100

158

1000

2568

10000

35660

100000

456574

1000000

5565709

10000000

65657060

100000000

756570557

1000000000

8565705523

 

由於Longint範圍的限制,無法計算再大的數值,但是可以通過使用int64或者extended來解決。

不過,我們有點不甘心,只計算位數難免用途不廣,其實,我們可以用來計算n!的前幾位!

什麼原理呢?

例如計算1000!時,我們得到的數為2567.6046080272230971441871361093945……

而我們從下表中可以看出點東西來:

n

 

5

0.698970004336019……

50

1.698970004336019……

500

2.698970004336019……

5000

3.698970004336019……

我們可以得知:一個數增大10倍,它的常用對數整數部分將增加1,而小數部分不變。

我們得到的更重要的結果是:一個數的常用對數的小數部分所對應的冪的各位數字(不考慮小數點)與原數各位數字相同。

所以,對於1000!來說:
2567.6046080272230971441871361093945……的小數部分為
0.6046080272230971441871361093945……,
它的10為底的冪為4.0235372577197005393836315765635……,
也就是說1000!這個2568位數的前幾位為40235372577197005393836315765635……。

Well!我們可以求出一個階乘的前幾位數了!

但是,不要高興太早!這裡的精度無法保證!

下面讓我們來看看精確度方面:

特別說明的是,這裡說的精確度是指n!與Stirling逼近公式的精確度,而不是位數計算的精確度,因為位數計算基本上是精確的!

n

逼近值

n!準確值

相對誤差

10

3598695.619

3628800

0.008365359

20

2.42279E+18

2.4329E+18

0.004175011

30

2.64517E+32

2.65253E+32

0.002781536

40

8.14E+47

8.16E+47

0.002085461

50

3.04E+64

3.04E+64

0.001668034

60

8.31E+81

8.32E+81

0.001389841

70

1.20E+100

1.20E+100

0.001191177

80

7.15E+118

7.16E+118

0.001042204

90

1.48E+138

1.49E+138

0.000926351

100

9.32E+157

9.33E+157

0.000833678

110

1.59E+178

1.59E+178

0.000757861

120

6.68E+198

6.69E+198

0.000694684

130

6.46E+219

6.47E+219

0.000641230

140

1.35E+241

1.35E+241

0.000595414

150

5.71E+262

5.71E+262

0.000555709

160

4.71E+284

4.71E+284

0.000520968

170

7.25E+306

7.26E+306

0.000490316

 

 

可以看到,本身它的相對誤差就很小,並且這個誤差是在逐漸減小。

當n=1000的時候,相對誤差只有0.00008333680287333986657587……!

這個誤差可以保證我們能夠取得階乘值的前3-4位準確值,但是,還是有點遺憾,感覺不是太好,我們最起碼需要7-8位的精確度,可是,我們能做到嗎?

可以!

還記得嗎?我們曾得到了一個更精確的算式:

(其中的θ是不定的,在(0,1)區間之內。)

 

從這個式子裡我們也可以看出準確值和逼近值之比就是

隨著n的增大,顯然這個值是隨n逐漸縮小的!並且是縮得很小,這也符合我們上面表格所體現的內容。

可是,這裡的θ是不固定的,要是固定的,我們就可以求出準確值。但是,有沒有想看看θ的變化趨勢呢?我們用上面算出的相對誤差反算θ

(計算公式為ln(相對誤差+1)×12×n=θ

n

θ

10

0.999667612

20

0.999916726

30

0.999962975

40

0.999979170

50

0.999986668

60

0.999990741

70

0.999993198

80

0.999994792

90

0.999995885

100

0.999996667

110

0.999997245

120

0.999997685

130

0.999998028

140

0.999998299

150

0.999998519

160

0.999998698

170

0.999998847

180

0.999998971

190

0.999999077

200

0.999999167

我們驚喜地發現θ竟然十分接近1,而且在逐漸地逼近1,這是個令人興奮地訊息!

實際上,即使是求1的階乘,θ也會達到0.9727376027,這是一個本身就是一個很“精確”的數字了!當n=1000時,θ將0.99999996665875876427498746773752,與1的差別只有0.000000033341241235725012532263(約等於3.33412×10-8)!

如果可以精確到這樣,我們就太高興了!

我們把θ用1帶入,然後求值,這次得到的結果出奇的好!

下表是經過θ=1帶入修正的各項指標:

n

修正後逼近值

n!準確值

二者比例

10

3.62881005142693E+06

3.62880000000000E+06

0.999997230103867

20

2.43290285233215E+18

2.43290200817664E+18

0.999999653025394

30

2.65252887092925E+32

2.65252859812191E+32

0.999999897151981

40

8.15915318654567E+47

8.15915283247897E+47

0.999999956604970

50

3.04140938775049E+64

3.04140932017133E+64

0.999999977780315

60

8.32098721974147E+81

8.32098711274139E+81

0.999999987140939

70

1.19785717669724E+100

1.19785716699698E+100

0.999999991901990

80

7.15694574345356E+118

7.15694570462638E+118

0.999999994574895

90

1.48571597014272E+138

1.48571596448176E+138

0.999999996189743

100

9.33262157031762E+157

9.33262154439441E+157

0.999999997222301

110

1.58824554483731E+178

1.58824554152274E+178

0.999999997913062

120

6.68950292420235E+198

6.68950291344912E+198

0.999999998392522

130

6.46685549739670E+219

6.46685548922047E+219

0.999999998735671

140

1.34620124893450E+241

1.34620124757175E+241

0.999999998987707

150

5.71338396114816E+262

5.71338395644585E+262

0.999999999176966

160

4.71472363918940E+284

4.71472363599206E+284

0.999999999321839

170

7.25741561941125E+306

7.25741561530799E+306

0.999999999434611

180

2.00896062594820E+00

2.00896062499134E+00

0.999999999523704

190

9.68032267917558E+00

9.68032267525524E+00

0.999999999595020

可以看到,這裡的差別已經很小了!

當n=1000,二者比例達到了0.999999999997221,這足以保證10位的準確度!

到這裡,我們就找到了一個精確度高並且速度快的演算法來計算階乘的前幾位!

我們求階乘前幾位的最後公式為


frac(n)為取小數部分

順便說一下,以上的資料都是用Free Pascal計算的,使用的是Extended格式,由此看來,Extended的精確度也是很高的!

用Pascal語言來書寫:

10**frac(0.5*Ln(2*n*3.1415926) / Ln(10) + n * Ln(n / Exp(1)) / Ln(10))*exp(1/12/n)

有個技巧是:在FP中,可以使用a**b來計算ab,可以不必再用exp(y*ln(x))。

 

筆者希望讀者仔細思考,如何書寫精度最高,速度最快(儘管速度已經是毫秒級了!)?還請注意資料溢位地情況。

還有,為了輸出前幾位,需要下點功夫處理一下輸出問題。筆者在這裡就簡略了。

 

別急,我們的“征途”還沒完,我們要繼續精確化!

 

下面是來自http://mathworld.wolfram.com/StirlingsSeries.html的資料,介紹Stirling’s Series,即“斯特林級數”,也就是Stirling逼近的原始式。比較抽象,讀者瀏覽一下即可。

筆者英語不佳,勉強翻譯了一下,便於大家閱讀。

Stirling’s Series

斯特林級數

Theasymptotic series for the
gamma function is given by

這個Γ函式的漸進級數如下(譯者注:Γ為γ的大寫,希臘字母,讀做“伽馬”)

The coefficient of can given explicitly by

的係數可以明確地給出

where is the number of permutations ofnwithk
permutation cycles all of which are ≥3(Comtet 1974, p. 267). Another formula for the s is given by the recurrence relation

上面的是一個以k為迴圈的n的變數,它是始終≥3的(Comtet
1974, p. 267)。另外的一個計算的關係如下

with , then

當,則

where is the
double factorial (Borwein and Corless 1999).

這裡的的意思是雙階乘(Borwein and Corless 1999).

(譯者注:把階乘的定義引申,定義N!! = N*(N-2)*(N-4)*…,若N為偶數,則乘至2,反之,則乘至1,而0!!
= 0。我們稱之為雙階乘(Double Factorial))

The series for is obtained by adding an additional factor ofx,

級數是由x+1的Γ函式獲得,

The expansion of is what is usually called Stirling’s series. It is given by the simple analytic expression

的展開式通常被叫做斯特林級數。它簡單地分析表示為,

where  is aBernoulli number. Interestingly, while the numerators in this expansion are the same as those of for the first several
hundred terms, they differ at n=574, 1185, 1240, 1269, 1376, 1906, 1910, … , with the corresponding ratios being 37, 103, 37, 59, 131, 37, 67, …

這裡地是伯努力數。有趣的是,當n每增加數百後,會出現重複,比如當n=574,
1185, 1240, 1269, 1376, 1906 , 1910 , … 時,對應的是37, 103, 37, 59, 131, 37, 67, …

 

對於以上內容,單就本文所討論的主題來說,沒有太大用途,許多人在用Stirling公式計算大數階乘的時候(注意,他們是直接計算階乘近似值,而筆者是通過常用對數反算近似值),常常不使用“Stirling逼近”而直接使用“Stirling級數展開式”,這樣主要是因為他們注意到了“Stirling逼近”簡單式的誤差過大,轉而使用10-100項的“Stirling級數展開式”。在本文中,筆者使用“Stirling逼近”的“精確式”,採用修正的方法求近似值,已經取得了不亞於使用100項的“Stirling級數展開式”的精確度,並且避免了階乘數值的溢位。

筆者在上文也說過,筆者看了臺灣蔡永裕先生的《談Stirling公式的改良》一文,感覺非常有水平,筆者有時很有疑問,臺灣地區的計算機學、數學等科學的發展水平還是比較高的!經常性的筆者搜尋數學、計算機學的內容都會搜到許多臺灣網站的內容,可能還有一個原因就是臺灣地區的計算機普及率較高,資料上網率較高。

這裡兩個網址:

臺灣“中央研究院”數學研究所(數學研究所)http://www.math.sinica.edu.tw/

《數學傳播》雜誌(《數學傳播》)http://www.math.sinica.edu.tw/media/default.jsp

讀者能看得順繁體字更好,如果看不大習慣,就用一些軟體來“轉換”一下。

再次言歸正傳,蔡永裕先生的原理與筆者基本相同,都是對Stirling逼近公式進行修正,蔡先生文中聯想到公式:

於是設e的漸近相等值E,將Stirling公式變為(筆者注:即≈)

反算得漸近於1,於是設

其中Fn為修正引數,則匯出

然後反算得資料表格。繼而欲求Fn的表示式,經過試驗,選擇了

得到了Stirling公式的修正版:

其常用對數形式為:

這樣一來,精確度提高了很多,當計算1000!時,蔡先生的演算法誤差僅有-2.971E-13,而如果使用原始Stirling式的誤差為-8.33299E-05,筆者之前的演算法誤差是1.118E-12,略遜於蔡式演算法,於是筆者思考了一下,使用二次修正的方法來實現精確化。

方法如下:

對於公式:

因為θ接近於1,則令,f(n)為修正函式。則

為了尋找f(n)的式子,從簡單的一次函式試起,我們先試一試f(n)/n,發現竟然是等差數列,在除以n後,得到的是一個常量!

n

θ

f(n)

f(n)/n

f(n)/n2

50

0.999986668189609

75008.56753

1500.171351

30.00342701

100

0.999996666761655

300008.5492

3000.085492

30.00085492

150

0.999998518536300

675008.1019

4500.054013

30.00036009

200

0.999999166673422

1200009.727

6000.048637

30.00024319

250

0.999999466670049

1875011.893

7500.047571

30.00019028

300

0.999999629630836

2700008.792

9000.029307

30.00009769

350

0.999999727877187

3674811.34

10499.46097

29.99845991

400

0.999999791661586

4799882.935

11999.70734

29.99926834

450

0.999999835405298

6075529.694

13501.1771

30.00261577

500

0.999999866704792

7502145.139

15004.29028

30.00858056

550

0.999999889796665

9074135.523

16498.42822

29.99714222

600

0.999999907419129

10801367.44

18002.27906

30.00379843

650

0.999999921104972

12675069.96

19500.10763

30.00016558

700

0.999999931990204

14703764.09

21005.37728

30.00768183

750

0.999999940767808

16882711.46

22510.28195

30.01370927

800

0.999999947926410

19203592.46

24004.49057

30.00561321

850

0.999999953865914

21675947.14

25501.11429

30.00131093

900

0.999999958850112

24301402.95

27001.55883

30.00173203

950

0.999999963096340

27097582.94

28523.77151

30.02502264

1000

0.999999966659766

29993790.67

29993.79067

29.99379067

 

顯然,它們都與30相近;而且,我們完全可以認為,這裡的誤差是由計算誤差引起的!

 

於是,我們得到f(n)=30n2關係式。

繼而,有

“前幾位”公式:

這樣一來

n

二次修正版逼近值
科學計數法系數

n!準確值
科學計數法系數

相對誤差

50

3.04140932016361

3.04140932017133

-2.5383029012E-12

100

9.33262154439367

9.33262154439441

-7.9380946261E-14

150

5.71338395644579

5.71338395644585

-1.0547118734E-14

200

7.88657867364788

7.88657867364790

-2.5535129566E-15

看!僅僅n=200的時候,誤差就小於!

不過,由於畢竟f(n)=30n2是有誤差的,所以誤差不會一直減少,而會波動,正如上面的求f(n)/n2的表格,它的值是波動的。

這是因為:我們不可能用固定的多項式來表示出對數型變化!

但我們把誤差降到了最小,當n=1000時,
逼近值4.0238726007709361277668E+2568與準確值4.0238726007709377354370E+2568的誤差僅有-3.9953306821471308041868495761274E-16

事實上,在高等數學裡,根據Taylor展式可以算得:

看似我們屬於“碰對了”,其實這就是“數學實驗”的魅力!

 

到此為止,我們可以說獲得了成功!!!

 

在程式設計方面,注意點很多,主要還是溢位的問題,比如1/(360*n*n*n)對於較大的可能溢位,而寫成1/360/n/n/n就可以避免。

exp(frac(0.5*Ln(2*n*3.14159265358979323)/Ln(10)+n*Ln(n/Exp(1))/Ln(10))*ln(10))*exp(1/12/n-1/360/n/n)

 

以上內容是筆者看書時偶爾思考到的,並認真得進行了思考和實驗,具體地試驗比較,這種做法叫做“數學實驗”。

《數學實驗》是在我國高等學校中新開設的一門課程。現在還處於試點和摸索階段,有許多不同的想法和作法。課程目的,是使學生掌握數學實驗的基本思想和方法,即不把數學看成先驗的邏輯體系,而是把它視為一門“實驗科學”,從問題出發,藉助計算機,通過學生親自設計和動手,體驗解決問題的過程,從實驗中去學習、探索和發現數學規律。既然是實驗課而不是理論課,最重要的就是要讓學生自己動手,自己藉助於計算機去“折騰”數學,在“折騰”的過程中去學習,去觀察,去探索,去發現,而不是由老師教他們多少內容。既不是由老師教理論,主要的也不是由老師去教計算機技術或教演算法。不著意追求內容的系統性、完整性。而著眼於激發學生自己動手和探索的興趣。

數學實驗可以包括兩部分主要內容:第一部分是基礎部分,圍繞高等數學的基本內容,讓學生充分利用計算機及軟體(比較有名的是Mathematica)的數值功能和圖形功能展示基本概念與結論,去體驗如何發現、總結和應用數學規律。另一部分是高階部分,以高等數學為中心向邊緣學科發散,可涉及到微分幾何,數值方法,數理統計,圖論與組合,微分方程,運籌與優化等,也可涉及到現代新興的學科和方向,如分形、混沌等。這部分的內容可以是新的,但不必強調完整性,教師介紹一點主要的思想,提出問題和任務,讓學生嘗試通過自己動手和觀察實驗結果去發現和總結其中的規律。即使總結不出來也沒有關係,留待將來再學,有興趣的可以自己去找參考書尋找答案。

筆者寫本文,一來總結自己所想,二來拋磚引玉,希望大家能在數學中尋找靈感,並優化使之適用於計算機程式設計,這才是演算法藝術

 

共11641字

寫於:2005年8月6日~8日

江蘇·贛榆

 

參考文獻:

《好玩的數學——不可思議的e》(陳任政著,科學出版社,2005);

《談Stirling公式的改良》(蔡永裕,臺灣亞洲聚合公司,刊自臺灣《數學傳播》20卷4期,民國85年12月-公元1996年12月);
pdf檔案下載:
http://www.math.sinica.edu.tw/math_media/pdf.php?m_file=ZDIwNC8yMDQwNA==

《談Stirling公式》(蔡聰明,載於數學傳播第十七卷第二期,作者當時任教於臺大數學系);
http://episte.math.ntu.edu.tw/articles/mm/mm_17_2_05/

清華大學線上學習--《組合數學》之(1.3節Stirling近似公式);
http://www.cs.cityu.edu.hk/~luoyan/mirror/tsinghua/combinemaths/1/1_3.htm

Stirling級數http://mathworld.wolfram.com/StirlingsSeries.html

入門篇之二

在《大數階乘之計算從入門到精通―入門篇之一》中,我們給出一個計算階乘的程式,它採用char型陣列存貯大數,1個元素表示1位十進位制數字,在計算時,一次乘法可計算一位數字和一個整數的乘積。該演算法具有簡單直觀的優點,但缺點也是明顯的,速度不快,佔用記憶體空間也較多,本文將給出一個改後的程式,有效的克服了這些缺點。

學過80×86彙編的人都知道,8086/8088的CPU可對兩個16位元的數相乘,其結果為32位元,80386及其後的32位CPU可對兩個32位元的數相乘,結果為64位元(以下寫作bit)。8086
CPU等16位CPU已完全淘汰,這是不去討論它。在32位c語言編譯器中,unsigned
long(DWORD)型變數可以表示一個32bit的整數,unsigned short(WORD)型變數可表示一個16bit的整數。兩個65535以內的數相乘,其結果完全可以存貯在一個unsigned
long變數中。另外,在好多32位編譯器中,也提供了64bit的整數型別(如在VC中,unsigned
__int64可表示一個64bit的整數,在GCC中,long
long可表示一個64位的整數)。同理兩個40億以內的數相乘,其結果可以用一個unsigned __int64
型的變數來儲存。讓一個具有一次可計算兩個32bit數乘法能力的CPU一次只計算1個1位10進位制數和一個整數的乘法,實在是一種浪費。下面我們提出兩種大數的表示法和運算方法。

第一種方法:用WORD型陣列表示大數,用DWORD型變數表示兩個WORD型變數的乘積。陣列的每個元素表示4位十進位制數。在運算時,從這個數的最後一個元素開始,依次乘以當前乘數並加上上次的進位,其和的低4位數依然存在原位置,高几位向前進位。當乘數i小於42萬時,其乘積加上進位可以用一個DWORD型變數表示,故這個程式可以計算上到42萬的階乘,當計算42萬的階乘時,佔用記憶體空間小於1.1兆位元組。至於速度,在計算1000/10000的階乘時,在迅馳1.7G電腦約需0.015/0.86秒。

 

#include "stdlib.h"

#include "stdio.h"

#include "math.h"

 

#define PI 3.1415926535897932384626433832795

#define RAD 10000

typedef unsigned long DWORD;

typedef unsigned short WORD;

//用stirling公式估算結果長度,稍微算得大一點

DWORD calcResultLen(DWORD n,DWORD rad)

{

    double r=0.5*log(2*PI) + (n+0.5)*log(n)-n;

   return (DWORD)(r/log(rad)+2);

}

 

void calcFac1(DWORD n)

{

    DWORD i,carry,prod,len;

    WORD *buff,*pHead,*pTail,*p;

if (n==0)        

      { printf("%d!=1",n); return; }

 //—計算並分配所需的儲存空間

 len=calcResultLen(n,RAD); 

 

     buff=(WORD*)malloc( sizeof(WORD)*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=(DWORD)(*p) * i +carry;

                     *p=(WORD)(prod % RAD);

                      carry=prod / RAD;

                  }

      while (carry>0)

                 {

                          pHead–;

                          *pHead=(WORD)(carry % RAD);

                          carry /= RAD;

                  }

    }

  //顯示計算結果

  printf("%d!=%d",n,*pHead); 

   for (p=pHead+1;p<=pTail;p++)

printf("%04d",*p);

     printf("/n");

   free(buff);//釋放分配的記憶體

}

int main(int argc, char* argv[])

{

DWORD n;

printf("please input n:");

scanf("%d",&n);

calcFac1(n);

return 0;

}

 

第二種方法:用DWORD型陣列表示大數,用unsigned __int64
表示兩個DWORD型變數的乘積。陣列的每個元素表示9位十進位制數。在運算時,從這個數的最後一個元素開始,依次乘以當前乘數i(i=1..n)並加上上次的進位,其和的低9位數依然存在原位置,高几位向前進位。從演算法上講,這個程式能夠計算到40億的階乘,在實際計算過程中,僅受限於記憶體的大小。至於速度,比前一個程式要慢一些,原因在於unsigned
__int64的除法較慢。我們將在下一篇文章給出解決方法,下面是採用該方法計算階乘的程式碼。

#define TEN9 1000000000

void calcFac2(DWORD n)

{

     DWORD *buff,*pHead,*pTail,*p;

     DWORD t,i,len;

     UINT64 carry,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 +carry;

                        *p=(DWORD)(prod % TEN9);

                        carry=prod / TEN9;

                  }

                  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);//釋放分配的記憶體

}

 入門篇之三彙編的威力

在上一篇文章《大數階乘之計算從入門到精通-入門篇之二》中,我們給出兩個計算大數階乘的程式,其中第2個程式由於用到64位整數的除法,速度反而更慢。在本文中,我們採用在C語言中嵌入彙編程式碼的方法,改進瓶頸部分,使計算速度提高到原先3倍多。

我們首先看一下計算大數階乘的核心程式碼(見下),可以看到,在每次迴圈中,需要計算1次64位的乘法,1次64位的加法,2次64位的除法(求餘視為除法)。我們知道,除法指令是慢於乘法指令的,對於64位的除法更是如此。在vc中,兩個64位數的除法是調aulldiv函式來實現的,兩個64位數的求餘運算是呼叫aullrem來實現的。通過分析這兩個函式的原始碼(彙編程式碼)得知,在做除法運算時,首先檢查除數,如果它小於2的32次方,則需兩次除法以得到一個64位商,如果除數大於等於232,運算將更為複雜。同樣在做求餘運算時,如果除數小於232,為了得到一個32位的餘數,也需2次除法。在我們這個例子中,除數為1000000000,小於232,因此,這段程式碼需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,*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+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",n,*pHead);

     for (p=pHead+1;p<=pTail;p++)

           printf("%09d",*p);

     printf("/n");

 free(buff);//釋放分配的記憶體

}

求N!的高精度演算法

 

本文是張一飛2001年寫的論文,原文可從http://oibh.kuye.cn/download/thesis/thesis2001_zhangyifei.zip處下載

 求N!的高精度演算法

本文中的演算法主要針對Pascal語言

這篇文章的內容

你瞭解高精度嗎?

你曾經使用過哪些資料結構?

你仔細思考過如何優化演算法嗎?

在這裡,你將看到怎樣成倍提速求N!的高精度演算法


Pascal中的標準整數型別

資料型別

值域

Shortint

-128~127

Byte

0~ 255

Integer

-32768~ 32768

Word

0~ 65535

LongInt

-2147483648~2147483647

Comp

-9.2e18~9.2e18

Comp雖然屬於實型,實際上是一個64位的整數

高精度演算法的基本思想

Pascal中的標準整數型別最多隻能處理在-263~263之間的整數。如果要支援更大的整數運算,就需要使用高精度

高精度演算法的基本思想,就是將無法直接處理的大整數,分割成若干可以直接處理的小整數段,把對大整數的處理轉化為對這些小整數段的處理


資料結構的選擇

每個小整數段保留儘量多的位

使用Comp型別

採用二進位制表示法


 

每個小整數段保留儘量多的位

一個例子:計算兩個15位數的和

Ø方法一

•分為15個小整數段,每段都是1位數,需要15次1位數加法

Ø方法二

•分為5個小整數段,每段都是3位數,需要5次3位數加法

Ø方法三

•Comp型別可以直接處理15位的整數,故1次加法就可以了

Ø比較

•用Integer計算1位數的加法和3位數的加法是一樣快的

•故方法二比方法一效率高

•雖然對Comp的操作要比Integer慢,但加法次數卻大大減少

•實踐證明,方法三比方法二更快


使用Comp型別

高精度運算中,每個小整數段可以用Comp型別表示

Comp有效數位為19~20位

求兩個高精度數的和,每個整數段可以保留17位

求高精度數與不超過m位整數的積,每個整數段可以保留18–m位

求兩個高精度數的積,每個整數段可以保留9位

如果每個小整數段保留k位十進位制數,實際上可以認為其只儲存了1位10k進位制數,簡稱為高進位制數,稱1位高進位制數為單精度數


採用二進位制表示法

採用二進位制表示,運算過程中時空效率都會有所提高,但題目一般需要以十進位制輸出結果,所以還要一個很耗時的進位制轉換過程。因此這種方法競賽中一般不採用,也不在本文討論之列.


演算法的優化

高精度乘法的複雜度分析

連乘的複雜度分析

設定快取

分解質因數求階乘

二分法求乘冪

分解質因數後的調整


高精度乘法的複雜度分析

計算n位高進位制數與m位高進位制數的積

Ø需要n*m次乘法

Ø積可能是n+m–1或n+m位高進位制數


連乘的複雜度分析(1)

一個例子:計算5*6*7*8

Ø方法一:順序連乘

•5*6=30,1*1=1次乘法

•30*7=210,2*1=2次乘法

•210*8=1680,3*1=3次乘法  共6次乘法

Ø方法二:非順序連乘

•5*6=30,1*1=1次乘法

•7*8=56,1*1= 1次乘法

•30*56=1680,2*2=4次乘法   共6次乘法


連乘的複雜度分析(2)

若“n位數*m位數=n+m位數”,則n個單精度數,無論以何種順序相乘,乘法次數一定為n(n-1)/2次

Ø證明:

•設F(n)表示乘法次數,則F(1)=0,滿足題設

•設k<n時,F(k)=k(k-1)/2,現在計算F(n)

•設最後一次乘法計算為“k位數*(n-k)位數”,則

•F(n)=F(k)+F(n-k)+k (n-k)=n(n-1)/2(與k的選擇無關)


設定快取(1)

一個例子:計算9*8*3*2

Ø方法一:順序連乘

•9*8=72,1*1=1次乘法

•72*3=216,2*1=2次乘法

•216*2=432,3*1=3次乘法

Ø方法二:非順序連乘

•9*8=72,1*1=1次乘法

•3*2=6,1*1=1次乘法

•72*6=432,2*1=2次乘法

共6次乘法

 


特點:n位數*m位數可能是n+m-1位數


設定快取(2)

考慮k+t個單精度數相乘a1*a2*…*ak*ak+1*…*ak+t

Ø設a1*a2*…*ak結果為m位高進位制數(假設已經算出)

Øak+1*…*ak+t結果為1位高進位制數

Ø若順序相乘,需要t次“m位數*1位數”,共mt次乘法

Ø可以先計算ak+1*…*ak+t,再一起乘,只需要m+t次乘法

在設定了快取的前提下,計算m個單精度數的積,如果結果為n位數,則乘法次數約為n(n–1)/2次,與m關係不大 

–設S=a1a2… am,S是n位高進位制數
–可以把乘法的過程近似看做,先將這m個數分為n組,每組的積仍然是一個單精度數,最後計算後面這n個數的積。時間主要集中在求最後n個數的積上,這時基本上滿足“n位數*m位數=n+m位數”,故乘法次數可近似的看做n(n-1)/2次


設定快取(3)

快取的大小

Ø設所選標準資料型別最大可以直接處理t位十進位制數

Ø設快取為k位十進位制數,每個小整數段儲存t–k位十進位制數

Ø設最後結果為n位十進位制數,則乘法次數約為

Øk/(n–k) ∑(i=1..n/k)i=(n+k)n/(2k(t–k)),其中k遠小於n

Ø要乘法次數最少,只需k (t–k)最大,這時k=t/2

Ø因此,快取的大小與每個小整數段大小一樣時,效率最高

Ø故在一般的連乘運算中,可以用Comp作為基本整數型別,每個小整數段為9位十進位制數,快取也是9位十進位制數


分解質因數求階乘

例:10!=28*34*52*7

Øn!分解質因數的複雜度遠小於nlogn,可以忽略不計

Ø與普通演算法相比,分解質因數後,雖然因子個數m變多了,但結果的位數n沒有變,只要使用了快取,乘法次數還是約為n(n-1)/2次

Ø因此,分解質因數不會變慢(這也可以通過實踐來說明)

Ø分解質因數之後,出現了大量求乘冪的運算,我們可以優化求乘冪的演算法。這樣,分解質因數的好處就體現出來了


二分法求乘冪

二分法求乘冪,即:

Øa2n+1=a2n*a

Øa2n=(an)2

Ø其中,a是單精度數

複雜度分析

Ø假定n位數與m位數的積是n+m位數

Ø設用二分法計算an需要F(n)次乘法

ØF(2n)=F(n)+n2,F(1)=0

Ø設n=2k,則有F(n)=F(2k)=∑(i=0..k–1)4i=(4k–1)/3=(n2–1)/3

與連乘的比較

Ø用連乘需要n(n-1)/2次乘法,二分法需要(n2–1)/3

Ø連乘比二分法耗時僅多50%

Ø採用二分法,複雜度沒有從n2降到nlogn


二分法求乘冪之優化平方演算法

怎樣優化

Ø(a+b)2=a2+2ab+b2

Ø例:123452=1232*10000+452+2*123*45*100

Ø把一個n位數分為一個t位數和一個n-t位數,再求平方

怎樣分

Ø設求n位數的平方需要F(n)次乘法

ØF(n)=F(t)+F(n-t)+t(n-t),F(1)=1

Ø用數學歸納法,可證明F(n)恆等於n(n+1)/2

Ø所以,無論怎樣分,效率都是一樣

Ø將n位數分為一個1位數和n–1位數,這樣處理比較方便


二分法求乘冪之複雜度分析

複雜度分析

Ø前面已經求出F(n)=n(n+1)/2,下面換一個角度來處理

ØS2=(∑(0≤i<n)ai10i)2=∑(0≤i<n)ai2102i+2∑(0≤i<j<n)aiaj10i+j

Ø一共做了n+C(n,2)=n(n+1)/2次乘法運算

Ø普通演算法需要n2次乘法,比改進後的慢1倍

改進求乘冪的演算法

Ø如果在用改進後的方法求平方,則用二分法求乘冪,需要(n+4)(n–1)/6次乘法,約是連乘演算法n(n–1)/2的三分之一


分解質因數後的調整(1)

為什麼要調整

Ø計算S=211310,可以先算211,再算310,最後求它們的積

Ø也可以根據S=211310=610*2,先算610,再乘以
2即可

Ø兩種演算法的效率是不同的


分解質因數後的調整(2)

什麼時候調整

Ø計算S=ax+kbx=(ab)xak

Ø當k<xlogab時,採用(ab)xak比較好,否則採用ax+kbx更快

Ø證明:

•可以先計算兩種演算法的乘法次數,再解不等式,就可以得到結論

•也可以換一個角度來分析。其實,兩種演算法主要差別在最後一步求積上。由於兩種方法,積的位數都是一樣的,所以兩個因數的差越大,乘法次數就越小

•∴當axbx–ak>ax+k–bx時,選用(ab)xak,反之,則採用ax+kbx

•∴axbx–ak>ax+k–bx

•∴(bx–ak)(ax+1)>0

•∴bx>ak

•這時k<xlogab


總結

內容小結

Ø用Comp作為每個小整數段的基本整數型別

Ø採用二進位制優化演算法

Ø高精度連乘時快取和快取的設定

Ø改進的求平方演算法

Ø二分法求乘冪

Ø分解質因數法求階乘以及分解質因數後的調整

應用

Ø高精度求乘冪(平方)

Ø高精度求連乘(階乘)

Ø高精度求排列組合


結束語

求N!的高精度演算法本身並不難,但我們仍然可以從多種角度對它進行優化。

其實,很多經典演算法都有優化的餘地。我們自己編寫的一些程式也不例外。只要用心去優化,說不準你就想出更好的演算法來了。

也許你認為本文中的優化毫無價值。確實是這樣,競賽中對高精度的要求很低,根本不需要優化。而我以高精度演算法為例,不過想談談如何優化一個演算法。我想說明的只有一點:演算法是可以優化的。