超长整数的基础运算 算法实现之乘、除篇

笔算乘法：

对于m位和n位的输入，传统的乘法需要m*n次基本的乘法，也即算法复杂度为O()。我们用纸和笔做乘法运算时，用乘数的每一位乘以被乘数的每一位并加上上一列的进位而产生一行适当移位的中间结果，然后再将各行中间结果相加即得到乘法的最终结果，例如10进制下计算189*34的过程如下表：

 笔算乘法的运算过程

本算法按照上述过程进行计算，但在计算机上最好是把内部的乘法和加法并行的执行，即计算每一行中间结果的同时将该行加到最终结果上，这样既可以省去不少步骤，也避免了中间结果的空间开销和内存管理操作，这个简单的优化能够一定程度上提高效率。本算法的伪代码如下表：

由于两个16位的乘积将超过16位所能表示的最大数，我们知道一个32位的数能够完整地保存两个16位数的乘积，上面的伪码中用一个32位来保存两个16位的数乘积再加上两个单字后的结果，下面讨论一下它的安全性（是否会产生溢出）：上面表达式计算的最大的结果为(2¹⁶ - 1)*(2¹⁶- 1) + (2¹⁶ - 1)，化简后为2³² - 2¹⁶，小于一个32位数的最大值，所以一个32位数能够保存该表达式的结果而不产生溢出，程序不会损失精度。

在乘法中有一类乘法较为特殊--自平方。针对这类乘法使用一般的计算法自然是可行的，但是可以采用更加快捷的方式，使用一次遍历整数每一位做乘法后即可解决一般乘法需要嵌套循环两次才能处理的结果，因为自平方的两个乘数是“对称”的，以A（a1a2a3a4）为例：每一位的乘法执行后都会有对应的位重复操作一次，比如a1*a4会执行两次，而且对应的乘积结果在相同的结果位置上累加。所以在此基础上自平方的算法效率要比一般乘法要高。

/*
乘法运算
*/
int mulHBInt(HBigInt *product, HBigInt *biA, HBigInt *biB){
	HBigInt biT;		//计算结果先保存在临时变量中
	register un_short *pWordA = biA->pBigInt;
	register un_short *pWordB = biB->pBigInt;
	register un_short *pPrd = NULL;
	long result = 0, i = 0, j = 0;

    //初始化临时大整数变量
	if((result = initHBInt(&biT,biA->length + biB->length)) != RETURN_OK_BINT)
		return result;
	
	biT.length = biA->length + biB->length;
	biT.sign = (biA->sign == biB->sign) ? 1: -1; //同号为正，异号为负
	pPrd = biT.pBigInt;
	
	for(i=0; i<biA->length; ++i) { // 按传统的纸笔乘法的方式计算乘积
		for(j=0; j<biB->length; ++j) {
			pPrd[i+j] = pPrd[i+j] + pWordA[i] * pWordB[j];
			// 如果超出单位长度能表示最大的值（本例中是2<<16）,则进行一次格式化处理
			if(pPrd[i+j] >> BIT_PRE_WORD) formatHBInt(&biT);
		}
	}
	//formatHBInt(&biT);

	trimHBInt(&biT);//去除高位无效的0

	extendHBInt(product,biT.length);
	assignHBInt(product,&biT);
	deleteHBInt(&biT); //清除临时变量

	return RETURN_OK_BINT;
}

此算法是模拟笔算除法的形式，变形而来。根据笔算除法的特性，我们知道最重要的也是较难的一步是，对商的估算。我们做除法运算，是根据被除数和除数最高几位来试商的，此算法也是使用这种方式。在试商时，为了提高速度，尝试使用的商不可能是从1开始直到i，i满足。本运算使用二分查找法来试商，从而提高了速度。

       除试商外，另外一个难点就是对除数进行位对齐。多数操作中，除数的位数小于被除数，所以在进行试商前需要对除数进行左移，即人为扩大除数的倍数（Bⁿ）。至于左移具体几位，需要对高位进行比较后确定，规则如下：

