用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

posted on 2011-01-19 13:47  龚敏敏  阅读(1825)  评论(0编辑  收藏  举报