[Bug c/30013] New: Multiple flaws in decimal floating-point arithmetic conversions fixed

["beebe at math dot utah dot edu" <gcc-bugzilla at gcc dot gnu dot org>]

Testing of my decimal floating-point arithmetic code with gcc-4.2 and
gcc-4.3 releases turned up several bugs in gcc's type conversions that
I've had to track down and fix.

They repair several problems:

(1) Loss of trailing digits by failure to provide a suitable precision
    in conversion formats; I use the values from the formula in David
    Matula's CACM 19(3) 716--723 June 1968 paper, expressed by this
    code snippet:

        ### matula(nbits): return number of decimal digits needed to ensure
        ### correct round-trip conversion between binary and decimal
        func matula(nbits) { return ceil(nbits/log2(10) + 1) }

(2) Failure to provide prototypes of library functions that are not
    necessarily supplied by system header files: if the gcc developers
    would routinely use the gcc -Werror-implicit-function-declaration
    option in builds, this problem could have, and should have, been
    caught before release.

(3) Incorrect data types for conversion functions.

(4) Failure to match sprintf() argument types with format items.

(5) Serious loss of precision for _Decimal128 by using double strtod()
    instead of long double strtold().

The test file given below can be used to check the conversions before
and after the patch.

Here are diffs suitable for input to patch:

103 airy> diff gcc/config/dfp-bit.c.org gcc/config/dfp-bit.c
516c516
<   sprintf (buf, BFP_FMT, x);
---
>   sprintf (buf, BFP_FMT, (double)x);

104 airy> diff gcc/config/dfp-bit.h.org gcc/config/dfp-bit.h
241c241
< #define BFP_FMT "%e"
---
> #define BFP_FMT "%.9e"
246c246
< #define BFP_FMT "%e"
---
> #define BFP_FMT "%.17e"
253,255c253,257
< #define BFP_FMT "%e"
< #define BFP_VIA_TYPE double
< #define STR_TO_BFP strtod
---
> #define BFP_FMT "%.36Le"
> #define BFP_VIA_TYPE long double
> /* strtold is declared in <stdlib.h> only for C99.  */
> extern long double strtold (const char *, char **);
> #define STR_TO_BFP strtold
422c424
< typedef float DFtype __attribute__ ((mode (DF)));
---
> typedef double DFtype __attribute__ ((mode (DF)));
424c426
< typedef float XFtype __attribute__ ((mode (XF)));
---
> typedef long double XFtype __attribute__ ((mode (XF)));

Here is the test file:

