java - (a * b) / c MulDiv 和处理中间乘法的溢出
问题描述
我需要做以下算术:
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
自己实现整个。我想要做的是比使用泛型更快地解决特定问题(a*b/c
在哪里a*b
可以溢出)BigInteger
。
Edit2:如果它可以以一种聪明的方式完成,那将是理想的,根本不会溢出,评论中出现了一些提示,但我仍在寻找一个正确的。
更新: 我尝试将 BigInteger 代码移植到我的特定需求,而不创建对象,在第一次迭代中,与使用 BigInteger(在我的开发 PC 上)相比,我的速度提高了约 46%。
Long.MIN_VALUE
然后我尝试了一些修改过的@David Eisenstat 解决方案,与 BigInteger 相比,它给了我约 56%(我从to运行 100_000_000_000 个随机输入Long.MAX_VALUE
)减少了运行时间(超过 2 倍)(与我改编的 BigInteger 算法相比约为 18%) .
在优化和测试上会有更多的迭代,但在这一点上,我认为我必须接受这个答案作为最好的答案。
解决方案
我一直在修补一种方法,(1)乘以21 位肢体上的学校算法(2)继续进行长除法a
,残差的不寻常表示使用 a来存储高阶位和 a存储低位。我不知道它是否可以与标准长除法竞争,但为了您的享受,b
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;
}
推荐阅读
- testing - ETL 测试中的关键错误是什么,我们如何找到它?
- java - 如何在 Java 中删除 JSON 多余的行?
- swift - 缩小全屏图像的大小
- r - 如何修复二项式 GLM 模型中的截距?
- python - Scrapy SitemapSpider 从站点地图获取更多信息
- vba - Word VBA:找到红色字体并给出章节名称
- microsoft-translator - 通过 python 代码进行 Microsoft 自定义翻译培训
- reactjs - 访问上一个状态
- twig - What is the purpose of using element['#object'] in twig templates in drupal 8. How does this work?
- ms-access - 参数表的最佳设计是什么?