0

gmp変数のマシン精度を取得しようとしています。

そのために、固定精度で gmp の精度を計算するために、wikipediaのコードを適合させました。

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 標準) も試しましたが、結果は正しいです...

4

1 に答える 1

3

ウィキペディアの式と一致する結果が得られ、値は精度によって異なります。

ただし、値 (および実効精度) は、四肢境界(1)を超える場合にのみ変化します。私にとって、これは 64 の倍数を意味します。

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

計算されたマシン イプシロンは

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 プログラムの出力はmantissa,'e',exponent、値が次の形式になっています。

0.mantissa * 10^exponent

(1) GMP は、浮動小数点数を指数 (基数 2) と仮数 (および符号) のペアとして表します。仮数は、符号なし整数の配列である limbs として維持されます。私にとって、四肢は64ビットで、32ビットシステムでは通常32ビット（iirc）です。そのため、目的の精度が(k-1)*LIMB_BITS(排他的) と (排他的)の間にある場合k*LIMB_BITS、仮数の配列にはリムが含まれk、それらすべてが使用されるため、有効な精度はk*LIMB_BITSビットになります。したがって、四肢の数が変化した場合にのみイプシロンが変化します。

于 2012-06-12T11:57:10.197 に答える