Bruce Evans wrote:
> On Wed, 5 Mar 2008, Colin Percival wrote:
>> David Schultz wrote:
>>> On Wed, Mar 05, 2008, Colin Percival wrote:
>>>> Yes and no. Double rounding isn't allowed;
>>>
>>> Yes, it is.
>>
>> No, it isn't. The following code
>
> Yes it is :-). Extra precision for intermediate values implies double
> rounding for the final results, and extra precision for intermediate
> values is a fundamental part of C (it used to be only that float
> expressions were (de-facto standard required to be) evaluated in double
> precision, but now it affects double expressions too, and C99 has vast
> complications to support i387'-ish extra precision.
Let me clarify: Double rounding isn't allowed by IEEE754. To quote from
the October 5th draft of IEEE 754R (which I happen to have in front of me
right now):
5.1 Overview
All conforming implementations of this standard shall provide the operations
listed in this clause for all supported floating-point formats available for
arithmetic. Each of the computational operations that return a numeric result
specified by this standard shall be performed as if it first produced an
intermediate result correct to infinite precision and with unbounded range,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
and then rounded that intermediate result to fit in the destination’s format
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
This is clarified further in clause 10:
10.1 Expression evaluation rules
Clause 5 of this standard specifies the result of a single arithmetic
operation going to an explicit destination. Every operation has an explicit
or implicit destination. One rounding occurs to fit the exact result into a
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
destination format. That result is reproducible in that the same operation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
applied to the same operands under the same attributes produces the same
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
result on all conforming implementations in all languages.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Programming language standards might define syntax for expressions that
combine one or more operations of this standard, producing a result to fit an
explicit or implicit final destination. When a variable with a declared format
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
is a final destination, as in format conversion to a variable, that declared
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
format of that variable governs its rounding. The format of an implicit
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
destination, or of an explicit destination without a declared format, is
defined by language standard expression evaluation rules.
In other words, in the code
double x, y, z;
z = x * y + 1.0;
it is perfectly legal for (x * y) to be computed with extra precision; but in
the code
double x, y, z;
z = x * y;
z = z - 1.0;
both the multiplication and the subtraction must be rounded *once* to double
precision.
> In C, with FLT_EVAL_METHOD = 2 (evaluate all FP expressions in the range
> and precision of the long double type) [...]
As someone wrote earlier, the authors of C99 were a bit confused. I can only
assume that they intended to be IEEE754-compliant and the rules concerning
FLT_EVAL_METHOD apply to *implicit* destinations only.
>>> 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.
>
> Nah, it's quite feasible, especially for transcendental functions since
> transcendental numbers in general can be approximated more closely by
> rationals than algebraic numbers.
In fact, one of the common ways to prove that a real number is transcendental
is to show that it can be approximated more closely by rationals than any
algebraic function can be.
> The evaluation just needs to be
> done in extra precision
How much extra precision do you need? If the exact value of a transcendental
function is extremely close to halfway between two consecutive floating-point
values, how do you decide which way to round it?
Computing transcendental functions to within (0.5+x)ulp for small positive
values of x is certainly feasible, but computing provably correctly rounded
transcendental functions is very very hard.
>> Fortunately, library functions aren't required to have any particular
>> error bounds.
>
> We require error bounds of < 1 ulp and aim for < 0.5[0]1 ulps.
What standard states that those bounds are required (or can be relied upon
when proving that code is correct)? Intel FPUs over the years have been
variously documented to compute transcendental functions to within 1.0, 2.5,
or 3.5 ulp -- are you saying that the 8087 doesn't comply with the standard
which it inspired?
Colin Percival