從零開始學算法：高精度計算

注：轉載請註明：http://www.cnblogs.com/ECJTUACM-873284962/ 

前言:因爲計算機運算是有模運算,數據範圍的表示有必定限制,如整型int(C++中int 與long相同)表達範圍是(-2^31~2^31-1),unsigned long(無符號整數)是(0~2^32-1),都約爲幾十億.若是採用實數型,則能保存最大的double只能提供15~16位的有效數字,即只能精確表達數百萬億的數.所以,在計算位數超過十幾位的數時,不能採用現有類型,只能本身編程計算.
高精度計算通用方法:高精度計算時通常用一個數組來存儲一個數,數組的一個元素對應於數的一位(固然,在之後的學習中爲了加快計算速度,也可用數組的一個元素表示數的多位數字,暫時不講),表示時,因爲數計算時可能要進位,所以爲了方便,將數由低位到高位依次存在數組下標對應由低到高位置上,另外,咱們申請數組大小時,通常考慮了最大的狀況,在不少狀況下,表示有富餘,即高位有不少0,可能形成無效的運算和判斷,所以,咱們通常將數組的第0個下標對應位置來存儲該數的位數.如數:3485(三千四百八十五)，表達在數組a[10]上狀況是:
下標　　0  　 1　   2   　3   　　4 　  5    6    7    8    9  
內容　　4   　5  　 8   　4   　　3  　 0    0    0    0    0
說明：位數   個位  十位  百位　千位
具體在計算加減乘除時方法就是小學時採用的列豎式方法．
注：高精度計算時通常用正數，對於負數，經過處理符號位的修正．
一．高精度數的存儲
１．如對數採用的字符串輸入

 1 #include <iostream>
 2 #include <cstring>
 3 using namespace std;
 4 const int N=100;//最多100位
 5 int main()
 6 {
 7 int a[N+1],i;
 8 string s1;
 9 cin>>s1;//數s1
10 memset(a,0,sizeof(a)); //數組清0
11 a[0]=s1.length(); //位數
12 for(i=1;i<=a[0];i++) a[i]=s1[a[0]-i]-'0';//將字符轉爲數字並倒序存儲．
13 return 0;
14 }

2.直接讀入

 1 #include <iostream>
 2 using namespace std;
 3 const int N=100;//最多100位
 4 int main()
 5 {
 6 int a[N+1],i,s,key;
 7 cin>>key;//數key
 8 memset(a,0,sizeof(a)); //數組清0
 9 i=0;//第0位
10 while(key)  //當key大於0
11 {
12   a[++i]=key%10;//取第i位的數
13   key=key/10;
14 }
15 a[0]=i; //共i位數
16 return 0;
17 }

3.直接初始化(用a[]存儲)
初始化爲0: memset(a,0,sizeof(a));
初始化爲1: memset(a,0,sizeof(a));a[0]=1;a[1]=1;

如下程序都只寫函數，不寫完整程序，全部高精度數存儲都知足上述約定。
二.高精度數比較

1 int compare(int a[],int b[])   //比較a和b的大小關係，若a>b則爲1，a<b則爲-1,a=b則爲0
2 {int i;
3 if (a[0]>b[0]) return 1;//a的位數大於b則a比b大
4 if (a[0]<b[0]) return -1;//a的位數小於b則a比b小
5 for(i=a[0];i>0;i--)  //從高位到低位比較
6      {if (a[i]>b[i]) return 1;
7       if (a[i]<b[i]) return -1;}
8 return 0;//各位都相等則兩數相等。
9 }

3、高精度加法

 1 int plus(int a[],int b[]) //計算a=a+b
 2 {int i,k;
 3 k=a[0]>b[0]?a[0]:b[0]; //k是a和b中位數最大的一個的位數
 4 for(i=1;i<=k;i++)
 5     {a[i+1]+=(a[i]+b[i])/10;  //如有進位，則先進位
 6     a[i]=(a[i]+b[i])%10;}  //計算當前位數字,注意：這條語句與上一條不能交換。
 7 if(a[k+1]>0) a[0]=k+1;  //修正新的a的位數（a+b最多隻能的一個進位）
 8                else a[0]=k;
 9 return 0;
10 }

 

4、高精度減法

 

 1 int gminus(int a[],int b[]);//計算a=a-b，返加符號位0:正數 1:負數
 2 { int flag,i
 3   flag=compare(a,b); //調用比較函數判斷大小
 4 if (falg==0)//相等
 5   {memset(a,0,sizeof(a));return 0;}  //若a=b，則a=0,也可在return前加一句a[0]=1,表示是 1位數0
 6 if(flag==1) //大於  
 7   {  for(i=1;i<=a[0];i++)
 8       {  if(a[i]<b[i]){ a[i+1]--;a[i]+=10;} //若不夠減則向上借一位
 9         a[i]=a[i]-b[i];}
10      while(a[a[0]]==0) a[0]--; //修正a的位數
11     return 0;}
12 if (flag==-1)//小於  則用a=b-a,返回-1
13     { for(i=1;i<=b[0];i++)       {  if(b[i]<a[i]){ b[i+1]--;b[i]+=10;} //若不夠減則向上借一位
14         a[i]=b[i]-a[i];}
15       a[0]=b[0];
16      while(a[a[0]]==0) a[0]--; //修正a的位數
17     return -1;}
18 }

