(a * b) / c MulDiv 和处理中间乘法溢出
(a * b) / c MulDiv and dealing with overflow from intermediate multiplication
我需要做以下运算：
long a,b,c;
long result = a*b/c;
虽然结果是 gua运行teed 适合 long，但乘法不是，所以它可能会溢出。
我试着一步一步地做（先乘后除），同时通过将 a*b 的中间结果拆分成一个最大为 4 的 int 数组（很像 BigInteger正在使用其 int[] mag 变量）。
在这里我被部门卡住了。我无法理解进行精确除法所需的位移位。我只需要商（不需要余数）。
假设的方法是：
public static long divide(int[] dividend, long divisor)
此外，我不考虑使用 BigInteger，因为这部分代码需要快速（我想坚持使用基元和基元数组）。
如有任何帮助，我们将不胜感激！
编辑：
我并不是要自己实现整个 BigInteger。我想做的是比使用通用 BigInteger.
更快地解决特定问题（a*b/c，其中 a*b 可能溢出）
Edit2：如果能以 聪明 的方式完成，那将是理想的，完全不会溢出，评论中出现了一些提示，但我仍在寻找一个是正确的。
更新：
我尝试将 BigInteger 代码移植到我的特定需求，而不创建对象，并且在第一次迭代中，与使用 BigInteger（在我的开发电脑上）相比，我的速度提高了约 46%。
然后我尝试了一点修改的@David Eisenstat 解决方案，它给了我 ~56%（我 运行 100_000_000_000 运行dom 输入从 Long.MIN_VALUE 到 Long.MAX_VALUE) 与 BigInteger 相比减少了 运行 倍（超过 2 倍）（与我采用的 BigInteger 算法相比大约减少了 18%）。
还会有更多的优化和测试迭代，但在这一点上，我想我必须接受这个答案是最好的。
将a/c和b/c分成整数和小数（余数）部分，那么你有：
a*b/c
= c * a/c * b/c
= c * (x/c + y/c) * (z/c + w/c)
= xz/c + xw/c + yz/c + yw/c where x and z are multiples of c
因此，您可以轻松计算前三个因子而不会溢出。根据我的经验，这通常足以涵盖典型的溢出情况。但是，如果您的 除数 太大，以至于 (a % c) * (b % c) 溢出，此方法仍然失败。如果这对你来说是一个典型的问题，你可能想看看其他方法（例如，将 a 和 b 中的最大者以及 c 除以 2，直到你不再有溢出，但是如何做到这一点而不引入额外的错误，因为这个过程中的偏差是非常重要的——你可能需要在一个单独的变量中保持 运行 的错误分数）
无论如何，上面的代码：
long a,b,c;
long bMod = (b % c)
long result = a * (b / c) + (a / c) * bMod + ((a % c) * bMod) / c;
如果速度是一个大问题（我假设至少在某种程度上是这样，因为你问这个），你可能需要考虑将 a/c 和 b/c 存储在变量并通过乘法计算 mod，例如将 (a % c) 替换为 (a - aDiv * c)——这允许您将每次调用的 4 个分区改为 2 个。
您可以使用最大公约数 (gcd) 来提供帮助。
a * b / c = (a / gcd(a,c)) * (b / (c / gcd(a,c)))
编辑：OP 让我解释上面的等式。基本上，我们有：
a = (a / gcd(a,c)) * gcd(a,c)
c = (c / gcd(a,c)) * gcd(a,c)
Let's say x=gcd(a,c) for brevity, and rewrite this.
a*b/c = (a/x) * x * b
--------------
(c/x) * x
Next, we cancel
a*b/c = (a/x) * b
----------
(c/x)
您可以更进一步。让 y = gcd(b, c/x)
a*b/c = (a/x) * (b/y) * y
------------------
((c/x)/y) * y
a*b/c = (a/x) * (b/y)
------------
(c/(xy))
这是获取 gcd 的代码。
static long gcd(long a, long b)
{
if (b == 0)
return a;
return gcd(b, a % b);
}
我一直在修补一种方法，即 (1) 将 a 和 b 与 21 位肢体上的学校算法相乘 (2) 继续进行长除以 c，残差 a*b - c*q 的不寻常表示使用 double 存储高位，使用 long 存储低位。我不知道它是否可以与标准长除法竞争，但为了您的享受，
public class MulDiv {
public static void main(String[] args) {
java.util.Random r = new java.util.Random();
for (long i = 0; true; i++) {
if (i % 1000000 == 0) {
System.err.println(i);
}
long a = r.nextLong() >> (r.nextInt(8) * 8);
long b = r.nextLong() >> (r.nextInt(8) * 8);
long c = r.nextLong() >> (r.nextInt(8) * 8);
if (c == 0) {
continue;
}
long x = mulDiv(a, b, c);
java.math.BigInteger aa = java.math.BigInteger.valueOf(a);
java.math.BigInteger bb = java.math.BigInteger.valueOf(b);
java.math.BigInteger cc = java.math.BigInteger.valueOf(c);
java.math.BigInteger xx = aa.multiply(bb).divide(cc);
if (java.math.BigInteger.valueOf(xx.longValue()).equals(xx) && x != xx.longValue()) {
System.out.printf("a=%d b=%d c=%d: %d != %s\n", a, b, c, x, xx);
}
}
}
// Returns truncate(a b/c), subject to the precondition that the result is
// defined and can be represented as a long.
private static long mulDiv(long a, long b, long c) {
// Decompose a.
long a2 = a >> 42;
long a10 = a - (a2 << 42);
long a1 = a10 >> 21;
long a0 = a10 - (a1 << 21);
assert a == (((a2 << 21) + a1) << 21) + a0;
// Decompose b.
long b2 = b >> 42;
long b10 = b - (b2 << 42);
long b1 = b10 >> 21;
long b0 = b10 - (b1 << 21);
assert b == (((b2 << 21) + b1) << 21) + b0;
// Compute a b.
long ab4 = a2 * b2;
long ab3 = a2 * b1 + a1 * b2;
long ab2 = a2 * b0 + a1 * b1 + a0 * b2;
long ab1 = a1 * b0 + a0 * b1;
long ab0 = a0 * b0;
// Compute a b/c.
DivBy d = new DivBy(c);
d.shift21Add(ab4);
d.shift21Add(ab3);
d.shift21Add(ab2);
d.shift21Add(ab1);
d.shift21Add(ab0);
return d.getQuotient();
}
}
public strictfp class DivBy {
// Initializes n <- 0.
public DivBy(long d) {
di = d;
df = (double) d;
oneOverD = 1.0 / df;
}
// Updates n <- 2^21 n + i. Assumes |i| <= 3 (2^42).
public void shift21Add(long i) {
// Update the quotient and remainder.
q <<= 21;
ri = (ri << 21) + i;
rf = rf * (double) (1 << 21) + (double) i;
reduce();
}
// Returns truncate(n/d).
public long getQuotient() {
while (rf != (double) ri) {
reduce();
}
// Round toward zero.
if (q > 0) {
if ((di > 0 && ri < 0) || (di < 0 && ri > 0)) {
return q - 1;
}
} else if (q < 0) {
if ((di > 0 && ri > 0) || (di < 0 && ri < 0)) {
return q + 1;
}
}
return q;
}
private void reduce() {
// x is approximately r/d.
long x = Math.round(rf * oneOverD);
q += x;
ri -= di * x;
rf = repairLowOrderBits(rf - df * (double) x, ri);
}
private static double repairLowOrderBits(double f, long i) {
int e = Math.getExponent(f);
if (e < 64) {
return (double) i;
}
long rawBits = Double.doubleToRawLongBits(f);
long lowOrderBits = (rawBits >> 63) ^ (rawBits << (e - 52));
return f + (double) (i - lowOrderBits);
}
private final long di;
private final double df;
private final double oneOverD;
private long q = 0;
private long ri = 0;
private double rf = 0;
}
让我想多了。
我希望简单的案例更快：让 double 来处理。
Newton-Raphson 可能是其他人更好的选择。
/** Multiplies both <code>factor</code>s
* and divides by <code>divisor</code>.
* @return <code>Long.MIN_VALUE</code> if result out of range,<br/>
* else <code>factorA * factor1 / divisor</code> */
public static long
mulDiv(long factorA, long factor1, long divisor) {
final double dd = divisor,
product = (double)factorA * factor1,
a1_d = product / dd;
if (a1_d < -TOO_LARGE || TOO_LARGE < a1_d)
return tooLarge();
if (-ONE_ < a1_d && a1_d < ONE_)
return 0;
if (-EXACT < product && product < EXACT)
return (long) a1_d;
long pLo = factorA * factor1, //diff,
pHi = high64(factorA, factor1);
if (a1_d < -LONG_MAX_ || LONG_MAX_ < a1_d) {
long maxdHi = divisor >> 1;
if (maxdHi < pHi
|| maxdHi == pHi
&& Long.compareUnsigned((divisor << Long.SIZE-1),
pLo) <= 0)
return tooLarge();
}
final double high_dd = TWO_POWER64/dd;
long quotient = (long) a1_d,
loPP = quotient * divisor,
hiPP = high64(quotient, divisor);
long remHi = pHi - hiPP, // xxx overflow/carry
remLo = pLo - loPP;
if (Long.compareUnsigned(pLo, remLo) < 0)
remHi -= 1;
double fudge = remHi * high_dd;
if (remLo < 0)
fudge += high_dd;
fudge += remLo/dd;
long //fHi = (long)fudge/TWO_POWER64,
fLo = (long) Math.floor(fudge); //*round
quotient += fLo;
loPP = quotient * divisor;
hiPP = high64(quotient, divisor);
remHi = pHi - hiPP; // should be 0?!
remLo = pLo - loPP;
if (Long.compareUnsigned(pLo, remLo) < 0)
remHi -= 1;
if (0 == remHi && 0 <= remLo && remLo < divisor)
return quotient;
fudge = remHi * high_dd;
if (remLo < 0)
fudge += high_dd;
fudge += remLo/dd;
fLo = (long) Math.floor(fudge);
return quotient + fLo;
}
/** max <code>double</code> trusted to represent
* a value in the range of <code>long</code> */
static final double
LONG_MAX_ = Double.valueOf(Long.MAX_VALUE - 0xFFF);
/** max <code>double</code> trusted to represent a value below 1 */
static final double
ONE_ = Double.longBitsToDouble(
Double.doubleToRawLongBits(1) - 4);
/** max <code>double</code> trusted to represent a value exactly */
static final double
EXACT = Long.MAX_VALUE >> 12;
static final double
TWO_POWER64 = Double.valueOf(1L<<32)*Double.valueOf(1L<<32);
static long tooLarge() {
// throw new RuntimeException("result too large for long");
return Long.MIN_VALUE;
}
static final long ONES_32 = ~(~0L << 32);
static long high64(long factorA, long factor1) {
long loA = factorA & ONES_32,
hiA = factorA >>> 32,
lo1 = factor1 & ONES_32,
hi1 = factor1 >>> 32;
return ((loA * lo1 >>> 32)
+loA * hi1 + hiA * lo1 >>> 32)
+ hiA * hi1;
}
（我从 IDE 中重新安排了一些代码，使 mulDiv() 位于顶部。
由于懒惰，我有一个用于符号处理的包装器 - 可能会在地狱冻结之前尝试并正确地完成它。
对于计时，迫切需要输入模型：
如何使得每个可能的结果都是同样可能的？）
也许不聪明，但有线性结果时间
#define MUL_DIV_TYPE unsigned int
#define BITS_PER_TYPE (sizeof(MUL_DIV_TYPE)*8)
#define TOP_BIT_TYPE (1<<(BITS_PER_TYPE-1))
//
// result = ( a * b ) / c, without intermediate overflow.
//
MUL_DIV_TYPE mul_div( MUL_DIV_TYPE a, MUL_DIV_TYPE b, MUL_DIV_TYPE c ) {
MUL_DIV_TYPE st, sb; // product sum top and bottom
MUL_DIV_TYPE d, e; // division result
MUL_DIV_TYPE i, // bit counter
j; // overflow check
st = 0;
sb = 0;
d = 0;
e = 0;
for( i = 0; i < BITS_PER_TYPE; i++ ) {
//
// Shift sum left to make space
// for next partial sum
//
st <<= 1;
if( sb & TOP_BIT_TYPE ) st |= 1;
sb <<= 1;
//
// Add a to s if top bit on b
// is set.
//
if( b & TOP_BIT_TYPE ) {
j = sb;
sb += a;
if( sb < j ) st++;
}
//
// Division.
//
d <<= 1;
if( st >= c ) {
d |= 1;
st -= c;
e++;
}
else {
if( e ) e++;
}
//
// Shift b up by one bit.
//
b <<= 1;
}
//
// Roll in missing bits.
//
for( i = e; i < BITS_PER_TYPE; i++ ) {
//
// Shift across product sum
//
st <<= 1;
if( sb & TOP_BIT_TYPE ) st |= 1;
sb <<= 1;
//
// Division, continued.
//
d <<= 1;
if( st >= c ) {
d |= 1;
st -= c;
}
}
return( d ); // remainder should be in st
}
我需要做以下运算：
long a,b,c;
long result = a*b/c;
虽然结果是 gua运行teed 适合 long，但乘法不是，所以它可能会溢出。
我试着一步一步地做（先乘后除），同时通过将 a*b 的中间结果拆分成一个最大为 4 的 int 数组（很像 BigInteger正在使用其 int[] mag 变量）。
在这里我被部门卡住了。我无法理解进行精确除法所需的位移位。我只需要商（不需要余数）。
假设的方法是：
public static long divide(int[] dividend, long divisor)
此外，我不考虑使用 BigInteger，因为这部分代码需要快速（我想坚持使用基元和基元数组）。
如有任何帮助，我们将不胜感激！
编辑：
我并不是要自己实现整个 BigInteger。我想做的是比使用通用 BigInteger.
a*b/c，其中 a*b 可能溢出）
Edit2：如果能以 聪明 的方式完成，那将是理想的，完全不会溢出，评论中出现了一些提示，但我仍在寻找一个是正确的。
更新： 我尝试将 BigInteger 代码移植到我的特定需求，而不创建对象，并且在第一次迭代中，与使用 BigInteger（在我的开发电脑上）相比，我的速度提高了约 46%。
然后我尝试了一点修改的@David Eisenstat 解决方案，它给了我 ~56%（我 运行 100_000_000_000 运行dom 输入从 Long.MIN_VALUE 到 Long.MAX_VALUE) 与 BigInteger 相比减少了 运行 倍（超过 2 倍）（与我采用的 BigInteger 算法相比大约减少了 18%）。
还会有更多的优化和测试迭代，但在这一点上，我想我必须接受这个答案是最好的。
将a/c和b/c分成整数和小数（余数）部分，那么你有：
a*b/c
= c * a/c * b/c
= c * (x/c + y/c) * (z/c + w/c)
= xz/c + xw/c + yz/c + yw/c where x and z are multiples of c
因此，您可以轻松计算前三个因子而不会溢出。根据我的经验，这通常足以涵盖典型的溢出情况。但是，如果您的 除数 太大，以至于 (a % c) * (b % c) 溢出，此方法仍然失败。如果这对你来说是一个典型的问题，你可能想看看其他方法（例如，将 a 和 b 中的最大者以及 c 除以 2，直到你不再有溢出，但是如何做到这一点而不引入额外的错误，因为这个过程中的偏差是非常重要的——你可能需要在一个单独的变量中保持 运行 的错误分数）
无论如何，上面的代码：
long a,b,c;
long bMod = (b % c)
long result = a * (b / c) + (a / c) * bMod + ((a % c) * bMod) / c;
如果速度是一个大问题（我假设至少在某种程度上是这样，因为你问这个），你可能需要考虑将 a/c 和 b/c 存储在变量并通过乘法计算 mod，例如将 (a % c) 替换为 (a - aDiv * c)——这允许您将每次调用的 4 个分区改为 2 个。
您可以使用最大公约数 (gcd) 来提供帮助。
a * b / c = (a / gcd(a,c)) * (b / (c / gcd(a,c)))
编辑：OP 让我解释上面的等式。基本上，我们有：
a = (a / gcd(a,c)) * gcd(a,c)
c = (c / gcd(a,c)) * gcd(a,c)
Let's say x=gcd(a,c) for brevity, and rewrite this.
a*b/c = (a/x) * x * b
--------------
(c/x) * x
Next, we cancel
a*b/c = (a/x) * b
----------
(c/x)
您可以更进一步。让 y = gcd(b, c/x)
a*b/c = (a/x) * (b/y) * y
------------------
((c/x)/y) * y
a*b/c = (a/x) * (b/y)
------------
(c/(xy))
这是获取 gcd 的代码。
static long gcd(long a, long b)
{
if (b == 0)
return a;
return gcd(b, a % b);
}
我一直在修补一种方法，即 (1) 将 a 和 b 与 21 位肢体上的学校算法相乘 (2) 继续进行长除以 c，残差 a*b - c*q 的不寻常表示使用 double 存储高位，使用 long 存储低位。我不知道它是否可以与标准长除法竞争，但为了您的享受，
public class MulDiv {
public static void main(String[] args) {
java.util.Random r = new java.util.Random();
for (long i = 0; true; i++) {
if (i % 1000000 == 0) {
System.err.println(i);
}
long a = r.nextLong() >> (r.nextInt(8) * 8);
long b = r.nextLong() >> (r.nextInt(8) * 8);
long c = r.nextLong() >> (r.nextInt(8) * 8);
if (c == 0) {
continue;
}
long x = mulDiv(a, b, c);
java.math.BigInteger aa = java.math.BigInteger.valueOf(a);
java.math.BigInteger bb = java.math.BigInteger.valueOf(b);
java.math.BigInteger cc = java.math.BigInteger.valueOf(c);
java.math.BigInteger xx = aa.multiply(bb).divide(cc);
if (java.math.BigInteger.valueOf(xx.longValue()).equals(xx) && x != xx.longValue()) {
System.out.printf("a=%d b=%d c=%d: %d != %s\n", a, b, c, x, xx);
}
}
}
// Returns truncate(a b/c), subject to the precondition that the result is
// defined and can be represented as a long.
private static long mulDiv(long a, long b, long c) {
// Decompose a.
long a2 = a >> 42;
long a10 = a - (a2 << 42);
long a1 = a10 >> 21;
long a0 = a10 - (a1 << 21);
assert a == (((a2 << 21) + a1) << 21) + a0;
// Decompose b.
long b2 = b >> 42;
long b10 = b - (b2 << 42);
long b1 = b10 >> 21;
long b0 = b10 - (b1 << 21);
assert b == (((b2 << 21) + b1) << 21) + b0;
// Compute a b.
long ab4 = a2 * b2;
long ab3 = a2 * b1 + a1 * b2;
long ab2 = a2 * b0 + a1 * b1 + a0 * b2;
long ab1 = a1 * b0 + a0 * b1;
long ab0 = a0 * b0;
// Compute a b/c.
DivBy d = new DivBy(c);
d.shift21Add(ab4);
d.shift21Add(ab3);
d.shift21Add(ab2);
d.shift21Add(ab1);
d.shift21Add(ab0);
return d.getQuotient();
}
}
public strictfp class DivBy {
// Initializes n <- 0.
public DivBy(long d) {
di = d;
df = (double) d;
oneOverD = 1.0 / df;
}
// Updates n <- 2^21 n + i. Assumes |i| <= 3 (2^42).
public void shift21Add(long i) {
// Update the quotient and remainder.
q <<= 21;
ri = (ri << 21) + i;
rf = rf * (double) (1 << 21) + (double) i;
reduce();
}
// Returns truncate(n/d).
public long getQuotient() {
while (rf != (double) ri) {
reduce();
}
// Round toward zero.
if (q > 0) {
if ((di > 0 && ri < 0) || (di < 0 && ri > 0)) {
return q - 1;
}
} else if (q < 0) {
if ((di > 0 && ri > 0) || (di < 0 && ri < 0)) {
return q + 1;
}
}
return q;
}
private void reduce() {
// x is approximately r/d.
long x = Math.round(rf * oneOverD);
q += x;
ri -= di * x;
rf = repairLowOrderBits(rf - df * (double) x, ri);
}
private static double repairLowOrderBits(double f, long i) {
int e = Math.getExponent(f);
if (e < 64) {
return (double) i;
}
long rawBits = Double.doubleToRawLongBits(f);
long lowOrderBits = (rawBits >> 63) ^ (rawBits << (e - 52));
return f + (double) (i - lowOrderBits);
}
private final long di;
private final double df;
private final double oneOverD;
private long q = 0;
private long ri = 0;
private double rf = 0;
}
我希望简单的案例更快：让 double 来处理。
Newton-Raphson 可能是其他人更好的选择。
/** Multiplies both <code>factor</code>s
* and divides by <code>divisor</code>.
* @return <code>Long.MIN_VALUE</code> if result out of range,<br/>
* else <code>factorA * factor1 / divisor</code> */
public static long
mulDiv(long factorA, long factor1, long divisor) {
final double dd = divisor,
product = (double)factorA * factor1,
a1_d = product / dd;
if (a1_d < -TOO_LARGE || TOO_LARGE < a1_d)
return tooLarge();
if (-ONE_ < a1_d && a1_d < ONE_)
return 0;
if (-EXACT < product && product < EXACT)
return (long) a1_d;
long pLo = factorA * factor1, //diff,
pHi = high64(factorA, factor1);
if (a1_d < -LONG_MAX_ || LONG_MAX_ < a1_d) {
long maxdHi = divisor >> 1;
if (maxdHi < pHi
|| maxdHi == pHi
&& Long.compareUnsigned((divisor << Long.SIZE-1),
pLo) <= 0)
return tooLarge();
}
final double high_dd = TWO_POWER64/dd;
long quotient = (long) a1_d,
loPP = quotient * divisor,
hiPP = high64(quotient, divisor);
long remHi = pHi - hiPP, // xxx overflow/carry
remLo = pLo - loPP;
if (Long.compareUnsigned(pLo, remLo) < 0)
remHi -= 1;
double fudge = remHi * high_dd;
if (remLo < 0)
fudge += high_dd;
fudge += remLo/dd;
long //fHi = (long)fudge/TWO_POWER64,
fLo = (long) Math.floor(fudge); //*round
quotient += fLo;
loPP = quotient * divisor;
hiPP = high64(quotient, divisor);
remHi = pHi - hiPP; // should be 0?!
remLo = pLo - loPP;
if (Long.compareUnsigned(pLo, remLo) < 0)
remHi -= 1;
if (0 == remHi && 0 <= remLo && remLo < divisor)
return quotient;
fudge = remHi * high_dd;
if (remLo < 0)
fudge += high_dd;
fudge += remLo/dd;
fLo = (long) Math.floor(fudge);
return quotient + fLo;
}
/** max <code>double</code> trusted to represent
* a value in the range of <code>long</code> */
static final double
LONG_MAX_ = Double.valueOf(Long.MAX_VALUE - 0xFFF);
/** max <code>double</code> trusted to represent a value below 1 */
static final double
ONE_ = Double.longBitsToDouble(
Double.doubleToRawLongBits(1) - 4);
/** max <code>double</code> trusted to represent a value exactly */
static final double
EXACT = Long.MAX_VALUE >> 12;
static final double
TWO_POWER64 = Double.valueOf(1L<<32)*Double.valueOf(1L<<32);
static long tooLarge() {
// throw new RuntimeException("result too large for long");
return Long.MIN_VALUE;
}
static final long ONES_32 = ~(~0L << 32);
static long high64(long factorA, long factor1) {
long loA = factorA & ONES_32,
hiA = factorA >>> 32,
lo1 = factor1 & ONES_32,
hi1 = factor1 >>> 32;
return ((loA * lo1 >>> 32)
+loA * hi1 + hiA * lo1 >>> 32)
+ hiA * hi1;
}
（我从 IDE 中重新安排了一些代码，使 mulDiv() 位于顶部。
由于懒惰，我有一个用于符号处理的包装器 - 可能会在地狱冻结之前尝试并正确地完成它。
对于计时，迫切需要输入模型：
如何使得每个可能的结果都是同样可能的？）
也许不聪明，但有线性结果时间
#define MUL_DIV_TYPE unsigned int
#define BITS_PER_TYPE (sizeof(MUL_DIV_TYPE)*8)
#define TOP_BIT_TYPE (1<<(BITS_PER_TYPE-1))
//
// result = ( a * b ) / c, without intermediate overflow.
//
MUL_DIV_TYPE mul_div( MUL_DIV_TYPE a, MUL_DIV_TYPE b, MUL_DIV_TYPE c ) {
MUL_DIV_TYPE st, sb; // product sum top and bottom
MUL_DIV_TYPE d, e; // division result
MUL_DIV_TYPE i, // bit counter
j; // overflow check
st = 0;
sb = 0;
d = 0;
e = 0;
for( i = 0; i < BITS_PER_TYPE; i++ ) {
//
// Shift sum left to make space
// for next partial sum
//
st <<= 1;
if( sb & TOP_BIT_TYPE ) st |= 1;
sb <<= 1;
//
// Add a to s if top bit on b
// is set.
//
if( b & TOP_BIT_TYPE ) {
j = sb;
sb += a;
if( sb < j ) st++;
}
//
// Division.
//
d <<= 1;
if( st >= c ) {
d |= 1;
st -= c;
e++;
}
else {
if( e ) e++;
}
//
// Shift b up by one bit.
//
b <<= 1;
}
//
// Roll in missing bits.
//
for( i = e; i < BITS_PER_TYPE; i++ ) {
//
// Shift across product sum
//
st <<= 1;
if( sb & TOP_BIT_TYPE ) st |= 1;
sb <<= 1;
//
// Division, continued.
//
d <<= 1;
if( st >= c ) {
d |= 1;
st -= c;
}
}
return( d ); // remainder should be in st
}