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

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

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

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

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

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

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

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

8、優(yōu)。這一計(jì)劃的進(jìn)展并不快，曾一度停止。目前，csdn上blog數(shù)量正在迅速地增加，我也萌生了寫blog的計(jì)劃，借此機(jī)會(huì)，對大數(shù)階乘之計(jì)算作一個(gè)整理，用文字和代碼詳述我的各個(gè)版本的算法和實(shí)現(xiàn)，同時(shí)也可能分析一些我在網(wǎng)上看到的別人寫的程序，當(dāng)然在這一過程中，我會(huì)繼續(xù)編寫未完成的版本或改寫以前己經(jīng)實(shí)現(xiàn)的版本，爭取使我公開的第一份代碼都是精品，這一過程可能是漫長的，但是我會(huì)盡力做下去。 菜鳥篇程序1，一個(gè)最直接的計(jì)算階乘的程序#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ù)雜了一些，使用了遞歸，一個(gè)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; 程序點(diǎn)評，這兩個(gè)程序在計(jì)算12以內(nèi)的數(shù)是正確，但當(dāng)n>12,程序的計(jì)算結(jié)果就完全錯(cuò)誤了，單從算法上講，程序并沒有錯(cuò)，可是這個(gè)程序到底錯(cuò)在什么地方呢？看來程序作者并沒有意識(shí)到，一個(gè)long型整數(shù)能夠表示的范圍是很有限的。當(dāng)n>=13時(shí)，計(jì)算結(jié)果溢出，在C語言，整數(shù)相乘時(shí)發(fā)生溢出時(shí)不會(huì)

11、產(chǎn)生任何異常，也不會(huì)給出任何警告。既然整數(shù)的范圍有限，那么能否用范圍更大的數(shù)據(jù)類型來做運(yùn)算呢？這個(gè)主意是不錯(cuò)，那么到底選擇那種數(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)行這個(gè)程序，將運(yùn)算結(jié)果并和windows計(jì)算器對比后發(fā)現(xiàn)，當(dāng)于在170以內(nèi)時(shí)，結(jié)果在誤差范圍內(nèi)是正確。但當(dāng)N>=171,結(jié)果就不能正確顯示了。這是為什么呢？和程序1類似，數(shù)據(jù)發(fā)生了溢出，即運(yùn)算結(jié)果超出的數(shù)據(jù)類型能夠表示的范圍?？磥鞢語言提供的數(shù)據(jù)類型不能滿足計(jì)算大數(shù)階乘的需要，為此只有兩個(gè)辦法。1.找一個(gè)能表示和處理大數(shù)的運(yùn)算的類庫。2.自己實(shí)現(xiàn)大數(shù)的存儲(chǔ)和運(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)實(shí)生活中，并不需要這么高的精度，比如銀河系的直徑有10萬光年，如果用原子核的直徑來度量，31位十進(jìn)制數(shù)就可使得誤差不超過一個(gè)原子核。 2.大數(shù)的表示： 2.1定點(diǎn)數(shù)和浮點(diǎn)數(shù)我們知道，在計(jì)算機(jī)中，數(shù)是存貯在內(nèi)存(RAM)中的。在內(nèi)存中存儲(chǔ)一個(gè)數(shù)有兩類格式，定點(diǎn)數(shù)和浮點(diǎn)數(shù)。定點(diǎn)數(shù)可以精確地表示一個(gè)整數(shù)，但數(shù)的范圍相對較小，如一個(gè)32比特的無符號(hào)整數(shù)可表示04294967295之間的數(shù)，可精確到910位數(shù)字(這里的數(shù)字指10進(jìn)制數(shù)字，如無特別指出，數(shù)字一律指10進(jìn)制數(shù)字),而一個(gè)8字節(jié)的無符號(hào)整數(shù)

14、則能精確到19位數(shù)字。浮點(diǎn)數(shù)能表示更大的范圍，但精度較低。當(dāng)表示的整數(shù)很大的，則可能存在誤差。一個(gè)8字節(jié)的雙精度浮點(diǎn)數(shù)可表示2.22*10-308到 1.79*10308之間的數(shù),可精確到1516位數(shù)字.2.2日常生活中的數(shù)的表示：對于這里提到的大數(shù)，上文提到的兩種表示法都不能滿足需求。為此，必需設(shè)計(jì)一種表示法來存儲(chǔ)大數(shù)。我們以日常生活中的十進(jìn)制數(shù)為例，看看是如何表示的。如一個(gè)數(shù)N被寫成"12345",則這個(gè)數(shù)可以用一個(gè)數(shù)組a來表示，a0=1,a1=2,a2=3,a3=4,a4=5,這時(shí)數(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ì)算機(jī)語言表示：在日常生活中，我們使用的阿拉伯?dāng)?shù)字只有09共10個(gè)，按照書寫習(xí)慣，一個(gè)字符表示1位數(shù)字。計(jì)算機(jī)中，我們常用的最小數(shù)據(jù)存儲(chǔ)單位是字節(jié)，C語言稱之為char,多個(gè)字節(jié)可表示一個(gè)更大的存儲(chǔ)單位。習(xí)慣上，兩個(gè)相鄰字節(jié)組合起來稱作一個(gè)短整數(shù)，在32位的C語言編譯器中稱之為short,匯編語語言一般記作word,4個(gè)相鄰的字節(jié)組合起來稱為一個(gè)長整數(shù)，在32位的C語言編譯器中稱之為long,匯編語言一般記作DWORD。在計(jì)算機(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)點(diǎn)。后者則具有容易顯示的優(yōu)點(diǎn)。我們試舉例說明： 例1：若一個(gè)大數(shù)用一個(gè)長為len的short型數(shù)組A來表示，并采用權(quán)從大到小的順序依次存放，數(shù)N表示為A0 * 65536(len-1)+A1 * 65536(len-2)+.Alen-1 * 655360,這時(shí)65536稱為基，其進(jìn)制2的16次方。 例2：若一個(gè)大數(shù)用一個(gè)長為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ù)組的每個(gè)元素可表示4位數(shù)字。一般地，這時(shí)數(shù)組的每一個(gè)元素為小于10000的數(shù)。類似的，可以用long型數(shù)組，基為2324294967296來表示一個(gè)大數(shù); 當(dāng)然可以用long型組，基為1000000000來表示，這種表示法，數(shù)組的每個(gè)元素可表示9位數(shù)字。當(dāng)然，也可以用char型數(shù)組，基為10。最后一種表示法，在新手寫的計(jì)算大數(shù)階乘程序最為常見，但計(jì)算速度卻是最慢的。使用更大的基，可以充分發(fā)揮CPU的計(jì)算能力，計(jì)算量將更少，計(jì)算速度更快，占用的存儲(chǔ)空間也更少。 2.4 大尾序和小尾序，我們在書寫一個(gè)數(shù)時(shí)，總是先寫權(quán)較大的數(shù)字，

