關(guān)于大數(shù)階乘的程序編寫

1、大數(shù)階乘序大數(shù)階乘的計算是一個有趣的話題，從中學(xué)生到大學(xué)教授，許多人都投入到這個問題的探索和研究之中，并發(fā)表了他們自己的研究成果。如果你用階乘作關(guān)鍵字在google上搜索，會找到許多此類文章，另外，如果你使用google學(xué)術(shù)搜索，也能找到一些計算大數(shù)階乘的學(xué)術(shù)論文。但這些文章和論文的深度有限，并沒有給出一個高速的算法和程序。 我和許多對大數(shù)階乘感興趣的人一樣，很早就開始編制大數(shù)階乘的程序。從2000年開始寫第一個大數(shù)階乘程序算起，到現(xiàn)在大約己有6-7年的時光，期間我寫了多個版本的階乘計算器，在階乘計算器的算法探討和程序的編寫和優(yōu)化上，我花費(fèi)了很大的時間和精力，品嘗了這一過程中的種種甘苦，我曾因

2、一個新算法的實現(xiàn)而帶來速度的提升而興奮，也曾因費(fèi)了九牛二虎之力但速度反而不及舊的版本而懊惱，也有過因解一個bug而通宵敖夜的情形。我寫的大數(shù)階乘的一些代碼片段散見于互聯(lián)網(wǎng)絡(luò)，而算法和構(gòu)想則常?？M繞在我的腦海。自以為，我對大數(shù)階乘計算器的算法探索在深度和廣度上均居于先進(jìn)水平。我常常想，應(yīng)該寫一個關(guān)于大數(shù)階乘計算的系列文章，一為整理自己的勞動成果，更重要的是可以讓同行分享我的知識和經(jīng)驗，也算為IT界做一點兒貢獻(xiàn)吧。 我的第一個大數(shù)階乘計算器始于2000年，那年夏天，我買了一臺電腦，開始在家專心學(xué)習(xí)VC，同時寫了我的第一個VC程序，一個仿制windows界面的計算器。該計算器的特點是高精度和高速度，

3、它可以將四則運(yùn)算的結(jié)果精確到6萬位以內(nèi)，將三角、對數(shù)和指數(shù)函數(shù)的結(jié)果精確到300位以內(nèi)，也可以計算開方和階乘等。當(dāng)時，我碰巧看到一個叫做實用計器的軟件。值得稱頌的是，該計算器的作者姜邊竟是一個高中生，他的這個計算器功能強(qiáng)大，獲得了各方很高的評價。該計算的功能之一是可以計算9000以內(nèi)階乘的精確值，且速度很快。在佩服之余，也激起我寫一個更好更快的程序的決心，經(jīng)過幾次改進(jìn)，終于使我的計算器在做階乘的精確計算時 (以9000!為例)，可比實用計算器快10倍。而當(dāng)精確到30多位有效數(shù)字時，可比windows自帶的計算器快上7500倍。其后的2001年1月，我在csdn上看到一個貼子，題目是“有誰可以用

4、四行代碼求出1000000的階乘”，我回復(fù)了這個貼子，給出一個相對簡潔的代碼，這是我在網(wǎng)上公布的第一個大數(shù)階乘的程序。這一時期，可以看作是我寫階乘計算器的第一個時期。 我寫階乘計算器的第二個時期始于2003年9月，在那時我寫了一組專門計算階乘的程序，按運(yùn)行速度來分，分為三個級別的版本，初級版、中級版和高級版。初級版本的算法許多人都能想到，中級版則采用大數(shù)乘以大數(shù)的硬乘法，高級版本在計算大數(shù)乘法時引入分治法。期間在csdn社區(qū)發(fā)了兩個貼子，“擂臺賽：計算n!(階乘)的精確值，速度最快者2000分送上”和“擂臺賽：計算n!(階乘)的精確值，速度最快者2000分送上（續(xù)）”。其高級算的版本完成于20

