bn_div.c

实现大整数的除法运算

/************************************************************
  Copyright (C), 2004, Aerospace Information. Co., Ltd.
  FileName: bn_div.cpp
  Author: LSX      Version :  1.0        Date:2004/9/16
  Description:     大数除法运算，根据OpenSSL源码改编      
  Version:         1.0
  Function List:
    1. bn_div      两个大整数做除法运算
  History:         
      <author>    <time>     <version >   <desc>
	    LSX     2004/9/16       1.0      加入注释
************************************************************/

#include <stdio.h>
//#include <openssl/bn.h>
#include "bn_lcl.h"

/* The old slow way */
#if 0

/**************************************************************************************
  Function:     BN_div(BIGNUM *dv, BIGNUM *rem, const BIGNUM *m, const BIGNUM *d,
	                   BN_CTX *ctx)
  Description:  两个大整数做除法运算
  Calls:        bn_check_top, BN_is_zero, BN_ucmp, BN_copy, BN_num_bits, BN_zero, 
                BN_wexpand, BN_lshift, BN_lshift1, BN_usub, BN_rshift1,
				BN_CTX_start, BN_CTX_get, BN_CTX_end
  Called by:    
  Input:        两个大整数m和d, d不为0
  Output:       dv = m / d
                rem = m % d
  Return:       1, if 正确相除
                0, otherwise
  Others:         
**************************************************************************************/
int BN_div(BIGNUM *dv, BIGNUM *rem, const BIGNUM *m, const BIGNUM *d, BN_CTX *ctx)
{
	int i,nm,nd;
	int ret = 0;
	BIGNUM *D;

	bn_check_top(m);
	bn_check_top(d);
	if (BN_is_zero(d))
	{
		fprintf(stderr,"BN_div by zero error!\n");
		return(0);
	}

	if (BN_ucmp(m,d) < 0)  /*m<d, we can do directly*/
	{
		if (rem != NULL)
		{ 
			if (BN_copy(rem,m) == NULL) 
			{
				return(0); 
			}
		}
		if (dv != NULL) BN_zero(dv);
		return(1);
	}

	BN_CTX_start(ctx);
	D = BN_CTX_get(ctx);
	if (dv == NULL) 
	{
		dv = BN_CTX_get(ctx);
	}
	if (rem == NULL) 
	{
		rem = BN_CTX_get(ctx);
	}
	if (D == NULL || dv == NULL || rem == NULL)
	{
		goto end;
	}

	nd=BN_num_bits(d);
	nm=BN_num_bits(m);
	if (BN_copy(D,d) == NULL) 
	{
		goto end;
	}
	if (BN_copy(rem,m) == NULL)
	{
		goto end;
	}

	/* The next 2 are needed so we can do a dv->d[0]|=1 later
	 * since BN_lshift1 will only work once there is a value :-) */
	BN_zero(dv);
	bn_wexpand(dv,1);
	dv->top=1;

	if (!BN_lshift(D,D,nm-nd))  /*使得除数与被除数的bits一样*/ 
	{
		goto end;
	}
	
	for (i=nm-nd; i>=0; i--)
	{
		if (!BN_lshift1(dv,dv)) 
		{
			goto end;
		}
		if (BN_ucmp(rem,D) >= 0)
		{
			dv->d[0]|=1;
			if (!BN_usub(rem,rem,D))
			{
				goto end;
			}
		}

		if (!BN_rshift1(D,D)) 
		{
			goto end;
		}
	}
	rem->neg=BN_is_zero(rem)?0:m->neg;
	dv->neg=m->neg^d->neg;
	ret = 1;
 end:
	BN_CTX_end(ctx);
	return(ret);
}

#else

