Chinaunix首页 | 论坛 | 博客
  • 博客访问: 236475
  • 博文数量: 32
  • 博客积分: 424
  • 博客等级: 准尉
  • 技术积分: 465
  • 用 户 组: 普通用户
  • 注册时间: 2012-03-02 10:23
文章分类

全部博文(32)

文章存档

2012年(32)

分类: C/C++

2012-03-12 10:38:28

代码下载:
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);

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

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

  9.     int x_len=l_len;
  10.     int y_len=r_len;

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

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

  15.     int carry;
  16.     carry=units_lshift(y_val,right->b_val,r_len,m); /*扩大y,通过左移来实现*/
  17.     assert(carry==0);
  18.     carry=units_lshift(x_val,left->b_val,l_len,m);  /*扩大x,通过左移实现*/

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


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


  27.     int i;
  28.     b_unit y_t=y_val[y_len-1];
  29.     b_unit y_t_1=y_val[y_len-2];
  30.     b_unit q,r;

  31.     for(i=x_len-1;i>=y_len;i--)  /*步骤3*/
  32.     {
  33.         b_unit x_i=x_val[i];
  34.         b_unit x_i_1=x_val[i-1];
  35.         b_unit x_i_2=x_val[i-2];

  36.         twob_unit xx=x_i<<SHIFT|x_i_1;
  37.         assert(x_i<=y_t);
  38.         /*步骤3.1,对q[i-t-1]估值*/
  39.         if(x_i==y_t)   /*情况1*/
  40.         {
  41.             q=BASE-1;   
  42.             r=(b_unit)(xx-(twob_unit)q*y_t);  /*计算r的值*/
  43.         }
  44.         else   /*情况2*/
  45.         {
  46.             q=(b_unit)(xx/y_t);   
  47.             r=(b_unit)(xx-(twob_unit)q*y_t); /*计算r的值*/
  48.         }
  49.         /*步骤3.2*/
  50.         while((twob_unit)q*y_t_1>((twob_unit)r<<SHIFT)+x_i_2) /*精确q值*/
  51.         {
  52.             q--;
  53.             r+=y_t;
  54.             if(r>=BASE)
  55.                 break;
  56.         }



  57.         twob_unit y_mul=0;
  58.         b_unit borrow=0;
  59.         int j;
  60.         /*步骤3.3   x=x-q[i-t-1]*y*b^(i-t-1)  */
  61.         for(j=0;j<y_len;j++)
  62.         {
  63.             y_mul=(twob_unit)y_val[j]*q+y_mul+borrow;
  64.             borrow=x_val[i-y_len+j]-(b_unit)(y_mul&MASK);
  65.             x_val[i-y_len+j]=borrow&MASK;
  66.             y_mul>>=SHIFT;
  67.             borrow>>=SHIFT;
  68.         }
  69.         y_mul+=borrow;

  70.         b_unit carry=0;


  71.  
  72.         assert(x_i-y_mul==-1||x_i-y_mul==0);
  73.         /*步骤3.4,前面对q估值,可以存在被q比真实的值大1,但这种情况很少发生*/
  74.         /*即把多减的,加回x*/
  75.         if((sb_unit)(x_i-y_mul)<0)
  76.         {
  77.             q--;
  78.             for(j=0;j<y_len;j++)
  79.             {
  80.                 carry=x_val[i-y_len+j]+y_val[j]+carry;
  81.                 x_val[i-y_len+j]=carry&MASK;
  82.                 carry>>=SHIFT;
  83.             }
  84.             assert(carry==1);
  85.         }


  86.         x_val[i]=0;
  87.         *c_qa--=q;

  88.     }

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

  92.     bg_format_len(y);
  93.     bg_format_len(qa);

  94.     *prem=y;
  95.     *pdiv=qa;




  96.     bg_free(x);

  97.     return 1;
  98. }
上面代码只是在除数长度为大于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);

  5.     assert(y_len==1);

  6.     BGInteger* x=bg_malloc(x_len+1);
  7.     BGInteger* y=bg_malloc(y_len);


  8.     b_unit* x_val=x->b_val;
  9.     b_unit* y_val=y->b_val;

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

  12.     if(x->b_val[x_len-1]>=y->b_val[y_len-1])  /*保证x[x_len-1]
  13.     {
  14.         x_len++;
  15.     }
  16.     int k_len=x_len-y_len;
  17.     BGInteger* qa=bg_malloc(k_len);   /*为商分配空间*/
  18.     b_unit* c_qa=qa->b_val+k_len-1;

  19.     int i;
  20.     b_unit q;
  21.     b_unit r;
  22.     b_unit yt=y_val[y_len-1];

  23.     for(i=x_len-1;i>=y_len;i--) /*一位数的除法比较简单,就像在草稿纸上做除法一样*/
  24.     {
  25.         b_unit x_i=x_val[i];
  26.         b_unit x_i_1=x_val[i-1];
  27.         twob_unit xx=x_i<<SHIFT|x_i_1;
  28.         q=xx/yt;
  29.         r=xx-(twob_unit)yt*q;

  30.         *c_qa--=q;
  31.         x_val[i-1]=r;
  32.     }
  33.     memcpy(y_val,x_val,y_len*sizeof(*y_val));

  34.     bg_format_len(qa);
  35.     bg_format_len(y);

  36.     bg_free(x);

  37.     *pdiv=qa;
  38.     *prem=y;
  39.     return 1;
  40. }




阅读(3218) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~