5、03年11月。此后，郭先強(qiáng)于2004年5月10日也發(fā)表了系列貼子，“擂臺：超大整數(shù)高精度快速算法”、“擂臺：超大整數(shù)高精度快速算法(續(xù))”和“擂臺：超大整數(shù)高精度快速算法(續(xù)2)”, 該貼重點展示了大數(shù)階乘計算器的速度。這個貼子一經(jīng)發(fā)表即引起了熱列的討論，除了我和郭先強(qiáng)先生外，郭雄輝也寫了同樣功能的程序，那時，大家都在持續(xù)改進(jìn)自己的程序，看誰的程序更快。初期，郭先強(qiáng)的稍稍領(lǐng)先，中途郭子將apfloat的代碼應(yīng)用到階乘計算器中，使得他的程序勝出，后期（2004年8月后）在我將程序作了進(jìn)一步的改進(jìn)后，其速度又稍勝于他們。在這個貼子中，arya提到一個開放源碼的程序，它的大數(shù)乘法采用FNTCRT(快

6、速數(shù)論變換中國剩余定理)。郭雄輝率先采用apflot來計算大數(shù)階乘，后來郭先強(qiáng)和我也參于到apfloat的學(xué)習(xí)和改進(jìn)過程中。在這點上，郭先強(qiáng)做得非常好，他在apfloat的基礎(chǔ)上，成功地優(yōu)化和改時算法，并應(yīng)用到大數(shù)階乘計算器上，同時他也將FNT算法應(yīng)用到他的超大整數(shù)高精度快速算法庫中，并在2004.10.18正式推出V3.0.2.1版。此后，我在2004年9月9日也完成一個采用FNT算法的版本，但卻不及郭先強(qiáng)的階乘計算器快。那時，郭先強(qiáng)的程序是我們所知的運(yùn)算速度最快的階乘計算器，其速度超過久負(fù)盛名的數(shù)學(xué)軟件Mathematica和Maple,以及通用高精度算法庫GMP。我寫階乘計算器的第三個時

7、間約開始于2006年，在2005年8月收到北大劉楚雄老師的一封e-mail，他提到了他寫的一個程序在計算階乘時比我們的更快。這使我非常吃驚，在詢問后得知，其核心部分使用的是ooura寫的FFT函數(shù)。ooura的FFT代碼完全公開，是世界上運(yùn)行的最快的FFT程序之一，從這點上，再次看到了我們和世界先進(jìn)水平的差距。佩服之余，我決定深入學(xué)習(xí)FFT算法，看看能否寫出和ooura速度相當(dāng)或者更快的程序，同時一個更大計劃開始形成，即寫一組包括更多算法的階乘計算器，包括使用FFT算法的終極版和使用無窮級數(shù)的stirling公式來計算部分精度的極速版，除此之外，我將重寫和優(yōu)化以前的版本，力爭使速度更快，代碼更

8、優(yōu)。這一計劃的進(jìn)展并不快，曾一度停止。目前，csdn上blog數(shù)量正在迅速地增加，我也萌生了寫blog的計劃，借此機(jī)會，對大數(shù)階乘之計算作一個整理，用文字和代碼詳述我的各個版本的算法和實現(xiàn)，同時也可能分析一些我在網(wǎng)上看到的別人寫的程序，當(dāng)然在這一過程中，我會繼續(xù)編寫未完成的版本或改寫以前己經(jīng)實現(xiàn)的版本，爭取使我公開的第一份代碼都是精品，這一過程可能是漫長的，但是我會盡力做下去。 菜鳥篇程序1，一個最直接的計算階乘的程序#include "stdio.h"#include "stdlib.h" int main(int argc, char* argv)

9、long i,n,p; printf("n=?"); scanf("%d",&n); p=1; for (i=1;i<=n;i+) p*=i; printf("%d!=%dn",n,p); return 0; 程序2，稍微復(fù)雜了一些，使用了遞歸，一個c+初學(xué)者寫的程序#include <iostream.h> long int fac(int n); void main() int n; cout<<"input a positive integer:" cin>>

10、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èi)的數(shù)是正確，但當(dāng)n>12,程序的計算結(jié)果就完全錯誤了，單從算法上講，程序并沒有錯，可是這個程序到底錯在什么地方呢？看來程序作者并沒有意識到，一個long型整數(shù)能夠表示的范圍是很有限的。當(dāng)n>=13時，計算結(jié)果溢出，在C語言，整數(shù)相乘時發(fā)生溢出時不會

