大数阶乘(一)
序
大數階乘的計算是一個有趣的話題，從中學生到大學教授，許多人都投入到這個問題的探索和研究之中，并發(fā)表了他們自己的研究成果。如果你用階乘作關鍵字在google上搜索，會找到許多此類文章，另外，如果你使用google學術搜索，也能找到一些計算大數階乘的學術論文。但這些文章和論文的深度有限，并沒有給出一個高速的算法和程序。
?
我和許多對大數階乘感興趣的人一樣，很早就開始編制大數階乘的程序。從2000年開始寫第一個大數階乘程序算起，到現在大約己有6-7年的時光，期間我寫了多個版本的階乘計算器，在階乘計算器的算法探討和程序的編寫和優(yōu)化上，我花費了很大的時間和精力，品嘗了這一過程中的種種甘苦，我曾因一個新算法的實現而帶來速度的提升而興奮，也曾因費了九牛二虎之力但速度反而不及舊的版本而懊惱，也有過因解一個bug而通宵敖夜的情形。我寫的大數階乘的一些代碼片段散見于互聯網絡，而算法和構想則常常縈繞在我的腦海。自以為，我對大數階乘計算器的算法探索在深度和廣度上均居于先進水平。我常常想，應該寫一個關于大數階乘計算的系列文章，一為整理自己的勞動成果，更重要的是可以讓同行分享我的知識和經驗，也算為IT界做一點兒貢獻吧。
我的第一個大數階乘計算器始于2000年，那年夏天，我買了一臺電腦，開始在家專心學習VC，同時寫了我的第一個VC程序，一個仿制windows界面的計算器。該計算器的特點是高精度和高速度，它可以將四則運算的結果精確到6萬位以內，將三角、對數和指數函數的結果精確到300位以內，也可以計算開方和階乘等。當時，我碰巧看到一個叫做實用計器的軟件。值得稱頌的是，該計算器的作者姜邊竟是一個高中生，他的這個計算器功能強大，獲得了各方很高的評價。該計算的功能之一是可以計算9000以內階乘的精確值，且速度很快。在佩服之余，也激起我寫一個更好更快的程序的決心，經過幾次改進，終于使我的計算器在做階乘的精確計算時 (以9000!為例)，可比實用計算器快10倍。而當精確到30多位有效數字時，可比windows自帶的計算器快上7500倍。其后的2001年1月，我在csdn上看到一個貼子，題目是“有誰可以用四行代碼求出1000000的階乘”，我回復了這個貼子，給出一個相對簡潔的代碼，這是我在網上公布的第一個大數階乘的程序。這一時期，可以看作是我寫階乘計算器的第一個時期。
?
我寫階乘計算器的第二個時期始于2003年9月，在那時我寫了一組專門計算階乘的程序，按運行速度來分，分為三個級別的版本，初級版、中級版和高級版。初級版本的算法許多人都能想到，中級版則采用大數乘以大數的硬乘法，高級版本在計算大數乘法時引入分治法。期間在csdn社區(qū)發(fā)了兩個貼子，“擂臺賽：計算n!(階乘)的精確值，速度最快者2000分送上”和“擂臺賽：計算n!(階乘)的精確值，速度最快者2000分送上（續(xù)）”。其高級算的版本完成于2003年11月。此后，郭先強于2004年5月10日也發(fā)表了系列貼子，“擂臺：超大整數高精度快速算法”、“擂臺：超大整數高精度快速算法(續(xù))”和“擂臺：超大整數高精度快速算法(續(xù)2)”, 該貼重點展示了大數階乘計算器的速度。這個貼子一經發(fā)表即引起了熱列的討論，除了我和郭先強先生外，郭雄輝也寫了同樣功能的程序，那時，大家都在持續(xù)改進自己的程序，看誰的程序更快。初期，郭先強的稍稍領先，中途郭子將apfloat的代碼應用到階乘計算器中，使得他的程序勝出，后期（2004年8月后）在我將程序作了進一步的改進后，其速度又稍勝于他們。在這個貼子中，arya提到一個開放源碼的程序，它的大數乘法采用FNT＋CRT(快速數論變換＋中國剩余定理)。郭雄輝率先采用apflot來計算大數階乘，后來郭先強和我也參于到apfloat的學習和改進過程中。在這點上，郭先強做得非常好，他在apfloat的基礎上，成功地優(yōu)化和改時算法，并應用到大數階乘計算器上，同時他也將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公式來計算部分精度的極速版，除此之外，我將重寫和優(yōu)化以前的版本，力爭使速度更快，代碼更優(yōu)。這一計劃的進展并不快，曾一度停止。
目前，csdn上blog數量正在迅速地增加，我也萌生了寫blog的計劃，借此機會，對大數階乘之計算作一個整理，用文字和代碼詳述我的各個版本的算法和實現，同時也可能分析一些我在網上看到的別人寫的程序，當然在這一過程中，我會繼續(xù)編寫未完成的版本或改寫以前己經實現的版本，爭取使我公開的第一份代碼都是精品，這一過程可能是漫長的，但是我會盡力做下去。
菜鳥篇
程序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語言，整數相乘時發(fā)生溢出時不會產生任何異常，也不會給出任何警告。既然整數的范圍有限，那么能否用范圍更大的數據類型來做運算呢？這個主意是不錯，那么到底選擇那種數據類型呢？有人想到了double類型，將程序1中l(wèi)ong型換成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計算器對比后發(fā)現，當于在170以內時，結果在誤差范圍內是正確。但當N>=171,結果就不能正確顯示了。這是為什么呢？和程序1類似，數據發(fā)生了溢出，即運算結果超出的數據類型能夠表示的范圍。看來C語言提供的數據類型不能滿足計算大數階乘的需要，為此只有兩個辦法。1.找一個能表示和處理大數的運算的類庫。2.自己實現大數的存儲和運算問題。方法1不在本文的討論的范圍內。本系列的后續(xù)文章將圍繞方法2來展開。
?
大數的表示
?
?
1.大數，這里提到的大數指有效數字非常多的數，它可能包含少則幾十、幾百位十進制數，多則幾百萬或者更多位十進制數。有效數字這么多的數只具有數學意義，在現實生活中，并不需要這么高的精度，比如銀河系的直徑有10萬光年，如果用原子核的直徑來度量，31位十進制數就可使得誤差不超過一個原子核。
?
2.大數的表示：
2.1定點數和浮點數
我們知道，在計算機中，數是存貯在內存(RAM)中的。在內存中存儲一個數有兩類格式，定點數和浮點數。定點數可以精確地表示一個整數，但數的范圍相對較小，如一個32比特的無符號整數可表示0－4294967295之間的數，可精確到9－10位數字(這里的數字指10進制數字，如無特別指出，數字一律指10進制數字),而一個8字節(jié)的無符號整數則能精確到19位數字。浮點數能表示更大的范圍，但精度較低。當表示的整數很大的，則可能存在誤差。一個8字節(jié)的雙精度浮點數可表示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位數字。計算機中，我們常用的最小數據存儲單位是字節(jié)，C語言稱之為char,多個字節(jié)可表示一個更大的存儲單位。習慣上，兩個相鄰字節(jié)組合起來稱作一個短整數，在32位的C語言編譯器中稱之為short,匯編語語言一般記作word,4個相鄰的字節(jié)組合起來稱為一個長整數，在32位的C語言編譯器中稱之為long,匯編語言一般記作DWORD。在計算機中，按照權的不同，數的表示可分為兩種，2進制和10進制，嚴格說來，應該是2^k進制和10^K進制，前者具占用空間少，運算速度快的優(yōu)點。后者則具有容易顯示的優(yōu)點。我們試舉例說明：
例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。最后一種表示法，在新手寫的計算大數階乘程序最為常見，但計算速度卻是最慢的。使用更大的基，可以充分發(fā)揮CPU的計算能力，計算量將更少，計算速度更快，占用的存儲空間也更少。
?
2.4 大尾序和小尾序，我們在書寫一個數時，總是先寫權較大的數字，后寫權較小的數字，但計算機中的數并不總是按這個的順序存放。小尾（Little Endian）就是低位字節(jié)排放在內存的低端，高位字節(jié)排放在內存的高端。例如對于一個4字節(jié)的整數0x12345678,將在內存中按照如下順序排放, Intel處理器大多數使用小尾(Little Endian)字節(jié)序。
Address[0]: 0x78
Address[1]: 0x56
Address[2]: 0x34
Address[3]:0x12
大尾（Big Endian）就是高位字節(jié)排放在內存的低端，低位字節(jié)排放在內存的高端。例如對于一個4字節(jié)的整數0x12345678,將在內存中按照如下順序排放, Motorola處理器大多數使用大尾(Big Endian)字節(jié)序。
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 大數的鏈式存儲法
如果我們搜索大數階乘的源代碼，就會發(fā)現，有許多程序采用鏈表存儲大數。盡管這種存儲方式能夠表示大數，也不需要事先知道一個特定的數有多少位有效數字，可以在運算過程中自動擴展鏈表長度。但是，如果基于運算速度和內存的考慮，強烈不建議采用這種存儲方式，因為：
?
?
1. 這種存儲方式的內存利用率很低。基于大數乘法的計算和顯示，一般需要定義雙鏈表，假如我們用1個char表示1位十進制數，則可以這樣定義鏈表的節(jié)點：
struct _node
{
struct _node* pre;
struct _node* next;
char n;
};
當編譯器采用默認設置，在通常的32位編譯器，這個結構體將占用12字節(jié)。但這并不等于說，分配具有1000個節(jié)點的鏈表需要1000*12字節(jié)。不要忘記，操作系統或者庫函數在從內存池中分配和釋放內存時，也需要維護一個鏈表。實驗表明，在VC編譯的程序，一個節(jié)點總的內存占用量為 sizeof(struct _node) 向上取16的倍數再加8字節(jié)。也就是說，采用這種方式表示n位十進制數需要 n*24字節(jié)，而采用1個char型數組僅需要n字節(jié)。
?
?
2采用鏈表方式表示大數的運行速度很慢.
2.1如果一個大數需要n個節(jié)點，需要調用n次malloc(C)或new(C++)函數，采用動態(tài)數組則不要用調用這么多次malloc.
?
2.2 存取數組表示的大數比鏈表表示的大數具有更高的cache命中率。數組的各個元素的地址是連續(xù)的，而鏈表的各個節(jié)點在內存中的地址是不連續(xù)的，而且具有更大的數據量。因此前者的cache的命中率高于后者，從而導致運行速度高于后者。
?
2.3對數組的順序訪問也比鏈表快，如p1表示數組當前元素的地址，則計算數組的下一個地址時一般用p1++,而對鏈表來說則可能是p2=p2->next,毫無疑問，前者的執(zhí)行速度更快。
?
?
近似計算之一
?
在＜階乘之計算從入門到精通－菜鳥篇＞中提到，使用double型數來計算階乘，當n>170,計算結果就超過double數的最大范圍而發(fā)生了溢出，故當n＞170時，就不能用這個方法來計算階乘了，果真如此嗎？No,只要肯動腦筋,辦法總是有的。
通過windows計算器，我們知道，171！＝1.2410180702176678234248405241031e+309，雖然這個數不能直接用double型的數來表示，但我們可以用別的方法來表示。通過觀察這個數，我們發(fā)現，這個數的表示法為科學計算法，它用兩部分組成，一是尾數部分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;
double MAX_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個字節(jié)(64bit)，我們把最低位記作bit0,最高位記作bit63,則一個浮點數各個部分定義為:第一部分尾數：bit0至bit51，共計52bit，第二部分階碼:bit52-bit62,共計11bit，第三部分符號位:bit63，0:表示正數，1表示負數。如一個數為0.xxxx * 2^ exp,則exp表示指數部分，范圍為－1023到1024，實際存儲時采用移碼的表示法，即將exp的值加上0x3ff，使其變?yōu)橐粋€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) //繼續(xù)相乘可能溢出，調整之
{
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;
}
?
程序優(yōu)化的威力：
程序4是一個很好的程序，它可以計算到1千萬（當n更大時，p->n2可能溢出）的階乘，但是從運行速度上講，它仍有潛力可挖，在采用了兩種優(yōu)化技術后，（見程序5），速度竟提高5倍，甚至超出筆者的預計。
第一種優(yōu)化技術，將頻繁調用的函數定義成inline函數，inline是C++關鍵字，如果使用純C的編譯器，可采用MACRO來代替。
第二種優(yōu)化技術，將循環(huán)體內的代碼展開，由一個循環(huán)步只做一次乘法，改為一個循環(huán)步做32次乘法。
實際速度：計算10000000！，程序4需要0.11667秒，程序5只需要0.02145秒。測試環(huán)境為迅馳1.7G,256M RAM。關于程序優(yōu)化方面的內容，將在后續(xù)文章專門講述。下面是被優(yōu)化后的部分代碼，其余代碼不變。
?
程序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次方
轉自：liangbch@263.net
轉載于:https://blog.51cto.com/charlie555/783098
總結
- 上一篇： Java web开发——Servlet详
- 下一篇： 单片机快速开平方的算法