Home My Page Projects MPFR
Summary Activity Forums Tracker Lists Tasks Docs News SCM Files

[#11300] mpfr_set_ld problem with IEEE quad

Date:
2010-10-17 20:07
Priority:
3
State:
Closed
Submitted by:
Jakub Jelinek (jakub)
Assigned to:
Paul Zimmermann (zimmerma)
Category:
general
Resolution:
Fixed
Target Version:
3.1
 
Summary:
mpfr_set_ld problem with IEEE quad

Detailed description
https://bugzilla.redhat.com/show_bug.cgi?id=643505

On s390/s390x (and probably any other arch with IEEE quad long double):


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

int
main (void)
{
mpfr_set_default_prec (__LDBL_MANT_DIG__);
mpfr_set_emin (-__LDBL_MAX_EXP__-__LDBL_MANT_DIG__+4);
mpfr_set_emax (__LDBL_MAX_EXP__);
mpfr_t xx, yy;
long double l = 0x1.23456789abcdef0123456789abcdp-913L;
mpfr_init_set_ld (xx, l, GMP_RNDN);
mpfr_printf ("%Ra\n%La\n%La\n", xx, l, mpfr_get_ld (xx, GMP_RNDN));
mpfr_clear (xx);
l = 0x1.23456789abcdef0123456789abcdp-914L;
mpfr_init_set_ld (yy, l, GMP_RNDN);
mpfr_printf ("%Ra\n%La\n%La\n", yy, l, mpfr_get_ld (yy, GMP_RNDN));
mpfr_clear (yy);
return 0;
}

prints correct number in the first case, but when it is divided by two, it
prints (and similarly retrieves with mpfr_get_ld) complete garbage:

0x9.1a2b3c4d5e6f78091a2b3c4d5e68p-916
0x1.23456789abcdef0123456789abcdp-913
0x1.23456789abcdef0123456789abcdp-913
-0xb.fffffffffffffffffffffffffff8p-1028
0x1.23456789abcdef0123456789abcdp-914
-0x1.7fffffffffffffffffffffffffffp-1025

I think the bug can be quite obviously seen in set_ld.c, the shift_exp -= NNNN
and x /= div1N stuff is only safe if
inexact = mpfr_set_d (u, (double) x, GMP_RNDZ);
hasn't been done in that function, because otherwise all of what has been
accumulated into t will be wrongly multiplied by 2^-1024 etc.
What happens when trying to mpfr_set_ld on
0x1.23456789abcdef0123456789abcdp-914L
is that we have:
1st iteration x = 8.215640181713713163755129115618857e-276
2nd iteration x = -9.9763047253779979586613445625595544e-293
3rd iteration x = -4.1720134847010025932941863449982576e-309
div10 = 5.5626846462680034577255817933310102e-309
As 3rd iteration x is smaller in absolute value than div10, we divide it by
0x1p-1024 (i.e. multiply by 0x1p1024 to get something like 0.75 in x) and set
shift_exp to -1024. But we don't adjust t in any way, so instead of the whole
value being roughly 8.215640181713713163755129115618857e-276 it will be
roughly 8.215640181713713163755129115618857e-276 * 0x1p-1024 plus roughly 0.75
* 0x1p-1024.

Shifting up (untested patch attached to above bugzilla bugreport) might work for IEEE quad, but might cause issues with IBM double double long double format where the low double may be many thousands bits away from the upper double.

Followup

Message
Date: 2010-10-20 16:44
Sender: Jakub Jelinek

On s390x-linux the tester program successfully passed 1 billion of long double fmal tests, so I think this can be closed.
Thanks for fixing it.
Date: 2010-10-19 16:47
Sender: Jakub Jelinek

Thanks, the patch looks good and on a second thought, this shouldn't be a problem for double double format, because even if the upper double is close to DBL_MAX and the smaller double close to __DBL_DENORM_MIN__, in double double div11 is zero.

As for testing, it will take some time before I get access to a s390x box again, hopefully by tomorrow evening I'll be able to report back.
Date: 2010-10-19 15:01
Sender: Vincent Lefèvre

Paul,
long double is quad on Sparc Solaris (not x86) and on HP-UX.
Date: 2010-10-19 14:28
Sender: Paul Zimmermann

Dear Jakub,
I have (hopefully) fixed that bug in revision 7222.
It is hard to check for us since we have no access to a
machine with quad double, moreover on our machines the
special code for little endian double extended format is
used, thus we do not exercise the "generic" code which had the
bug. By manually forcing the generic code, I was able to
construct another example which fails with little endian
double extended, and my patch makes it work, thus I hope it
will work for you too.

Once you confirm it works, I will close this bug.

Paul

PS: by the way, what does "configure" say on the following line on this processor? Does it recognize the quad format?

checking format of `long double' floating point... IEEE extended, little endian

Attached Files:

Changes:

Field Old Value Date By
Target Version2.42011-05-02 14:46vlefevre
close_date2011-05-02 14:462011-05-02 14:46vlefevre
ResolutionAccepted2010-10-20 20:45zimmerma
close_date2010-10-20 20:452010-10-20 20:45zimmerma
status_idOpen2010-10-20 20:45zimmerma
close_dateNone2010-10-19 15:13zimmerma
status_idASSIGNED2010-10-19 15:13zimmerma
close_date2010-10-19 15:012010-10-19 15:01vlefevre
ResolutionNone2010-10-19 14:28zimmerma
close_date2010-10-19 14:282010-10-19 14:28zimmerma
assigned_tonone2010-10-19 14:28zimmerma
status_idOpen2010-10-19 14:28zimmerma