in commons-math-legacy-core/src/main/java/org/apache/commons/math4/legacy/core/dfp/Dfp.java [1682:1888]
public Dfp divide(Dfp divisor) {
int[] dividend; // current status of the dividend
int[] quotient; // quotient
int[] remainder; // remainder
int qd; // current quotient digit we're working with
int nsqd; // number of significant quotient digits we have
int trial = 0; // trial quotient digit
int minadj; // minimum adjustment
boolean trialgood; // Flag to indicate a good trail digit
int md = 0; // most sig digit in result
int excp; // exceptions
// make sure we don't mix number with different precision
if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
final Dfp result = newInstance(getZero());
result.nans = QNAN;
return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
}
Dfp result = newInstance(getZero());
/* handle special cases */
if (nans != FINITE || divisor.nans != FINITE) {
if (isNaN()) {
return this;
}
if (divisor.isNaN()) {
return divisor;
}
if (nans == INFINITE && divisor.nans == FINITE) {
result = newInstance(this);
result.sign = (byte) (sign * divisor.sign);
return result;
}
if (divisor.nans == INFINITE && nans == FINITE) {
result = newInstance(getZero());
result.sign = (byte) (sign * divisor.sign);
return result;
}
if (divisor.nans == INFINITE && nans == INFINITE) {
field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
result = newInstance(getZero());
result.nans = QNAN;
result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
return result;
}
}
/* Test for divide by zero */
if (divisor.mant[mant.length - 1] == 0) {
field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
result = newInstance(getZero());
result.sign = (byte) (sign * divisor.sign);
result.nans = INFINITE;
result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
return result;
}
dividend = new int[mant.length + 1]; // one extra digit needed
quotient = new int[mant.length + 2]; // two extra digits needed 1 for overflow, 1 for rounding
remainder = new int[mant.length + 1]; // one extra digit needed
/* Initialize our most significant digits to zero */
dividend[mant.length] = 0;
quotient[mant.length] = 0;
quotient[mant.length + 1] = 0;
remainder[mant.length] = 0;
/* copy our mantissa into the dividend, initialize the quotient while we are at it */
for (int i = 0; i < mant.length; i++) {
dividend[i] = mant[i];
quotient[i] = 0;
remainder[i] = 0;
}
/* outer loop. Once per quotient digit */
nsqd = 0;
for (qd = mant.length + 1; qd >= 0; qd--) {
/* Determine outer limits of our quotient digit */
// r = most sig 2 digits of dividend
final int divMsb = dividend[mant.length] * RADIX + dividend[mant.length - 1];
int min = divMsb / (divisor.mant[mant.length - 1] + 1);
int max = (divMsb + 1) / divisor.mant[mant.length - 1];
trialgood = false;
while (!trialgood) {
// try the mean
trial = (min + max) / 2;
/* Multiply by divisor and store as remainder */
int rh = 0;
for (int i = 0; i < mant.length + 1; i++) {
int dm = (i < mant.length) ? divisor.mant[i] : 0;
final int r = (dm * trial) + rh;
rh = r / RADIX;
remainder[i] = r - rh * RADIX;
}
/* subtract the remainder from the dividend */
rh = 1; // carry in to aid the subtraction
for (int i = 0; i < mant.length + 1; i++) {
final int r = ((RADIX - 1) - remainder[i]) + dividend[i] + rh;
rh = r / RADIX;
remainder[i] = r - rh * RADIX;
}
/* Lets analyze what we have here */
if (rh == 0) {
// trial is too big -- negative remainder
max = trial - 1;
continue;
}
/* find out how far off the remainder is telling us we are */
minadj = (remainder[mant.length] * RADIX) + remainder[mant.length - 1];
minadj /= divisor.mant[mant.length - 1] + 1;
if (minadj >= 2) {
min = trial + minadj; // update the minimum
continue;
}
/* May have a good one here, check more thoroughly. Basically its a good
* one if it is less than the divisor */
trialgood = false; // assume false
for (int i = mant.length - 1; i >= 0; i--) {
if (divisor.mant[i] > remainder[i]) {
trialgood = true;
break;
}
if (divisor.mant[i] < remainder[i]) {
break;
}
}
if (remainder[mant.length] != 0) {
trialgood = false;
}
if (!trialgood) {
min = trial + 1;
}
}
/* Great we have a digit! */
quotient[qd] = trial;
if (trial != 0 || nsqd != 0) {
nsqd++;
}
if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
// We have enough for this mode
break;
}
if (nsqd > mant.length) {
// We have enough digits
break;
}
/* move the remainder into the dividend while left shifting */
dividend[0] = 0;
System.arraycopy(remainder, 0, dividend, 1, mant.length);
}
/* Find the most sig digit */
md = mant.length; // default
for (int i = mant.length + 1; i >= 0; i--) {
if (quotient[i] != 0) {
md = i;
break;
}
}
/* Copy the digits into the result */
for (int i = 0; i < mant.length; i++) {
result.mant[mant.length - i - 1] = quotient[md - i];
}
/* Fixup the exponent. */
result.exp = exp - divisor.exp + md - mant.length;
result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
if (result.mant[mant.length - 1] == 0) { // if result is zero, set exp to zero
result.exp = 0;
}
if (md > (mant.length - 1)) {
excp = result.round(quotient[md - mant.length]);
} else {
excp = result.round(0);
}
if (excp != 0) {
result = dotrap(excp, DIVIDE_TRAP, divisor, result);
}
return result;
}