5、高精度乘法1(高精度乘單精度數，單精度數是指一般的整型數)

1 int multi1(int a[],long  key) //a=a*key,key是單精度數  
2 {int i,k;
3 if (key==0){memset(a,0,sizeof(a));a[0]=1;return 0;} //單獨處理key=0
4 for(i=1;i<=a[0];i++)a[i]=a[i]*key;//先每位乘起來
5 for(i=1;i<=a[0];i++){a[i+1]+=a[i]/10;a[i]%=10;} //進位
6 //注意上一語句退出時i=a[0]+1
7 while(a[i]>0) {a[i+1]=a[i]/10;a[i]=a[i]%10;i++;a[0]++];}  //繼續處理超過原a[0]位數的進位,修正a的位數
8 return 0;
9 }

6、 高精度除以低精度；
算法：按照從高位到低位的順序，逐位相除。在除到第j位時，該位在接受了來自第j+1位的餘數後與除數相除，若是最高位爲零，則商的長度減一。源程序以下：

 1 #include  <stdio.h>
 2 #define   N  500
 3 main()
 4 {
 5   int  a[N] = {0}, c[N] = {0};
 6   int  i, k, d, b;
 7   char  a1[N];  
 8   printf("Input 除數:");
 9   scanf("%d", &b);
10   printf("Input 被除數:");
11   scanf("%s", a1);
12   k = strlen(a1);
13   for(i = 0; i < k; i++)  a[i] = a1[k - i - 1] - '0';
14   d = 0;
15   for(i = k - 1; i >= 0 ; i--)
16   {
17      d = d * 10 + a[i];
18      c[i] = d / b;
19      d = d % b;      
20   }   
21   while(c[k - 1] == 0 && k > 1)  k--;  
22   printf("商=");
23   for(i = k - 1; i >= 0; i--)  printf("%d", c[i]);
24   printf("\n餘數=%d", d);   
25 }     

7、高精度乘以高精度（要求用盡量少的存儲單元）；
算法：用數組保存兩個高精度數，而後逐位相乘，注意考慮進位和總位數。源程序以下：

 1 #include  <stdio.h>
 2 main()
 3 {
 4   int  a[240] = {0}, b[240] = {0}, c[480] = {0};
 5   int  i, j, ka, kb, k;
 6   char  a1[240], b1[240];
 7   gets(a1);   
 8   ka = strlen(a1);
 9   gets(b1);   
10   kb = strlen(b1);
11   k = ka + kb;
12   for(i = 0; i < ka; i++)  a[i] = a1[ka-i-1] - '0';
13   for(i = 0; i < kb; i++)  b[i] = b1[kb-i-1] - '0';
14   for(i = 0; i < ka; i++)
15     for(j = 0; j < kb; j++)
16     {
17       c[i + j] = c[i + j] + a[i] * b[j];
18       c[i + j +1] = c[i + j +1] + c[i + j]/10;
19       c[i + j] = c[i + j] % 10;
20     }
21   if(!c[k])  k--;
22   for(i = k-1; i >= 0; i--)  printf("%d", c[i]);        
23 }     

