Subject: Re: accuracy of "long double"
To: Neil Booth <>
From: Matthias Drochner <>
List: port-amd64
Date: 08/20/2007 15:17:52
This is a multipart MIME message.

Content-Type: text/plain; charset=us-ascii said:
> b) Exponents have a wider range, so overflows or underflows that
>    would happen on doubles don't necessarily happen on long
>    doubles.

Yes, I had figured that too, see my second reply to the PR.

> d) Despite many claims to the contrary that you will read on the
>    net, and on GCC lists, there is no double-rounding problem.

The famous "paranoia" test program finds a number of rounding
flaws eg on Linux. I've extracted the failing tests for
multiplication - see attachment.
The tests succeed on NetBSD. If I set the FPCW to its initial
value, ie 64-bit rounding (by enabling the "#if 0" part),
the tests fail. With -mfpmath=sse, the tests succeed.
So doubles are still not rounded correctly if the FPU is used
in extended precision mode.

best regards

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. 

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

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

#define FLOAT double

FLOAT Zero = 0.0, One = 1.0, OneAndHalf = 1.5;
FLOAT T, X, Y, Y1, Y2, Z;
FLOAT U2 = EPSILON; /* probed in "paranoia" */


	Y2 = One + U2;
	X = OneAndHalf - U2;
	Z = (X - U2) * Y2;
	Z = Z - X;
	return (Z == Zero);


	Y1 = One - U2;
	X = OneAndHalf - U2;
	Y = OneAndHalf + U2;
	T = Y * Y1;
	T = T - X;
	return (T <= Zero);

#if 0
        const int fpcw = 0x37f;
        asm("fldcw %0"::"m"(fpcw));
	if (!mult1())
		printf("* rounding error 1\n");
	if (!mult2())
		printf("* rounding error 2\n");
	return 0;