% cat dec-cast-bug.c
/***********************************************************************
[23-Nov-2006]
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void
test(long double x)
{
    float f;
    double d;
    long double y;
    _Decimal32 d32;
    _Decimal64 d64;
    _Decimal128 d128;

    d32 = (_Decimal32)x;
    (void)printf(" (_Decimal32)(%.21Lg) = 0x%08x\n", x,
                 ((unsigned int *)(&d32))[0]);

    d64 = (_Decimal64)x;
    (void)printf(" (_Decimal64)(%.21Lg) = 0x%08x_%08x\n", x,
                 ((unsigned int *)(&d64))[1],
                 ((unsigned int *)(&d64))[0]);

    d128 = (_Decimal128)x;
    (void)printf("(_Decimal128)(%.21Lg) = 0x%08x_%08x_%08x_%08x\n", x,
                 ((unsigned int *)(&d128))[3],
                 ((unsigned int *)(&d128))[2],
                 ((unsigned int *)(&d128))[1],
                 ((unsigned int *)(&d128))[0]);

    (void)printf("\n");

    f = (float)x;
    (void)printf("      (float)(%.21Lg) = %.9g             = %a\n", x,
(double)f, (double)f);

    d = (double)x;
    (void)printf("     (double)(%.21Lg) = %.17g     = %a\n", x, (double)d,
(double)d);

    y = (long double)x;
    (void)printf("(long double)(%.21Lg) = %.21Lg = %La\n", x, (long double)y,
(long double)y);

    (void)printf("\n");

    y = (long double)(float)x;
    (void)printf("      (long double)(float)(%.21Lg) = %.21Lg = %La\n", x, y,
y);

    y = (long double)(double)x;
    (void)printf("     (long double)(double)(%.21Lg) = %.21Lg = %La\n", x, y,
y);

    y = (long double)(long double)x;
    (void)printf("(long double)(long double)(%.21Lg) = %.21Lg = %La\n", x, y,
y);

    (void)printf("\n");

    f = (float)(_Decimal32)x;
    (void)printf(" (float)(_Decimal32)(%.21Lg) = %.9g = %a\n", x, (double)f,
(double)f);

    f = (float)(_Decimal64)x;
    (void)printf(" (float)(_Decimal64)(%.21Lg) = %.9g = %a\n", x, (double)f,
(double)f);

    f = (float)(_Decimal128)x;
    (void)printf("(float)(_Decimal128)(%.21Lg) = %.9g = %a\n", x, (double)f,
(double)f);

    (void)printf("\n");

    d = (double)(_Decimal32)x;
    (void)printf(" (double)(_Decimal32)(%.21Lg) = %.17g = %a\n", x, (double)d,
(double)d);

    d = (double)(_Decimal64)x;
    (void)printf(" (double)(_Decimal64)(%.21Lg) = %.17g = %a\n", x, (double)d,
(double)d);

    d = (double)(_Decimal128)x;
    (void)printf("(double)(_Decimal128)(%.21Lg) = %.17g = %a\n", x, (double)d,
(double)d);

    (void)printf("\n");

    y = (long double)(_Decimal32)x;
    (void)printf(" (long double)(_Decimal32)(%.21Lg) = %.21Lg = %La\n", x,
(long double)y, (long double)y);

    y = (long double)(_Decimal64)x;
    (void)printf(" (long double)(_Decimal64)(%.21Lg) = %.21Lg = %La\n", x,
(long double)y, (long double)y);

    y = (long double)(_Decimal128)x;
    (void)printf("(long double)(_Decimal128)(%.21Lg) = %.21Lg = %La\n", x,
(long double)y, (long double)y);

    (void)printf("\n");

    y = (long double)(float)(_Decimal32)x;
    (void)printf(" (long double)(float)(_Decimal32)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    y = (long double)(float)(_Decimal64)x;
    (void)printf(" (long double)(float)(_Decimal64)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    y = (long double)(float)(_Decimal128)x;
    (void)printf("(long double)(float)(_Decimal128)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    (void)printf("\n");

    y = (long double)(double)(_Decimal32)x;
    (void)printf(" (long double)(double)(_Decimal32)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    y = (long double)(double)(_Decimal64)x;
    (void)printf(" (long double)(double)(_Decimal64)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    y = (long double)(double)(_Decimal128)x;
    (void)printf("(long double)(double)(_Decimal128)(%.21Lg) = %.21Lg = %La\n",
x, y, y);

    (void)printf("\n");

    y = (long double)(_Decimal32)x;
    (void)printf(" (long double)(_Decimal32)(%.21Lg) = %.21Lg = %La\n", x, y,
y);

    y = (long double)(_Decimal64)x;
    (void)printf(" (long double)(_Decimal64)(%.21Lg) = %.21Lg = %La\n", x, y,
y);

    y = (long double)(_Decimal128)x;
    (void)printf("(long double)(_Decimal128)(%.21Lg) = %.21Lg = %La\n", x, y,
y);
}

int
main(int argc, char* argv[])
{
    test(3.141592653589793238462643383279502884L);

    return (EXIT_SUCCESS);
}


-- 
           Summary: Multiple flaws in decimal floating-point arithmetic
                    conversions fixed
           Product: gcc
           Version: 4.3.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: c
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: beebe at math dot utah dot edu
 GCC build triplet: x86_64-unknown-linux-gnu
  GCC host triplet: x86_64-unknown-linux-gnu
GCC target triplet: x86_64-unknown-linux-gnu


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=30013