长整数－－除法与取模

5666阅读　0评论2012-03-12　NosicLin
分类：C/C++

代码下载：
git clone git://git.code.sf.net/p/cutility/code cutility-code
长整数总共有两个文件一个为big_integer.h，一个为big_integer.c

在看这一篇文章这前你需要了解长整数的数据结构

(一)简介

在长整数所有运算中，除法算得上是最难的一部份。除法的算法和前面介绍的乘法的步骤相似，类似于在草稿纸上计算多位数的除法。

(二)除法与取模

(1)两长整数绝对值的除法

算法描述为：
说明：
  1. 一个n位的数除以一个t位的数，商的长度不会大于n-t+1位
  2. 除数y的长度大于或等于2，除数x的长度大于或等于除数的长度
  3. 除数x与被除数y必须为正数
对于数值，BGInteger是用一个数组来存储。数组中每一个元素有16位，其中低15位用于存储数据。这种数据结构并不能很好适应上面的算法：

问题1：
步骤3.1是对q[i-t-1]估值，步骤3.2用一个循环来最终精确q[i-t-1]的值。在但步骤3.2中面临着这样一个问题，BGInteger数组中每一位有效数据为15位，但循环条件(q[i-t-1]*(y[t]*b+y[t-1]))>x[i]*b*b+x[i-1]*b+x[i-2])中，大于符号左边与右边的算式都需要至少45位才能完整的保存结果，保证不会产生数据溢出。但是在32位机上，不能直接保存与处理45位数据。所以步骤3.2不能直接使用，需要一定的转换。

转换过程如下：
其中 r 的值如下：
因为q[i-t-1]在步骤3.1的估值为：


所以r和q[i-t-1]的长度不会大于15位。

r是关于q[i-t-1]的函数，当q[i-t-1]变化时，r也会随之变化：
所以步骤3.2改为：
问题2：
在对商进行估值时，可能遇到这样的情况，假设十进制：
x[i]=0，x[i-1]=9；y[t]=1，y[t-1]=9。步骤3.1对于q进行估值为9，但是真正的值不会超过5，这样就会造成在步骤3.2至少循环4次，才能精确q的值。在BGInteger中以15位为一个单元来运算，很可能出现在步骤3.2循环上千至万次，大大降低了运算的效率。造成这种情况出现在原因主要是因为 y[t]*b-y[t-1]的值大小，例如上面y[t]=1，y[t-1]=9，b=10；y[t]*b-y[t-1]=1*10-9=1。
解决这种情况的一种方法为：除数与被除数扩大相同的倍数，让y[t]尽可能的大。
例如，对上面x与y同时扩大5倍得：x[i]=4 ，x[i-1]=5；y[t]=9，y[t-1]=5。对q估值为45/9＝5。非常接近精确值。

代码实现：

