大整数加法代码

Question

我正在尝试使用以下代码在 CUDA 中实现大整数加法

__global__ void add(unsigned *A, unsigned *B, unsigned *C/*output*/, int radix){

    int id = blockIdx.x * blockDim.x + threadIdx.x; 
    A[id ] = A[id] + B[id]; 

    C[id ] =  A[id]/radix;
    __syncthreads();
    A[id] =  A[id]%radix + ((id>0)?C[id -1]:0);
    __syncthreads();
    C[id] = A[id];
}

但它不能正常工作，而且我现在不知道如何处理额外的进位。谢谢

Answer 1

TL;DR 构建一个进位先行加法器，其中每个单独的加法器添加模基数，而不是模 2

添加需要进位

您模型中的问题是您的进位有波动。参见 Rippling carry adders。

如果您使用的是 FPGA，那将不是问题，因为它们有专门的逻辑来快速完成 (carry chains, they're cool)。但是，唉，你在 GPU 上！

也就是说，对于给定的 id，您只知道输入进位（因此您是否要对 A[id]+B[id] 或 A[id]+B[id]+1 求和）当所有总和都较小 id 值已经计算出来。事实上，最初，你只知道第一个进位。

    A[3]+B[3] + ?     A[2]+B[2] + ?     A[1]+B[1] + ?     A[0]+B[0] + 0       
        |                 |                 |                 |
        v                 v                 v                 v
       C[3]              C[2]              C[1]              C[0]

表征进位输出

而且每个和还有进位输出，图中没有。因此，您必须将这个更大方案中的加法视为具有 3 个输入和 2 个输出的函数：(C, c_out) = add(A, B, c_in)

为了不等待 O(n) 完成总和（其中 n 是您的总和被分割成的项目数），您可以预先计算每个 id 处的所有可能结果。这不是很大的工作量，因为 A 和 B 没有改变，只有进位。所以你有 2 个可能的输出：(c_out0, C) = add(A, B, 0) 和 (c_out1, C') = add(A, B, 1).

现在有了所有这些结果，我们基本上需要实现一个 carry lookahead unit。

为此，我们需要计算出每个和的进位输出的函数 P 和 G :

Pa.k.a。以下所有定义
- 传播
- "if a carry comes in, then a carry will go out of this sum"
- c_out1 && !c_out0
- A + B == radix-1
Ga.k.a。以下所有定义
- 生成
- "whatever carry comes in, a carry will go out of this sum"
- c_out1 && c_out0
- c_out0
- A + B >= radix

换句话说，c_out = G or (P and c_in)。所以现在我们有了一个算法的开始，它可以很容易地告诉我们每个 id 的进位输出直接作为其进位输入的函数：

在每个 id，计算 C[id] = A[id]+B[id]+0
获取G[id] = C[id] > radix -1
得到P[id] = C[id] == radix-1

对数树

现在我们可以在 O(log(n)) 内完成，尽管树状的东西在 GPU 上很讨厌，但仍然比等待更短。实际上，从彼此相邻的 2 个加法中，我们可以得到一组 G 和一组 P :

对于id和id+1：

step = 2
if id % step == 0, do steps 6 through 10, otherwise, do nothing
group_P = P[id] and P[id+step/2]
group_G = (P[id+step/2] and G[id]) or G[id+step/2]
c_in[id+step/2] = G[id] or (P[id] and c_in[id])
step = step * 2
if step < n, go to 5

最后（每次用较少 ids 对树的每个级别重复步骤 5-10 之后），所有内容都将用 Ps 和 Gs 是你计算出来的，c_in[0] 是 0。在维基百科页面上有用于按 4 而不是 2 分组的公式，这将使您在 O(log_4(n)) 而不是 O(log_2(n)).

算法结束：

在每个id，得到c_in[id]
return (C[id]+c_in[id]) % radix

利用硬件

我们在这最后一部分真正做的是用逻辑模拟进位先行加法器的电路。然而，我们已经在硬件中有做类似事情的加法器（根据定义）。

让我们用基于 2 的基数来替换我们基于基数的 P 和 G 的定义，就像我们硬件内部的逻辑一样，模拟 2 位的总和 a 和 b 在每个阶段：如果 P = a ^ b （异或）和 G = a & b （逻辑与）。换句话说，a = P or G 和 b = G。因此，如果我们创建一个 intP 整数和一个 intG 整数，其中每个位分别是 P 和 G 我们从每个 ids 总和（限制us 到 64 和），那么加法 (intP | intG) + intG 与我们精心设计的逻辑方案具有完全相同的进位传播。

我猜这些整数的归约仍然是对数运算，但这是可以预料的。

有趣的是，和的每一位都是其进位输入的函数。实际上，和的每一位最终都是 3 位 a+b+c_in % 2 的函数。

如果在那位P == 1，那么a + b == 1，因此a+b+c_in % 2 == !c_in
否则，a+b 是 0 或 2，并且 a+b+c_in % 2 == c_in

因此我们可以简单地形成整数（或者更确切地说是位数组）int_cin = ((P|G)+G) ^ P，其中 ^ 为 xor。

因此我们的算法有一个替代结尾，替换第 4 步及之后的步骤：

在每个 id，将 P 和 G 移动 id ：P = P << id 和 G = G << id
进行 OR 归约得到 intG 和 intP，它们是 [=11= 的所有 P 和 G 的 OR ] 0..63
计算（一次）int_cin = ((P|G)+G) ^ P
在每个 id，得到 `c_in = int_cin & (1 << id) ? 1 : 0;
return(C[id]+c_in) % radix

PS ：另外，如果 radix 很大，请注意数组中的整数溢出。如果不是，那么整个事情就没有意义了，我想...

PPS ：在备用结尾中，如果您有超过 64 个项目，请用 P 和 G 来表征它们，就好像 radix 是 2^64，并在更高层次上重新运行相同的步骤（减少，得到c_in）然后回到较低层次应用7和P+G+carry in from higher level

大整数加法代码

Big integer addition code

c

math

cuda

biginteger

添加需要进位

表征进位输出

对数树

利用硬件