长整数－－除法与取模

代码下载：

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。非常接近精确值。

代码实现：

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

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

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

    int x_len=l_len;
    int y_len=r_len;

    b_unit* x_val=x->b_val;
    b_unit* y_val=y->b_val;

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

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

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


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


    int i;
    b_unit y_t=y_val[y_len-1];
    b_unit y_t_1=y_val[y_len-2];
    b_unit q,r;

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

        twob_unit xx=x_i<<SHIFT|x_i_1;
        assert(x_i<=y_t);
        /*步骤3.1,对q[i-t-1]估值*/
        if(x_i==y_t)   /*情况1*/
        {
            q=BASE-1;
            r=(b_unit)(xx-(twob_unit)q*y_t);  /*计算r的值*/
        }
        else   /*情况2*/
        {
            q=(b_unit)(xx/y_t);
            r=(b_unit)(xx-(twob_unit)q*y_t); /*计算r的值*/
        }
        /*步骤3.2*/
        while((twob_unit)q*y_t_1>((twob_unit)r<<SHIFT)+x_i_2) /*精确q值*/
        {
            q--;
            r+=y_t;
            if(r>=BASE)
                break;
        }



        twob_unit y_mul=0;
        b_unit borrow=0;
        int j;
        /*步骤3.3   x=x-q[i-t-1]*y*b^(i-t-1)  */
        for(j=0;j<y_len;j++)
        {
            y_mul=(twob_unit)y_val[j]*q+y_mul+borrow;
            borrow=x_val[i-y_len+j]-(b_unit)(y_mul&MASK);
            x_val[i-y_len+j]=borrow&MASK;
            y_mul>>=SHIFT;
            borrow>>=SHIFT;
        }
        y_mul+=borrow;

        b_unit carry=0;



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


        x_val[i]=0;
        *c_qa--=q;

    }

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

    bg_format_len(y);
    bg_format_len(qa);

    *prem=y;
    *pdiv=qa;




    bg_free(x);

    return 1;
}

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

代码为：

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

    assert(y_len==1);

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


    b_unit* x_val=x->b_val;
    b_unit* y_val=y->b_val;

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

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

    int i;
    b_unit q;
    b_unit r;
    b_unit yt=y_val[y_len-1];

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

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

    bg_format_len(qa);
    bg_format_len(y);

    bg_free(x);

    *pdiv=qa;
    *prem=y;
    return 1;
}