大数除法运算，利用乘法优化

在大数除法过程中，如果能有意识的利用乘法运算，应当能够提升计算速度的，因为乘法速度比较快。

比如10000进制表示法的大数。先用比除数多一个数位的基本单位去除以除数，得到基本单位的商d和余数r，那么如果这个数位上是n，用n去乘商和余数，可以得到近似的商nd和余数nr。

假如除数是一个数位的x，那么先用两位数{1,0}，实际是10000，去除x得商d和余数r，那么数{n, 0}除x，会得到商nd和余数nr。这样就得到了大部分的结果了。

用数组来模拟一下就是：

	for(i=0; i<n-1; i++) {
		wc=0;
jump:
		y=a[i]*d;
		z=a[i]*r;
		d[i]+=y;
		a[i+1]+=z;
		if (a[i+1]>=unit) {
			a[i]=a[i+1]/unit;
			a[i+1]%=unit;
			print ++wc;
			goto jump;
		}
	}

下面在y_tab解释程序中尝试一下。为了简单假设为10进制。被除数通过共享数组output[]传去，除数通过参数传进去。

output[]={9, 9, -1};

func div(x) {
	for(i=0; output[i]>=0 &&i<199; i++) {
		a[i]=output[i];
	}
	a[i]=-1; n=i;
	...
}

x是个位数。注意如果变为10000进制，个位数的意思是0~9999。

求出x的表b[i]和 c[i]，b[i]存的是x*2^n，c[i]存的是指数：

	z = 1 *unit;
	b[i=0]=x;
	c[0]=c=1;
	while(z>=x) {
		x+=x;
		b[++i]=x;
		c+=c;
		c[i]=c;
	}
	x=b[0];
	ci=i;

进而求得d和r：

	--i;y=0;
	while(i>=0) {
		if(z>=b[i]) {
			z-=b[i];
			y+=c[i];
		}
		--i;
	}

	d=y; r=z;

然后套用这个结果，得到商d[i]和余数a[i+1]：

	for(i=0; i<n-1; i++) {
		wc=0;
jump:
		y=a[i]*d;
		z=a[i]*r;
		d[i]+=y;
		a[i+1]+=z;
		if (a[i+1]>=unit) {
			a[i]=a[i+1]/unit;
			a[i+1]%=unit;
			print ++wc;
			goto jump;
		}
	}

最后处理一下尾数：

	z=a[i];j=i;
	i=ci;
	--i;y=0;
	while(i>=0) {
		if(z>=b[i]) {
			z-=b[i];
			y+=c[i];
		}
		--i;
	}
	i=j;
	d[i-1]+= y;
	w= z;

并报告结果：

	print x,"\b:", "\b";
	for(i=0; i<n-1; i++) print d[i], "\b";
	print "...", w;

把这些代码连起来就是：

func div(x)
{
	unit=10;

	for(i=0; output[i]>=0 &&i<199; i++) {
		a[i]=output[i];
	}
	a[i]=-1; n=i;

	z = 1 *unit;
	b[i=0]=x;
	c[0]=c=1;
	while(z>=x) {
		x+=x;
		b[++i]=x;
		c+=c;
		c[i]=c;
	}
	x=b[0];
	ci=i;

	--i;y=0;
	while(i>=0) {
		if(z>=b[i]) {
			z-=b[i];
			y+=c[i];
		}
		--i;
	}

	d=y; r=z;
	for(i=0; i<n-1; i++) {
		wc=0;
jump:
		y=a[i]*d;
		z=a[i]*r;
		d[i]+=y;
		a[i+1]+=z;
		if (a[i+1]>=unit) {
			a[i]=a[i+1]/unit;
			a[i+1]%=unit;
			print ++wc;
			goto jump;
		}
	}

	z=a[i];j=i;
	i=ci;
	--i;y=0;
	while(i>=0) {
		if(z>=b[i]) {
			z-=b[i];
			y+=c[i];
		}
		--i;
	}
	i=j;
	d[i-1]+= y;
	w= z;

	print x,"\b:", "\b";
	for(i=0; i<n-1; i++) print d[i], "\b";
	print "...", w;	
}

最后需要处理d[]的进位，处理后d[]传到output[]数组。这部分是商。余数通过return 返回。这里先省略。

看一下运行结果吧：

for(i=2; i<10; i++) div(i);
2: 49 ... 1
1
3: 33 ... 0
1
2
4: 24 ... 3
5: 19 ... 4
1
2
6: 16 ... 3
1
2
7: 14 ... 1
1
2
8: 12 ... 3
1
9: 11 ... 0

如果unit取10000。output[]={9,9,-1}，现在的意思是90009。重新修改上面div函数，修改unit=10000。最后输入：

div(2345);
2345: 38 ... 899

这样也行。得到商38余数899。除数多位数的情况这里不展示了。感兴趣的可以自己改。