用2个float模拟double

在有些设备上只有float没有double，比如前几代GPU、部分移动设备等。当非得用到double精度的时候该怎么办？

我记得去年在某个地方见到过用2个float模拟double的作法，经过一番玩命地搜索，得来全不费功夫，就在CUDA SDK的Mandelbrot例子里找到了2个float模拟double乘法的函数。甚至，GTX280上的double也是类似的方法模拟出来的，所 以慢的惊人，只有float八分之一的速度。

先show一下模拟乘法的函数dsmul：

//  This function multiplies DS numbers A and B to yield the DS product C.
__device__ inline  void  dsmul( float   & c0,  float   & c1,
     const   float  a0,  const   float  a1,  const   float  b0,  const   float  b1)
{
     //  This splits dsa(1) and dsb(1) into high-order and low-order words.
     float  cona  =  a0  *   8193.0f ;
     float  conb  =  b0  *   8193.0f ;
     float  sa1  =  cona  -  (cona  -  a0);
     float  sb1  =  conb  -  (conb  -  b0);
     float  sa2  =  a0  -  sa1;
     float  sb2  =  b0  -  sb1;

     //  Multilply a0 * b0 using Dekker's method.
     float  c11  =  a0  *  b0;
     float  c21  =  (((sa1  *  sb1  -  c11)  +  sa1  *  sb2)  +  sa2  *  sb1)  +  sa2  *  sb2;

     //  Compute a0 * b1 + a1 * b0 (only high-order word is needed).
     float  c2  =  a0  *  b1  +  a1  *  b0;

     //  Compute (c11, c21) + c2 using Knuth's trick, also adding low-order product.
     float  t1  =  c11  +  c2;
     float  e  =  t1  -  c11;
     float  t2  =  ((c2  -  e)  +  (c11  -  (t1  -  e)))  +  c21  +  a1  *  b1;

     //  The result is t1 + t2, after normalization.
    c0  =  e  =  t1  +  t2;
    c1  =  t2  -  (e  -  t1);
}  //  dsmul

在NV的论坛上还找到了高人写的一系列运算函数。虽然是CUDA的，但要改成别的语言和环境也是轻而易举的事情：dsmath.h