11、產(chǎn)生任何異常，也不會給出任何警告。既然整數(shù)的范圍有限，那么能否用范圍更大的數(shù)據(jù)類型來做運(yùn)算呢？這個主意是不錯，那么到底選擇那種數(shù)據(jù)類型呢？有人想到了double類型，將程序1中l(wèi)ong型換成double類型，結(jié)果如下:#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

12、; printf("%lf!=%.16gn",n,p); return 0; 運(yùn)行這個程序，將運(yùn)算結(jié)果并和windows計算器對比后發(fā)現(xiàn)，當(dāng)于在170以內(nèi)時，結(jié)果在誤差范圍內(nèi)是正確。但當(dāng)N>=171,結(jié)果就不能正確顯示了。這是為什么呢？和程序1類似，數(shù)據(jù)發(fā)生了溢出，即運(yùn)算結(jié)果超出的數(shù)據(jù)類型能夠表示的范圍?？磥鞢語言提供的數(shù)據(jù)類型不能滿足計算大數(shù)階乘的需要，為此只有兩個辦法。1.找一個能表示和處理大數(shù)的運(yùn)算的類庫。2.自己實現(xiàn)大數(shù)的存儲和運(yùn)算問題。方法1不在本文的討論的范圍內(nèi)。本系列的后續(xù)文章將圍繞方法2來展開。大數(shù)的表示 1.大數(shù)，這里提到的大數(shù)指有效數(shù)字非常多的數(shù)，

13、它可能包含少則幾十、幾百位十進(jìn)制數(shù)，多則幾百萬或者更多位十進(jìn)制數(shù)。有效數(shù)字這么多的數(shù)只具有數(shù)學(xué)意義，在現(xiàn)實生活中，并不需要這么高的精度，比如銀河系的直徑有10萬光年，如果用原子核的直徑來度量，31位十進(jìn)制數(shù)就可使得誤差不超過一個原子核。 2.大數(shù)的表示： 2.1定點數(shù)和浮點數(shù)我們知道，在計算機(jī)中，數(shù)是存貯在內(nèi)存(RAM)中的。在內(nèi)存中存儲一個數(shù)有兩類格式，定點數(shù)和浮點數(shù)。定點數(shù)可以精確地表示一個整數(shù)，但數(shù)的范圍相對較小，如一個32比特的無符號整數(shù)可表示04294967295之間的數(shù)，可精確到910位數(shù)字(這里的數(shù)字指10進(jìn)制數(shù)字，如無特別指出，數(shù)字一律指10進(jìn)制數(shù)字),而一個8字節(jié)的無符號整數(shù)

14、則能精確到19位數(shù)字。浮點數(shù)能表示更大的范圍，但精度較低。當(dāng)表示的整數(shù)很大的，則可能存在誤差。一個8字節(jié)的雙精度浮點數(shù)可表示2.22*10-308到 1.79*10308之間的數(shù),可精確到1516位數(shù)字.2.2日常生活中的數(shù)的表示：對于這里提到的大數(shù)，上文提到的兩種表示法都不能滿足需求。為此，必需設(shè)計一種表示法來存儲大數(shù)。我們以日常生活中的十進(jìn)制數(shù)為例，看看是如何表示的。如一個數(shù)N被寫成"12345",則這個數(shù)可以用一個數(shù)組a來表示，a0=1,a1=2,a2=3,a3=4,a4=5,這時數(shù)N=a4*100+a3*101+a2*102+a1*103+a0*104,(104表示

15、10的4次方，下同),10i可以叫做權(quán)，在日常生活中，a0被稱作萬位，也說是說它的權(quán)是10000，類似的，a1被稱作千位，它的權(quán)是1000。 2.3大數(shù)在計算機(jī)語言表示：在日常生活中，我們使用的阿拉伯?dāng)?shù)字只有09共10個，按照書寫習(xí)慣，一個字符表示1位數(shù)字。計算機(jī)中，我們常用的最小數(shù)據(jù)存儲單位是字節(jié)，C語言稱之為char,多個字節(jié)可表示一個更大的存儲單位。習(xí)慣上，兩個相鄰字節(jié)組合起來稱作一個短整數(shù)，在32位的C語言編譯器中稱之為short,匯編語語言一般記作word,4個相鄰的字節(jié)組合起來稱為一個長整數(shù)，在32位的C語言編譯器中稱之為long,匯編語言一般記作DWORD。在計算機(jī)中，按照權(quán)的不

