This code divides two high precision fixed point numbers. It is part of a larger package providing high precision fixed point support. Other than minor formatting for HTML presentation, this code is exactly as found in the original source code.
This code is made available by First Data Corporation for the restricted purpose of demonstrating C code written by Michael Lee Finney. Neither this code nor extracts thereof may be used in programs outside of First Data Corporation without their express permission.
/* Copyright (c) 1993..1996 by First Data Corporation. All rights reserved. */
/*
* Divide two fixed point numbers:
*
* q <- x / y
* r <- x % y
*
* None of the fixed point numbers may be synonymous. False is
* returned if either division by zero is attempted, or if the scaling
* factor of x is less than the scaling factor of y.
*/
bool _extlink fixPointDivide(
FixPoint * x,
FixPoint * y,
FixPoint * q,
FixPoint * r)
{
signed4 lenV = y->length;
FixPoint p;
signed4 * pu;
signed4 * lu;
signed4 * pv;
signed4 * lv;
signed4 * sv;
signed4 * tv;
signed4 * pq;
signed4 * lq;
signed4 esq;
signed4 esr;
signed4 tq[2];
signed4 tr[2];
signed4 nrml;
signed4 tmp;
signed4 zc;
signed4 lenU;
signed4 lenQ;
bool negV;
bool negQ;
bool negR;
/*
* If the divisor is zero, or the scaling constraints are violated
* then exit immediately, returning an error condition. The division
* will succeed under all other circumstances.
*/
unless (lenV and x->scale >= y->scale)
return false;
/*
* If the dividend is less than the divisor in magnitude then the
* result will always have the quotient as zero and the remainder
* equal to the dividend.
*/
when (fixPointAbsCompare(x, y) < 0)
{
fixPointZero(q);
fixPointCopy(r, x);
return true;
}
lenU = x->length;
negR = (bool)(lenU < 0);
negV = (bool)(lenV < 0);
negQ = (bool)(negR ^ negV);
lenV = labs(lenV);
lenU = labs(lenU);
/*
* If the divisor is a single digit then use the single digit
* division routine.
*/
when (lenV == 1)
{
tmp = fixPointDivideByInteger(q, x, negV ? -(y->value[0]) : y->value[0]);
when (q->length)
q->scale = x->scale - y->scale;
unless (tmp)
{
fixPointZero(r);
return true;
}
if (tmp < 0)
r->length = -1, r->value[0] = -tmp;
else
r->length = 1, r->value[0] = tmp;
r->scale = x->scale;
return true;
}
/*
* In order to make good estimates of individual quotient digits,
* we require that:
*
* v[lenV-1] >= R / 2
*
* where v[lenV-1] is the leading digit of the divisor. If this does
* not hold, then we multiply both dividend and divisor by:
*
* R / (v[lenV-1] + 1).
*
* This may increase the length of u by one, but the length of v will
* not increase. We also require that:
*
* u[lenU-1] < v[lenV-1]
*
* as an initial constraint, so that we always know that the leading lenV
* digits of u are less than v. This constraint will be maintained
* throughout the division. If this constraint is initially violated we
* add a leading zero digit to u. This happens exactly in those cases
* where the length of u was not increased by the normalization step
* above, so the length of u is always increased by one. This may cause
* the leading digit of the quotient to be zero, so an adjustment in the
* length of the quotient may be necessary when the division is complete.
* However, no more than a single leading digit of the quotient may be
* zero.
*/
pv = y->value, lv = pv + lenV;
if (lv[-1] < fixPointRadix / 2)
{
nrml = fixPointRadix / (lv[-1] + 1);
fixPointMultiplyByInteger(r, x, nrml);
fixPointMultiplyByInteger(&p, y, nrml);
lenU = labs(r->length);
pv = p.value, lv = pv + lenV;
}
else
{
nrml = 0;
fixPointCopy(r, x);
}
pu = r->value, lu = pu + lenU;
unless (lu[-1] < lv[-1])
*lu++ = 0, lenU++;
lenQ = lenU - lenV;
pq = q->value, lq = pq + lenQ;
tv = lv--;
while (pq < lq)
{
lq--; lu--;
/*
* Each pass through this loop computes one digit of the quotient.
*
* At this point the leading digits of the partial remainder are:
*
* lu[0], lu[-1], lu[-2], . . ., lu[1-lenV], lu[-lenV] = pu[0]
*
* and the digits of the divisor are:
*
* lv[0], lv[-1], lv[-2], . . ., lv[1-lenV]
*
* and the next digit of the quotient is to placed at lq[0]. We know
* that the divisor is at least two digits because the single digit
* case was handled separately. The length of the partial remainder
* is at least three digits. This can be seen by recognizing that the
* current length of the partial remainder is given by:
*
* lenU + ((lq - pq + 1) - lenQ)
*
* Note that initial length is lenU, and that on the first pass through
* this loop: lq - pq + 1 = lenQ, so the length can be rewritten as:
*
* lenU - (lenU - lenV) + (lq - pq + 1) = lenV + lq - pq + 1
*
* Since: lq >= pq, we have: lq - pq >= 0, and since: lenV >= 2
*
* lenV + lq - pq + 1 >= lenV + 1 >= 3.
*
* Further, the constraint that the leading lenV digits of the partial
* remainder are always less than the divisor means that at each step
* of this algorithm lu[0] will become zero.
*/
if (lu[0] < lv[0])
{
/*
* Estimate the next digit of the quotient based on the
* leading digits of u and v:
*
* esq <- lu[0].lu[-1] / lv[0]
* esr <- lu[0].lu[-1] % lv[0].
*/
fixPointDigDivide(lu[0], lu[-1], lv[0], &esq, &esr);
}
else
{
esq = fixPointRadix - 1;
esr = lu[-1] + lv[0];
}
/*
* The current estimate is known to be in the range:
*
* esq - 2 <= q <= esq
*
* so using the next digit, we will improve the estimate so that
* it is either correct or one too large. The probability of it
* being one too large, on the average, is on the order of:
*
* 2 / b.
*
* To accomplish this, as long as:
*
* esq * lv[-1] > esr * R + lu[-2]
*
* we will reduce the estimate by one.
*/
/*
* Calculate: esq * lv[-1] as: tq[0].tq[1]
*/
tmp = lv[-1];
fixPointDigMult(&tq[0], &tq[1], esq, tmp);
/*
* Calculate: esr * R + lu[-2] as: tr[0].tr[1]
*/
tr[0] = esr;
tr[1] = lu[-2];
while (tq[0] > tr[0] or (tq[0] == tr[0] and tq[1] > tr[1]))
{
/*
* We have: esq * lv[-1] > esr * R + lu[-2]
* so we must reduce eq by one, and adjust tq and tr
* correspondingly.
*/
esq--;
tmp = lv[-1];
fixPointDigSubtractCarry(&tmp, &tq[1], tq[1], tmp = 0;);
fixPointDigSubtractCarry(&tmp, &tq[0], tq[0], tmp = 0;);
tr[0] += lv[0];
}
/*
* We have an estimate for the next digit of the quotient.
* Now we must subtract the divisor multiplied by that digit
* from the partial remainder. If we have a carry then our
* estimate was one too large, so reduce it by one and add
* back the divisor to correct the partial remainder.
*/
when (esq)
{
pu = lu - lenV;
sv = pv;
if (esq > 1)
{
DigitVectorMsub(&zc, pu, sv, tv, esq);
}
else
{
DigitVectorSubtract(&zc, pu, sv, tv);
}
when (zc)
{
pu = lu - lenV;
sv = pv;
DigitVectorAdd(pu, sv, tv);
esq--;
}
}
lq[0] = esq;
}
/*
* Determine the length, sign and scaling of the quotient.
*/
unless (lq[lenQ-1])
lenQ--;
q->length = negQ ? -lenQ : lenQ;
q->scale = x->scale - y->scale;
/*
* Determine the length, sign and scaling of the remainder.
*/
r->length = 1;
until (*lu)
lu--;
r->length = lu - r->value + 1;
if (r->length)
{
when (negR)
fixPointNegate(r);
r->scale = x->scale;
}
else
r->scale = 0;
when (nrml)
fixPointDivideByInteger(r, r, nrml);
return true;
}