18、后寫權(quán)較小的數(shù)字，但計(jì)算機(jī)中的數(shù)并不總是按這個(gè)的順序存放。小尾（Little Endian）就是低位字節(jié)排放在內(nèi)存的低端，高位字節(jié)排放在內(nèi)存的高端。例如對于一個(gè)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)存的高端。例如對于一個(gè)4字節(jié)的整數(shù)0x12345678,將在內(nèi)存中按照如下順序排放, Motorola處理器大多數(shù)使用

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

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

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

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

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

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

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

26、使其表示的數(shù)的大小不變。這樣由于尾數(shù)很小，再乘以一個(gè)數(shù)就不會(huì)溢出了，下面給出完整的代碼。 程序3.#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能夠計(jì)算的最大的n值，如果你想計(jì)算更大的數(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); /計(jì)算n的階乘 printfResult(&

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

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

31、t數(shù)，最高位清0。2.“>>4”,將其右移4位以清除最低位的4bit尾數(shù)，變成一個(gè)11bit的數(shù)（最高位5位為零）。3(rank-0x3ff)”, 減去0x3ff還原成真實(shí)的指數(shù)部分。以下為完整的代碼。 程序4：#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能夠計(jì)算的最大的n值，如果你想計(jì)算更大的數(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浮點(diǎn)數(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浮點(diǎn)數(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);/求計(jì)算結(jié)果的常用對數(shù)int logxN=(int)(floor(logx); /logx的整數(shù)部分sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);/轉(zhuǎn)化為科學(xué)計(jì)算法形式的字符串 int main(int argc, char* argv)struct bigNum r;char buff

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

36、nline函數(shù)，inline是C+關(guān)鍵字，如果使用純C的編譯器，可采用MACRO來代替。第二種優(yōu)化技術(shù)，將循環(huán)體內(nèi)的代碼展開，由一個(gè)循環(huán)步只做一次乘法，改為一個(gè)循環(huán)步做32次乘法。實(shí)際速度：計(jì)算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浮點(diǎn)數(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浮點(diǎn)數(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次方近似計(jì)算之二 在階乘之計(jì)算從入門到精通近似計(jì)算之一中，我們采用兩個(gè)數(shù)來表示中間結(jié)果，使得計(jì)算的范圍擴(kuò)大到1千萬，并可0.02秒內(nèi)完成10000000！的計(jì)算。在保證接近16位有效數(shù)字的前提下，有沒有更快的算法呢。很幸運(yùn)，有一個(gè)叫做Stirling的公式，它可以快速計(jì)算出一個(gè)較大的數(shù)的階乘，而且數(shù)越大，精度越高。有 :/mathworld.wolfram 查找Stirlings Series可找到2個(gè)利用斯特林公式求階乘的公式，一個(gè)是普通形式，一個(gè)是對數(shù)形式。前一個(gè)公式包含一個(gè)無窮級(jí)數(shù)和的形式，包含級(jí)數(shù)前5項(xiàng)的公式為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ù)的底。后一個(gè)公式也為一個(gè)無窮級(jí)數(shù)和的形式，一般稱這個(gè)級(jí)數(shù)為斯特林(Stirling)級(jí)數(shù),包括級(jí)數(shù)前3項(xiàng)的n!的對數(shù)公式為：ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n3 + 1/1260/n5)-,下面我們分別用這兩個(gè)公式計(jì)算n!.完整的代碼如下：#include "stdafx.h"#include "math.h"#d

44、efine PI 3.1415926535897932384626433832795#define E 2.7182818284590452353602874713527struct bigNum double n; /科學(xué)計(jì)數(shù)法表示的尾數(shù)部分 int e; /科學(xué)計(jì)數(shù)法表示的指數(shù)部分,若一個(gè)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, /級(jí)數(shù)的和 item; /級(jí)數(shù)的每一項(xiàng) 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, /級(jí)數(shù)的和 item; /級(jí)數(shù)的每一項(xiàng) 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); /用第一種方法，計(jì)算n的階乘 printfResult(&r,buff); /將結(jié)果轉(zhuǎn)化一個(gè)字符串 printf("%d!=%s by method 1n",n,buff); calcFac2(&r,n); /用第二種方法，計(jì)算n的階乘 printfResult(&r,buff); /將結(jié)果轉(zhuǎn)化一個(gè)字符串 printf("%d!=%s by method 2n",n,buff); return 0;程序運(yùn)行時(shí)間的測量 算

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