将一个 int64 拆分为两个 int32，执行数学运算，然后重新加入

Question

我在具有 64 位整数限制的硬件限制下工作。不支持浮点数。我正在处理需要乘除的非常大的整数。乘法时，我遇到了 64 位溢出。我正在用python制作一个解决方案的原型。这就是我的功能：

upper = x >> 32 #x is cast as int64 before being passed to this function
lower = x & 0x00000000FFFFFFFF

temp_upper = upper * y // z #Dividing first is not an option, as this is not the actual equation I am working with. This is just to make sure in my testing I overflow unless I do the splitting.
temp_lower = lower * y // z

return temp_upper << 32 | lower

这在一定程度上可行，但我最终会损失很多精度（我的结果有时会相差几百万）。从表面上看，这似乎是因为分裂而发生的。如果足够的话，它会将上部向右移动。然后当我将它移回原位时，我有一个零间隙。

不幸的是，这个主题很难用谷歌搜索，因为任何带有上/下的东西都会带来关于向上/向下舍入的结果。任何关于拆分 int 的操作都会返回将它们拆分为 char 数组的结果。任何关于 int 算术的东西都会用整数数学提出基本代数。也许我只是不擅长谷歌搜索。但是你们能给我一些关于如何做到这一点的指示吗？

像这样拆分只是我正在尝试的事情，它不一定是解决方案。我需要做的就是暂时超过 64 位整数限制。最终结果将在 64bit 以下（除法部分之后）。我记得在大学里学习如何将其拆分，然后进行数学计算并重新组合。但不幸的是，正如我所说，我无法在网上找到有关如何对其进行实际数学运算的任何内容。

最后，我的人数有时很少。所以我不能砍掉正确的部分。我需要结果基本上等同于我使用 int128 之类的东西。

我想另一种看待这个问题的方式是这样。由于我对拆分 int64 没有任何问题，我们可以忘记那部分。所以我们可以假装两个 int64 被喂给我，一个在上，一个在下。我不能将它们组合起来，因为它们不适合单个 int64。所以我需要先将它们除以 Z。组合步骤很容易。如何进行除法？

谢谢。

Answer 1

据我了解，您想执行(x*y)//z。

您的号码 x,y,z 都适合 64 位，除了中间 x*y 需要 128 位。

你的问题确实和分裂有关：你有

h * y = qh * z + rh
l * y = ql * z + rl

h * y << 32 + l*y = (qh<<32 + ql) * z + (rh<<32 + rl)

但没有说(rh<<32 + rl) < z，在你的情况下l*y 的高位与h * y 的低位重叠，所以你得到了错误的商，可能有很多单位。

你应该做的第二个操作是：

rh<<32 + l * y = ql' * z + rl'

然后得到总商qh<<32 + ql'

当然，在计算左操作数时必须注意避免溢出...

由于您只拆分 x*y 的一个操作数，我假设中间结果总是适合 96 位。

如果这是正确的，那么您的问题是将 3 个 32 位肢体 x*y 除以 2 个 32 位肢体 z。
因此，它就像 Burnigel - Ziegler 分治算法一样。
算法可以这样分解：

• 以唐叶为例，得到乘法x*y的3个肢体a2,a1,a0
• 将z分成2个分支z1,z0
• 执行div32( (a2,a1,a0) , (z1,z0) )

这里有一些伪代码，只处理正操作数，不保证正确，但你对实现有个概念：

p = 1<<32;

function (a1,a0) = split(a)
    a1 = a >> 32;
    a0 = a - (a1 * p);

function (a2,a1,a0) = mul22(x,y)
    (x1,x0) = split(x) ;
    (y1,y0) = split(y) ;
    (h1,h0) = split(x1 * y1);
    assert(h1 == 0); -- assume that results fits on 96 bits
    (l1,l0) = split(x0 * y0);
    (m1,m0) = split((x1 - x0) * (y0 - y1));  -- karatsuba trick
    a0 = l0;
    (carry,a1) = split( l1 + l0 + h0 + m0 );
    a2 = l1 + m1 + h0 + carry;
    
function (q,r) = quorem(a,b)
    q = a // b;
    r = a - (b * q);

function (q1,q0,r0) = div21(a1,a0,b0)
   (q1,r1) = quorem(a1,b0);
   (q0,r0) = quorem( r1 * p + a0 , b0 );
   (q1,q0) = split( q1 * p + q0 );
   
function q = div32(a2,a1,a0,b1,b0)
    (q,r) = quorem(a2*p+a1,b1*p+b0);
    q = q * p;
    (a2,a1)=split(r);
    if a2<b1
        (q1,q0,r)=div21(a2,a1,b1);
        assert(q1==0); -- since a2<b1...
    else
        q0=p-1;
        r=(a2-b1)*p+a1+b1;
    (d1,d0) = split(q0*b0);
    r = (r-d1)*p + a0 - d0;
    while(r < 0)
        q = q - 1;
        r = r + b1*p + b0;
   
function t=muldiv(x,y,z)
    (a2,a1,a0) = mul22(x,y);
    (z1,z0) = split(z);
    if z1 == 0
        (q2,q1,r1)=div21(a2,a1,z0);
        assert(q2==0); -- otherwise result will not fit on 64 bits
        t = q1*p + ( ( r1*p + a0 )//z0);
    else
        t = div32(a2,a1,a0,z1,z0);