small numerical differences in floating point result between wintel and Sun/SPARC

[Dick Hendrickson]

Gordon Burditt wrote:
...



The IEEE standard requires that, for numbers in this range, the
result must correctly rounded to the target precision. Since they're
both the same type and KIND (hence, the same target precision),
and since there is but one value that's closest to one tenth that's
representable as an IEEE single precision number, they must
give the same result.

Fortran is a little more restrictive in this case. It
requires that all constants written in the same form
have the same value. So 0.1 always has the same value
within a program. The non-IEEE parts of Fortran don't
specify how accurate the representation must be.

Dick Hendrickson

 

[Dik T. Winter]

> Also I have an inherent mistrust of options like /Op. gcc certainly
> doesn't get everything right with its version of that option, maybe other
> compilers do but it can be hard on x86 especially and I wouldn't trust
> it without some verification.

Indeed. Is an expression like a * b + c indeed calculated with 53 bits
mantissa only?

 

[JS]

How can you switch FPU to 64 bit mode? Thanks.

 

[JS]

Sorry, I hasn't been able to create a small program that demonstrates
the problem.

 

[Gordon Burditt]

How can you switch FPU to 64 bit mode? Thanks.

In an implementation-specific manner, if it's possible at all. Some
implementations have fpsetprec(FP_PE). Other Intel x86-based
implementations might require escaping to assembly language.
I know of no guarantee that if you set the precision, that
the compiler won't shortly thereafter change it back to what
it wants.

Gordon L. Burditt

 

[Tim Prince]

Also I have an inherent mistrust of options like /Op. gcc certainly
doesn't get everything right with its version of that option, maybe other
compilers do but it can be hard on x86 especially and I wouldn't trust
it without some verification.

Indeed. Is an expression like a * b + c indeed calculated with 53 bits
mantissa only?[/QUOTE]
Yes, because 53-bit precision mode is set, not on account of /Op. floats
also will be done with effective promotion of intermediates to double,
unless SSE is selected.

 

[Tim Prince]

In an implementation-specific manner, if it's possible at all. Some
implementations have fpsetprec(FP_PE). Other Intel x86-based
implementations might require escaping to assembly language.
I know of no guarantee that if you set the precision, that
the compiler won't shortly thereafter change it back to what
it wants.

Every compiler I've seen preserves the precision mode. True, this requires
code which saves the FPU state and restores it after fiddling modes to
support int cast. Compilers which set 53-bit mode do so only in
initialization, prior to execution.

 

[Merrill & Michele]

Sorry, I hasn't been able to create a small program that demonstrates
the problem.

Does the Standard have anything to say about machine epsilon? MPJ

 

[Victor Bazarov]

Does the Standard have anything to say about machine epsilon? MPJ

Which Standard do you mean to ask about? You cross-posted to three
different language newsgroups. The C++ Standard says that the member
'epsilon' of 'std::numeric_limits' returns the difference between 1
and the least value greater than 1 that is representable. It does not
say specify the actual value returned from 'epsilon'. The C Standard
does give the maximum acceptable values of 'FLT_EPSILON', 'DBL_EPSILON'
and 'LDBL_EPSILON', but your implementation is free to provide its own,
smaller values. And I have no idea about the Fortran Standard.

V

 

[James Giles]

Merrill & Michele wrote: ....

Which Standard do you mean to ask about? You cross-posted to three
different language newsgroups. The C++ Standard says that the member
'epsilon' of 'std::numeric_limits' returns the difference between 1
and the least value greater than 1 that is representable. It does not
say specify the actual value returned from 'epsilon'. The C Standard
does give the maximum acceptable values of 'FLT_EPSILON', 'DBL_EPSILON'
and 'LDBL_EPSILON', but your implementation is free to provide its own,
smaller values. And I have no idea about the Fortran Standard.

Fortran has a generic intrinsic function called EPSILON. The
following code prints the values of the machine epsilon for single
and double precision:

Print *, epsilon(1.0), epsilon(1.0d0)
End

The return value is b^(1-p), where b is the base of the floating point
numbers (usually 2) and p is the number of base-b digits in the significand
(including any hidden normalization digit). For IEEE single, the result
is 2^(-23). For IEEE double, the result is 2^(-52). Note that the result
doesn't depend on the argument's value, only it's type attributes, so
EPSILON(0.0) returns the same result as EPSILON(1.0), since both
arguments are default precision reals. If the implementation were to
support additional real precisions (like quad, or double-extended)
the same inquiry function could be applied to them. Similarly, if
an implementation were to support the proposed IEEE decimal
floating point, the same inquiry function could be applied to those
as well. The result is in the same floating-point precision (and
base) as the argument.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

 

[CBFalconer]

.... snip ...

Does the Standard have anything to say about machine epsilon? MPJ

c:\dnld>grep -n epsilon \stds\n869.txt
28465: epsilon FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON

Now is that so hard you couldn't do it for yourself?