16、同，數(shù)的表示可分為兩種，2進(jìn)制和10進(jìn)制，嚴(yán)格說來，應(yīng)該是2k進(jìn)制和10K進(jìn)制，前者具占用空間少，運(yùn)算速度快的優(yōu)點。后者則具有容易顯示的優(yōu)點。我們試舉例說明： 例1：若一個大數(shù)用一個長為len的short型數(shù)組A來表示，并采用權(quán)從大到小的順序依次存放，數(shù)N表示為A0 * 65536(len-1)+A1 * 65536(len-2)+.Alen-1 * 655360,這時65536稱為基，其進(jìn)制2的16次方。 例2：若一個大數(shù)用一個長為len的short型數(shù)組A來表示并采用權(quán)從大到小的順序依次存放，數(shù)N=A0 * 10000(len-1)+A1 * 10000(len-2)+.Alen-1 *

17、100000,這里10000稱為基，其進(jìn)制為10000，即：104，數(shù)組的每個元素可表示4位數(shù)字。一般地，這時數(shù)組的每一個元素為小于10000的數(shù)。類似的，可以用long型數(shù)組，基為2324294967296來表示一個大數(shù); 當(dāng)然可以用long型組，基為1000000000來表示，這種表示法，數(shù)組的每個元素可表示9位數(shù)字。當(dāng)然，也可以用char型數(shù)組，基為10。最后一種表示法，在新手寫的計算大數(shù)階乘程序最為常見，但計算速度卻是最慢的。使用更大的基，可以充分發(fā)揮CPU的計算能力，計算量將更少，計算速度更快，占用的存儲空間也更少。 2.4 大尾序和小尾序，我們在書寫一個數(shù)時，總是先寫權(quán)較大的數(shù)字，

18、后寫權(quán)較小的數(shù)字，但計算機(jī)中的數(shù)并不總是按這個的順序存放。小尾（Little Endian）就是低位字節(jié)排放在內(nèi)存的低端，高位字節(jié)排放在內(nèi)存的高端。例如對于一個4字節(jié)的整數(shù)0x12345678,將在內(nèi)存中按照如下順序排放, Intel處理器大多數(shù)使用小尾(Little Endian)字節(jié)序。 Address0: 0x78 Address1: 0x56Address2: 0x34 Address3:0x12大尾（Big Endian）就是高位字節(jié)排放在內(nèi)存的低端，低位字節(jié)排放在內(nèi)存的高端。例如對于一個4字節(jié)的整數(shù)0x12345678,將在內(nèi)存中按照如下順序排放, Motorola處理器大多數(shù)使用

19、大尾(Big Endian)字節(jié)序。Address0: 0x12 Address1: 0x34Address2: 0x56 Address3:0x78 類似的，一個大數(shù)的各個元素的排列方式既可以采用低位在前的方式，也可以采用高位在前的方式，說不上那個更好，各有利弊吧。我習(xí)慣使用高位在前的方式。2.5 不完全精度的大數(shù)表示：盡管以上的表示法可準(zhǔn)確的表示一個整數(shù)，但有時可能只要求計算結(jié)果只精確到有限的幾位。如用 windows自帶的計算器計算1000的階乘時，只能得到大約32位的數(shù)字，換名話說，windows計算器的精度為32位。1000的階乘是一個整數(shù)，但我們只要它的前幾位有效數(shù)字，象windo

