用两个uint表示一个双精度，然后乘以？

Question

我正在使用CUDA (GPGPU编程)进行一些研究，由于新的硬件架构，与生俱来的双精度性能比单精度性能(高出24倍！)。我决定尝试使用两个uint来表示一个双精度。这样，我可以运行DP计算，而不会对性能造成重大影响。

例如，假设我们想要使用此方法表示double 10.12。

uint real = 10，decimal = 12；

因此，'real.decimal‘在视觉上代表我们的替身。

让我们称这个类型为“single”。

给定:单个a= 10.12，b= 20.24；

乘、除、加、减两个单数的有效算法是什么？

单个c=a* b；

请记住，要使其正常工作，不能进行任何DP计算或使用大于32位的数据类型。

Answer 1

如果您想替换双精度，您可以考虑使用"double single“表示，其中操作数是称为"head”和"tail“的单精度数字对，并且满足|tail| <= 0.5 * ulp (|head|)的规范化要求。CUDA的float2是保存所需浮点数对的合适类型。例如，x分量表示“尾”，y分量表示“头”。

请注意，这并不能提供与双精度完全相同的精度。该表示的精度被限制为49位(来自每个单精度尾数的24位加上尾部的符号位1位)与双精度的53位。这些操作不会提供IEEE舍入，而且由于乘法代码中临时量的溢出，范围也可能会有所减小。不正常的数量也可能不能准确地表示。

我不认为性能会比使用GPU的本地双精度操作好很多，因为寄存器压力会更高，并且每个“双单”操作所需的单精度操作的数量相当多。当然，您可以尝试测量这两种变体的性能。

下面是"double single“格式的加法示例。注释中注明了该算法的源代码。我相信引用的工作是基于D.Priest的工作，他在20世纪80年代末或90年代初对这个主题进行了研究。有关“双一”乘法的示例，可以查看CUDA附带的文件math_functions.h中的函数__internal_dsmul()。

关于CUDA中的64位整数操作的快速注释(正如一位评论者指出的那样，这是一个潜在的替代方案)。当前的GPU硬件对本机64位整数操作的支持非常有限，基本上提供了浮点类型之间的转换，以及在加载和存储中使用64位地址的索引。另外，64位整数操作是通过本机32位操作实现的，通过查看反汇编代码(cuobjdump --dump-sass)可以很容易地看出这一点。

/* Based on: Andrew Thall, Extended-Precision Floating-Point Numbers for GPU 
   Computation. Retrieved from http://andrewthall.org/papers/df64_qf128.pdf
   on 7/12/2011.
*/
__device__ float2 dsadd (float2 a, float2 b)
{
    float2 z;
    float t1, t2, t3, t4, t5, e;

    t1 = __fadd_rn (a.y, b.y);
    t2 = __fadd_rn (t1, -a.y);
    t3 = __fadd_rn (__fadd_rn (a.y, t2 - t1), __fadd_rn (b.y, -t2));
    t4 = __fadd_rn (a.x, b.x);
    t2 = __fadd_rn (t4, -a.x);
    t5 = __fadd_rn (__fadd_rn (a.x, t2 - t4), __fadd_rn (b.x, -t2));
    t3 = __fadd_rn (t3, t4);
    t4 = __fadd_rn (t1, t3);
    t3 = __fadd_rn (t1 - t4, t3);
    t3 = __fadd_rn (t3, t5);
    z.y = e = __fadd_rn (t4, t3);
    z.x = __fadd_rn (t4 - e, t3);
    return z;
}

Answer 2

@njuffa发布的答案几乎可以肯定是实现这一愿望的最聪明的方式。我不知道它会有多快，但由于它利用了设备上最高吞吐量的引擎(单精度浮点引擎)，它似乎是所有竞争对手中最好的。尽管如此，你确实要求一个定点表示，所以我想我应该继续并发布这篇文章，只是将这个问题作为一个想法的仓库。

这段代码有很多注意事项：

减法并不比大多数有符号算术(你提到的uint)更慢，所以alternatives

• doesn't is questionable

• didn't tested

• probably

2. no overflow/underflow checking
3. ignored请求绝对不使用任何64位表示或操作。当然，我没有使用过任何64位(DP)浮点运算，但我确实使用了64位整数运算。我的理由是(在某些情况下)我将固有的64位表示拆分为两个32位数量，执行一堆32位操作，然后重新组装是毫无意义的，因为如果我请求64位(整数)操作，机器(编译器)无论如何都会这样做。
4. 未优化。你也许可以用PTX做得更好，但它可能属于“擦亮大便”这一类。

