logo标签列表

CUDA下的大整数加法

ccudagpgputhrust
11

我正在GPU上开发一种加密算法，目前卡在实现大整数加法的算法上。通常，大整数是按照一串32位字来表示的。

举个例子，我们可以使用一个线程来相加两个32位字。为了简单起见，假设要相加的数字具有相同的长度，并且每个块中的线程数等于字数。那么：

__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扩展，但是大整数包似乎还没有实现。也许有人可以给我一个提示如何在CUDA上实现？

2
GPU可以直接处理长达64位（long long）的数据。关于128位数据的处理方法在这个SO问题/答案中有详细说明。 - Robert Crovella
4
加法在算术上不够复杂，不能有效地跨线程进行拆分（至少我暂时想不到）。对于乘法，你可以让每个线程计算一个列的部分 32x32 -> 64 位积的总和，然后在最后传播进位。你还可以通过将加法的结果作为单独的和向量和进位向量来延迟进位传播。很多事情都取决于确切的算法背景。 - njuffa
是的，我同意加法不像乘法那样计算密集。但是，如果我有一个长度为2048甚至更长的整数，在使用一个CUDA线程循环执行加法时将非常低效，因为这个操作对我来说看起来非常尴尬并行，除了进位传播。 - user1545642
我已经拼凑出一个CUDA内核，可以并行添加多达1024个64位无符号量，并且还能够并行处理这些问题的批次。从内核计算的角度来看，假设我们并行处理大量问题，它比我的天真CPU代码快约10倍。当你考虑数据复制时间时，它与CPU时间大致相同。没有PTX，只有C代码，所以我相信它可以运行得更快，但我不知道能提高多少速度。如果您想查看它，我会将其发布为答案。此外，我对其正确性不作任何声明。 - Robert Crovella
以下是CUDA Thrust扩展的链接，提供了大整数乘法的代码：http://cuda-thrust-extensions.googlecode.com/svn/trunk/big%20integer/ - user2023370
显示剩余3条评论
2个回答
8
你说得对，通过前缀和计算可以进行进位传递，但是定义二元函数以及证明其可结合性（并行前缀和所需）需要一些技巧。实际上，这个算法在理论上被用于 进位预测加法器。
假设我们有两个大整数 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);