20、ws計算器這樣，只能表示部分有效數(shù)字的表示法叫不完全精度，不完全精度不但占用空間省，更重要的是，在只要求計算結(jié)果為有限精度的情況下，可大大減少計算量。大數(shù)的不完全精度的表示法除了需要用數(shù)組存儲有數(shù)數(shù)字外，還需要一個數(shù)來表示第一個有效數(shù)字的權(quán)，10的階乘約等于4.023872600770937e+2567，則第一個有效數(shù)字的權(quán)是102567，這時我們把2567叫做階碼。在這個例子中，我們可以用一個長為16的char型數(shù)組和一個數(shù)來表示，前者表示各位有效數(shù)字，數(shù)組的各個元素依次為:4,0,2,3,8,7,2,6,0,0,7,7,0,9,3,7，后者表示階碼，值為2567。 2.6 大數(shù)的鏈?zhǔn)酱鎯Ψ?/p>

21、 如果我們搜索大數(shù)階乘的源代碼，就會發(fā)現(xiàn)，有許多程序采用鏈表存儲大數(shù)。盡管這種存儲方式能夠表示大數(shù)，也不需要事先知道一個特定的數(shù)有多少位有效數(shù)字，可以在運(yùn)算過程中自動擴(kuò)展鏈表長度。但是，如果基于運(yùn)算速度和內(nèi)存的考慮，強(qiáng)烈不建議采用這種存儲方式，因為：1. 這種存儲方式的內(nèi)存利用率很低?；诖髷?shù)乘法的計算和顯示，一般需要定義雙鏈表，假如我們用1個char表示1位十進(jìn)制數(shù)，則可以這樣定義鏈表的節(jié)點： struct _node struct _node* pre;struct _node* next;char n; ;當(dāng)編譯器采用默認(rèn)設(shè)置，在通常的32位編譯器，這個結(jié)構(gòu)體將占用12字節(jié)。但這并不等于

22、說，分配具有1000個節(jié)點的鏈表需要1000*12字節(jié)。不要忘記，操作系統(tǒng)或者庫函數(shù)在從內(nèi)存池中分配和釋放內(nèi)存時，也需要維護(hù)一個鏈表。實驗表明，在VC編譯的程序，一個節(jié)點總的內(nèi)存占用量為 sizeof(struct _node) 向上取16的倍數(shù)再加8字節(jié)。也就是說，采用這種方式表示n位十進(jìn)制數(shù)需要 n*24字節(jié)，而采用1個char型數(shù)組僅需要n字節(jié)。2采用鏈表方式表示大數(shù)的運(yùn)行速度很慢.2.1如果一個大數(shù)需要n個節(jié)點，需要調(diào)用n次malloc(C)或new(C+)函數(shù)，采用動態(tài)數(shù)組則不要用調(diào)用這么多次malloc.2.2 存取數(shù)組表示的大數(shù)比鏈表表示的大數(shù)具有更高的cache命中率。數(shù)組的各

23、個元素的地址是連續(xù)的，而鏈表的各個節(jié)點在內(nèi)存中的地址是不連續(xù)的，而且具有更大的數(shù)據(jù)量。因此前者的cache的命中率高于后者，從而導(dǎo)致運(yùn)行速度高于后者。2.3對數(shù)組的順序訪問也比鏈表快，如p1表示數(shù)組當(dāng)前元素的地址，則計算數(shù)組的下一個地址時一般用p1+,而對鏈表來說則可能是p2=p2->next,毫無疑問，前者的執(zhí)行速度更快。 近似計算之一 在階乘之計算從入門到精通菜鳥篇中提到，使用double型數(shù)來計算階乘，當(dāng)n>170,計算結(jié)果就超過double數(shù)的最大范圍而發(fā)生了溢出，故當(dāng)n170時，就不能用這個方法來計算階乘了，果真如此嗎？No,只要肯動腦筋,辦法總是有的。通過windows

24、計算器，我們知道，171！1.2410180702176678234248405241031e+309，雖然這個數(shù)不能直接用double型的數(shù)來表示，但我們可以用別的方法來表示。通過觀察這個數(shù)，我們發(fā)現(xiàn)，這個數(shù)的表示法為科學(xué)計算法，它用兩部分組成，一是尾數(shù)部分1.2410180702176678234248405241031，另一個指數(shù)部分309。不妨我們用兩個數(shù)來表示這個超大的數(shù)，用double型的數(shù)來表示尾數(shù)部分，用一個long型的數(shù)來表示指數(shù)部分。這會涉及兩個問題：其一是輸出，這好說，在輸出時將這兩個部分合起來就可以了。另一個就是計算部分了，這是難點所在（其實也不難）。下面我們分析一下，

