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

注：转载请注明：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
说明：位数   个位  十位  百位　千位
具体在计算加减乘除时方法就是小学时采用的列竖式方法．
注：高精度计算时一般用正数，对于负数，通过处理符号位的修正．
一．高精度数的存储
１．如对数采用的字符串输入

#include <iostream>
#include <cstring>
using namespace std;
const int N=100;//最多100位
int main()
{
int a[N+1],i;
string s1;
cin>>s1;//数s1
memset(a,0,sizeof(a)); //数组清0
a[0]=s1.length(); //位数
for(i=1;i<=a[0];i++) a[i]=s1[a[0]-i]-'0';//将字符转为数字并倒序存储．
return 0;
}

2.直接读入

#include <iostream>
using namespace std;
const int N=100;//最多100位
int main()
{
int a[N+1],i,s,key;
cin>>key;//数key
memset(a,0,sizeof(a)); //数组清0
i=0;//第0位
while(key)  //当key大于0
{
  a[++i]=key%10;//取第i位的数
  key=key/10;
}
a[0]=i; //共i位数
return 0;
}

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

以下程序都只写函数，不写完整程序，所有高精度数存储都满足上述约定。
二.高精度数比较

int compare(int a[],int b[])   //比较a和b的大小关系，若a>b则为1，a<b则为-1,a=b则为0
{int i;
if (a[0]>b[0]) return 1;//a的位数大于b则a比b大
if (a[0]<b[0]) return -1;//a的位数小于b则a比b小
for(i=a[0];i>0;i--)  //从高位到低位比较
     {if (a[i]>b[i]) return 1;
      if (a[i]<b[i]) return -1;}
return 0;//各位都相等则两数相等。
}

三、高精度加法

int plus(int a[],int b[]) //计算a=a+b
{int i,k;
k=a[0]>b[0]?a[0]:b[0]; //k是a和b中位数最大的一个的位数
for(i=1;i<=k;i++)
    {a[i+1]+=(a[i]+b[i])/10;  //若有进位，则先进位
    a[i]=(a[i]+b[i])%10;}  //计算当前位数字,注意：这条语句与上一条不能交换。
if(a[k+1]>0) a[0]=k+1;  //修正新的a的位数（a+b最多只能的一个进位）
               else a[0]=k;
return 0;
}

 

四、高精度减法

 

int gminus(int a[],int b[]);//计算a=a-b，返加符号位0:正数 1:负数
{ int flag,i
  flag=compare(a,b); //调用比较函数判断大小
if (falg==0)//相等
  {memset(a,0,sizeof(a));return 0;}  //若a=b，则a=0,也可在return前加一句a[0]=1,表示是 1位数0
if(flag==1) //大于  
  {  for(i=1;i<=a[0];i++)
      {  if(a[i]<b[i]){ a[i+1]--;a[i]+=10;} //若不够减则向上借一位
        a[i]=a[i]-b[i];}
     while(a[a[0]]==0) a[0]--; //修正a的位数
    return 0;}
if (flag==-1)//小于  则用a=b-a,返回-1
    { for(i=1;i<=b[0];i++)       {  if(b[i]<a[i]){ b[i+1]--;b[i]+=10;} //若不够减则向上借一位
        a[i]=b[i]-a[i];}
      a[0]=b[0];
     while(a[a[0]]==0) a[0]--; //修正a的位数
    return -1;}
}

五、高精度乘法1(高精度乘单精度数，单精度数是指通常的整型数)

int multi1(int a[],long  key) //a=a*key,key是单精度数  
{int i,k;
if (key==0){memset(a,0,sizeof(a));a[0]=1;return 0;} //单独处理key=0
for(i=1;i<=a[0];i++)a[i]=a[i]*key;//先每位乘起来
for(i=1;i<=a[0];i++){a[i+1]+=a[i]/10;a[i]%=10;} //进位
//注意上一语句退出时i=a[0]+1
while(a[i]>0) {a[i+1]=a[i]/10;a[i]=a[i]%10;i++;a[0]++];}  //继续处理超过原a[0]位数的进位,修正a的位数
return 0;
}

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

#include  <stdio.h>
#define   N  500
main()
{
  int  a[N] = {0}, c[N] = {0};
  int  i, k, d, b;
  char  a1[N];  
  printf("Input 除数:");
  scanf("%d", &b);
  printf("Input 被除数:");
  scanf("%s", a1);
  k = strlen(a1);
  for(i = 0; i < k; i++)  a[i] = a1[k - i - 1] - '0';
  d = 0;
  for(i = k - 1; i >= 0 ; i--)
  {
     d = d * 10 + a[i];
     c[i] = d / b;
     d = d % b;      
  }   
  while(c[k - 1] == 0 && k > 1)  k--;  
  printf("商=");
  for(i = k - 1; i >= 0; i--)  printf("%d", c[i]);
  printf("\n余数=%d", d);   
}     

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

