This is a multipart MIME message.

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


neil@daikokuya.co.uk 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
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_34042622989200
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
#define EPSILON DBL_EPSILON

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

int
mult1()
{

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

int
mult2()
{

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

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

--==_Exmh_34042622989200--