8、高精度除以高精度（要求用盡量少的存儲單元）；
算法：用計算機模擬手算除法，把除法試商轉化爲連減。

  1 #include  <stdio.h>
  2 #define   N  500
  3 int  bj(int a[], int b[], int k1, int k2)   /*比較大小函數*/
  4 {
  5    int i, t, flag;       /*flag做標誌位*/
  6    if(k1 < k2)  
  7      flag = 0;           /*被除數小於除數返回0*/
  8    else if(k1 > k2)  
  9           flag = 1;      /*被除數大於除數返回1*/
 10         else
 11           {              /*被除數和除數位數相等則逐位進行比較*/
 12             i = k1;
 13             t = 0;
 14             while(t == 0 && i > 0)
 15             {
 16               if(a[i] > b[i]) {t = 1; flag = 1;}
 17               else if(a[i] == b[i])  i--;
 18               else  {t = 1; flag = 0;}        
 19             }
 20             if(i == 0 && t == 0)  flag = 2;     /*被除數等於除數返回2*/
 21           }
 22   return flag;           
 23 }
 24 int  jf(int a[], int b[], int k1, int k2)       /*減法運算*/
 25 {
 26   int  i, k, d[N];
 27   for(i = 0; i < k2; i++)  d[i] = b[i];        /*把除數賦給數組d*/
 28   for(i = k2; i < N; i++)  d[i] = 0;          /*d數組無數據的高位置0*/
 29   k = k1 - k2 - 1;                            /*計算減法起始位置*/
 30   if(k < 0)  k = 0;
 31   if(k > 0)
 32   {
 33     for(i = k2 - 1; i >= 0; i--)  d[i + k] = d[i];  /*移動減數位數與被減數對齊*/
 34     for(i = 0; i < k; i++)  d[i] = 0;            /*移動後的其他位置0*/
 35   }  
 36   for(i = 0; i < k1; i++)
 37   {
 38     if(a[i] >= d[i])  a[i] -= d[i];
 39     else
 40     {
 41       a[i + 1] = a[i + 1] - 1;
 42       a[i] = 10 + a[i] - d[i];  
 43     }     
 44   }   
 45   return k;
 46 }
 47 main()
 48 {
 49   int  a[N] = {0}, b[N] = {0}, c[N] = {0}, d[N] = {0};
 50   int  i, ka, kb, m, t, t1, t2, k, x, kd, kk;
 51   char  a1[N], b1[N];  
 52   printf("Input 被除數:");
 53   scanf("%s", a1);
 54   ka = strlen(a1);
 55   for(i = 0; i < ka; i++)  a[i] = a1[ka - i -1] - '0';
 56   printf("Input 除數:");
 57   scanf("%s", b1);
 58   kb = strlen(b1);
 59   for(i = 0; i < kb; i++)  b[i] = b1[kb - i -1] - '0';
 60   kd = ka;    /*保存被除數位數  */
 61   t2 = bj(a, b, ka, kb);
 62   m = 0;
 63   do
 64   {
 65     while(a[ka - 1] == 0)  ka--;
 66     t = bj(a, b, ka, kb);   
 67     if(t >= 1)
 68     {
 69       k = jf(a, b, ka, kb);
 70       c[k]++;      
 71       if(k > m)  m = k;
 72       t1 = 0;
 73       for(i = k; i <= m; i++)
 74       {
 75         x = c[i] + t1;
 76         c[i] = x % 10;
 77         t1 = x / 10;     
 78       }
 79       if(t1 > 0)  {m++; c[m] = t1;  }     
 80     }   
 81   }while(t == 1);
 82   if(t2 == 0)  
 83   {
 84     printf("商=0");  
 85     printf("\n餘數=");
 86     for(i = kd - 1; i >= 0; i--)  printf("%d", a[i]);
 87     exit(1);  
 88   }
 89   if(t2 == 2)
 90   {
 91     printf("商 = 1");  
 92     printf("\n餘數 = 0");
 93     exit(1);  
 94   }
 95   kk = kd;
 96   while(!c[kd - 1])  kd--;
 97   printf("商 = ");
 98   for(i = kd - 1; i >= 0; i--)  printf("%d", c[i]);
 99   while(!a[kk])  kk--;
100   printf("\n餘數 = ");
101   if(kk < 0)  
102   {
103     printf("0");  
104     exit(1);
105   }
106   for(i = kk; i >= 0; i--)  printf("%d", a[i]);
107 } 

下面給出一些案例：

問題1. N！，要求精確到P位（0〈P〈1000〉。
算法：結果用數組a保存，開始時a[0]=1，依次乘以數組中各位，注意進位和數組長度的變化。源程序以下：

 1 #include   <stdio.h>
 2 #define    M   1000
 3 main()
 4 {
 5   int a[M], i, n, j, flag = 1;
 6   printf("n=");
 7   scanf("%d",&n);
 8   printf("n!=");
 9   a[0] = 1;
10   for(i = 1; i < M; i++) a[i] = 0;
11    for(j = 2; j <= n; j++)
12    {
13      for(i = 0; i < flag; i++) a[i] *= j;
14      for(i = 0; i < flag; i++)
15        if(a[i] >= 10)
16        {
17          a[i+1] += a[i]/10;
18          a[i] = a[i] % 10;
19          if(i == flag-1)  flag++;
20        }
21     }
22   for(j = flag - 1; j >= 0; j--)
23     printf("%d", a[j]);
24 } 

