大数运算之C 语言大数演算法

（一）
- 從字串讀取
- 加法
- 減法
- 乘法
- 比較大小
- 除法
(二)
- 定義 macro
- 加法
- 減法
- 乘法
- 比較
- 除法
(三)
- [0] 一次性進位問題 - 不用
- [1] 改善資料結構
- [2] 從字串讀入大數
- [3] 輸出大數結果
- [4] 大數加法
- [5] 大數減法
- [6] 大數乘法
- [7] 大數乘整數
- [7] 大數除整數
(四)
- [1] 算有幾位數
- [2] 精算值 - 暴力法
- [4] 精算值 - 質因數分解
- [5] 判斷尾數有幾個零
- [6] 其他假大數階乘問題

引文：

本文系一台湾计科大佬所写，由于本人实在是感觉写得非常好，
于是翻译工作也没有做（害怕和原文有出入）。

（一）

大數問題我認為是新手必練題型之一，但實在沒時間去 implement 做 sample code，把目前知道的作法先且粗略紀錄下來，以下探討「暫」以「無號大數」為標的，
下面的 code 憑印象之演示。

資料結構

一般可以用 link 去做大數，只要每多一 digit 就去新增一 node，但效能變得非常差，不論在釋放記憶體還是做 travel 時都很糟，一般還是以 array 居多。
array 方式是將每個位數視為一個整數存進去，假設要存的數字為 763421 ，則儲存方式為

arr[0]=1, arr[1]=2, arr[2]=4, arr[3]=3, arr[4]=6, arr[5]=7。

再來是該大數要陣列大小要多大、資料型態為何？一般在探討時都簡化問題，假定陣列大小是固定的，也就是先定義一 macro - MAX。

#define MAX 200
int  big[MAX];

這樣便可存下 200 個位數之數字，超過的話就把 MAX 再調大。但說明了這是「簡化問題」，實際上還是用 malloc/ new 把內容放在 heap 裡面做動態管理好點。

從字串讀取

假設一傳入字串 char str[] = “12345”，代表著

str[0]='1', str[1]='2', str[2]='3', str[3]='4', str[4]='5'

但實際上欲存的大數陣列是從高位存到低位，也就是

big[4]=1, big[3]=2, big[2]=3, big[1]=4, big[0]=5

於是必須反著讀 str，同時還考慮到字串存的為 ascii value (假定為 ascii value)，

‘0’ 存 48，‘1’ 存 49… 故存到大數裡面時再扣掉 offset value = ‘0’

#define MAX 200
int big[MAX];void read_from_string(int* big, char* str)
{int i=0, len=strlen(str);memset(big, 0, sizeof(int)*MAX);for(i=len-1; i>=0; --i)big[len - i - 1] = str[i] - '0';
}

memset 指令是用來把 big 陣列先全都清零，要用 for loop 慢慢跑也行。從整數 (正整數) 轉到大數的話就方便多了，觀念是將 x%10 之結果丟給大數 (取最後一位)，再將 x=x/10，如果最後 x 除到 0 時，就代表轉完了。

void read_from_int(int *big, int x)
{int i=0;memset(big, 0, sizeof(int)*MAX); // 將 big 先清 0while(x!=0){big[i++]=x%10;x/=10;}
}

顯示

顯示時一律從最高位開始顯示，但必須注意，若大數存的是 3 ，整個陣列會是 00…03，在 3 前面有 199 個 0，所以必須先把前導 0 全都拿掉後再行輸出。

void big_print(int *big)
{int i=MAX-1;for(i=MAX-1;i>=0 && big[i]==0; --i);while(i>=0) printf("%d", big[i]), --i;
}

但上面的 code 有問題，原因在於這種寫法，數值 0 就不會輸出，於是把去前導0之判別式小改一下就可用。

void big_print(int *big)
{int i=MAX-1;for(i=MAX-1;i>0 && big[i]==0; --i);while(i>=0) printf("%d", big[i]), --i;
}

這樣便可正常顯示大數。

加法

大數加法沒什麼技巧，和一般的小學的直式減法沒什麼兩樣，假設 rst = a+b

先算 rst[i] = a[i] + b[i] + carry ，其中 carry 代表 “上一筆運算之進位”。

若 rst[i] >=10 時，將 carry 設 1，再把 rst[i] 扣 10；

另一思維是，直接套公式 carry = rst[i] / 10 ，rst[i] = rst[i]%10，這裡採此方式。