因此，我将其作为一个框架提供给您，作为娱乐，也许您可以看到，与使用SP引擎相比，使用“真正的”定点格式的速度要慢得多。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 10000
#define nTPB 512
#define SCALE 1
#define MUL_ERR_LIM (0.0000001 * SCALE * SCALE)
#define ADD_ERR_LIM (0.0000001 * SCALE)

#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)


// fixed point number representation: 32 bit whole portion and 32 bit decimal
typedef union {
  struct {
    unsigned d;   // decimal portion
    unsigned w;   // whole number portion
  } part;
  unsigned long val;
} fxp;


__device__ __host__ inline fxp fxp_add(fxp a, fxp b){
  fxp temp;
  temp.val = a.val + b.val;
  return temp; 
}

__device__ __host__ inline fxp fxp_sub(fxp a, fxp b){
  fxp temp;
  temp.val = a.val - b.val;
  return temp;
}

__device__ __host__ inline unsigned fxp_whole(fxp a){
  return a.part.w;
}

__device__ __host__ inline unsigned fxp_decimal(fxp a){
  return a.part.d;
}

__device__ __host__ inline void fxp_set_whole(fxp *a, unsigned val){
  a->part.w = val;
}

__device__ __host__ inline void fxp_set_decimal(fxp *a, unsigned val){
  a->part.d = val;
}

__device__ __host__ inline fxp float_to_fxp(float val){
  fxp temp;
  temp.part.w = (unsigned) truncf(val);
  temp.part.d = (unsigned) rintf((val - truncf(val)) * 0x0000000100000000ul);
  return temp;
}

__device__ __host__ inline float fxp_to_float(fxp val){
  return val.part.w + (val.part.d/(float)0x0000000100000000ul);
}

__device__ __host__ inline fxp double_to_fxp(double val){
  fxp temp;
  temp.part.w = (unsigned) trunc(val);
  temp.part.d = (unsigned) rint((val - trunc(val)) * 0x0000000100000000ul);
  return temp;
}

__device__ __host__ inline double fxp_to_double(fxp val){
  return val.part.w + (val.part.d/(double)0x0000000100000000ul);
}

__device__ __host__ fxp fxp_mul(fxp a, fxp b){
  fxp temp;
  unsigned long ltemp = ((unsigned long)a.part.w * (unsigned long)b.part.d) + ((unsigned long)a.part.d * (unsigned long)b.part.w);
  unsigned utemp = (unsigned) (ltemp & 0x0FFFFFFFFul);
  temp.part.w = (unsigned) (ltemp >> 32);
  temp.part.d = (unsigned) (((unsigned long)a.part.d * (unsigned long)b.part.d) >> 32) + utemp;
  temp.part.w += (a.part.w * b.part.w) + ((temp.part.d < utemp) ? 1:0);
  return temp;
}

__global__ void fxp_test(float *a, float *b, float *s, float *p, double *da, double *db, double *ds, double *dp, unsigned n){

  unsigned idx = threadIdx.x + (blockDim.x * blockIdx.x);
  if (idx < n){
   float la = a[idx];
   float lb = b[idx];
   double lda = da[idx];
   double ldb = db[idx];
   s[idx] = fxp_to_float(fxp_add(float_to_fxp(la), float_to_fxp(lb)));
   p[idx] = fxp_to_float(fxp_mul(float_to_fxp(la), float_to_fxp(lb)));
   ds[idx] = fxp_to_double(fxp_add(double_to_fxp(lda), double_to_fxp(ldb)));
   dp[idx] = fxp_to_double(fxp_mul(double_to_fxp(lda), double_to_fxp(ldb)));
   }
}