/*
除法运算
a 被除数
b 除数
c 商
d 余数
*/
int divHBInt(HBigInt *a, HBigInt *b, HBigInt *c, HBigInt *d){
	HBigInt ta,tb,tc,temp,tc_1;  //临时变量
	long result=0,first=1,middle=0,last=0;
	un_short mark=0;
	long len_ta,len_tb,len,i;
	int re;

	if(0 == b->length || (1 == b->length && 0 == b->pBigInt[0])) return -1; //除数不能为0
	//除数为1，则商等于被除数，余数为0
	if(1 == b->length && 1 == b->pBigInt[0]) {
		assignHBInt(c,a);
		setZeroHBInt(d);
		hbi_add_int(d,0);
		return RETURN_OK_BINT;
	}
	
	re = compareHBInt(a,b);
	if(-1 == re) { //除数大于被除数，商为0，余数为被除数
		setZeroHBInt(c);
		assignHBInt(d,a);
		return RETURN_OK_BINT;
	}
	
	if(0 == re) { //除数等于被除数，商为1，余数为0
		setZeroHBInt(c);
		setZeroHBInt(d);
		hbi_add_int(c,1);
		hbi_add_int(d,0);
		return RETURN_OK_BINT;
	}

	//初始化临时变量
	initHBInt(&ta,INITIAL_BINT);
	initHBInt(&tb,INITIAL_BINT);
	initHBInt(&tc,INITIAL_BINT);
	initHBInt(&temp,INITIAL_BINT);
	initHBInt(&tc_1,INITIAL_BINT);

	len_ta=ta.length;
	len_tb=ta.length;
	len=len_ta-len_tb;
	i=0;

	assignHBInt(&ta,a);				//对临时变量赋值
	assignHBInt(&tb,b);
	hbi_add_int(&tc,1);				//初始商为1
	extendHBInt(&temp,ta.length);	//扩展临时变量长度
	extendHBInt(&tc_1,ta.length);

	//二分法试商计算
	while(compareHBInt(&ta,b) > 0 ) {  // 保证ta大于b时循环
		len_ta=ta.length;
		len_tb=tb.length;
		len=len_ta-len_tb;
		i=0;
		mark=0;
		
		if(ta.pBigInt[len_ta-1] <= tb.pBigInt[len_tb-1]) {
			while(i<len_tb){
				if(ta.pBigInt[len_ta-1-i] < tb.pBigInt[len_tb-1-i]) {
					mark=ta.pBigInt[len_ta-1];
					len--;
					i++;
					break;
				} else i++;
			}
		}

		Left_shift_word(&tb,len);	//对除数进行左移位，使其跟被除数有相同数位
		Left_shift_word(&tc,len);	//商同时左移位
		
		result=mark*CARRY_RADIX+ta.pBigInt[len_ta-1-i];
		first=1;
		last= result / tb.pBigInt[len_tb-1+len] + 1;	//二分查找法的上限
		middle=(last+first)/2 + 1;						//二分查找的中间值

		for(; (first < last) && (first != middle); ) {	//使用二分查找法来试商
			assignHBInt(&temp,&tb);
			hbi_mul_radix(&temp,middle);
			
			if(-1 == compareHBInt(&ta,&temp)) last=middle;
			else first=middle;
			
			middle=(last+first)/2;
		}

		assignHBInt(&temp,&tb);
		hbi_mul_radix(&temp,middle);
		subHBInt(&ta,&ta,&temp);	//被除数-除数*商
		
		hbi_mul_radix(&tc,middle);  //得到所试的商
		
		addHBInt(&tc_1,&tc_1,&tc);	//对商进行累加

		setZeroHBInt(&tc);          //临时保存商的值清0
		hbi_add_int(&tc,1);			//对商的临时变量赋1      

		setZeroHBInt(&temp);	     
		assignHBInt(&tb,b);
	}
	
	if(compareHBInt(&ta,b) == 0) {
		hbi_add_int(c,1);
		setZeroHBInt(d);
	} else {
		if(c!=NULL)	assignHBInt(c,&tc_1);	//得到除法运算的商
		if(d!=NULL) assignHBInt(d,&ta);		//除法运算的余数
	}
	
	// 回收临时变量空间
	deleteHBInt(&ta); 
	deleteHBInt(&tb);
	deleteHBInt(&tc);
	deleteHBInt(&tc_1);
	deleteHBInt(&temp);

	return RETURN_OK_BINT;
}