void big_add(int *rst, int *a, int *b)
{int i, sum, carry;for(carry=i=0; i<MAX; ++i){rst[i] = a[i] + b[i] + carry;carry = rst[i] / 10;rst[i]%=10;}
}

減法

加法是用 carry，減法是用 borrow，和上敘思考一樣，以直式減法概念下去做。

void big_sub(int *rst, int *a, int *b)
{int i, borrow;for(borrow=i=0; i<MAX; ++i) {rst[i] = a[i]-b[i]-borrow; // 扣上一次借位if(rst[i]<0) {borrow=1, rst[i]+=10; // 需借位} else {borrow=0; // 無需借位}}
}

乘法

照直式之乘法規則，a[i] * b[j] 最後結果是放到 rst[i+j] 裡面去。

正常是每乘完一次之後，去判斷 rst[i+j] 是否大於 10，如果大於 10 就進行進位調整，但其實不用這麼麻煩，

只要整個乘法都做完後，再對 rst 做一次性調整就好。

有個細節稍注意，如果 a[i] 為 0 的話，就直接做下一個，節省時間

void big_mul(int *rst, int *a, int *b)
{int i, j, carry;memset(rst, 0, sizeof(int)*MAX); // 先清0for(i=0; i<MAX; ++i) {if(a[i]==0) continue;for(j=0; i+j<MAX; ++j)rst[i+j]+= a[i]*b[j];}// 一次性調整for(carry=i=0; i<MAX; ++i) {rst[i]+=carry;carry = rst[i] / 10;rst[i] %= 10;}}

比較大小

比較大小時從最高位開始一位逐一比較，code 當參考，

a>b 傳 + ； a<b 傳 -；a=b 傳 0int big_compare(int *a, int *b)
{int i=MAX-1;while(i>0 && a[i]==b[i]) --i;return a[i]-b[i];
}

除法

除法很麻煩，同時要做除法前，建議再多做下面幾個函式

void big_addi(int *big, int x); // 大數+正整數
void big_subi(int *big, int x); // 大數-正整數
void big_muli(int *big, int x); // 大數*正整數
void big_divi(int *big, int x); // 大數/正整數

要完成 big_div，big_muli 一定要做，big_divi 看個人有沒有興趣做，

其他的用不用得到不一定，要看個人 code 怎麼寫。

欲求 rst = a/b (整數部份)，一種低效但簡單的方式是，先將 rst 歸零，之後 a 不停減 b，rst 不停加一，直到 b>a 為止，rst 即為所求。

另一種方式概念，先假設 a=123456， b=345，a 用了 6 位數， b 用了 3位數，

直接先把 b 移成 6 位數，得到移動了 delta= (dig of a) - (dig of b) = 3 位數。

a = 123456 ，b = 345000 [ 商=0，餘=123456 ]，b 右移一個 digita = 123456， b = 34500 [ 商=3，餘=19956，b 右移一個 digit，

依此類推下去，直到 delta 變負數為止。以此為概念，整個變化如下。


delta=3, a = 123456，b = 345000
Q = a/b = 0，a = a - Q*b = 123456，
b 除 10 (記憶體右搬1個元素)
rst[delta] = Q ，rst[3]=0。
delta=2, a = 123456，b = 34500 [ans=0]
Q = a/b = 3，a = a - Q*b = 19956，
b 除 10 (記憶體右搬1個元素)
rst[delta] = Q ，rst[2]=3。
delta=1, a = 19956， b = 3450 [ans=3]
Q = a/b = 5，a = a - Q*b = 2706，
b 除10 (記憶體右搬1個元素)
rst[delta] = Q ，rst[1]=5。
delta=0, a = 2706， b = 345 [ans=35]
Q = a/b = 7，a = a - Q*b = 291，
b 除10 (記憶體右搬1個元素)
rst[delta] = Q ，rst[0]=7。

所以最後得到之 ans=357，而 a=291 就是餘數，
驗證一下： 123456 = 345 * 357 + 291，無誤。

可想一下如果是 10^ 8 除以 1，原本慢慢扣要執行 10^8 次，

但這種方式在外回圈執行 8 次，即使用效率最差的線性查找法找商，

內回圈執行 10 次，差別變成 10^8 與 10*8 次的差別。

類似的方法不只一種，上面這想法是死 K 算盤本想到的，

而算盤本裡面還有其他三、四種除法演算法 (二進位)，

其中還包含了非回復式除法，有興趣可自行看。

上述下來，要將除法做得像樣，又要有額外幾個副函式

陣列左移
陣列右移
大數二分法

陣列左移、陣列右移其中一個可以用 memmove (筆者並沒很喜愛用此指令)，較習慣全用 for loop 慢慢 assigned；但問題較大的是二分法的實作，有心想做份好一點的大數，這裡將會是第一個瓶頸難題。

(二)

一個陣列元素裡面只用了 0~9，10 進制，這數值大小其實用不到 int big[MAX] 方式宣告，只要用 char big[MAX] 就行了，這麼一來也只是省空間，對速度沒助益。

由於 int 可表達範圍約 0~ 2.1 * 10 ^ 9 ，若考慮乘法問題，一個數最大可表示到 sqrt (2.1 * 10^9) 約 65536 左右，只用 10 進位方式實在是太慢了，於是改成「萬」進位，原因為 10000 為最接近 65536 之 10^n 數字。

先提，大數演算法效率如何，一半關鍵因素在於儲存進制之方式。

定義 macro

接下來的版本，建議先全用 #define 方式定義起來。

#define MAX 200
#define CAP 10000
#define CLEN 4
typedef int  btype;

上面有些定義是有爭議的，特別是 CLEN 部份，這裡定義 CLEN 是為了到時輸出用的，要使得 10^CLEN = CAP，但也有些人認為這到時候再用算的就行了，筆者認為既然常用，那就一次定出來就好。而 typedef 部份未必會做，但通常是做了較佳，因擴充性會較好，因 btype 代表的是該大數陣列要用怎樣的資料型態做儲存。

為了程式碼說明方便，這裡只有 typedef 那裡不採用，並採用的是萬進位 (即 CAP 為 10000)，

也假設主程式裡面宣告為 int a[MAX], b[MAX]; 作為大數之用。

從字串讀取

之前有提過，要從字串轉存成大數時，必須從後面開始存起，在萬進位也一樣，考慮 char str[] = “12345” 時，

big[1] = 1，big[0] = 2345。

一樣都是從最後面讀出來，只是多了一個 pwr 做累計，每次都乘 10 ，乘到 CAP 時就回到 1。

void big_read(int *big, const char* str)
{int len = strlen(str);int i, j, pwr;memset(big, 0, MAX*sizeof(big[0]));for(j=0, i=len-1; i>=0; ++j)for(pwr=1; pwr!=CAP && i>=0; pwr*=10, --i)big[j]+=(str[i]-'0')*pwr;
}

顯示

考慮大數為 big[2]=1, big[1]=2, big[2]=3，它所代表的不是 123，而是 100020003。

一開始先去前導 0 ，找到第一個非零之元素，直接輸出；之後的元素全都以 0 填滿其 CLEN 長度。

void big_print(int *big)
{int i=MAX-1;while(i>=0 && big[i--]==0);++i;printf("%d", big[i--]);while(i>=0) printf("%0*d", CLEN, big[i--]);
}

加法

仿直接加法沒什麼兩樣，只是把 10 進位換成 CAP 進位而已。

void big_add(int *rst, int *a, int *b)
{int i, sum, carry;for(carry=i=0; i<MAX; ++i){rst[i] = a[i] + b[i] + carry;carry = rst[i] / CAP;rst[i]%=CAP;}
}

減法

仿直接減法，繼續把 10 進位換成 CAP 進位。

void big_sub(int *rst, int *a, int *b)
{int i, borrow;for(borrow=i=0; i<MAX; ++i) {rst[i] = a[i]-b[i]-borrow; // 扣上一次借位if(rst[i]<0) {borrow=1, rst[i]+=CAP; // 需借位} else {borrow=0; // 無需借位}}
}

乘法

這裡再提一下一次進位要注意的地方。

上篇 [大數] C 語言大數演算法 for beginner 裡面，筆者用的是一次進位法則，

可以那麼用是因為原本是 10 進位，要在同一個 rst[i+j] 擠到爆的可能性非常小

(可以推，有興趣可自己推)，但這次已經是用萬進制，若最後 rst[i+j] 是

(9999 * 9999) 累加起來，可容許 rst[i+j] = a[i] * b[j] 約有 20 次這種情況，
超過這次數的話，到時若不先進位，結果肯定會爆。解決方式有兩種，

一種是將原本的萬進制，調成千進制；另一種是放棄一次進位，直接改逐次進位。

這裡程式碼乃以一次進位為主。

void big_mul(int *rst, int *a, int *b)
{int i, j, carry;memset(rst, 0, sizeof(int)*MAX); // 先清0for(i=0; i<MAX; ++i) {if(a[i]==0) continue;for(j=0; i+j<MAX; ++j)rst[i+j]+= a[i]*b[j];}// 一次性調整for(carry=i=0; i<MAX; ++i) {rst[i]+=carry;carry = rst[i] / CAP;rst[i] %= CAP;}
}

比較

這部份一模一樣 (愈寫愈心虛，一直參考上一篇)

int big_compare(int *a, int *b)
{int i=MAX-1;while(i>0 && a[i]==b[i]) --i;return a[i]-b[i];
}

除法

關鍵都在 implement 大數二分法。求單一商數時，

10 進制裡最大差別是 4/10；而萬進制裡最大差別變成了 13 / 10000；

這裡應看出效果到底差在哪了。Implement 等有時間補上。

另，取模運算，實際就是除法運算，只是最後傳回值是餘數不是商數。

其他改善

截至目前為止，以 10^n 進制而言，所有提過的地方有個地方可改善 - 紀錄使用元素個數

假設 rst = a* b，其中 a 佔 n1 位數，b 佔 n2 位數，n = max(n1, n2)，

得到 rst 結果為 x 位數，若增加一數值，紀錄目前這個大數已用了多少元素，

則速度會再快！可應用以下特性

(1) 加法最多只需執行 n 次，x 最大值為 n+1

(2) 減法最多只需執行 n 次，x 最大值為 n

(3) 乘法最多只需執行 min (2n, MAX) 次，x 最大值為 min (2n+1, MAX)

(4) 除法最多只需執行 n 次，x 最大值為 n

(5) 比較 a, b 大小，就直接從 n1, n2 開始比較

(6) 做輸出時直接從 n1 / n2 / n 開始輸出，節省去前導 0 之動作。

(7) size==0 && rst[0]==0 ，可判斷 0；同樣也可判斷 1。

  如此對乘法速度有所幫助，也可防除零問題。

上面 (1) ~ (4) 這些位數關係要推導我懶 (會推，但懶得推)，

不過用「觀查法」就可得到結論，考慮 10 進位模式，對於各運算考慮

9999 (+-*/ ) 99999999 (+-*/ ) 11 (+-*/ ) 99991 (+-*/ ) 1

之四種特例，應不難觀查出結果。

(三)

[0] 一次性進位問題 - 不用

這問題筆者細思過，我們以一般加法為例，一種是逐次進位，一種是一次性調整，兩種寫法如下所示。

// 逐次進位
void big_add(int * A , int * B , int * C)
{int i, j, c;for(c=i=0; i<SIZE; ++i) {C[i] = A[i] + B[i] + c;if(C[i]>=10) c=1, C[i]-=10;else c=0;}
}// 一次進位
void big_add(int * A , int * B , int * C)
{int i;for(i=0; i<SIZE; ++i)C[i] = A[i] + B[i];for(i=0; i<SIZE-1; ++i) {C[i+1] += C[i]/10;C[i] %= 10;}
}

在加、減法可能比較看不出來，但在乘法的時候，一次性進位缺點會不少。在萬進位情況下，使用一次性進位容易溢位；65536進位下更容易溢位，所以筆者打算全都用逐次進位之方式。但對乘法而言，似乎不多人會寫逐次進位，所以寫起來有三個 loop，其實一樣只要兩個 loop 就行了。

void big_mul(int * A , int * B , int * C)
{int i, j, ij, carry=0;memset( (void*)C, 0, BSIZE);for(i=0; i<SIZE; ++i){for(j=0; (ij=i+j) < SIZE; ++j){C[ij] += carry + A[i] * B[j];carry = C[ij] / 10;C[ij] %= 10;}}
}

上述逐次進位方式並不會造成任何 ov 可能性，接下來我們全都用逐次進位方式，不再一口氣進位。

[1] 改善資料結構

一般我們從字串「讀入」開始，是直接轉成大數，這其實也沒問題，但一般大數都是直接以 int arr[SIZE] 方式表達，這也合理。而一份簡單降低一點 BigO 常數的技巧是 : 紀錄用了多少位數。

ok，我不確定這是十分有效的，只是直覺它是可以有效的。為方便說明，我還是給靜態陣列大小，不做動態配置。但事實上這篇講完之後，其實要實作動態自動調整陣列大小，已不是問題，只是為避免內文冗長，所以用靜態說明。

在定義大數時候，我給了一份 struct

#define BIG_SIZE 20
#define BIG_CAP  10
#define BIG_LEN  1typedef struct tagBig{int arr[BIG_SIZE]; // 大數陣列int dig;           // 使用位數
}Big;

dig 代表這個陣列裡面實際使用的位數，試想一下，如果是 00001 + 000002 ，實際上是只要加一個位數就好 (他們兩個的 dig 都是 1)，但一般傳統給的是把整個 BIG_SIZE 都算完。換而言之，傳統之方式， BIG_SIZE 給的愈大，計算愈久，但表示範圍廣；BIG_SIZE 給的愈小，計算時間不會拉長，表示範圍長。

另外這份程式其實可以再加一個 int cur_size，做動態方式調整，這裡就不再敘述。而傳遞方式以 pass by pointer ( C++ 可用 pass by reference) 方式傳遞，否則會造成 struct 裡之 array 初始化時間過長問題。

接下來所有的作法都只要稍想一下就可以了。然後有四份 macro 先定義下來，方便使用。

#define MAX(a,b) ( ((a)>=(b)) ? (a) : (b) )
#define MIN(a,b) ( ((a)<=(b)) ? (a) : (b) )
#define BSIZE(n) ( sizeof(int) * n)
#define SWAP(TYPE, a, b) { TYPE t; t=a; a=b; }

[2] 從字串讀入大數

和之前程式碼雷同，不過為了省時間，並沒有調用 memset、填零之動作。但多了「去字串前導零」動作，另外也紀錄了使用了多少位數。

void big_from_string(Big * big , const char * str)
{int i, j, len, pwr;int * arr = big->arr;// 去前導0while(*str=='0' && *str!='\0') ++str;// 轉入大數len = strlen(str);arr[0] = 0;for(j=0, i=len-1; i>=0; ++j, arr[j]=0)for(pwr=1; pwr!=BIG_CAP && i>=0; pwr*=10, --i)arr[j]+=(str[i]-'0')*pwr;// 計算位數big->dig = j;
}

換而言之，即使 BIG_SIZE = 100，實際使用 dig = 10，但其實不會把其他 90 個大小填入 0 ，以期節省時間。另注意件事，假設 capacity = 10 (10進位)，當轉出來之數值為 0 時， dig = 0；1~9 時，dig=1；10~19 時，dig=2… 依此類推。這麼設計理由，是認為在做以下的運算會比較方便。

[3] 輸出大數結果

若 dig = 0 時直接輸出 0 ，否則從 dig 開始輸出即可。

// 輸出
void big_print(Big * big)
{int i = big->dig;//printf("(dig = %d)", i);if(!i) printf("0");else {int *arr = big->arr;printf("%d", arr[--i]); // 輸出第一位數字for(--i; i>=0; --i) // 輸出其他數字printf("%0*d", BIG_LEN, arr[i]);}
}

這裡給一份測試碼


int main()
{char s1[2000], s2[2000];Big A, B;    while(scanf("%s%s", s1,s2)==2){big_from_string(&A, s1); printf("\nA = "), big_print(&A);big_from_string(&B, s2); printf("\nB = "), big_print(&B);}return 0;
}

加上 stdio.h、stdlib.h、string.h，上面四段合起來可做一份測試了。

[4] 大數加法

(1) 結果位數 : 結果是 Max(A的位數 , B 的位數) + 1，最後的 +1 不一定會有，如果最後還有進位的話就有 +1 。

(2) 偷用 memcpy : 試想一下 123456 + 14 的時候，其實真正在計算的只有 3 位數，56+14，還有一個 carry + 4，其他的是用 memcpy 直接照抄。

(3) 特殊值判斷 : +0 的話就直接 memcpy 。

// C = A + B
void big_add(Big * BigA , Big * BigB , Big * BigC)
{int i, carry;int len_a = BigA->dig, len_b = BigB->dig;int *A = BigA->arr, * B = BigB->arr;int *C = BigC->arr;// 特殊值判斷if(len_a==0 && len_b==0) { // A = B = 0BigC->dig = 0;return ;} else if(len_a==0) { // A = 0, C = A + B = Bmemcpy((void*)C, (void*)B, BSIZE(len_b));BigC->dig = len_b;return ;} else if(len_b==0) { // B = 0, C = A + B = Amemcpy((void*)C, (void*)A, BSIZE(len_a));BigC->dig = len_a;return ;}// 保證 len_a >= len_bif(len_a < len_b) {SWAP(int , len_a, len_b);SWAP(int * , A, B);}A[len_a]=0, ++len_a;// 加到 len_bfor(carry=i=0; i<len_b; ++i) {C[i]  = A[i] + B[i] + carry;carry = C[i]/BIG_CAP;C[i]%=BIG_CAP;}// 若有進位，將連續之 BIG_CAP-1 全改 0if(carry)while(A[i]==BIG_CAP2)C[i]=0, ++i;C[i] = A[i] + carry, ++i; // 補上最後一個進位// 剩下之 A 全丟到 Cif(i<len_a)memcpy( (void*)(C+i), (void*)(A+i), BSIZE(len_a-i));// 調整 len_a, 寫回結果BigC->dig = len_a - !C[len_a-1];
}

[5] 大數減法

(1) 結果位數 : 結果是 0~Max(A的位數 , B 的位數)，這裡最後是用 polling 方式去查。

(2) 偷用 memcpy : 和上面一樣

(3) 特殊值判斷 : -0 的話就直接 memcpy 。

(4) 我們不考慮結果為負的時候之處理。

// C = A - B , 逐次借位
void big_sub(Big * BigA , Big * BigB , Big * BigC)
{int i, borrow;int len_a = BigA->dig, len_b = BigB->dig;int *A = BigA->arr, * B = BigB->arr;int *C = BigC->arr;// 特殊值判斷if(len_a==0 && len_b==0) { // A = B = 0BigC->dig = 0;return ;} else if(len_a==0 || len_b > len_a) { // A = 0, C = A - B = -B// 直接設 0 不處理 < 或以十補數方式處理 >BigC->dig = 0;return ;} else if(len_b==0) { // B = 0, C = A - B = Amemcpy((void*)C, (void*)A, BSIZE(len_a));BigC->dig = len_a;return ;}// 先做 len_b 個減法for(borrow=i=0; i<len_b; ++i) {C[i] = A[i] - B[i] - borrow;if(C[i] < 0) C[i]+=BIG_CAP, borrow=1;else borrow=0;}// A < B , 溢位處理if(len_a<=len_b && borrow){BigC->dig = 0;return ;}// 把借位處理完, 遇到 0 直接給 BIG_CAP-1if(borrow)while(A[i]==0)C[i]=BIG_CAP2, ++i;C[i]=A[i]-borrow, ++i; // 補上最後借位// 剩下的 A 全給 Cif(i < len_a)memcpy( (void*)(C+i), (void*)(A+i), BSIZE(len_a-i));//回頭查剩幾位, 寫回 Cfor(i=len_a-1 ; i>0 && C[i]==0 ; --i);if(i==0) BigC->dig = (C[0]!=0);else BigC->dig = i+1;
}

[6] 大數乘法

大數乘法這是效率低的版本 O(n^2)，另外有 O(n^log(3)) 、O(nlogn) 版本。

(1) 結果位數 : 正確結果應該是 (A位數 + B位數 -1 ) ~ (A 位數 + B 位數 + 1)，那個 -1 蠻多人略過、覺得不對，舉個例，２　＊　４　＝　８，就是有 -1 的情況。

(2) 這裡不能用 memcpy 加速。

(3) 特殊值判斷 : 有 0 設 0 ，有 1 設 A/B。

(4) 這裡沒處理位數不夠放之情況，但註解裡會提出要注意。

(5) 拿掉 carry 變數。

//
// 大數乘大數, 逐次進位, 無 ovvoid big_mul(Big * BigA , Big * BigB , Big * BigC)
{int i, j,ij;int len_a = BigA->dig, len_b = BigB->dig;int *A = BigA->arr, *B = BigB->arr, *C=BigC->arr;int len_c = len_a + len_b + 1;// 特殊值判斷if(!len_a || !len_b) { // A==0 , B==0, C = A * B = 0BigC->dig = 0;return ;} else if(len_a==1 && A[0]==1) { // A==1, C = A * B = Bmemcpy( (void*)C, (void*)B , BSIZE(len_b));BigC->dig = len_b;return ;} else if(len_b==1 && B[0]==1) { // B==1, C = A * 1 = Amemcpy( (void*)C, (void*)A , BSIZE(len_a));BigC->dig = len_a;return ;}// 做 arr-size , len_c 之防呆// 開始做乘法memset((void*)C, 0, BSIZE(len_c));for(i = 0; i < len_a ; ++i) {for( j=0 ; j<len_b ; ++j) {ij = i+j;C[ij] += A[i] * B[j];C[ij+1] += C[ij]/BIG_CAP;            C[ij]%=BIG_CAP;}}// 修正 len_cif(C[len_c-1]==0) --len_c;if(C[len_c-1]==0) --len_c;BigC->dig = len_c;
}

[7] 大數乘整數

大數乘整數習慣上還是用一次性進位。

void big_mul_int(Big * BigA , int numB , Big * BigC)
{int i, carry;int *A = BigA->arr, * C = BigC->arr;int len_a = BigA->dig;// 特殊值判斷if( !len_a || !numB) { // 設0BigC->dig = 0;return ;} /* 其他 A, numB = 1 之判斷 */// 開始做乘法for( carry = i = 0; i<len_a; ++i){C[i] = A[i] * numB + carry; // (!!!) 注意這行carry = C[i] / BIG_CAP;C[i] %= BIG_CAP;}// 再調整進位while(carry) {C[i] = carry % BIG_CAP;carry /= BIG_CAP;++i;}BigC->dig = i;
}

注意到上面的算法其實會有溢位的可能性，假設為 10 進位時，當 numB = 2147483647 就溢位了，會溢位就不叫大數運算了。但相對的，這種算法在針對 numB 是小量的，比如說在 BIG_CAP 以內，且 numB 也在 BIG_CAP 內的話，其實是不會有溢位問題的。

上面這段碼是拿來做大數除大數用的乘法，因在過程中的設計的確用這段碼就行了。要考慮不溢位的話其實很簡單，再將 numB 轉成一個小量的大數 (BIG_CAP = 10 的時候了不起陣列大小也才 10 而已) 即可。

void big_mul_int(Big * BigA , int numB , Big * BigC)
{int i, j, ij, carry;int *A = BigA->arr, * C = BigC->arr;int len_a = BigA->dig;int len_b, len_c;int B[20];// 特殊值判斷if( !len_a || !numB) { // 設0BigC->dig = 0;return ;} /* 其他 A, numB = 1 之判斷 */// 將整數 numB 轉到 array 裡去for(len_b=0; numB; ++len_b, numB/=BIG_CAP)B[len_b] = numB%BIG_CAP;len_c = len_a + len_b +1;// 開始做乘法memset((void*)C, 0, BSIZE(len_c));for(i = 0; i < len_a ; ++i) {for(carry = j=0 ; j<len_b ; ++j) {ij = i+j;C[ij] += A[i] * B[j];C[ij+1] += C[ij]/BIG_CAP;            C[ij]%=BIG_CAP;}}// 修正 len_cif(C[len_c-1]==0) --len_c;if(C[len_c-1]==0) --len_c;BigC->dig = len_c;
}

[7] 大數除整數

實際上考慮直式除法便行。

///
// 大數除整數, 無溢位問題
// BigA : 被除數, numB : 除數 , C : 商
// 傳回 : 成功傳回餘數(>=0), 失敗傳回負值
//
int big_div_int(Big * BigA , int numB, Big * BigC)
{int * A = BigA->arr, * C = BigC->arr;int i, len_a = BigA->dig;int r, len_c; //餘數// 特殊值判斷if(numB<=0) { // 除零錯誤, 除到負數BigC->dig=0;return -1;} else if(!len_a) {// 分子為 0BigC->dig = 0;return 0; // 餘數設 0} else if(numB==1) { // 分母為 1memcpy( (void*)BigC, (void*)BigA, len_a);BigC->dig = len_a;return 0; // 餘數為 0}// 開始做除法r=0;for(i=len_a-1; i>=0; --i) {r = r *BIG_CAP + A[i];C[i] = r / numB;r %= numB;}for(i=len_a-1; i>0 && C[i]==0; --i);BigC->dig = i + !(i==0 && C[0]==0);//     if(i==0 && C[0]==0) BigC->dig = 0;//     else BigC->dig=i+1;return r; // 傳回餘數
}

(四)

[1] 算有幾位數

嚴格講起，這是個假大數問題。

可用斯特靈公式逼近，但拿到的是一個近似值。要細算的話推導如下。

n! = n * (n-1) * (n-2) * .... * 2 * 1 ，算位數取 log10，log10(n!) = log10(  n * (n-1) * (n-2) * .... * 2 * 1  ) ，log10(n!) = log10(n) + log10(n-1) + log10(n-2) + .... + log10(2) + log10(1) ，

最後結果出來再取 ceil 。

雖浮點數有誤差，但有效是 53 位數，應付到幾十萬階乘是絕對足夠 (速度變慢倒是真的)

[2] 精算值 - 暴力法

一樣用大數去做，先放上暴力法

// --------------------------------------
// 精算 n! 之值, 傳回有幾位數 , slow
int fact(int n, int * Big, int size)
{int i, j, dig=1, carry=0;int mini, pwr, num;if(n<0) {Big[0]=0; return -1;}else if(n==0 || n==1) {Big[0]=1; return 1;}Big[0] = 1;for(i=2 ; i<=n; ++i){// 進行大數乘法if(size < dig) mini = size;else mini = dig;for(carry =j=0 ; j<mini; ++j) {Big[j] = carry + Big[j] * i;carry  = Big[j] / BIG_CAP;Big[j]%=BIG_CAP;}// 調進位while(carry && j<size){Big[j] = carry % BIG_CAP;carry /= BIG_CAP;++j;}dig=j;}// 計算有幾個位數dig = (j-1)*BIG_LEN;for(num=Big[j-1]; num ; num/=10) ++dig;return dig;
}

注意中間使用的 dig 計算，拿來加快乘法速度。試過，用萬進位可行。

[4] 精算值 - 質因數分解

10! = 10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2，把每個數做質因數分解，得到，
10! = 2^8 * 3^4 * 5^2  * 7，接下來就好辦了。

找一個階乘的質因式分解，其實有技巧，但需要較大的空間(一般題目應該都算夠用，目前看過找較大的是上萬)。一開始最好先建質數表，可以普通篩法或線性篩法都行 ( 幾十萬以內都很快)。篩法得到質數為 {2,3,5,7} 接著…

(1) 一開始拿質數 2 出來，

i = 1 , cnt = 0;
i = i * 2 = 2 , 10 / 2 = 5 , cnt+=5, cnt=5
i = i * 2 = 4 , 10 / 4 = 2 , cnt+=2, cnt=7
i = i * 2 = 8 , 10 / 8 = 1, cnt+=1 , cnt=8

(2) 再拿出質數 3 出來

i = 1 , cnt = 0;
i = i * 3 = 3 , 10 / 3 = 3 , cnt+=3, cnt=3
i = i * 3 = 9 , 10 / 9 = 1 , cnt+=1, cnt=4
i = i * 3 = 27 , 結束

其他拿質數 5、質數 7 都一樣。另外可以額外處理地方是在質數 7 。

因為是 10!，10/2 = 5，而 7 > 5，故可以判定，7 一定只有一個，就是自己，其它的絕不可能再分解下去 (很簡單的數論應用而已)。換比較通用的講法是：在求 n 階乘，做上述動作時，只要做到 n/2 便可，其它大於 n/2 的全都乘上去就行了。

fast-power 有個地方要注意，由於出來的結果可能會 overflow，故建議最好事先判斷有沒有可能會 ov。沒 ov 的話用大數做 fast_power 有點浪費，有 ov 就一定得用大數做 fast-power。

再來是看想不想進一步優化。再以 10! 為例。10! = 2^8 * 3^4 * 5^2 * 7^{1，考慮前兩項就好，(2}8) * (3^4) = (6^4) * (2^4)，效果略有所不同。拆法有很多種 (特別注意到這個例子的次方剛好是 8 4 2 1 )，考慮溢位、速度之間做取捨之衡量，要設計好不容易。其他的不再贅述。

[5] 判斷尾數有幾個零

這是個假大數問題。

要構成一個 0 ，一定要有 2/5 同時存在相乘，如 8 * 5 = 222 *5 ，只能取出 1 對 2/5 出來，所以只有一個零；又如 16 * 125 = 2^4 * 5^3，可以取出 3 對出來，所以尾數有 3 個零。

然後很明顯的一點是，[1,n] 裡面，可以被 2 分解的個數，一定比可以比 5 分解的個數還多。

再來是，像 25 可以再被多提一個 5 出來、125 又可以再多提一個 5 出來，所以最後就是在算有幾個數字可以被 5 / 25 / 125 / 625 / … 分解，總合就是答案。

#include <stdio.h>int fact_tail_zero(int n)
{int ret=0;while(n/=5) ret+=n;return ret;
}int main()
{int x;while(scanf("%d", &x)==1)printf("tail zero : %d\n", fact_tail_zero(x));return 0;
}

[6] 其他假大數階乘問題

(1) 求 N! 最右邊第一個非零數字 O(logn)

(2) N! 之二進位表示法中，最低位元 1 之位置 -> 其實就是在求 [1,N] 質因數含 2 的個數。