对于定点数学库,我需要一个快速的96位64位特定除法算法

Question

我目前正在编写一个快速的32.32定点数学库.我成功地使加法,减法和乘法工作正确,但我完全坚持分裂.

对那些不记得的人提醒一下:32.32定点数是一个具有32位整数部分和32位小数部分的数字.

我想出的最好的算法需要96位整数除法,这是编译器通常没有内置函数的东西.

无论如何,这里是:

G = 2^32

notation: x is the 64-bit fixed-point number, x1 is its low nibble and x2 is its high

G*(a/b) = ((a1 + a2*G) / (b1 + b2*G))*G      // Decompose this

G*(a/b) = (a1*G) / (b1*G + b2) + (a2*G*G) / (b1*G + b2)

如您所见,(a2*G*G)保证大于常规的64位整数.如果我的编译器实际上支持uint128_t,我只需执行以下操作:

((uint128_t)x << 32) / y)

好吧他们不是,我需要一个解决方案.谢谢您的帮助.

Answer 1

您可以将较大的除法分解为多个用较少位进行除法的块.正如另一张已经提到过的海报,该算法可以在Knuth的TAOCP中找到.

但是,没有必要买这本书!

黑客高兴网站上有一个代码用C实现了算法.它是用32位算术编写的64位无符号除法,因此你无法直接删除代码.要从64位到128位,你必须将所有类型,掩码和常数加宽两个,例如,short变为int,0xffff变为0xffffffffll等.

在这个简单易行的改变后,您应该可以进行128位分割.

代码在这里:http://www.hackersdelight.org/HDcode/divlu.c(由于行尾,可能会在Web浏览器中严重包装.如果是这样,只需保存代码并用记事本打开它).

由于您的最大值只需要96位,因此64位分区中的一个将始终返回零,因此您甚至可以简化代码.

哦 - 在我忘记之前:代码只适用于无符号值.要从有符号转换为无符号除法,您可以执行以下操作(伪代码样式):

fixpoint Divide (fixpoint a, fixpoint b)
{
  // check if the integers are of different sign:
  fixpoint sign_difference = a ^ b; 

  // do unsigned division:
  fixpoint x = unsigned_divide (abs(a), abs(b));

  // if the signs have been different: negate the result.
  if (sign_difference < 0)
  {
     x = -x;
  }

  return x;
}

该网站本身也值得一试:http : //www.hackersdelight.org/

希望能帮助到你.

顺便说一下 - 你正在做的好任务..你介意告诉我们你需要什么定点库吗？

顺便说一下 - 除法的普通移位和减法算法也可以.

如果您的目标是x86,则可以使用MMX或SSE内在函数实现它.该算法仅依赖于原始操作,因此它的执行速度也非常快.

Answer 2

更好的自我调整答案：
请原谅答案的 C# 主义，但以下内容应该适用于所有情况。可能有一个解决方案可以更快地找到正确的班次，但我必须比现在思考得更深入。但这应该是相当有效的：

int upshift = 32;
ulong mask = 0xFFFFFFFF00000000;
ulong mod = x % y;
while ((mod & mask) != 0)
{
     // Current upshift of the remainder would overflow... so adjust
     y >>= 1;
     mask <<= 1;
     upshift--;

     mod = x % y;
}
ulong div = ((x / y) << upshift) + (mod << upshift) / y;

简单但不安全的答案：如果余数在高 32 位中设置了任何位，则
此计算可能会导致余数向上移位时溢出x % y，从而导致不正确的答案。

((x / y) << 32) + ((x % y) << 32) / y

第一部分使用整数除法并给出答案的高位（将它们向上移）。

第二部分从高位除法的剩余部分（无法再除的位）计算低位，向上移位然后除法。