25、用什么方法可以保證不會溢出呢？我們考慮170！，這個數(shù)約等于7.257415e+306，可以用double型來表示，但當(dāng)這個數(shù)乘以171就溢出了。我們看看這個等式：7.257415e+3067.257415e+306 * 100 (注1)（如用兩個數(shù)來表示，則尾數(shù)部分7.257415e+306，指數(shù)部分0)(7.257415e+306 / 10300 )* (100*10300)=（7.257415e6）*（10 300）（如用兩個數(shù)來表示，則尾數(shù)部分7.257415e+6，指數(shù)部分300） 依照類似的方法，在計算過程中，當(dāng)尾數(shù)很大時，我們可以重新調(diào)整尾數(shù)和指數(shù)，縮小尾數(shù)，同時相應(yīng)地增大指數(shù)，

26、使其表示的數(shù)的大小不變。這樣由于尾數(shù)很小，再乘以一個數(shù)就不會溢出了，下面給出完整的代碼。 程序3.#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能夠計算的最大的n值，如果你想計算更大的數(shù)對數(shù)，可將其改為更大的值#define MAX_MANTISSA (1e308/MAX_N) /最大尾數(shù) struct bigNum double n1; /表示尾數(shù)部分 int n2; /表示指數(shù)部分 ; void calcFac(struct bigNum *p,int n) int i;

27、 double MAX_POW10_LOG=(floor(log10(1e308/MAX_N); /最大尾數(shù)的常用對數(shù)的整數(shù)部分, 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(str

28、uct 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 buff32; int n; printf("n=?"); scanf("%d",&n); calcFac(&r,n); /計算n的階乘 printfResult(&

29、amp;r,buff); /將結(jié)果轉(zhuǎn)化一個字符串 printf("%d!=%sn",n,buff); return 0; 以上代碼中的數(shù)的表示方式為：數(shù)的值等于尾數(shù)乘以 10 指數(shù)部分，當(dāng)然我們也可以表示為：尾數(shù) 乘以 2 指數(shù)部分，這們將會帶來這樣的好處，在調(diào)整尾數(shù)部分和指數(shù)部分時，不用除法，可以依據(jù)浮點數(shù)的格式直讀取階碼和修改階碼(上文提到的指數(shù)部分的標(biāo)準(zhǔn)稱呼)，同時也可在一定程序上減少誤差。為了更好的理解下面的代碼，有必要了解一下浮點數(shù)的格式。浮點數(shù)主要分為32bit單精度和64bit雙精度兩種。本方只討論64bit雙精度(double型)浮點數(shù)的格式，一個doubl

30、e型浮點數(shù)包括8個字節(jié)(64bit)，我們把最低位記作bit0,最高位記作bit63,則一個浮點數(shù)各個部分定義為:第一部分尾數(shù)：bit0至bit51，共計52bit，第二部分階碼:bit52-bit62,共計11bit，第三部分符號位:bit63，0:表示正數(shù)，1表示負(fù)數(shù)。如一個數(shù)為0.xxxx * 2 exp,則exp表示指數(shù)部分，范圍為1023到1024，實際存儲時采用移碼的表示法，即將exp的值加上0x3ff，使其變?yōu)橐粋€0到2047范圍內(nèi)的一個值。函數(shù)GetExpBase2 中各語句含義如下：1.“(*pWord & 0x7fff)”，得到一個bit48-bit63這個16bi

31、t數(shù)，最高位清0。2.“>>4”,將其右移4位以清除最低位的4bit尾數(shù)，變成一個11bit的數(shù)（最高位5位為零）。3(rank-0x3ff)”, 減去0x3ff還原成真實的指數(shù)部分。以下為完整的代碼。 程序4：#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能夠計算的最大的n值，如果你想計算更大的數(shù)對數(shù)，可將其改為更大的值#define MAX_MANTISSA (1e308/MAX_N) /最大尾數(shù)typedef unsigned short WORD; str

32、uct bigNum double n1; /表示尾數(shù)部分int n2; /表示階碼部分;short GetExpBase2(double a) / 獲得 a 的階碼 / 按照IEEE 754浮點數(shù)格式，取得階碼，僅僅適用于Intel 系列 cpu WORD *pWord=(WORD *)(&a)+3; WORD rank = ( (*pWord & 0x7fff) >>4 ); return (short)(rank-0x3ff); double GetMantissa(double a) / 獲得 a 的 尾數(shù)/ 按照IEEE 754浮點數(shù)格式，取得尾數(shù)，僅僅適

33、用于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)/繼續(xù)相乘可能溢出，調(diào)整之 p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1

34、); p->n1 *=(double)i; void printfResult(struct bigNum *p,char buff)double logx=log10(p->n1)+ p->n2 * log10(2);/求計算結(jié)果的常用對數(shù)int logxN=(int)(floor(logx); /logx的整數(shù)部分sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);/轉(zhuǎn)化為科學(xué)計算法形式的字符串 int main(int argc, char* argv)struct bigNum r;char buff

