From owner-freebsd-numerics@FreeBSD.ORG Sat Aug 10 03:31:44 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id A50E5619 for ; Sat, 10 Aug 2013 03:31:44 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail104.syd.optusnet.com.au (mail104.syd.optusnet.com.au [211.29.132.246]) by mx1.freebsd.org (Postfix) with ESMTP id 4FCF524B3 for ; Sat, 10 Aug 2013 03:31:40 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail104.syd.optusnet.com.au (Postfix) with ESMTPS id B2FDE421D30; Sat, 10 Aug 2013 13:31:31 +1000 (EST) Date: Sat, 10 Aug 2013 13:31:30 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Zimmermann Paul Subject: Re: help requested with long double issue on 32-bit FreeBSD In-Reply-To: Message-ID: <20130810122319.T1129@besplex.bde.org> References: MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1613949660-1376105490=:1129" X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=Yos2GeoX c=1 sm=1 tr=0 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 a=PO7r1zJSAAAA:8 a=_JvqIC1MqocA:10 a=JzwRw_2MAAAA:8 a=dsHV8SD-vYkA:10 a=nlC_4_pT8q9DhB4Ho9EA:9 a=cz2ZRIgtxKwA:10 a=wJWlkF7cXJYA:10 a=Zeg_oI3sAAAA:8 a=LxBWjLRv-lQ41HUyHYIA:9 a=xxJZeSBEJMgmUUgp:21 a=_nIdGiT0i-FIxw2P:21 a=45ClL6m2LaAA:10 Cc: freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 10 Aug 2013 03:31:44 -0000 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-1613949660-1376105490=:1129 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Fri, 9 Aug 2013, Zimmermann Paul wrote: > [this mail was sent on the MPFR list, but Steve Kargl suggested to post i= t > also here] > > we need help to investigate an error we get with the development version = of > MPFR on i686-freebsd. The error occurs on the hydra.nixos.org continuous > integration platform. Unfortunately we don't have direct access to the > corresponding computer, thus isolating the issue is not easy. > > The error we get can be seen at http://hydra.nixos.org/build/5666171/log/= raw. > Starting from the long double value in 3, and subtracting e=3D0.25, we ge= t the > long double value in 4, which is obviously wrong in the last bits: > > 3: d=3D2.632048527199790448306230411512629530e-01 > 4: d=3D1.320485271997904469509776959057489876e-02 > ^^^^^^^^^^^^^^^^^^^^ > > We believe this is a compiler bug (GCC 4.2.1 was used in this test). > > If you have access to a 32-bit computer under FreeBSD, please could you t= ry > to compile the development version and run the "tset_ld" test, after doin= g > export GMP_CHECK_RANDOMIZE=3D1376689137? > > If you can reproduce the error, please tell Vincent Lef=C3=A8vre and myse= lf > (Vincent.Lefevre at ens-lyon.fr and Paul.Zimmermann at inria.fr); we'll h= elp > you to isolate the problem (e.g., trying different optimization levels). > > If you can't reproduce (on a 32-bit FreeBSD machine), please tell us too; > this might indicate the bug was fixed in a later gcc version. I didn't run the tests, but it was easy to see that this is just the result of FreeBSD's default rounding precision in 32-bit mode being 53 bits and not switching to 64 bits in the tests. Long doubles having 64 mantissa bits obviously cannot work right with a rounding precision of 53, so something must switch to 64 bits in any program or part of a program that uses long doubles and actually needs them. Also, gcc supports the default rounding precision to a fault. It rounds even literal constants to 53 bits at compile time, so it is difficult to even initialize d to the above values. The following programs demonstrates several variations of the problem. % #include %=20 % unsigned char cd[10] =3D { % 0x05, 0x7d, 0x5f, 0x29, 0x55, 0xc9, 0xc2, 0x86, 0xfd, 0x3f, % }; %=20 % main() % { % =09long double d =3D 2.632048527199790448306230411512629530e-01L; % =09long double e =3D 1.320485271997904469509776959057489876e-02L; %=20 % #ifdef FIX_RUNTIME_ROUNDING_PRECISION % =09fpsetprec(FP_PE); % #endif % #ifdef FIX_COMPILE_TIME_ROUNDING_PRECISION % =09d =3D *(long double *)&cd; % #endif % =09printf("%La\n", d); % =09printf("%.37Le\n", d); % =09printf("%La\n", e); % =09printf("%.37Le\n", e); % =09printf("%La\n", d - 0.25); % =09printf("%.37Le\n", d - 0.25); % } This is easiest to test using CC=3Dclang, since clang is not bug-for-bug compatible with gcc4-2.1 and doesn't round even long double constant expressions to 53 bits at compile time. So it gives your 3:d directly. gcc needs the magic FIX_COMPILE_TIME_ROUNDING_PRECISION to give same 3:d. (Even hex constants are rounded to 53 bits by gcc.) Then the result of d - 0.25 is the same as your 4:d in the default rounding precision, but FIX_RUNTIME_ROUNDING_PRECISION fixes this. The bug is also easy to see by printing the values in hex on a system without the bug. On FreeBSD-amd64: % 0x1.0d8592aa52befa0ap-2 % 2.6320485271997904483062304115126295301e-01 % 0x1.b0b2554a57df4p-7 % 1.3204852719979044695097769590574898757e-02 % 0x1.b0b2554a57df414p-7 % 1.3204852719979044830623041151262953008e-02 The intermediate result 4:d is the final result rounded to 53 bits. Note that this was printed by a program compiled by clang. clang presumably rounding 4:d to 64 bits, but it ends up with only 53 nonzero bits since it was produced by a program that rounded it to that many and the decimal format has more than enough precision to recover the 53-bit value when rounded to 64 bits. It is interesting that the decimal digits for d and e line up but the binary digits don't. I prefer a different %a format that lines up the binary digits more often, but this is unfortunately not supported by printf(). This is hard to fix. Programs using long double precision should switch the rounding precision at runtime using fpsetprec() as above. But naive programmers won't know to do this, and apparently even people at INRIA don't know this :-). Newer long double functions in FreeBSD libm do the switch at runtime using macros that hide most of the complications including avoiding switching the rounding precision if it is already 64 bits (switching the precision is very slow, but testing to avoid it when unnecessary only takes large code when it is in fact unnecessary). Compile-time rounding of constant expressions is even harder to fix. FreeBSD libm used to mostly use the method of writing long double constants as the sum of 2 double constants, with a volatile hack to force the evaluation at runtime (in 64-bit precision iff someone remembered to switch the rounding precision). Newer parts mostly use macro that hides some of the details of initializing a union containing the necessary values in bits. Bruce --0-1613949660-1376105490=:1129--