NetBSD Problem Report #36770

From martin@duskware.de  Fri Aug 10 21:40:50 2007
Return-Path: <martin@duskware.de>
Received: from mail.netbsd.org (mail.netbsd.org [204.152.190.11])
	by narn.NetBSD.org (Postfix) with ESMTP id CFA4063B8EA
	for <gnats-bugs@gnats.netbsd.org>; Fri, 10 Aug 2007 21:40:49 +0000 (UTC)
Message-Id: <20070810203346.809BB63B8EA@narn.NetBSD.org>
Date: Fri, 10 Aug 2007 20:33:46 +0000 (UTC)
From: M.Drochner@fz-juelich.de
Reply-To: M.Drochner@fz-juelich.de
To: netbsd-bugs-owner@NetBSD.org
Subject: "long double" arithmetics doesn't work on i386
X-Send-Pr-Version: www-1.0

>Number:         36770
>Category:       port-i386
>Synopsis:       "long double" arithmetics doesn't work on i386
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    port-i386-maintainer
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Fri Aug 10 21:45:00 +0000 2007
>Last-Modified:  Thu Aug 16 21:10:02 +0000 2007
>Originator:     Matthias Drochner
>Release:        current/i386
>Organization:
>Environment:
gcc version 4.1.3 20070620 prerelease (NetBSD nb1 20070620)
>Description:
This seems to be a toolchain issue.
Appearently, gcc fails to emit instuctions which set the i387
floating point control word to 64-bit rounding where "long double"
data types are used. This makes that 53-bit rounding is used which
is the default on NetBSD (for good reasons afaiu).
So one doesn't get additional precision from the use of "long double".

>How-To-Repeat:
Try the program below (taken from an sqlite3 selftest).
"xxx" implements a poor man's atof(); while it doesn't do proper
rounding the use of "long double" internally should make that
the "double" result is correct anyway (as is the result of atof(),
for comparision).
Setting the FPCW, as done in the "#ifdef FIX" part, makes that
the result is correct. This should be done by gcc automatically.


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

double
xxx(const char *s)
{
        long double r = 0.0;

#ifdef FIX
        const int fpcw = 0x37f;
        asm("fldcw %0"::"m"(fpcw));
#endif

        while (isdigit(*s)) {
                r = r * 10.0 + (*s - '0');
                s++;
        }
        return r;
}

main()
{
        const char in[] = "9223372036854774800";
        double f = xxx(in);
        double g = atof(in);
        printf("%s %.20g %.20g\n", in, f, g);
}

>Fix:
For now I can only tell that setting TARGET_96_ROUND_53_LONG_DOUBLE
in the gcc configuration doesn't help.

>Audit-Trail:
From: Matthias Drochner <M.Drochner@fz-juelich.de>
To: gnats-bugs@NetBSD.org
Cc: port-i386-maintainer@NetBSD.org, gnats-admin@NetBSD.org,
	netbsd-bugs@NetBSD.org
Subject: Re: port-i386/36770: "long double" arithmetics doesn't work on i386

Date: Mon, 13 Aug 2007 15:48:59 +0200

 M.Drochner@fz-juelich.de said:
 > gcc fails to emit instuctions which set the i387 floating point
 > control word to 64-bit rounding where "long double" data types are
 > used

 After some research I came to the conclusion that gcc doesn't
 intend to do so. (It might kill performance anyway.)
 It does instead live with the fact that the FPU does only
 provide IEEE "double" accuracy, even with "long double"
 variables. Since the standards (ISO C and IEEE754) do not
 require that "long double" is any better that "double", this
 is legal.

 > setting TARGET_96_ROUND_53_LONG_DOUBLE
 > in the gcc configuration doesn't help

 What this does is to dumb down compile time calculations
 on "long double" types to use 53-bit accuracy (while using
 the 80/96-bit format where data are put to .data).
 So this should be set for NetBSD/i386 (and amd64) for
 consistency between compile-time and runtime calculations,
 otherwise we might see different behaviour depending on
 eg optimization levels.

 > Try the program below (taken from an sqlite3 selftest).

 sqlite3's assumption that use of "long double" will fix
 the flaws of the conversion algorithm is not portable.
 I've filed a bug report there (#2570).

 This further means for NetBSD/i386:
 -Use of "long double" is just a waste of memory, and memory
  bandwidth.
 +It is almost as easy to implement expl(), logl() and all the other
  "long double" variants of math functions on i386 as it is for
  platforms where "long double" is identical to "double": just
  use wrappers which convert arguments, no need to implement
  algorithms which provide higher accuracies.

 best regards
 Matthias


 Forschungszentrum Juelich GmbH
 52425 Juelich

 Sitz der Gesellschaft: Juelich
 Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
 Vorsitzende des Aufsichtsrats: MinDirig'in Baerbel Brumme-Bothe
 Vorstand: Prof. Dr. Achim Bachem (Vorsitzender), Dr. Ulrich Krafft (stellv. 
 Vorsitzender)