点击(此处)折叠或打开

  1. static int abs_divrem(BGInteger* left,BGInteger* right,BGInteger** pdiv,BGInteger** prem)
  2. {
  3.     /*取绝对值*/
  4.     int l_len=abs(left->b_len);
  5.     int r_len=abs(right->b_len);

  7.     BUG_ON(l_len<r_len||r_len<2,"Error Length");

  9.     BGInteger* x=bg_malloc(l_len+1);  /*为被除数x分配空间*/
  10.     BGInteger* y=bg_malloc(r_len);    /*为除数y分配空间*/

  12.     int x_len=l_len;
  13.     int y_len=r_len;

  15.     b_unit* x_val=x->b_val;
  16.     b_unit* y_val=y->b_val;

  18.     int d=bg_bit_in_unit(right->b_val[r_len-1]);/*计算除数y需要扩大的倍数*/
  19.     int m=SHIFT-d;                              /*通过判数最高位有几位为0来确定*/

  21.     int carry;
  22.     carry=units_lshift(y_val,right->b_val,r_len,m); /*扩大y，通过左移来实现*/
  23.     assert(carry==0);
  24.     carry=units_lshift(x_val,left->b_val,l_len,m);  /*扩大x，通过左移实现*/

  26.     if(carry!=0||x_val[x_len-1]>=y_val[y_len-1])    /*判断x的高位是否比y小*/
  27.     {                                               /*或者左移是否有溢出*/
  28.         x_val[x_len]=carry;    /*保证x的高位比y的高位小，这样就可以跳过步骤2*/
  29.         x_len++;
  30.     }


  33.     int k=x_len-y_len;
  34.     BGInteger* qa =bg_malloc(k); /*为商分配空间，因为x[x_len]
  35.     b_unit* c_qa=qa->b_val+k-1;  /*超过k位，k=x_len-y_len*/


  38.     int i;
  39.     b_unit y_t=y_val[y_len-1];
  40.     b_unit y_t_1=y_val[y_len-2];
  41.     b_unit q,r;

  43.     for(i=x_len-1;i>=y_len;i--)  /*步骤3*/
  44.     {
  45.         b_unit x_i=x_val[i];
  46.         b_unit x_i_1=x_val[i-1];
  47.         b_unit x_i_2=x_val[i-2];

  49.         twob_unit xx=x_i<<SHIFT|x_i_1;
  50.         assert(x_i<=y_t);
  51.         /*步骤3.1,对q[i-t-1]估值*/
  52.         if(x_i==y_t)   /*情况1*/
  53.         {
  54.             q=BASE-1;   
  55.             r=(b_unit)(xx-(twob_unit)q*y_t);  /*计算r的值*/
  56.         }
  57.         else   /*情况2*/
  58.         {
  59.             q=(b_unit)(xx/y_t);   
  60.             r=(b_unit)(xx-(twob_unit)q*y_t); /*计算r的值*/
  61.         }
  62.         /*步骤3.2*/
  63.         while((twob_unit)q*y_t_1>((twob_unit)r<<SHIFT)+x_i_2) /*精确q值*/
  64.         {
  65.             q--;
  66.             r+=y_t;
  67.             if(r>=BASE)
  68.                 break;
  69.         }



  73.         twob_unit y_mul=0;
  74.         b_unit borrow=0;
  75.         int j;
  76.         /*步骤3.3   x=x-q[i-t-1]*y*b^(i-t-1)  */
  77.         for(j=0;j<y_len;j++)
  78.         {
  79.             y_mul=(twob_unit)y_val[j]*q+y_mul+borrow;
  80.             borrow=x_val[i-y_len+j]-(b_unit)(y_mul&MASK);
  81.             x_val[i-y_len+j]=borrow&MASK;
  82.             y_mul>>=SHIFT;
  83.             borrow>>=SHIFT;
  84.         }
  85.         y_mul+=borrow;

  87.         b_unit carry=0;


  90.  
  91.         assert(x_i-y_mul==-1||x_i-y_mul==0);
  92.         /*步骤3.4，前面对q估值，可以存在被q比真实的值大1，但这种情况很少发生*/
  93.         /*即把多减的，加回x*/
  94.         if((sb_unit)(x_i-y_mul)<0)
  95.         {
  96.             q--;
  97.             for(j=0;j<y_len;j++)
  98.             {
  99.                 carry=x_val[i-y_len+j]+y_val[j]+carry;
  100.                 x_val[i-y_len+j]=carry&MASK;
  101.                 carry>>=SHIFT;
  102.             }
  103.             assert(carry==1);
  104.         }


  107.         x_val[i]=0;
  108.         *c_qa--=q;

  110.     }

  112.     /*因为除数与被除数都同时扩大了2^m倍，所以余数也同时会扩大2^m倍*/
  113.     /*余数左移m位，得到真正的余数*/
  114.     units_rshift(y_val,x_val,y_len,m);

  116.     bg_format_len(y);
  117.     bg_format_len(qa);

  119.     *prem=y;
  120.     *pdiv=qa;




  125.     bg_free(x);

  127.     return 1;
  128. }
上面代码只是在除数长度为大于2的情况下才能使用，如果长度为1，得还需要格外的实现
代码为：

点击(此处)折叠或打开

  1. static int abs_divrem1(BGInteger* left,BGInteger* right,BGInteger** pdiv,BGInteger** prem)
  2. {
  3.     int x_len=abs(left->b_len);
  4.     int y_len=abs(right->b_len);

  6.     assert(y_len==1);

  8.     BGInteger* x=bg_malloc(x_len+1);
  9.     BGInteger* y=bg_malloc(y_len);


  12.     b_unit* x_val=x->b_val;
  13.     b_unit* y_val=y->b_val;

  15.     memcpy(x_val,left->b_val,x_len*sizeof(*x_val)); /*为x分配空间*/
  16.     memcpy(y_val,right->b_val,y_len*sizeof(*y_val)); /*为y分配空间*/

  18.     if(x->b_val[x_len-1]>=y->b_val[y_len-1])  /*保证x[x_len-1]
  19.     {
  20.         x_len++;
  21.     }
  22.     int k_len=x_len-y_len;
  23.     BGInteger* qa=bg_malloc(k_len);   /*为商分配空间*/
  24.     b_unit* c_qa=qa->b_val+k_len-1;

  26.     int i;
  27.     b_unit q;
  28.     b_unit r;
  29.     b_unit yt=y_val[y_len-1];

  31.     for(i=x_len-1;i>=y_len;i--) /*一位数的除法比较简单，就像在草稿纸上做除法一样*/
  32.     {
  33.         b_unit x_i=x_val[i];
  34.         b_unit x_i_1=x_val[i-1];
  35.         twob_unit xx=x_i<<SHIFT|x_i_1;
  36.         q=xx/yt;
  37.         r=xx-(twob_unit)yt*q;

  39.         *c_qa--=q;
  40.         x_val[i-1]=r;
  41.     }
  42.     memcpy(y_val,x_val,y_len*sizeof(*y_val));

  44.     bg_format_len(qa);
  45.     bg_format_len(y);

  47.     bg_free(x);

  49.     *pdiv=qa;
  50.     *prem=y;
  51.     return 1;
  52. }




上一篇：长整数－－乘法
下一篇：长整数－－左移与右移