shithub: femtolisp

ref: 4e18c3987bb3f2940be37db6ccb3d7986902e3a1
dir: /3rd/mp/mpdiv.c/

#include "platform.h"

// division ala knuth, seminumerical algorithms, pp 237-238
// the numbers are stored backwards to what knuth expects so j
// counts down rather than up.

void
mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder)
{
	int j, s, vn, sign, qsign, rsign;
	mpdigit qd, *up, *vp, *qp;
	mpint *u, *v, *t;

	assert(quotient != remainder);
	assert(divisor->flags & MPnorm);

	// divide by zero
	if(divisor->top == 0){
		fprintf(stderr, "mpdiv: by zero\n");
		exit(2);
	}

	// division by one or small powers of two
	if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){
		int64_t r = 0;
		if(dividend->top > 0)
			r = (int64_t)dividend->sign * (dividend->p[0] & divisor->p[0]-1);
		if(quotient != nil){
			sign = divisor->sign;
			for(s = 0; ((divisor->p[0] >> s) & 1) == 0; s++)
				;
			mpright(dividend, s, quotient);
			if(sign < 0)
				quotient->sign ^= (-mpmagcmp(quotient, mpzero) >> 31) << 1;
		}
		if(remainder != nil){
			remainder->flags |= dividend->flags & MPtimesafe;
			vtomp(r, remainder);
		}
		return;
	}
	assert((dividend->flags & MPtimesafe) == 0);

	// quick check
	if(mpmagcmp(dividend, divisor) < 0){
		if(remainder != nil)
			mpassign(dividend, remainder);
		if(quotient != nil)
			mpassign(mpzero, quotient);
		return;
	}

	qsign = divisor->sign * dividend->sign;
	rsign = dividend->sign;

	// D1: shift until divisor, v, has hi bit set (needed to make trial
	//     divisor accurate)
	qd = divisor->p[divisor->top-1];
	for(s = 0; (qd & mpdighi) == 0; s++)
		qd <<= 1;
	u = mpnew((dividend->top+2)*Dbits + s);
	if(s == 0 && divisor != quotient && divisor != remainder) {
		mpassign(dividend, u);
		v = divisor;
	} else {
		mpleft(dividend, s, u);
		v = mpnew(divisor->top*Dbits);
		mpleft(divisor, s, v);
	}
	up = u->p+u->top-1;
	vp = v->p+v->top-1;
	vn = v->top;

	// D1a: make sure high digit of dividend is less than high digit of divisor
	if(*up >= *vp){
		*++up = 0;
		u->top++;
	}

	// storage for multiplies
	t = mpnew(4*Dbits);

	qp = nil;
	if(quotient != nil){
		mpbits(quotient, (u->top - v->top)*Dbits);
		quotient->top = u->top - v->top;
		qp = quotient->p+quotient->top-1;
	}

	// D2, D7: loop on length of dividend
	for(j = u->top; j > vn; j--){

		// D3: calculate trial divisor
		mpdigdiv(up-1, *vp, &qd);

		// D3a: rule out trial divisors 2 greater than real divisor
		if(vn > 1) for(;;){
			memset(t->p, 0, 3*Dbytes);	// mpvecdigmuladd adds to what's there
			mpvecdigmuladd(vp-1, 2, qd, t->p);
			if(mpveccmp(t->p, 3, up-2, 3) > 0)
				qd--;
			else
				break;
		}

		// D4: u -= v*qd << j*Dbits
		sign = mpvecdigmulsub(v->p, vn, qd, up-vn);
		if(sign < 0){

			// D6: trial divisor was too high, add back borrowed
			//     value and decrease divisor
			mpvecadd(up-vn, vn+1, v->p, vn, up-vn);
			qd--;
		}

		// D5: save quotient digit
		if(qp != nil)
			*qp-- = qd;

		// push top of u down one
		u->top--;
		*up-- = 0;
	}
	if(qp != nil){
		assert((quotient->flags & MPtimesafe) == 0);
		mpnorm(quotient);
		if(quotient->top != 0)
			quotient->sign = qsign;
	}

	if(remainder != nil){
		assert((remainder->flags & MPtimesafe) == 0);
		mpright(u, s, remainder);	// u is the remainder shifted
		if(remainder->top != 0)
			remainder->sign = rsign;
	}

	mpfree(t);
	mpfree(u);
	if(v != divisor)
		mpfree(v);
}