/* BN_div computes  dv := num / divisor,  rounding towards zero, and sets up
 * rm  such that  dv*divisor + rm = num  holds.
 * Thus:
 *     dv->neg == num->neg ^ divisor->neg  (unless the result is zero)
 *     rm->neg == num->neg                 (unless the remainder is zero)
 * If 'dv' or 'rm' is NULL, the respective value is not returned.
 */

 /**************************************************************************************
  Function:     BN_div(BIGNUM *dv, BIGNUM *rem, const BIGNUM *m, const BIGNUM *d,
	                   BN_CTX *ctx)
  Description:  两个大整数做除法运算
  Calls:        bn_check_top, BN_is_zero, BN_ucmp, BN_copy, BN_num_bits, BN_zero, 
                BN_wexpand, BN_lshift, BN_lshift1, BN_usub, BN_rshift1,
				BN_CTX_start, BN_CTX_get, BN_CTX_end
  Called by:    
  Input:        两个大整数m和d, d不为0
  Output:       dv = m / d
                rem = m % d
  Return:       1, if 正确相除
                0, otherwise
  Others:         
**************************************************************************************/
int BN_div(BIGNUM *dv, BIGNUM *rm, const BIGNUM *num, const BIGNUM *divisor, BN_CTX *ctx)
{
	int norm_shift,i,j,loop;
	BIGNUM *tmp,wnum,*snum,*sdiv,*res;
	BN_ULONG *resp,*wnump;
	BN_ULONG d0,d1;
	int num_n,div_n;

	bn_check_top(num);
	bn_check_top(divisor);

	if (BN_is_zero(divisor))
	{
		fprintf(stderr,"BN_div by zero error!\n");
		return(0);
	}

	if (BN_ucmp(num,divisor) < 0)
	{
		if (rm != NULL)
		{
			if (BN_copy(rm,num) == NULL) 
			{
				return(0); 
			}
		}
		if (dv != NULL) 
		{
			BN_zero(dv);
		}
		return(1);
	}

	BN_CTX_start(ctx);
	tmp=BN_CTX_get(ctx);
	snum=BN_CTX_get(ctx);
	sdiv=BN_CTX_get(ctx);
	if (dv == NULL)
	{
		res=BN_CTX_get(ctx);
	}
	else
	{
		res=dv;
	}
	if (sdiv == NULL || res == NULL)
	{
		goto err;
	}
	tmp->neg=0;

	/* First we normalise the numbers */
	norm_shift=BN_BITS2-((BN_num_bits(divisor))%BN_BITS2);
	/* add by ZQS,no need normalist! */
	if (norm_shift==BN_BITS2)
	{
		norm_shift=0;
	}
	if (!(BN_lshift(sdiv,divisor,norm_shift))) 
	{
		goto err;
	}
	sdiv->neg=0;
	norm_shift+=BN_BITS2;
	if (!(BN_lshift(snum,num,norm_shift))) 
	{
		goto err;
	}
	snum->neg=0;
	div_n=sdiv->top;
	num_n=snum->top;
	loop=num_n-div_n;

	/* Lets setup a 'window' into snum
	 * This is the part that corresponds to the current
	 * 'area' being divided */
	BN_init(&wnum);
	wnum.d=	 &(snum->d[loop]);
	wnum.top= div_n;
	wnum.dmax= snum->dmax+1; /* a bit of a lie */

	/* Get the top 2 words of sdiv */
	/* i=sdiv->top; */
	d0=sdiv->d[div_n-1];
	d1=(div_n == 1)?0:sdiv->d[div_n-2];

	/* pointer to the 'top' of snum */
	wnump= &(snum->d[num_n-1]);

	/* Setup to 'res' */
	res->neg= (num->neg^divisor->neg);
	if (!bn_wexpand(res,(loop+1))) 
	{
		goto err;
	}
	res->top=loop;
	resp= &(res->d[loop-1]);

	/* space for temp */
	if (!bn_wexpand(tmp,(div_n+1))) 
	{
		goto err;
	}

	if (BN_ucmp(&wnum,sdiv) >= 0)
	{
		if (!BN_usub(&wnum,&wnum,sdiv))
		{
			goto err;
		}
		*resp=1;
		res->d[res->top-1]=1;
	}
	else
	{
		res->top--;
	}
	if (res->top == 0)
	{
		res->neg = 0;
	}
	resp--;

	for (i=0; i<loop-1; i++)
	{
		BN_ULONG q,l0;
		BN_ULONG n0,n1,rem=0;

		n0=wnump[0];
		n1=wnump[-1];
		if (n0 == d0)
		{
			q=BN_MASK2;
		}
		else  /* n0 < d0 */
		{
			BN_ULLONG t2;

			q=(BN_ULONG)(((((BN_ULLONG)n0)<<BN_BITS2)|n1)/d0);
			/*
			 * rem doesn't have to be BN_ULLONG. The least we
			 * know it's less that d0, isn't it?
			 */
			rem=(n1-q*d0)&BN_MASK2;
			t2=(BN_ULLONG)d1*q;

			for (;;)
			{
				if (t2 <= ((((BN_ULLONG)rem)<<BN_BITS2)|wnump[-2]))
				{
					break;
				}
				q--;
				rem += d0;
				if (rem < d0) 
				{
					break;  /* don't let rem overflow */
				}
				t2 -= d1;
			}
		}

		l0=bn_mul_words(tmp->d,sdiv->d,div_n,q);
		wnum.d--; wnum.top++;
		tmp->d[div_n]=l0;
		for (j=div_n+1; j>0; j--)
		{
			if (tmp->d[j-1])
			{
				break;
			}
		}
		tmp->top=j;

		j=wnum.top;
		if (!BN_sub(&wnum,&wnum,tmp)) 
		{
			goto err;
		}

		snum->top=snum->top+wnum.top-j;

		if (wnum.neg)
		{
			q--;
			j=wnum.top;
			if (!BN_add(&wnum,&wnum,sdiv)) 
			{
				goto err;
			}
			snum->top+=wnum.top-j;
		}
		*(resp--)=q;
		wnump--;
	}
	if (rm != NULL)
	{
		/* Keep a copy of the neg flag in num because if rm==num
		 * BN_rshift() will overwrite it.
		 */
		int neg = num->neg;
		BN_rshift(rm,snum,norm_shift);
		if (!BN_is_zero(rm))
		{
			rm->neg = neg;
		}
	}
	BN_CTX_end(ctx);
	return(1);
err:
	BN_CTX_end(ctx);
	return(0);
}

#endif