c - 扩展范围的浮点型

Question

我需要在 C 中表示非常小阶的浮点数（例如 0.6745 × 2 ^-3000）。这种支持必须独立于平台（适用于 CPU 和 GPU-CUDA）。有效数字的大长度不是必需的。

我不能使用高精度库（GMP、MPFR 等），因为它们不能在 GPU 上运行。另一方面，CUDA 不支持该long double类型。有什么解决办法吗？是否有可能以某种方式实现自定义浮点类型？

score 1 · Accepted Answer

您可以在对数空间中工作，即将每个数字表示为 e ^x其中 x 是您的标准浮点类型：

可以使用log-sum-exp 技巧来执行加法/减法（以及更普遍的求和），即
- e ^x +e ^y = e ^x (1+e ^y-x ) = e ^{x + log(1+exp(yx))}
乘法/除法变成加法/减法
- e ^x × e ^x = e ^x+y
提升权力相当简单
- (e ^x )^(e ^x ) = e ^{x exp(y)}

score 0 · Accepted Answer

我写了简单的解决方案（使用这项工作）：

#include <math.h>
#include <stdint.h>

#define DOUBLE_PRECISION 53

/*  DOUBLE PRECISION FLOATING-POINT TYPE WITH EXTENDED EXPONENT     */

typedef struct Real {
    double sig;     //significand
    long exp;       //binary exponent
} real;

/*  UNION FOR DIVISION DOUBLE BY 2^POW                              */

union DubleIntUnion
{
    double      dvalue;
    uint64_t    ivalue;
};


/*  PLACE SIGNIFICAND OF REAL NUMBER IN RANGE [1, 2)            */

inline real adjust(real x){
    real y;
    y.exp = x.exp;
    y.sig = x.sig;
    if(y.sig == 0){
        y.exp = 0;
    } else if (fabs(y.sig) >= 2.0){
        y.exp = y.exp + 1;
        y.sig = y.sig / 2;
    } else if(fabs(y.sig) < 1){
        y.exp = y.exp - 1;
        y.sig = y.sig * 2;
    }
    return y;
}

/*  PLACE SIGNIFICAND OF REAL NUMBER IN RANGE [1, 2) FOR TINY NUMBER    */
/*  FOR EXAMPLE, AFTER SUBTRATION OR WHEN SET REAL FROM DOUBLE          */

inline real adjusttiny(real x){
    real y;
    y.exp = x.exp;
    y.sig = x.sig;
    while(1){
        x.exp = y.exp;
        x.sig = y.sig;
        y = adjust(x);
        if(x.exp == y.exp && x.sig == y.sig)
            break;
    }
    return y;
}

real set(double x){
    real y;
    real z;
    y.sig = x;
    y.exp = 0;
    return adjusttiny(y);
};

real set(real x){
    real y;
    y.exp = x.exp;
    y.sig = x.sig;
    return y;
};

/*  ARITHMETIC OPERATIONS   */

//divide x by 2^pow. Assert that x.exp - pow > e_min
inline double div2pow(const double x, const int pow)
{
    DubleIntUnion diu;
    diu.dvalue = x;
    diu.ivalue -= (uint64_t)pow << 52;      // subtract pow from exponent
    return diu.dvalue;
}

//summation
inline real sum(real x, real y){            
    real sum;
    int dexp = abs(x.exp - y.exp);

    if (x.exp > y.exp){
        sum.exp = x.exp;
        if(dexp <= DOUBLE_PRECISION){           
            sum.sig = div2pow(y.sig, dexp);     // divide y by 2^(x.exp - y.exp)
            sum.sig = sum.sig + x.sig;
        } else sum.sig = x.sig;
    } else if (y.exp > x.exp){
        sum.exp = y.exp;
        if(dexp <= DOUBLE_PRECISION){           
            sum.sig = div2pow(x.sig, dexp);     // divide x by 2^(y.exp - x.exp)
            sum.sig = sum.sig + y.sig;
        } else
            sum.sig = y.sig;
    } else {
        sum.exp = x.exp;
        sum.sig = x.sig + y.sig;
    }
    return adjust(sum);
}

//subtraction
inline real sub(real x, real y){            
    real sub;
    int dexp = abs(x.exp - y.exp);

    if (x.exp > y.exp){
        sub.exp = x.exp;
        if(dexp <= DOUBLE_PRECISION){           
            sub.sig = div2pow(y.sig, dexp); // divide y by 2^(x.exp - y.exp)
            sub.sig = x.sig - sub.sig;
        } else sub.sig = x.sig;
    } else if (y.exp > x.exp){
        sub.exp = y.exp;
        if(dexp <= DOUBLE_PRECISION){           
            sub.sig = div2pow(x.sig, dexp); // divide x by 2^(y.exp - x.exp)
            sub.sig = sub.sig - y.sig; 
        } else sub.sig = -y.sig;
    } else {
        sub.exp = x.exp;
        sub.sig = x.sig - y.sig;
    }
    return adjusttiny(sub);
}

//multiplication
inline real mul(real x, real y){            
    real product;
    product.exp = x.exp + y.exp;
    product.sig = x.sig * y.sig;
    return adjust(product);
}

//division
inline real div(real x, real y){            
    real quotient;
    quotient.exp = x.exp - y.exp;
    quotient.sig = x.sig / y.sig;
    return adjust(quotient);
}

乍一看工作正常。也许我错过了一些东西或者可以加速实施？

我如何实现功能floor和ceil这些数字？

score 0 · Accepted Answer

如果您需要非常大的指数，那么对称级别索引算法可能适合您的需求。然而，准确度更难预测，因此您可能需要更高的 LI（水平指数）值精度来进行补偿。提高精度的一种常见方法是双双运算，这在 CUDA 上也有很多使用。

CUDA 上还有许多多精度库，例如CUMP

更多信息：

c - 扩展范围的浮点型

3 回答 3

Related

Reference