35、32;int n;printf("n=?"); scanf("%d",&n);calcFac(&r,n); /計算n的階乘printfResult(&r,buff); /將結(jié)果轉(zhuǎn)化一個字符串printf("%d!=%sn",n,buff);return 0; 程序優(yōu)化的威力：程序4是一個很好的程序，它可以計算到1千萬（當(dāng)n更大時，p->n2可能溢出）的階乘，但是從運(yùn)行速度上講，它仍有潛力可挖，在采用了兩種優(yōu)化技術(shù)后，（見程序5），速度竟提高5倍，甚至超出筆者的預(yù)計。第一種優(yōu)化技術(shù)，將頻繁調(diào)用的函數(shù)定義成i

36、nline函數(shù)，inline是C+關(guān)鍵字，如果使用純C的編譯器，可采用MACRO來代替。第二種優(yōu)化技術(shù)，將循環(huán)體內(nèi)的代碼展開，由一個循環(huán)步只做一次乘法，改為一個循環(huán)步做32次乘法。實際速度：計算10000000！，程序4需要0.11667秒，程序5只需要0.02145秒。測試環(huán)境為迅馳1.7G,256M RAM。關(guān)于程序優(yōu)化方面的內(nèi)容，將在后續(xù)文章專門講述。下面是被優(yōu)化后的部分代碼，其余代碼不變。 程序5的部分代碼： inline short GetExpBase2(double a) / 獲得 a 的階碼/ 按照IEEE 754浮點數(shù)格式，取得階碼，僅僅適用于Intel 系列 cpu WOR

37、D *pWord=(WORD *)(&a)+3; WORD rank = ( (*pWord & 0x7fff) >>4 ); return (short)(rank-0x3ff);inline double GetMantissa(double a) / 獲得 a 的 尾數(shù) / 按照IEEE 754浮點數(shù)格式，取得尾數(shù)，僅僅適用于Intel 系列 cpu WORD *pWord=(WORD *)(&a)+3; *pWord &= 0x800f; /清除階碼 *pWord |= 0x3ff0; /重置階碼 return a; void calcFac

