|
From: Irek S. <isz...@gm...> - 2013-03-06 10:12:02
|
fwd to valgrind-developers because this is an issue for developers Irek ---------- Forwarded message ---------- From: Irek Szczesniak <isz...@gm...> Date: Wed, Mar 6, 2013 at 12:22 AM Subject: [Valgrind-users] Floating point DATA CORRUPTION with 80bit long long double!! (Re: floating point print error in ada) To: Philippe Waroquiers <phi...@sk...> Cc: Bob Rossi <bo...@br...>, val...@li... On Tue, Mar 5, 2013 at 11:56 PM, Philippe Waroquiers <phi...@sk...> wrote: > On Tue, 2013-03-05 at 18:54 +0100, Lionel Cons wrote: >> (1) in https://bugs.kde.org/show_bug.cgi?id=197915#c9 is a joke: >> >Julian Seward 2010-07-12 15:58:25 UTC >> > As per comment #0, adding support for 80-bit floats is low priority, because (1) AIUI the majority of floating point code is portable >> > and restricts itself to 64-bit values, >> >> The majority of _consumer_ software uses using double (aka 64bit >> float), but the majority of _scientific_ software (for example the >> whole NIH bioinformtics software stack or 99.9% of CERNs simulation >> software) is relying on long long double aka 80bit or 128bit floats >> (depending on platform, AMD64 uses 80bits). valgrind is useless for >> such software. > I am not too sure about the proportion of consumer software > versus scientific software. Assuming there is more consumer software, > the note (1) above is not such a joke at the end :). > > Reading the gcc manual, wouldn't it be a good idea to have the > scientific software to be rewritten (or at least compilable) > so as to use sse ? > gcc manual tells for `sse' > The resulting code should be considerably faster in the > majority of cases and avoid the numerical instability What does 'numerical instability' mean? Unless 387 has bugs I don't know about I think this refers to the general limitations of floating point being implemented using a fixed-width datatype - which wider-than-64bit-double datatypes can partially work around (at least the amount of work to compensate using algorithms like the Kahan summation can be greatly reduced, for example using a Kahan derivative for 64bit floating point requires a 370% overhead to compensate while using a 80bit floating point datatype reduces this to less than 120%). Remember SSE, together with MMX, fall into the category of multimedia accelerators (i.e. designed for throughput and low precision calculations) and are not designed for high precision calculations. This is not an option. Otherwise you have to compensate using the Kahan summation algorithm or similar, which puts the performance win ad absurdum. > problems of 387 code, but may break some existing code that > expects temporaries to be 80bit. It's not the temporaries which have to be 80bit. I think Lionel is referring to the fact that scientific software usually requires the maximum precision available. 64bit double usually will NOT work, at least not with extra work. > > This is the default choice for the x86-64 compiler. > > And as a bonus, you can run it under Valgrind on x86/amd64 :). > If this scientific code is fully portable, Where does the claim come from that this is more portable? If the platform doesn't have a long long double which is wider than double than this is very unfortunate and requires a LOT of extra calculations, but if the platform (like AMD64) supports a wider datatype then it is of course used at all costs. > then you could also decide > to run it under Valgrind e.g. on ppc32/ppc64 systems. > > Note: at my work, we are using Ada/gnat/gcc on x86/amd64. > The application code is compiled with sse. > We contemplated recompiling and/or changing the > Ada runtime to use sse only and fully avoid the 80 bits. > However, as at the end, we found very little impact (at least for > our apps) of running the 80 bits runtime on Valgrind, we have kept > the default gnat runtime (which uses 80 bits floats here and there). > As long as these 80 bits computation are "ok" if computation is > in reality done with 64 bits float, then not much impact. > YMMV. Uh. The impact here is DATA CORRUPTION, caused by doing calculations to be meant using 80bit datatypes being done with 64bit datatypes and the remaining bits filled with garbage. If this can't be solved or is not going to be solved then valgrind should ABORT instead of causing a silent corruption of data. Irek ------------------------------------------------------------------------------ Symantec Endpoint Protection 12 positioned as A LEADER in The Forrester Wave(TM): Endpoint Security, Q1 2013 and "remains a good choice" in the endpoint security space. For insight on selecting the right partner to tackle endpoint security challenges, access the full report. http://p.sf.net/sfu/symantec-dev2dev _______________________________________________ Valgrind-users mailing list Val...@li... https://lists.sourceforge.net/lists/listinfo/valgrind-users -- Irek |
|
From: Roland M. <rol...@nr...> - 2013-03-06 16:43:05
|
On Wed, Mar 6, 2013 at 10:19 AM, <pa...@fr...> wrote:
[snip]
> Is this really true? I would expect that Valgrind does all calculations at 64bit double precision and converts. This will result in truncation and/or under/overflow.
>
> My feeling is that unless the usually small numerical differences change your control flow, then just ignore the differences. The aim of testing with Valgrind isn't to validate numerical results, it is to validate memory use or performance or threading.
Erm... the problem is that this can cause major issues in numerical
simulations... ;-((
Below is a small example (yes yes... I know it's not 100% portable...
it should use at least |LDBL_DECIMAL_DIG| ... but I'm not going to
twist my mind for that unless there is demand for it) for AMD64/64bit:
-- snip --
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#include <stdio.h>
int main(int ac, char *av[])
{
long double i;
unsigned long long numiterations=0ULL;
puts("# start.");
for(i=1.L ; i < 1.00000000001L ; i=nextafterl(i, 5.L))
{
numiterations++;
}
printf("# done (after %llu/0x%llx iterations).\n",
numiterations, numiterations);
return EXIT_SUCCESS;
}
-- snip --
(this is more or less reduced from real-world simulation code which
uses |nextafterl()|&co. to iterate over some variations of the input
values to weed-out stuff like singularities etc.)
The expected output looks like this:
-- snip --
# start.
# done (after 92233720/0x57f5ff8 iterations)
-- snip --
On valgrind it just hangs (e.g. if you pass in a |double| casted to
|long double| into |nextafterl()|'s first argument and then cast the
|long double| result to |double| then the result is usually identical
to the input value... as result the algorithm is stuck in an endless
loop).
----
Bye,
Roland
--
__ . . __
(o.\ \/ /.o) rol...@nr...
\__\/\/__/ MPEG specialist, C&&JAVA&&Sun&&Unix programmer
/O /==\ O\ TEL +49 641 3992797
(;O/ \/ \O;)
|
|
From: John R. <jr...@bi...> - 2013-03-06 17:27:21
|
> long double i; > for(i=1.L ; i < 1.00000000001L ; i=nextafterl(i, 5.L)) > On valgrind it just hangs (e.g. if you pass in a |double| casted to > |long double| into |nextafterl()|'s first argument and then cast the > |long double| result to |double| then the result is usually identical > to the input value... as result the algorithm is stuck in an endless > loop). Noticing the unexpected i == nextafter(i, ...) is one of the checks performed by higher-quality software when verifying the execution environment before launching into lengthy computations that depend on that environment. -- |
|
From: Roland M. <rol...@nr...> - 2013-03-06 17:36:13
|
On Wed, Mar 6, 2013 at 6:28 PM, John Reiser <jr...@bi...> wrote:
>> long double i;
>
>> for(i=1.L ; i < 1.00000000001L ; i=nextafterl(i, 5.L))
>
>> On valgrind it just hangs (e.g. if you pass in a |double| casted to
>> |long double| into |nextafterl()|'s first argument and then cast the
>> |long double| result to |double| then the result is usually identical
>> to the input value... as result the algorithm is stuck in an endless
>> loop).
>
> Noticing the unexpected i == nextafter(i, ...) is one of the checks
> performed by higher-quality software when verifying the execution environment
> before launching into lengthy computations that depend on that environment.
Huh ? My example never used "==" anywhere (first guess is that the
mail body encoding used "==" to represent a single '=' character) ...
and using '"==" to compare a floating-point value against another
value is likely not going to generate the results you may expect...
... which is more or less a deja-vu since we had that topic in the
AT&T AST/ksh93 development list:
---------- Forwarded message ----------
From: Roland Mainz <rol...@nr...>
Date: Fri, Feb 8, 2013 at 12:48 PM
Subject: [ast-users] Floating-point arithmetic, rounding and testing
for equality... / was: Re: Fwd: Floating point oddities with pow()?
Accuracy problem?
Cc: ast-users <ast...@re...>, ast...@re...
[snip]
> What we're facing here seems that different implementations of a libm
> function may result different results at the last digits, right? Is
> there a function to test how accurate the libm functions are, and is
> there a libm function which tests for equality in a specified
> precision?
Testing whether two _floating_point_-values are equal is a tricky
business (e.g. the C99/ksh93 operator "==" is usually of little use...
;-( ) because of the rounding errors or (AFAIK in your case when the
i387 version of IEEE 754-2008 floating-point values are used) extra
precision caused by the issue that |long double| is 80bit on x86/AMD64
but some machine instructions have extra bits available during
calculation (like the new FMA instructions vs. a "manual"
multiply–accumulate operation... or |pow()| vs. manual "power of" ...)
or the registers have more bits than they store in memory.
Another issue may be (or often is) the base10<---->base2 conversion
which causes rounding errors when strings are converted to IEEE
754-2008 floating-point values and back. As "fix" C99 ([1]) added the
printf "%a" format to represent floating-point values in a hexadecimal
floating-point representation which can be used (with your example) to
"visualise" the issue that the difference is usually off by one or two
bits at the trailing end.
Keeping the email short (due lack of time... I may try a more detailed
answer when I have time) ... the standards only define |isgreater()|,
|isgreaterequal()|, |islessequal()|, |islessgreater()|,
|isunordered()| but intentionally no |isequal()| because the same
operation using different implementations is usually unlikely to
result in _exactly_ the same result.
A quick&&dirty implementation for a |isequal()| (maybe a better term
may be |isnearequal()|) may look like this (ksh93 syntax):
-- snip --
# test whether values in variables "a" and "b" are
# "equal" with a precision of 0.000001
if (( fabs(a - b) < 0.000001 )) ; then
-- snip --
(yes yes... I'm using the wrong mathematical terms (e.g. "precision" ?
"resolution" ? "distance" ?) in this case... please correct me... ;-(
)
[1]=ksh93 supports C99 "hexfloat" via it's printf builtin (note that
you have to use $ float varname ; ... ; printf "%a\n" varname # and
NOT $ float varname ; ... ; printf "%a\n" ${varname} # because the
${varname} means that the internal IEEE 754-2008 floating-point value
is converted to a base10 floating-point string value first and THEN
passed to the printf builtin utility) and via $ typeset -lX varname #
----
Bye,
Roland
--
__ . . __
(o.\ \/ /.o) rol...@nr...
\__\/\/__/ MPEG specialist, C&&JAVA&&Sun&&Unix programmer
/O /==\ O\ TEL +49 641 3992797
(;O/ \/ \O;)
_______________________________________________
ast-users mailing list
ast...@li...
http://lists.research.att.com/mailman/listinfo/ast-users
--
__ . . __
(o.\ \/ /.o) rol...@nr...
\__\/\/__/ MPEG specialist, C&&JAVA&&Sun&&Unix programmer
/O /==\ O\ TEL +49 641 3992797
(;O/ \/ \O;)
--
__ . . __
(o.\ \/ /.o) rol...@nr...
\__\/\/__/ MPEG specialist, C&&JAVA&&Sun&&Unix programmer
/O /==\ O\ TEL +49 641 3992797
(;O/ \/ \O;)
|