int main(){

  fxp a,b,c;
  float x,y,z, xp, yp, zp;
  double dx, dy, dz, dxp, dyp, dzp;
  if (sizeof(unsigned) != 4) {printf("unsigned type size error: %d\n", sizeof(unsigned)); return 1;}
  if (sizeof(unsigned long) != 8) {printf("unsigned long type error: %d\n", sizeof(unsigned long)); return 1;}
// test host side
  x = 76.705116;
  y = 1.891480;

  a = float_to_fxp(x);
  b = float_to_fxp(y);
// test conversions
  xp = fxp_to_float(a);
  yp = fxp_to_float(b);
  printf("con: a = %f, a should be %f, b = %f, b should be %f\n", xp, x, yp, y);
// test multiply
  c = fxp_mul(a, b);
  z = x*y;
  zp = fxp_to_float(c);
  printf("mul: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);
//test add
  c = fxp_add(a, b);
  z = x+y;
  zp = fxp_to_float(c);
  printf("add: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);
//test subtract
  c = fxp_sub(a, b);
  z = x-y;
  zp = fxp_to_float(c);
  printf("sub: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);

// now test doubles
  dx = 6.7877;
  dy = 5.2444;

  a = double_to_fxp(dx);
  b = double_to_fxp(dy);
// test conversions
  dxp = fxp_to_double(a);
  dyp = fxp_to_double(b);
  printf("dbl con: a = %f, a should be %f, b = %f, b should be %f\n", dxp, dx, dyp, dy);
// test multiply
  c = fxp_mul(a, b);
  dz = dx*dy;
  dzp = fxp_to_double(c);
  printf("double mul: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);
//test add
  c = fxp_add(a, b);
  dz = dx+dy;
  dzp = fxp_to_double(c);
  printf("double add: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);
//test subtract
  c = fxp_sub(a, b);
  dz = dx-dy;
  dzp = fxp_to_double(c);
  printf("double sub: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);

// test device side
  float *h_a, *d_a, *h_b, *d_b, *h_s, *d_s, *h_p, *d_p;
  double *h_da, *d_da, *h_db, *d_db, *h_ds, *d_ds, *h_dp, *d_dp;

  if ((h_a=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_b=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_s=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_p=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_da=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_db=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_ds=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
  if ((h_dp=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}

  cudaMalloc((void **)&d_a, N*sizeof(float));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_b, N*sizeof(float));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_s, N*sizeof(float));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_p, N*sizeof(float));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_da, N*sizeof(double));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_db, N*sizeof(double));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_ds, N*sizeof(double));
  cudaCheckErrors("cudamalloc fail");
  cudaMalloc((void **)&d_dp, N*sizeof(double));
  cudaCheckErrors("cudamalloc fail");

  for (unsigned i = 0; i < N; i++){
    h_a[i] = (float)drand48() * SCALE;
    h_b[i] = (float)drand48() * SCALE;
    h_da[i] = drand48()*SCALE;
    h_db[i] = drand48()*SCALE;
    }

  cudaMemcpy(d_a, h_a, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(d_b, h_b, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(d_da, h_da, N*sizeof(double), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(d_db, h_db, N*sizeof(double), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudamemcpy fail");

  fxp_test<<<(N+nTPB-1)/nTPB, nTPB>>>(d_a, d_b, d_s, d_p, d_da, d_db, d_ds, d_dp, N);
  cudaCheckErrors("kernel fail");


  cudaMemcpy(h_s, d_s, N*sizeof(float), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(h_p, d_p, N*sizeof(float), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(h_ds, d_ds, N*sizeof(double), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudamemcpy fail");
  cudaMemcpy(h_dp, d_dp, N*sizeof(double), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudamemcpy fail");

  for (unsigned i=0; i<N; i++){
    if (fabsf(h_s[i] - (h_a[i] + h_b[i])) > ADD_ERR_LIM) {printf("float add mismatch at %d, a: %f b: %f gpu: %f  cpu: %f\n", i, h_a[i], h_b[i],  h_s[i], (h_a[i] + h_b[i])); return 1;}
    if (fabsf(h_p[i] - (h_a[i] * h_b[i])) > MUL_ERR_LIM) {printf("float mul mismatch at %d, a: %f b: %f gpu: %f  cpu: %f\n", i, h_a[i], h_b[i],  h_p[i], (h_a[i] * h_b[i])); return 1;}
    if (fabs(h_ds[i] - (h_da[i] + h_db[i])) > ADD_ERR_LIM) {printf("double add mismatch at %d, a: %f b: %f gpu: %f  cpu: %f\n", i, h_da[i], h_db[i], h_ds[i], (h_da[i] + h_db[i])); return 1;}
    if (fabs(h_dp[i] - (h_da[i] * h_db[i])) > MUL_ERR_LIM) {printf("double mul mismatch at %d, a: %f b: %f gpu: %f  cpu: %f\n", i, h_da[i], h_db[i], h_dp[i], (h_da[i] * h_db[i])); return 1;}
    }
  printf("GPU results match!\n");
  return 0;
} 

另外，分子动力学研究代码AMBER运行快速的on GPUs，并实现了几种不同的算术模型来实现其速度，包括SP和DP计算as well as an SPFP format的混合，其中一些工作使用定点表示，而SP用于其他工作。然而，我不知道它的细节。但显然，它可以取得良好的效果。