#include  <stdio.h>
main()
{
  int  a[240] = {0}, b[240] = {0}, c[480] = {0};
  int  i, j, ka, kb, k;
  char  a1[240], b1[240];
  gets(a1);   
  ka = strlen(a1);
  gets(b1);   
  kb = strlen(b1);
  k = ka + kb;
  for(i = 0; i < ka; i++)  a[i] = a1[ka-i-1] - '0';
  for(i = 0; i < kb; i++)  b[i] = b1[kb-i-1] - '0';
  for(i = 0; i < ka; i++)
    for(j = 0; j < kb; j++)
    {
      c[i + j] = c[i + j] + a[i] * b[j];
      c[i + j +1] = c[i + j +1] + c[i + j]/10;
      c[i + j] = c[i + j] % 10;
    }
  if(!c[k])  k--;
  for(i = k-1; i >= 0; i--)  printf("%d", c[i]);        
}     

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

#include  <stdio.h>
#define   N  500
int  bj(int a[], int b[], int k1, int k2)   /*比较大小函数*/
{
   int i, t, flag;       /*flag作标志位*/
   if(k1 < k2)  
     flag = 0;           /*被除数小于除数返回0*/
   else if(k1 > k2)  
          flag = 1;      /*被除数大于除数返回1*/
        else
          {              /*被除数和除数位数相等则逐位进行比较*/
            i = k1;
            t = 0;
            while(t == 0 && i > 0)
            {
              if(a[i] > b[i]) {t = 1; flag = 1;}
              else if(a[i] == b[i])  i--;
              else  {t = 1; flag = 0;}        
            }
            if(i == 0 && t == 0)  flag = 2;     /*被除数等于除数返回2*/
          }
  return flag;           
}
int  jf(int a[], int b[], int k1, int k2)       /*减法运算*/
{
  int  i, k, d[N];
  for(i = 0; i < k2; i++)  d[i] = b[i];        /*把除数赋给数组d*/
  for(i = k2; i < N; i++)  d[i] = 0;          /*d数组无数据的高位置0*/
  k = k1 - k2 - 1;                            /*计算减法起始位置*/
  if(k < 0)  k = 0;
  if(k > 0)
  {
    for(i = k2 - 1; i >= 0; i--)  d[i + k] = d[i];  /*移动减数位数与被减数对齐*/
    for(i = 0; i < k; i++)  d[i] = 0;            /*移动后的其余位置0*/
  }  
  for(i = 0; i < k1; i++)
  {
    if(a[i] >= d[i])  a[i] -= d[i];
    else
    {
      a[i + 1] = a[i + 1] - 1;
      a[i] = 10 + a[i] - d[i];  
    }     
  }   
  return k;
}
main()
{
  int  a[N] = {0}, b[N] = {0}, c[N] = {0}, d[N] = {0};
  int  i, ka, kb, m, t, t1, t2, k, x, kd, kk;
  char  a1[N], b1[N];  
  printf("Input 被除数:");
  scanf("%s", a1);
  ka = strlen(a1);
  for(i = 0; i < ka; i++)  a[i] = a1[ka - i -1] - '0';
  printf("Input 除数:");
  scanf("%s", b1);
  kb = strlen(b1);
  for(i = 0; i < kb; i++)  b[i] = b1[kb - i -1] - '0';
  kd = ka;    /*保存被除数位数  */
  t2 = bj(a, b, ka, kb);
  m = 0;
  do
  {
    while(a[ka - 1] == 0)  ka--;
    t = bj(a, b, ka, kb);   
    if(t >= 1)
    {
      k = jf(a, b, ka, kb);
      c[k]++;      
      if(k > m)  m = k;
      t1 = 0;
      for(i = k; i <= m; i++)
      {
        x = c[i] + t1;
        c[i] = x % 10;
        t1 = x / 10;     
      }
      if(t1 > 0)  {m++; c[m] = t1;  }     
    }   
  }while(t == 1);
  if(t2 == 0)  
  {
    printf("商=0");  
    printf("\n余数=");
    for(i = kd - 1; i >= 0; i--)  printf("%d", a[i]);
    exit(1);  
  }
  if(t2 == 2)
  {
    printf("商 = 1");  
    printf("\n余数 = 0");
    exit(1);  
  }
  kk = kd;
  while(!c[kd - 1])  kd--;
  printf("商 = ");
  for(i = kd - 1; i >= 0; i--)  printf("%d", c[i]);
  while(!a[kk])  kk--;
  printf("\n余数 = ");
  if(kk < 0)  
  {
    printf("0");  
    exit(1);
  }
  for(i = kk; i >= 0; i--)  printf("%d", a[i]);
} 

