CODEWITHSUNDEEP
floating-pointprecisionmpfr
pmor
pmor

Reputation: 6276

How to compute DECIMAL_DIG from value returned by mpfr_get_prec()?

Context: C standard has xxx_DECIMAL_DIG macros, which can be used in printf to maintain precision of floating-point value. Example for double:

printf("%.*g\n", DBL_DECIMAL_DIG, x);

A simple question: how to compute the similar DECIMAL_DIG from value returned by mpfr_get_prec?

Example:

#include <stdio.h>
#include <mpfr.h>

#define PREC 128

int prec2decimal_dig(int prec)
{
    return prec / 3;    // example, not accurate
}

int main (void)
{
    mpfr_t x;

    mpfr_init2(x, PREC);
    mpfr_set_d(x, 2.0 / 3.0, MPFR_RNDN);
    mpfr_printf("%.*Rf\n", prec2decimal_dig(PREC), x);
    return 0;
}

$ gcc t78a.c -lmpfr -lgmp -oa && ./a
0.666666666666666629659232512494781985878944

What would be the precise implementation of prec2decimal_dig?

Extra: base-2: why DECIMAL_DIG - DIG == 2 or 3?.

Upvotes: 1

Views: 83

Answers (2)

vinc17
vinc17

Reputation: 3466

In practice, the native floating-point numbers are almost always in radix 2 (i.e. FLT_RADIX = 2 and b = 2 in the description of DECIMAL_DIG in the C standard), like MPFR. In this case, you can use mpfr_get_str_ndigits with radix b = 10 (note that b in the C standard is the input radix, i.e. 2, and b in the MPFR manual is the output radix, i.e. 10 in this context). This function is described as:

Return the minimal integer m such that any number of p bits, when output with m digits in radix b with rounding to nearest, can be recovered exactly when read again, still with rounding to nearest. More precisely, we have m = 1 + ceil(p × log(2)/log(b)), with p replaced by p − 1 if b is a power of 2.

Note that if FLT_RADIX = 10 on your platform, DECIMAL_DIG is just p. And I don't know any other platform with another radix.

Upvotes: 2

chux
chux

Reputation: 153456

XXX_DECIMAL_DIG is derived from:

  • 'b', the base. Examples: 2, 10, ...

  • 'p', the number of digits (in base b) encoded, both explicitly or implicitly.

When b is a power of 10:

XXX_DECIMAL_DIG = p * log10(b)

Otherwise

XXX_DECIMAL_DIG = ceiling(1 + (p * log10(b)))

Use 28/93 for a fraction just larger than log102.

Example for common double

DBL_DECIMAL_DIG = (int) ceil(1 + (53 * log10(2));
DBL_DECIMAL_DIG = (int) ceil(1 + (53 * 28.0/93));
DBL_DECIMAL_DIG = 1 + ((53 * 28 + 92)/93);
DBL_DECIMAL_DIG = 17;

#define LOG10_2_NUM 28
#define LOG10_2_DEN 93
int prec2decimal_dig(int prec) {
  return 1 + ((prec*28 + 92)/93);
  // or 
  return (prec*LOG10_2_NUM + LOG10_2_DEN*2 - 1)/LOG10_2_DEN;
}

Upvotes: 1

Related Questions