Вычислить точность машины на произвольной точности gmp

Question

Я пытаюсь получить точность машины на 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 standard), и результат правильный ...

Answer 1

Я получаю результаты, соответствующие wikipedia's formula, и значения зависят от точности.

Однако значение - и эффективная точность - изменяются только при пересечении границы конечности ⁽¹⁾. Для меня это означает, что кратные 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 

⁽¹⁾ GMP представляет собой число с плавающей точкой в ​​виде пары экспоненты (для основания 2) и мантиссы (и знак). Мантисса поддерживается как массив целых чисел без знака, конечности. Для меня конечности 64 бит, на 32-битных системах обычно 32 бита (iirc). Поэтому, когда желаемая точность находится между (k-1)*LIMB_BITS (исключая) и k*LIMB_BITS (включительно), массив для мантиссы содержит k конечности, и все они используются, поэтому эффективная точность равна k*LIMB_BITS бит. Поэтому эпсилон изменяется только при изменении количества конечностей.