下面给出一些案例：

问题1. N！，要求精确到P位（0〈P〈1000〉。
算法：结果用数组a保存，开始时a[0]=1，依次乘以数组中各位，注意进位和数组长度的变化。源程序如下：

#include   <stdio.h>
#define    M   1000
main()
{
  int a[M], i, n, j, flag = 1;
  printf("n=");
  scanf("%d",&n);
  printf("n!=");
  a[0] = 1;
  for(i = 1; i < M; i++) a[i] = 0;
   for(j = 2; j <= n; j++)
   {
     for(i = 0; i < flag; i++) a[i] *= j;
     for(i = 0; i < flag; i++)
       if(a[i] >= 10)
       {
         a[i+1] += a[i]/10;
         a[i] = a[i] % 10;
         if(i == flag-1)  flag++;
       }
    }
  for(j = flag - 1; j >= 0; j--)
    printf("%d", a[j]);
} 

问题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。对于个位单独考虑，每次左移一位。源程序如下：

#include <stdio.h>
#include <math.h>
#define  MAX  100000
main()
{
   int p;
   int i, j;
   scanf("%d", &p);
   printf("%d\n", (int)(p * log10(2.0)) + 1);
   long  store[110] = {0};
   store[0] = 1;
   int left = p % 10;
   p /= 10;
    for(i = 1; i <= p; i++)
    {
      for(j = 0; j <= 100; j++)
        store[j] <<= 10;
      for(j = 0; j <= 100; j++)
      {
        if(store[j] >= MAX)
        {
          store[j + 1] += store[j] / MAX;
          store[j] %= MAX;
        }
      }
    }
    for(i = 1; i <= left; i++)
    {
      for(j = 0; j <= 100; j++)
        store[j] <<= 1;
      for(j = 0; j <= 100; j++)
      {
        if(store[j] >= MAX)
        {
          store[j + 1] += store[j] / MAX;
          store[j] %= MAX;
        }
      }
    }
    store[0] -= 1;
    for(i = 1; i < 100; i++)
    {
      if(store[i - 1] < 0)
      {
         store[i] -= 1;
         store[i - 1] += MAX;
      }
      else
        break;
    }
    for(i = 99; i >= 0; i--)
    {
      printf("%05d", store[i]);
      if((100 - i) % 10 == 0)
          printf("\n");
    }
} 

问题3. 有一个正整数N（N可能达到120位），它是由若干个不大于65535的正整数相乘而得到的。请把这个数分解成素数因子（质因子）的乘积。
输入：输入文件只有一行为N的值。
输出：（1）素数因子由小到大分行输出；
（2）每一行输出一个素数因子和该素数因子的个数，用一个空格分开；
（3）如果正整数N的分解中有一个以上的大于65535的素数，请按照（1）、（2）的要求输出分解中的小于65535的素数后，在下一行输出
“DATA  ERROR！”。
算法：先将2到65535之间的所有素数保存在数组中，用这个数去除数组中的每一个数，得到一个质因数就打印出来。源程序如下：

#include  <stdio.h>
#include  <math.h>
int length, temp[120];
int sushu(int a[])
{
  int i, j, k = 0, m;
  for(i = 2; i <= 65537; i++)
  {
    m = sqrt(i);
    for(j = 2; j <= m; j++)
      if(i % j == 0)  break;
    if(j > m)
    {
      a[k] = i;
      k++;   
    }         
  }
 return k;     
}
int divide(int a[], int k)
{
  int i, d = 0;
  for(i = length - 1; i >= 0; i--)
  {
     d = d * 10 + a[i];
     temp[i] = d / k;
     d = d % k;      
  }   
  if(!d)
    {
      while(temp[length - 1] == 0 && length > 1)  length--;
      for(i = 0; i < length; i++)
      {
        a[i] = temp[i];
        temp[i] = 0;      
      }
      for(i = length; i < 120; i++) a[i] = 0;   
    }
    else  
      for(i = 0; i < length; i++)  temp[i] = 0;
  return d;     
}
main()
{
  int i, k, s, d;       /*s计数器; d余数*/
  int a[6600], b[120] = {0}, c[120] = {0};
  char b1[120];
  gets(b1);
 length = strlen(b1);
  for(i = 0; i < length; i++)  b[i] = b1[length - i - 1] - '0';
  k = sushu(a);
  for(i = 0; i < k; i++)
  {
    s = 0;
    d = divide(b, a[i]);
    while(!d)
    {
      s++;
      d = divide(b, a[i]);         
    }
    if(i == k - 1)  
 
    {  
      printf("Data Error!");
      break;
    }