用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
__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