問題2. 麥森數
【問題描述】形如2P-1的素數稱爲麥森數，這時P必定也是個素數。但反過來不必定，即若是P是個素數，2P-1不必定也是素數。到1998年末，人們已找到了37個麥森數。最大的一個是P=3021377，它有909526位。麥森數有許多重要應用，它與徹底數密切相關。
任務：從文件中輸入P（1000<P<3100000），計算2P-1的位數和最後500位數字（用十進制高精度數表示）
【輸入格式】
文件中只包含一個整數P（1000<P<3100000）
【輸出格式】
第一行：十進制高精度數2P-1的位數。
第2-11行：十進制高精度數2P-1的最後500位數字。（每行輸出50位，共輸出10行，不足500位時高位補0）
沒必要驗證2P-1與P是否爲素數。
【輸入樣例】
1279
【輸出樣例】
386
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000104079321946643990819252403273640855
38615262247266704805319112350403608059673360298012
23944173232418484242161395428100779138356624832346
49081399066056773207629241295093892203457731833496
61583550472959420547689811211693677147548478866962
50138443826029173234888531116082853841658502825560
46662248318909188018470682222031405210266984354887
32958028878050869736186900714720710555703168729087
算法：2的冪能夠轉化成左移運算，爲了提升運算速度，可每次左移10位，即每次乘210。對於個位單獨考慮，每次左移一位。源程序以下：

 1 #include <stdio.h>
 2 #include <math.h>
 3 #define  MAX  100000
 4 main()
 5 {
 6    int p;
 7    int i, j;
 8    scanf("%d", &p);
 9    printf("%d\n", (int)(p * log10(2.0)) + 1);
10    long  store[110] = {0};
11    store[0] = 1;
12    int left = p % 10;
13    p /= 10;
14     for(i = 1; i <= p; i++)
15     {
16       for(j = 0; j <= 100; j++)
17         store[j] <<= 10;
18       for(j = 0; j <= 100; j++)
19       {
20         if(store[j] >= MAX)
21         {
22           store[j + 1] += store[j] / MAX;
23           store[j] %= MAX;
24         }
25       }
26     }
27     for(i = 1; i <= left; i++)
28     {
29       for(j = 0; j <= 100; j++)
30         store[j] <<= 1;
31       for(j = 0; j <= 100; j++)
32       {
33         if(store[j] >= MAX)
34         {
35           store[j + 1] += store[j] / MAX;
36           store[j] %= MAX;
37         }
38       }
39     }
40     store[0] -= 1;
41     for(i = 1; i < 100; i++)
42     {
43       if(store[i - 1] < 0)
44       {
45          store[i] -= 1;
46          store[i - 1] += MAX;
47       }
48       else
49         break;
50     }
51     for(i = 99; i >= 0; i--)
52     {
53       printf("%05d", store[i]);
54       if((100 - i) % 10 == 0)
55           printf("\n");
56     }
57 } 

問題3. 有一個正整數N（N可能達到120位），它是由若干個不大於65535的正整數相乘而獲得的。請把這個數分解成素數因子（質因子）的乘積。
輸入：輸入文件只有一行爲N的值。
輸出：（1）素數因子由小到大分行輸出；
（2）每一行輸出一個素數因子和該素數因子的個數，用一個空格分開；
（3）若是正整數N的分解中有一個以上的大於65535的素數，請按照（1）、（2）的要求輸出分解中的小於65535的素數後，在下一行輸出
「DATA  ERROR！」。
算法：先將2到65535之間的全部素數保存在數組中，用這個數去除數組中的每個數，獲得一個質因數就打印出來。源程序以下：

 1 #include  <stdio.h>
 2 #include  <math.h>
 3 int length, temp[120];
 4 int sushu(int a[])
 5 {
 6   int i, j, k = 0, m;
 7   for(i = 2; i <= 65537; i++)
 8   {
 9     m = sqrt(i);
10     for(j = 2; j <= m; j++)
11       if(i % j == 0)  break;
12     if(j > m)
13     {
14       a[k] = i;
15       k++;   
16     }         
17   }
18  return k;     
19 }
20 int divide(int a[], int k)
21 {
22   int i, d = 0;
23   for(i = length - 1; i >= 0; i--)
24   {
25      d = d * 10 + a[i];
26      temp[i] = d / k;
27      d = d % k;      
28   }   
29   if(!d)
30     {
31       while(temp[length - 1] == 0 && length > 1)  length--;
32       for(i = 0; i < length; i++)
33       {
34         a[i] = temp[i];
35         temp[i] = 0;      
36       }
37       for(i = length; i < 120; i++) a[i] = 0;   
38     }
39     else  
40       for(i = 0; i < length; i++)  temp[i] = 0;
41   return d;     
42 }
43 main()
44 {
45   int i, k, s, d;       /*s計數器; d餘數*/
46   int a[6600], b[120] = {0}, c[120] = {0};
47   char b1[120];
48   gets(b1);
49  length = strlen(b1);
50   for(i = 0; i < length; i++)  b[i] = b1[length - i - 1] - '0';
51   k = sushu(a);
52   for(i = 0; i < k; i++)
53   {
54     s = 0;
55     d = divide(b, a[i]);
56     while(!d)
57     {
58       s++;
59       d = divide(b, a[i]);         
60     }
61     if(i == k - 1)  
62  
63     {  
64       printf("Data Error!");
65       break;
66     }