38、(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

39、.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 *=

40、(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

41、+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-

42、>n1 *=(double)(i); 注1：100，表示10的0次方近似計算之二 在階乘之計算從入門到精通近似計算之一中，我們采用兩個數(shù)來表示中間結(jié)果，使得計算的范圍擴(kuò)大到1千萬，并可0.02秒內(nèi)完成10000000！的計算。在保證接近16位有效數(shù)字的前提下，有沒有更快的算法呢。很幸運(yùn)，有一個叫做Stirling的公式，它可以快速計算出一個較大的數(shù)的階乘，而且數(shù)越大，精度越高。有 :/mathworld.wolfram 查找Stirlings Series可找到2個利用斯特林公式求階乘的公式，一個是普通形式，一個是對數(shù)形式。前一個公式包含一個無窮級數(shù)和的形式，包含級數(shù)前5項的公式為n!=

43、sqrt(2*PI*n)*(n/e)n*(1+ 1/12/n+ 1/288/n2 139/51840/n3-571/2488320/n4+),這里的PI為圓周率，e為自然對數(shù)的底。后一個公式也為一個無窮級數(shù)和的形式，一般稱這個級數(shù)為斯特林(Stirling)級數(shù),包括級數(shù)前3項的n!的對數(shù)公式為：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"#d

44、efine PI 3.1415926535897932384626433832795#define E 2.7182818284590452353602874713527struct bigNum double n; /科學(xué)計數(shù)法表示的尾數(shù)部分 int e; /科學(xué)計數(shù)法表示的指數(shù)部分,若一個bigNum為x,這x=n * 10e;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 printfRe

45、sult(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/n2 -139/51840/n3-571/2488320/n4+)void calcFac1(struct bigNum *p,int n) double logx, s, /級數(shù)的和 item; /級數(shù)的每一項 int i; / xn=

46、 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/n2 -139/51840/n3-571/2488320/n4+) for (item=1.0,s=0.0,i=0;i<sizeof(a1)/sizeof(double);i+) s+= item * a1i; item /= (double)n; /item= 1/(

47、ni) p->n *=s;/ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n3 + 1/1260/n5 -.)void calcFac2(struct bigNum *p,int n) double logR; double s, /級數(shù)的和 item; /級數(shù)的每一項 int i; logR=0.5*log(2.0*PI)+(double)n+0.5)*log(n)-(double)n; /s= (1/12/n -1/360/n3 + 1/1260/n5) for (item=1/(double)n,s=0.0,i=0;i<

48、sizeof(a2)/sizeof(double);i+) s+= item * a2i; item /= (double)(n)* (double)n; /item= 1/(n(2i+1) logR+=s; /根據(jù)換底公式，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 buff32; int n; printf("n=?&qu

49、ot;); scanf("%d",&n); calcFac1(&r,n); /用第一種方法，計算n的階乘 printfResult(&r,buff); /將結(jié)果轉(zhuǎn)化一個字符串 printf("%d!=%s by method 1n",n,buff); calcFac2(&r,n); /用第二種方法，計算n的階乘 printfResult(&r,buff); /將結(jié)果轉(zhuǎn)化一個字符串 printf("%d!=%s by method 2n",n,buff); return 0;程序運(yùn)行時間的測量 算

50、法的好壞有好多評價指標(biāo)，其中一個重要的指標(biāo)是時間復(fù)雜度。如果兩個程序完成一個同樣的任務(wù)，即功能相同，處理的數(shù)據(jù)相同，那么運(yùn)行時間較短者為優(yōu)。操作系統(tǒng)和庫函數(shù)一般都提供了對時間測量的函數(shù)，這么函數(shù)一般都會返回一個代表當(dāng)前時間的數(shù)值，通過在運(yùn)行某個程序或某段代碼之前調(diào)用一次時間函數(shù)來得到一個數(shù)值，在程序或者代碼段運(yùn)行之后再調(diào)用一次時間函數(shù)來得到另一個數(shù)值，將后者減去前者即為程序的運(yùn)行時間。在windwos平臺(指windwow95及以后的版本，下同),常用的用于測量時間的函數(shù)或方法有三種：1.API函數(shù)GetTickCount或C函數(shù)clock,2.API函數(shù)QueryPerformanceCounter，3:匯編指令RDSTC1API函數(shù)GetTickCount：函數(shù)原形：DWORD GetTickCount(VOID)；該函數(shù)取回從電腦開機(jī)至現(xiàn)在的毫秒數(shù)，即每個時間單位為1毫秒。他的分辨率比較低，常用在測量用時較長程序，如果你的程序用時為100毫秒以上，可以使用這個函數(shù).另一個和GetTickCount類似的函數(shù)是clock，該函數(shù)的回的時間的單位的是CLOCKS_PER_SEC，在windows95/2000操作系統(tǒng)，該值是1000,也就是說，在windows平臺，這兩個函數(shù)的功能幾乎完全相同。2.API函數(shù)QueryPerformanceCounter：函數(shù)原形：BOOL Q