用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