从零开始学算法：高精度计算

前言:由于计算机运算是有模运算,数据范围的表示有一定限制,如整型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 2 #include 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;//数s110 memset(a,0,sizeof(a)); //数组清011 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 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大于011 {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，ab[0]) return 1;//a的位数大于b则a比b大4 if (a[0]0;i--) //从高位到低位比较6 {if (a[i]>b[i]) return 1;7 if (a[i]

三、高精度加法

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 }

四、高精度减法

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]

五、高精度乘法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=04 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]+17 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 }

六、 高精度除以低精度； 算法：按照从高位到低位的顺序，逐位相除。在除到第j位时，该位在接受了来自第j+1位的余数后与除数相除，如果最高位为零，则商的长度减一。源程序如下：

1 #include 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 }

七、高精度乘以高精度（要求用尽可能少的存储单元）； 算法：用数组保存两个高精度数，然后逐位相乘，注意考虑进位和总位数。源程序如下：

1 #include 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 }

八、高精度除以高精度（要求用尽可能少的存储单元）； 算法：用计算机模拟手算除法，把除法试商转化为连减。

1 #include 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 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（1000

【输出格式】

第一行：十进制高精度数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 2 #include 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 else49 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 2 #include 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 main44 {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 }