c - CUDA 大整数加法

Question

我一直在 GPU 上开发加密算法，目前坚持使用一种算法来执行大整数加法。大整数通常表示为一堆 32 位字。

例如，我们可以用一个线程来添加两个32位的字。为简单起见，假设 要添加的数字具有相同的长度和每个 block 的线程数 == 字数。然后:

__global__ void add_kernel(int *C, const int *A, const int *B) {
     int x = A[threadIdx.x];
     int y = B[threadIdx.x];
     int z = x + y;
     int carry = (z < x);
     /** do carry propagation in parallel somehow ? */
     ............

     z = z + newcarry; // update the resulting words after carry propagation
     C[threadIdx.x] = z;
 }

我很确定有一种方法可以通过一些棘手的减少程序来进行进位传播，但我无法弄清楚..

我看过CUDA thrust extensions但是大整数包似乎还没有实现。 也许有人可以给我提示如何在 CUDA 上做到这一点？

Answer 1

你是对的，进位传播可以通过前缀和计算来完成，但是为这个操作定义二元函数并证明它是关联的(并行前缀和需要)有点棘手。事实上，该算法(理论上)用于 Carry-lookahead adder .

假设我们有两个大整数 a[0..n-1] 和 b[0..n-1]。 然后我们计算 (i = 0..n-1):

s[i] = a[i] + b[i]l;
carryin[i] = (s[i] < a[i]);

我们定义了两个函数:

generate[i] = carryin[i];
propagate[i] = (s[i] == 0xffffffff);

具有相当直观的含义:generate[i] == 1 表示进位产生于 位置 i while propagate[i] == 1 意味着进位将从位置传播 (i - 1) 到 (i + 1)。我们的目标是计算用于更新结果和 s[0..n-1] 的函数 carryout[0..n-1]。 carryout 可以递归计算如下:

carryout[i] = generate[i] OR (propagate[i] AND carryout[i-1])
carryout[0] = 0

这里 carryout[i] == 1 如果进位是在位置 i 生成的，或者它有时更早生成并传播到位置 i。最后，我们更新结果总和:

s[i] = s[i] + carryout[i-1];  for i = 1..n-1
carry = carryout[n-1];

现在可以非常简单地证明进位函数确实是二元关联的，因此可以应用并行前缀和计算。为了在 CUDA 上实现这一点，我们可以将标志“生成”和“传播”合并到一个变量中，因为它们是互斥的，即:

cy[i] = (s[i] == -1u ? -1u : 0) | carryin[i];

换句话说，

cy[i] = 0xffffffff  if propagate[i]
cy[i] = 1           if generate[i]
cy[u] = 0           otherwise

然后，可以验证以下公式计算进位函数的前缀和:

cy[i] = max((int)cy[i], (int)cy[k]) & cy[i];

对于所有 k < i。下面的示例代码显示了 2048 字整数的大加法。这里我使用了 512 个线程的 CUDA block :

// add & output carry flag
#define UADDO(c, a, b) \ 
     asm volatile("add.cc.u32 %0, %1, %2;" : "=r"(c) : "r"(a) , "r"(b));
// add with carry & output carry flag
#define UADDC(c, a, b) \ 
     asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(c) : "r"(a) , "r"(b));

#define WS 32

__global__ void bignum_add(unsigned *g_R, const unsigned *g_A,const unsigned *g_B) {

extern __shared__ unsigned shared[];
unsigned *r = shared; 

const unsigned N_THIDS = 512;
unsigned thid = threadIdx.x, thid_in_warp = thid & WS-1;
unsigned ofs, cf;

uint4 a = ((const uint4 *)g_A)[thid],
      b = ((const uint4 *)g_B)[thid];

UADDO(a.x, a.x, b.x) // adding 128-bit chunks with carry flag
UADDC(a.y, a.y, b.y)
UADDC(a.z, a.z, b.z)
UADDC(a.w, a.w, b.w)
UADDC(cf, 0, 0) // save carry-out

// memory consumption: 49 * N_THIDS / 64
// use "alternating" data layout for each pair of warps
volatile short *scan = (volatile short *)(r + 16 + thid_in_warp +
        49 * (thid / 64)) + ((thid / 32) & 1);

scan[-32] = -1; // put identity element
if(a.x == -1u && a.x == a.y && a.x == a.z && a.x == a.w)
    // this indicates that carry will propagate through the number
    cf = -1u;

// "Hillis-and-Steele-style" reduction 
scan[0] = cf;
cf = max((int)cf, (int)scan[-2]) & cf;
scan[0] = cf;
cf = max((int)cf, (int)scan[-4]) & cf;
scan[0] = cf;
cf = max((int)cf, (int)scan[-8]) & cf;
scan[0] = cf;
cf = max((int)cf, (int)scan[-16]) & cf;
scan[0] = cf;
cf = max((int)cf, (int)scan[-32]) & cf;
scan[0] = cf;

int *postscan = (int *)r + 16 + 49 * (N_THIDS / 64);
if(thid_in_warp == WS - 1) // scan leading carry-outs once again
    postscan[thid >> 5] = cf;

__syncthreads();

if(thid < N_THIDS / 32) {
    volatile int *t = (volatile int *)postscan + thid;
    t[-8] = -1; // load identity symbol
    cf = t[0];
    cf = max((int)cf, (int)t[-1]) & cf;
    t[0] = cf;
    cf = max((int)cf, (int)t[-2]) & cf;
    t[0] = cf;
    cf = max((int)cf, (int)t[-4]) & cf;
    t[0] = cf;
}
__syncthreads();

cf = scan[0];
int ps = postscan[(int)((thid >> 5) - 1)]; // postscan[-1] equals to -1
scan[0] = max((int)cf, ps) & cf; // update carry flags within warps
cf = scan[-2];

if(thid_in_warp == 0)
    cf = ps;
if((int)cf < 0)
    cf = 0;

UADDO(a.x, a.x, cf) // propagate carry flag if needed
UADDC(a.y, a.y, 0)
UADDC(a.z, a.z, 0)
UADDC(a.w, a.w, 0)
((uint4 *)g_R)[thid] = a;
}

请注意，宏 UADDO/UADDC 可能不再是必需的，因为 CUDA 4.0 具有相应的内在函数(但我不完全确定)。

另请注意，虽然并行缩减非常快，但如果您需要连续添加几个大整数，最好使用一些冗余表示(上面的评论中建议)，即首先累加结果在 64 位字中进行加法运算，然后在“一次扫描”中在最后执行一次进位传播。