Date: Thu, 6 Mar 2008 02:57:50 -0500 From: David Schultz <das@FreeBSD.ORG> To: Colin Percival <cperciva@FreeBSD.ORG> Cc: src-committers@FreeBSD.ORG, Bruce Evans <bde@FreeBSD.ORG>, cvs-all@FreeBSD.ORG, cvs-src@FreeBSD.ORG Subject: Re: cvs commit: src/sys/i386/include _types.h Message-ID: <20080306075750.GA48995@zim.MIT.EDU> In-Reply-To: <47CF9586.70707@freebsd.org> References: <200803051121.m25BLE03035426@repoman.freebsd.org> <20080305182531.GS68971@server.vk2pj.dyndns.org> <20080306021222.GA46783@zim.MIT.EDU> <47CF5D19.3090100@freebsd.org> <20080306033246.GA47280@zim.MIT.EDU> <47CF7EBF.6000009@freebsd.org> <20080306063452.GB48339@zim.MIT.EDU> <47CF9586.70707@freebsd.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, Mar 05, 2008, Colin Percival wrote:
> David Schultz wrote:
> > On Wed, Mar 05, 2008, Colin Percival wrote:
> >> Setting the i387 FPU to 53-bit precision gives standards-compliant
> >> behaviour whether people are using "double" or "long double".
> >
> > Not quite.
>
> How is it not standards-compliant?
First, because there are other issues, like constant folding and
dynamic precision and rounding, that gcc doesn't get right.
Second, because of the issues I mentioned with library functions,
and because of ill-defined meanings of constants such as
LDBL_MANT_DIG, etc.
> >> Yes and no. Double rounding isn't allowed;
> >
> > Yes, it is.
>
> No, it isn't. The following code
[...]
> should never output "BUGGY MATH", since x + y should round up to
> 1 + 2^(-52); but if you set the i387 to use extended precision, it
> gets the wrong answer.
True, but this is idealistic, and standards do not guarantee what
you want because it is too hard to do efficiently in practice.
Even the IEEE 754R draft basically says that double rounding is
allowed as long as the language makes it clear what is going on:
The last operation of many expressions is an assignment to an
explicit final destination variable. As a part of expression
evaluation rules, language standards shall specify when the next
to last operation is performed by rounding at most once to the
format of the explicit final destination, and when by rounding as
many as two times, first to an implicit intermediate format, and
then to the explicit final destination format.
See the description of C99's FLT_EVAL_METHOD macro, which is part
of C99's answer to the above.
If you asked Kahan, he'd probably tell you that both strict and
non-strict evaluation should be options. He was vehemently opposed
to the idea that Java have only a strict mode (similar to the one
you propose).
> >>> The downside is that this breaks long doubles.
> >> Except that it doesn't. Using 53-bit precision for long doubles is
> >> a perfectly valid and standards-compliant option.
> >
> > It's impossible to implement efficient library functions that
> > produce correct long double results
>
> You could have stopped this sentence here -- for all practical purposes,
> correctly rounded trigonometric functions are not feasible.
First, I said nothing about trigonometric anything, and second,
your claim here isn't true anyway. FreeBSD provides trig functions
that are correctly rounded most of the time (which is what most
people call "correctly rounded"). Both the IBM Accurate
Mathematical Library and Sun's crlibm provide trig functions that
are correctly rounded all of the time.
> > when the FPU is set to 64-bit
> > mode and avoid double rounding and related problems when the FPU
> > is set to 53-bit mode.
>
> Fortunately, library functions aren't required to have any particular
> error bounds. It's always nice if they are correct to within 0.6ulp
> instead of 0.8ulp or suchlike -- but nobody should be proving code
> correctness based on assumed properties of library transcendental
> functions. People should, and do, prove code correct based on the
> assumption that double precision arithmetic behaves like double precision
> arithmetic, however (myself included).
No, you misunderstand my point. I'm not talking about a 1-ulp
error, I'm talking about billions of ulps. In some long double
functions, we actually have to do things like the following to get
results that are even sane at all:
#ifdef __i386__
if (fpgetprec() != FP_PE)
return (double_precision_foo(x));
#endif
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20080306075750.GA48995>