Hallo Waldek,
thanks for the quick answer. I see your argumentation - but in that case
I would have expected a value that is too small but not zero. Let me
argue about this ...
For example for an 64 bit IEEE real the result should be something like
2.22045E-16.
If - on an Intel processor - all variables are held in the 80 bit
floating point registers we would get something like 1.08420E-19 which
also is well representable in 64 bits even when truncated.
See the attach "C" based example when setting useregister to 1 or 0 and
recompiling it even with
gcc -O3 ... (please note am not a "C" programmer - so be patient with my
coding there ;-)
In Modula you can provoke that with a code like
PROCEDURE MachEps() : LONGREAL;
VAR eps;
BEGIN
eps:=0.5;
REPEAT
eps:=0.5*eps;
UNTIL ((1.0+eps) = 1.0);
RETURN 2.0*eps;
END MachEps;
where it is much more likely that you will get the 1.084E-19 result than
in my previous example.
But my point is that the optimizer "thinks" to be very clever and sets
eps to 0.0 because that is the only value where $1+\eps = 1$ in strict
mathematical logic - but nonsense in floating point arithmetic. So you
are right - one should not trust on the 2.22044E-16 result (which I
didn't ;-) but 0.0 is definitely wrong !
I detected that because some algorithms finalized too late (e.g.
calculation Taylor series expansions) and wasted computer time ... or
even showed other underflow errors.
In general the observation points (in my view) to some
"over-optimizations" which might occur in other case as well (e.g. test
"IF (x # x) THEN" to detect an NAN according to IEEE 754 or using guard
digits for problematic summations (attached an example procedure using
this technique - I need to investigate if GM2 is "problematic" here as
well ).
In general I would like to have at least a compiler hint such as
"eliminating redundant code" to give the programmer a chance to see
something might be not as expected in such cases.
Point me to my misunderstanding if I am mistaken - I am willing to learn.
Gruß
Michael
PS: Using MachEps:=LowLong.ulp(1.0) is the better choice - but that was
not present when I coded parts of my stuff in the area before ISO M2 ;-(
Am 01.12.2015 um 03:13 schrieb Waldek Hebisch:
> Michael Riedl wrote:
>> Hallo Gaius,
>>
>> another curious result.
>>
>> MODULE TestMachEps;
>>
>> FROM STextIO IMPORT WriteLn;
>> FROM SLongIO IMPORT WriteFloat;
>>
>> PROCEDURE MachEps() : LONGREAL;
>>
>> VAR eps,EinsPlusEps : LONGREAL;
>> BEGIN
>> eps:=3D0.5; EinsPlusEps:=3D1.5;
>> WHILE (EinsPlusEps # 1.0) DO
>> eps:=3D0.5*eps;
>> EinsPlusEps:=3D1.0+eps;
>> END;
>> RETURN 2.0*eps;
>> END MachEps;
>>
>> VAR eps : LONGREAL;
>>
>> BEGIN
>> eps :=3D MachEps();
>>
>> WriteLn; WriteLn;
>> WriteString(" MachEps =3D "); WriteFloat(eps,20,10);
>> WriteLn; WriteLn;
>> END TestMachEps.
>>
>> Compile attached with -O0 and with -O2. In the latter the machine=20
>> precision equals to zero ...
>> By the way - same behavior with REAL.
> This is classic example of nonportable program which may work
> or not depending on optimization settings/machine. In short:
> in the loop compiler may use higher precision (normal on
> 32-bit Intel processors). At return higher precision results
> gets rounded to lower precision. Due to underflow this gives
> 0. You should get result you expect using SSE math, maybe
> Gfortran uses this by default...
>
TestMachEps.c
SumVek.mod