根据 gmp 任意精度计算机器精度

Question

我正在尝试获取 gmp 上的机器精度变量。

为此，我改编了 wikipedia 中的代码，以固定精度计算 gmp 的精度:

int main( int argc, char **argv )
{
    long int precision = 100;

mpf_set_default_prec(precision); // in principle redundant, but who cares

    mpf_t machEps, one, temp; // "one" is the "1" in gmp. tmp is to comparison.
    mpf_init2(machEps, precision); 
    mpf_set_ui(machEps, 1); // eps = 1
    mpf_init_set(one,machEps); // ensure "one" has the same precision as machEps
    mpf_init_set(temp,machEps); // ensure "temp" has the same precision as machEps

    do {
        mpf_div_ui(machEps,machEps,2); // eps = eps/2
        mpf_div_ui(temp,machEps,2);    // temp = eps/2
        mpf_add(temp,temp,one);        // temp += 1
    }
    while ( mpf_cmp(temp,one)); // temp == 1
    /// print the result...
    char *t = new char[400];
    mp_exp_t expprt;
    mpf_get_str(NULL, &expprt, 10, 10, machEps);

    sprintf(t, "%se%ld", mpf_get_str(NULL, &expprt, 10, mpf_get_default_prec(), machEps), expprt);

    printf( "Calculated Machine epsilon: %s\n", t);
    return 0;
}

但是结果与维基百科的公式并不一致，也不随我设置的精度而变化。我缺少什么？我也尝试过 double 和 float (c 标准)，结果是正确的......

Answer 1

我得到的结果与 wikipedia's formula 一致，值取决于精度。

但是，该值和有效精度仅在跨越肢体边界时发生变化⁽¹⁾。对我来说，这意味着 64 的倍数，所以对于

(k-1)*64 < precision <= k*64

计算出的机器 epsilon 为

0.5^(k*64)

一些结果:

$ ./a.out 192
Calculated Machine epsilon: 15930919111324522770288803977677118055911045551926187860739e-57
$ ./a.out 193
Calculated Machine epsilon: 8636168555094444625386351862800399571116000364436281385023703470168591803162427e-77

比较:

Prelude> 0.5^192
1.5930919111324523e-58
Prelude> 0.5^256
8.636168555094445e-78

GMP 程序的输出采用尾数,'e',指数 形式，其中值为

0.mantissa * 10^exponent

⁽¹⁾ GMP 将 float 表示为一对指数(基数为 2)和尾数(和符号)。尾数被维护为无符号整数数组，即四肢。对我来说，肢体是 64 位的，在 32 位系统上通常是 32 位 (iirc)。因此，当所需精度介于 (k-1)*LIMB_BITS(不包括)和 k*LIMB_BITS(包括)之间时，尾数数组包含 k 肢体，并且全部使用，因此有效精度为 k*LIMB_BITS 位。因此，epsilon 仅当肢体数量发生变化时才会发生变化。