From: Matthias Drochner <M.Drochner@fz-juelich.de>
To: M.Drochner@fz-juelich.de
Cc: gnats-bugs@NetBSD.org, port-i386-maintainer@NetBSD.org,
	gnats-admin@NetBSD.org, netbsd-bugs@NetBSD.org
Subject: Re: port-i386/36770: "long double" arithmetics doesn't work on i386

Date: Thu, 16 Aug 2007 23:06:47 +0200

 This is a multipart MIME message.

 --==_Exmh_24243496863520
 Content-Type: text/plain; charset=us-ascii


 Some more findings:

 M.Drochner@fz-juelich.de said:
 > -Use of "long double" is just a waste of memory, and memory
 >  bandwidth. 

 Almost, but not quite: while the exponent is effectively
 limited to 53 bits, the mantissa has still 15 bits which
 is 4 more than double. So there is no gain in accuracy
 but a little more dynamic range.

 > +It is almost as easy to implement expl(), logl() and all the other
 >  "long double" variants

 So, in the light of this, it is not _that_ easy because the numeric
 ranges are different. The algorithms can still be reused I think.

 Another fact is that the LDBL_* macros of <float.h> don't
 reflect reality. Eg LDBL_MAX is already infinity. The long double
 support in gdtoa also suffers from this.
 See the appended program. On an i386/current box, it gives:
 1.189731495357231765021264e+4932
 inf
 1.189731495357231765021264e+4932
 inf
 1.189731495357231599977347e+4932
 1.189731495357231632999029e+4932
 1.189731495357231599977347e+4932
 1.189731495357231632999029e+4932

 which shows that both the compiler and gdtoa don't know about the
 53-bit mantissa -- LDBL_MAX gets converted to a number which will
 cause an overflow, and the other number has some excess bits in
 the mantissa which will get lost if touched by the FPU later.

 If I set the TARGET_96_ROUND_53_LONG_DOUBLE as said in my previous mail,
 it gives:
 inf
 inf
 1.189731495357231765021264e+4932
 inf
 1.189731495357231632999029e+4932
 1.189731495357231632999029e+4932
 1.189731495357231599977347e+4932
 1.189731495357231632999029e+4932

 which shows that the constants generated at compile time are
 consistent with reality now, but gdtoa still assumes mantissa
 bits not usable for further calculations.

 Adapting the LDBL_ constants to the 53-bit mantissa should fix
 this, but some review is needed because at some points just
 mantissa lengths are compared to decide whether a "long double"
 is different from a "double".

 best regards
 Matthias





 -----------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------
 Forschungszentrum Juelich GmbH
 52425 Juelich

 Sitz der Gesellschaft: Juelich
 Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
 Vorsitzende des Aufsichtsrats: MinDirig'in Baerbel Brumme-Bothe
 Vorstand: Prof. Dr. Achim Bachem (Vorsitzender), Dr. Ulrich Krafft (stellv. 
 Vorsitzender)
 -----------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------

 --==_Exmh_24243496863520
 Content-Type: text/plain ; name="stest.c"; charset=us-ascii
 Content-Description: stest.c
 Content-Disposition: attachment; filename="stest.c"

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

 #define TRUE_LDBL_MAX 1.1897314953572316E+4932L
 long double null;

 void
 pr(long double d)
 {

 	printf("%.25Lg\n", d);
 	d += null;
 	printf("%.25Lg\n", d);
 }

 #define pr2(x) pr(x);pr(strtold(#x, 0))
 #define pr3(x) pr2(x)

 int
 main()
 {

 	pr3(LDBL_MAX);
 	pr3(TRUE_LDBL_MAX);

 	return 0;
 }

 --==_Exmh_24243496863520--

Home	PR Database Search