（a * b）/ c MulDiv并处理中间乘法的溢出

小编典典

（a * b）/ c MulDiv并处理中间乘法的溢出

algorithm

我需要执行以下算法：

long a,b,c;
long result = a*b/c;

虽然保证结果适合long，但乘法不正确，因此可能会溢出。

我尝试逐步进行操作（先相乘然后相除），同时通过将中间结果拆分a*b为最大4个大小的int数组来处理溢出（就像BigInteger使用其int[] mag变量一样）。

在这里，我陷入了分裂。我无法理解进行精确除法所需的按位移位。我需要的只是商（不需要余数）。

假设的方法是：

public static long divide(int[] dividend, long divisor)

另外，我不考虑使用，BigInteger因为这部分代码需要很快（我想坚持使用基本体和基本体数组）。

任何帮助将非常感激！

编辑：我不是要BigInteger自己实现整个。我正在尝试做的是比使用generic更快地解决特定问题（a*b/c在哪里a*b可能溢出）BigInteger。

Edit2：如果能够以一种 巧妙的 方式完成它，那将是理想的，因为根本不会溢出，评论中浮现了一些技巧，但是我仍在寻找正确的技巧。

更新：
我尝试将BigInteger代码移植到我的特定需求，而不创建对象，并且在第一次迭代中，与使用BigInteger（在我的开发PC上）相比，我的速度提高了约46％。

然后我尝试了一下修改@大卫Eisenstat的解决方案，这给了我〜56％（我跑从100_000_000_000随机输入Long.MIN_VALUE到Long.MAX_VALUE）减少运行时间（超过2X）比较的BigInteger（即〜18％相比，我的适应BigInteger的算法中）
。

关于优化和测试的迭代将会更多，但是在这一点上，我认为我必须接受最佳答案。

阅读 327

2020-07-28

共1个答案

小编典典

我一直在尝试一种方法（1）进行乘法运算，a并b使用21位四肢上的school算法（2）进行长除法运算，并使用a
来存储高阶位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;
}

2020-07-28