具有相当直观的含义：生成器[i] == 1表示在位置i处生成进位，而传播[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

如果在位置i生成进位（carry），则 carryout[i] == 1。或者它有时会在较早的位置生成并传播到位置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 块：

// 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位字中累加加法结果，然后在“一次扫描”中执行一次进位传播来完成所有运算。
我试图编译这个，但是在这一行上出现了一个错误： volatile short *scan = (volatile short *)(r + 16 + thid_in_warp + (49 * (thid / 64)) + ((thid / 32) & 1); 似乎缺少一个闭括号。我尝试在分号之前添加一个闭括号。你可以检查一下吗？（修复后，在尝试使用它来添加两个2048x32位无符号整数时，我遇到了启动失败。我可能有自己代码中的错误。） - Robert Crovella
哎呀，你说得对，应该是(r + 16 + thid_in_warp + 49 * (thid / 64)) + ((thid / 32) & 1)。我已经修复了。请确保为算法分配足够的共享内存以使其正确运行，大约需要(49 * 512 / 64) + 32个字。 - user1545642
还有一件事：你是在编译64位吗？我刚意识到负偏移可能会非常麻烦，因为它们默认情况下不是符号扩展的，即在这里：postscan [((thid >> 5) - 1)]。 - user1545642
你的代码修改没问题，现在我可以成功编译并运行了，不会崩溃了。是的，我正在编译64位版本。然而，我的一个测试案例有一个大整数(2048x32bit)，所有位都设置为1 (即0xF...FFF)，另一个大整数= 0x1，在无符号整数偏移量为128的位置上失败了，CPU 的结果是零，而GPU的结果是0xFFFFFFFF。前128个无符号整数的结果确实与 CPU 的结果匹配（全为零）。我认为这可能是你的代码中一个warp分裂成另一个warp的地方。当然，这可能仍然是我的问题。 - Robert Crovella
谢谢您检查代码。这意味着一个warp未被传播。明天我会在我的机器上检查这个例子。另一个预防措施是，使用UADDO/UADDC宏可能并不总是安全的，因为编译器可以重新排序指令并弄乱进位标志。 - user1545642
显示剩余2条评论
4
我想将我的答案与@asm一起发布，以便这个SO问题可以成为一个思路库。与@asm类似，我也检测并存储进位条件以及“carry-through”条件，即当中间字结果全部为1（0xF...FFF）时，如果进位传播到该字，则会“carry-through”到下一个字。
我没有在代码中使用任何PTX或asm，所以我选择使用64位无符号整数而不是32位，以实现2048x32位的能力，使用1024个线程。
与@asm的代码相比，我并行进位传播方案有更大的差异。我构建了一个位打包数组（“carry”），其中每个位表示来自1024个线程中每个独立中间64位加法生成的进位条件。我还构建了一个位打包数组（“carry_through”），其中每个位表示单个64位中间结果的carry_through条件。对于1024个线程，这相当于每个位打包数组16x64位共享内存，因此总共使用的共享内存量为64 + 3个32位数量。使用这些位打包数组，我执行以下操作以生成组合传播进位指示器：
carry = carry | (carry_through ^ ((carry & carry_through) + carry_through);

（请注意，进位向左移动了一位：carry[i]表示a[i-1]+b[i-1]的结果产生了一个进位）
以下是详细说明：
  1. 将carry和carry_through进行按位与运算，生成可能出现进位的地方，这些地方会与一个或多个进位条件序列相互作用
  2. 将步骤1的结果加上carry_through，生成一个有变化的结果，它代表着所有受到进位影响的单词
  3. 将carry_through和步骤2的结果进行异或操作，显示受影响的结果，其中包含1位
  4. 将步骤3的结果和普通的进位指示器进行按位或运算，得到组合进位条件，然后用它来更新所有中间结果。
请注意，步骤2中的加法需要进行另一种多词加法（对于由超过64个单词组成的大整数）。 我相信这种算法是有效的，并且已经通过了我投入的测试用例。
以下是我实现这个算法的示例代码：
// parallel add of large integers
// requires CC 2.0 or higher
// compile with:
// nvcc -O3 -arch=sm_20 -o paradd2 paradd2.cu
#include <stdio.h>
#include <stdlib.h>

#define MAXSIZE 1024 // the number of 64 bit quantities that can be added
#define LLBITS 64  // the number of bits in a long long
#define BSIZE ((MAXSIZE + LLBITS -1)/LLBITS) // MAXSIZE when packed into bits
#define nTPB MAXSIZE

// define either GPU or GPUCOPY, not both -- for timing
#define GPU
//#define GPUCOPY

#define LOOPCNT 1000

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

// perform c = a + b, for unsigned integers of psize*64 bits.
// all work done in a single threadblock.
// multiple threadblocks are handling multiple separate addition problems
// least significant word is at a[0], etc.

__global__ void paradd(const unsigned size, const unsigned psize, unsigned long long *c, const unsigned long long *a, const unsigned long long *b){

  __shared__ unsigned long long carry_through[BSIZE];
  __shared__ unsigned long long carry[BSIZE+1];
  __shared__ volatile unsigned mcarry;
  __shared__ volatile unsigned mcarry_through;

  unsigned idx = threadIdx.x + (psize * blockIdx.x);
  if ((threadIdx.x < psize) && (idx < size)){
    // handle 64 bit unsigned add first
    unsigned long long cr1 = a[idx];
    unsigned long long lc = cr1 + b[idx];
    // handle carry
    if (threadIdx.x < BSIZE){
      carry[threadIdx.x] = 0;
      carry_through[threadIdx.x] = 0;
      }
    if (threadIdx.x == 0){
      mcarry = 0;
      mcarry_through = 0;
      }
    __syncthreads();
    if (lc < cr1){
      if ((threadIdx.x%LLBITS) != (LLBITS-1))  
        atomicAdd(&(carry[threadIdx.x/LLBITS]), (2ull<<(threadIdx.x%LLBITS)));
      else atomicAdd(&(carry[(threadIdx.x/LLBITS)+1]), 1);
      }
    // handle carry-through
    if (lc == 0xFFFFFFFFFFFFFFFFull) 
      atomicAdd(&(carry_through[threadIdx.x/LLBITS]), (1ull<<(threadIdx.x%LLBITS))); 
    __syncthreads();
    if (threadIdx.x < ((psize + LLBITS-1)/LLBITS)){
      // only 1 warp executing within this if statement
      unsigned long long cr3 = carry_through[threadIdx.x];
      cr1 = carry[threadIdx.x] & cr3;
      // start of sub-add
      unsigned long long cr2 = cr3 + cr1;
      if (cr2 < cr1) atomicAdd((unsigned *)&mcarry, (2u<<(threadIdx.x)));
      if (cr2 == 0xFFFFFFFFFFFFFFFFull) atomicAdd((unsigned *)&mcarry_through, (1u<<threadIdx.x));
      if (threadIdx.x == 0) {
        unsigned cr4 = mcarry & mcarry_through;
        cr4 += mcarry_through;
        mcarry |= (mcarry_through ^ cr4); 
        }
      if (mcarry & (1u<<threadIdx.x)) cr2++;
      // end of sub-add
      carry[threadIdx.x] |= (cr2 ^ cr3);
      }
    __syncthreads();
    if (carry[threadIdx.x/LLBITS] & (1ull<<(threadIdx.x%LLBITS))) lc++;
    c[idx] = lc;
  }
}

int main() {

  unsigned long long *h_a, *h_b, *h_c, *d_a, *d_b, *d_c, *c;
  unsigned at_once = 256;   // valid range = 1 .. 65535
  unsigned prob_size = MAXSIZE ; // valid range = 1 .. MAXSIZE
  unsigned dsize = at_once * prob_size;
  cudaEvent_t t_start_gpu, t_start_cpu, t_end_gpu, t_end_cpu;
  float et_gpu, et_cpu, tot_gpu, tot_cpu;
  tot_gpu = 0;
  tot_cpu = 0;


  if (sizeof(unsigned long long) != (LLBITS/8)) {printf("Word Size Error\n"); return 1;}
  if ((c = (unsigned long long *)malloc(dsize * sizeof(unsigned long long)))  == 0) {printf("Malloc Fail\n"); return 1;}

  cudaHostAlloc((void **)&h_a, dsize * sizeof(unsigned long long), cudaHostAllocDefault);
  cudaCheckErrors("cudaHostAlloc1 fail");
  cudaHostAlloc((void **)&h_b, dsize * sizeof(unsigned long long), cudaHostAllocDefault);
  cudaCheckErrors("cudaHostAlloc2 fail");
  cudaHostAlloc((void **)&h_c, dsize * sizeof(unsigned long long), cudaHostAllocDefault);
  cudaCheckErrors("cudaHostAlloc3 fail");

  cudaMalloc((void **)&d_a, dsize * sizeof(unsigned long long));
  cudaCheckErrors("cudaMalloc1 fail");
  cudaMalloc((void **)&d_b, dsize * sizeof(unsigned long long));
  cudaCheckErrors("cudaMalloc2 fail");
  cudaMalloc((void **)&d_c, dsize * sizeof(unsigned long long));
  cudaCheckErrors("cudaMalloc3 fail");
  cudaMemset(d_c, 0, dsize*sizeof(unsigned long long));

  cudaEventCreate(&t_start_gpu);
  cudaEventCreate(&t_end_gpu);
  cudaEventCreate(&t_start_cpu);
  cudaEventCreate(&t_end_cpu);

  for (unsigned loops = 0; loops <LOOPCNT; loops++){
  //create some test cases
  if (loops == 0){
  for (int j=0; j<at_once; j++)
  for (int k=0; k<prob_size; k++){
    int i= (j*prob_size) + k;
    h_a[i] = 0xFFFFFFFFFFFFFFFFull;
    h_b[i] = 0;
    }
    h_a[prob_size-1] = 0;
    h_b[prob_size-1] = 1;
    h_b[0] = 1;
  }
  else if (loops == 1){
  for (int i=0; i<dsize; i++){
    h_a[i] = 0xFFFFFFFFFFFFFFFFull;
    h_b[i] = 0;
    }
    h_b[0] = 1;
  }
  else if (loops == 2){
  for (int i=0; i<dsize; i++){
    h_a[i] = 0xFFFFFFFFFFFFFFFEull;
    h_b[i] = 2;
    }
    h_b[0] = 1;
  }
  else {
  for (int i = 0; i<dsize; i++){
    h_a[i] = (((unsigned long long)lrand48())<<33) + (unsigned long long)lrand48();
    h_b[i] = (((unsigned long long)lrand48())<<33) + (unsigned long long)lrand48();
    }
  }
#ifdef GPUCOPY
  cudaEventRecord(t_start_gpu, 0);
#endif
  cudaMemcpy(d_a, h_a, dsize*sizeof(unsigned long long), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudaMemcpy1 fail");
  cudaMemcpy(d_b, h_b, dsize*sizeof(unsigned long long), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudaMemcpy2 fail");
#ifdef GPU
  cudaEventRecord(t_start_gpu, 0);
#endif
  paradd<<<at_once, nTPB>>>(dsize, prob_size, d_c, d_a, d_b);
  cudaCheckErrors("Kernel Fail");
#ifdef GPU
  cudaEventRecord(t_end_gpu, 0);
#endif
  cudaMemcpy(h_c, d_c, dsize*sizeof(unsigned long long), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudaMemcpy3 fail");
#ifdef GPUCOPY
  cudaEventRecord(t_end_gpu, 0);
#endif
  cudaEventSynchronize(t_end_gpu);
  cudaEventElapsedTime(&et_gpu, t_start_gpu, t_end_gpu);
  tot_gpu += et_gpu;
  cudaEventRecord(t_start_cpu, 0);
  //also compute result on CPU for comparison
  for (int j=0; j<at_once; j++) {
    unsigned rc=0;
    for (int n=0; n<prob_size; n++){
      unsigned i = (j*prob_size) + n;
      c[i] = h_a[i] + h_b[i];
      if (c[i] < h_a[i]) {
        c[i] += rc;
        rc=1;}
      else {
        if ((c[i] += rc) != 0) rc=0;
        }
      if (c[i] != h_c[i]) {printf("Results mismatch at offset %d, GPU = 0x%lX, CPU = 0x%lX\n", i, h_c[i], c[i]); return 1;}
      }
    }
  cudaEventRecord(t_end_cpu, 0);
  cudaEventSynchronize(t_end_cpu);
  cudaEventElapsedTime(&et_cpu, t_start_cpu, t_end_cpu);
  tot_cpu += et_cpu;
  if ((loops%(LOOPCNT/10)) == 0) printf("*\n");
  }
  printf("\nResults Match!\n");
  printf("Average GPU time = %fms\n", (tot_gpu/LOOPCNT));
  printf("Average CPU time = %fms\n", (tot_cpu/LOOPCNT));

  return 0;
}
实际上，我相信我的进位传递可以进一步减少为：carry = carry | (carry_through ^ (carry + carry_through)); - Robert Crovella
这非常有用。您能否给出关于您的计算机(CPU、GPU、操作系统等)CPU和GPU平均使用时间的数据？ - user3891236
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接