Post by principiantePost by Michael Mair[snip: Demonstration of floating point excess precision "problem",
gcc options to get rid of it]
Post by principiantePost by Michael MairNonetheless, your original program and your "corrections" were both
not apt to expose the excess precision weaknesses as there were
still C issues.
Maybe there is a misunderstanding here. I think that the excess
precision is perfectly fine, the weakness stands in the interaction
between C and the excess-precision floating point engine.
Well, the problem is not exactly a C problem as the C semantics
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Post by principiantePost by Michael Mairwork fine; the problem is that gcc breaks the C semantics on
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Post by principiantePost by Michael Mairpurpose to cope with a broken architecture.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Post by principiantePost by Michael MairC89 does not tell you anything at all about floats and doubles
apart from <FLT/DBL/LDBL>_<EPSILON/MAX/MIN> (and the derived
stuff) and we only demand 6/10/10 digits precision and exponents
ranging from 10^-37 to 10^37.
This is completely acceptable. Apart from that, the trick with
the volatile-pointer evaluation works portably.
In C99, you can find out the rounding mode, the evaluation mode
and -- if available (with gcc never) -- can be sure that IEC 60559
(IEEE 754) support is available if __STD_IEC_559__ is defined,
have access to the floating point environment and so on.
Post by principianteApparently, there is no way to ensure portability; if run on
a processor without excess precision, my example with the
double type would give the expected result.
Yes. Which is why you enforce this behaviour by either making
your whole environment slow when compiling or enforcing the
correct precision with a cast (or the working replacement)
where it is necessary.
Infact the cast doens't work. This was quite a surprise for me.
"if (c == a/b)" and "if (c == (float)(a/b))" give exactly
the same code.
Yes, this is what I was referring to in the "underlined" part
above. Only some version of the *((volatile float *)&variable)
will work with gcc on all platforms because volatile enforces
an evaluation of the _memory_ where 'variable' is stored which
will enforce truncation/rounding to float.
Rather ugly.
Post by principiantePost by Michael MairPortability is mostly not possible across different architectures
if you do not use a language enforcing the use of IEEE floats
and doubles.
Post by principianteSee for instance
http://cch.loria.fr/documentation/IEEE754/ACM/addendum.html
which is an Addendum (by David Goldberg) of the well-known
paper What Every Computer Scientist Should Know About Floating-Point
Arithmetic, by Kahan.
I rather like the version on the Sun server where they have
another addendum suggesting to take up the C99 idea of defining
integer types like int8_t, int16_t, ..., int_least8_t, ... int_fast8_t,
... and so on and use it on floating point variables.
It's the same addendum, indeed.
You are right -- this completely escaped me as the "good parts"
(compared with http://docs.sun.com/source/806-3568/ncg_goldberg.html
[just search for 'C99']) are missing...
Post by principiantePost by Michael MairPost by principianteTo ensure portability, it seems that the only solution is to force
the compiler to store each result of a single computation in
the right precision memory slot.
Once again: Portability is usually not possible. Even if you have
the same base and the same mantissa/exponent bits, the architecture
may round differently for representation or operations.
Apart from that, C demands that all operations with floats are
evaluated in at least double precision which can give you problems with
IEEE 754 conformance.
I thought it was not anymore. On page 198 of the K&R, Second Edition,
it's written that "There are two changes here. First, arithmetic on
float operands may be done in single precision, rather then double;
the first edition specified that all floating arithmetic was double
precision".
Weeeell, this is what you can read out of it. But in order to do
anything with it, you have to be able to _find_out_ and _change_ the
evaluation method.
Most implementations evaluate either in double or in long double
(you can find this out by checking FLT_EVAL_METHOD in <float.h>
if you go for C99:
[from my gcc's <float.h>]
-1 indeterminate
0 evaluate all operations and constants just to the range and
precision of the type
1 evaluate operations and constants of type float and double
to the range and precision of the double type, evaluate
long double operations and constants to the range and
precision of the long double type
2 evaluate all operations and constants to the range and
precision of the long double type
)
And there are still some situations where you are effectively
working with doubles. I remember having "fun" with the arguments
to a variadic function some years ago :-/
Post by principiantePost by Michael MairUsually, all people pretend that the involved floating point numbers
were IEEE 754 conforming and that there are no differences in
rounding.
BTW: The respective parts of the GSL documentation also are good to
read to get started with floating point arithmetic problems and
possible solutions.
I have seen that there is a new Floating Point standard from IEEE
that is being discussed. Are they also addressing these kinds of
problems?
http://grouper.ieee.org/groups/754/revision.html
I have not had a close look at it yet but AFAIK they only provide
the standard -- it is up to languages or the people making the
hardware to provide support and demand/provide conforming
"implementations". Usually, the programming languages miss out on
this as they have to stay compatible. This is why I find having
a "stdfloat.h" with types flt32_t, flt64_t, flt80_t, flt128_t, ...
and flt32_IEC559_t, and so on such a charming idea.
As there is nearly no real world demand, it probably will never
happen in the PC area -- fast evaluation and give or take relative
errors in the magnitude of sqrt(DBL_EPSILON) at the end of the
computation are prefered (especially if your theoretical model
is coarse enough to easily outdo this error by several orders of
magnitude which is the case in many applications)...
Post by principianteThanks for the very interesting discussion,
You're welcome -- apart from that, I learned something too :-)
Before, I was not aware that gcc also gives you trouble with
floats, as (float) typecasts worked just fine for me in the
situations I have considered up to now...
Cheers
Michael
--
E-Mail: Mine is a gmx dot de address.