ffpe-trap=precision and I/O of real numbers

[Daniel Feenberg]

Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:

program test
write(*,*) 'a'
write(*,*) 50.
stop
end

and run it under gdb with a breakpoint at the first executable statement::

Breakpoint 1, test () at test.for:2
2 write(*,*) 'a'
(gdb) n
a
3 write(*,*) 50.
( gdb) n

Program received signal SIGFPE, Arithmetic exception.
0x00007ffff7ba5a21 in ?? () from /lib64/libgfortran.so.3

It appears that the attempt to write a small integer triggers an arithmetic exception. Which seems very unlikely.
.Wouldn't that be worthy of mention in the docs.?

Daniel Feenberg

[baf]

Something is broken in your installation/version of gfortran. Compiling
your program with version 11.2 I get the following output:

a
50.0000000

[Neil]

The issue seems to occur with gfortran 9.3.0:

laptop3:~> cat temp.f90
program test
write(*,*)'a'
write(*,*)50.
end

laptop3:~> gfortran --version
GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

laptop3:~> gfortran -ffpe-trap=precision temp.f90

laptop3:~> ./a.out
a

Program received signal SIGFPE: Floating-point exception - erroneous
arithmetic operation.

Backtrace for this error:
#0 0x7f6ada3dcd5a
#1 0x7f6ada3dbef5
#2 0x7f6ada21020f
#3 0x7f6ada62d385
#4 0x7f6ada62eb98
#5 0x7f6ada62ef35
#6 0x7f6ada62ff76
#7 0x5638fc62b25e
#8 0x5638fc62b2b0
#9 0x7f6ada1f10b2
#10 0x5638fc62b0dd
#11 0xffffffffffffffff
Floating point exception (core dumped)

[Cyrmag]

On 10/24/2021 4:13 PM, Daniel Feenberg wrote:
> Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
>
> program test
> write(*,*) 'a'
> write(*,*) 50.
> stop
> end
>

<--- CUT --->
>
> Daniel Feenberg
>

I see the bug with gcc version 10.2.0 (GCC) under Cygwin-64 bit on
Windows 10-64 bit.

-- CyrMag

[jfh]

Hmmm. My system is ubuntu with gfortran 9.3 and it shows the bug. But man gfortran does not include precision in the list of possible ffpe-trap=list options. Is that option legal? If so, what is it supposed to do? If not, should it give a compile-time error?

[JRR]

On 10/24/21 11:13 PM, Daniel Feenberg wrote:
> Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
>
> program test
> write(*,*) 'a'
> write(*,*) 50.
> stop
> end
>

> It appears that the attempt to write a small integer triggers an arithmetic exception. Which seems very unlikely.
> .Wouldn't that be worthy of mention in the docs.?
>
> Daniel Feenberg
>
>
>

I checked this for both gfortran 7.5 on Ubuntu 18 and 9.3 on Ubuntu 20,
which shows this feature. As was pointed out, -ffpe-trap=precision
doesn't exist (in the documentation).
gfortran
-ffpe-summary=invalid,zero,overflow,underflow,denormal,inexact prog.f90
shows
a
50.0000000
Note: The following floating-point exceptions are signalling:
IEEE_INEXACT_FLAG
So, correctly it seems to state that 50.000 cannot be represented
exactly. For some reason, the argument -ffpe-trap=precision
seems to switch on the option inexact. That to me seems to be indeed
some sort of bug for parsing the options of the ffpe-trap flag.
gfortran -ffpe-trap=invalid,zero,overflow,underflow,denormal prec.f90
works without any problems (leaving out the inexact setting)
Cheers,
JRR
--
Juergen Reuter
Theoretical Particle Physics
Deutsches Elektronen-Synchrotron
Hamburg, Germany
---------------------------------
invalid is desy .and. com is de

[Cyrmag]

On 10/25/2021 7:36 AM, JRR wrote:
> So, correctly it seems to state that 50.000 cannot be represented
> exactly.

That is not true; the exact representation of 50.0 as a 32-bit real is
Z'42480000'. Note that there are 19 zero bits at the least significant
end, so there is plenty of leeway for representing more complicated
numbers than 50.0 using a 23 bit significand.

There is something going wrong, causing a trap to be signaled.

-- CyrMag

[JRR]

Sorry, you are right, I misread the docu of gfortran. 'inexact' does not
mean that 50. cannot be represented exactly by e.g. 32bit-real, but that
there is a loss of precision during the performed operation. In that
case the real 50. loses precision when used within the write (*,*)
procedure. If you use write (*, "(ES10,3)") no trap is signaled.

[steve kargl]

If you compile your code, including the one in this thread, with
-ffpe-trap=inexact, it will likely die. For any (non-trivial) code,
which is doing floating point arithmetic, the code will raise FE_INEXACT.

For the particular problem with "write(*,*) 50.", the runtime library is
doing floating arithmetic to determine if it should write "50.00000" or
"0.5E+02" or a few other checks.

If you (generic you here) want to find the location, then compile your
code with "-g" and use gdb.

% cat > a.f90
write(*,*) 50.
end
% gfcx -o z -g a.f90
% ./z
50.0000000
% gfcx -o z -g -ffpe-trap=inexact a.f90
% ./z
Floating exception (core dumped)
% gdb ./z z.core

#0 0x0000000200a73dd9 in get_float_string (dtp=0x7fffffffe4b0,
f=<optimized out>, source=<optimized out>, kind=4, comp_d=1,
buffer=0x7fffffffde20 "", precision=<optimized out>, size=<optimized out>,
result=<optimized out>, res_len=<optimized out>)
at ../../../gccx/libgfortran/io/write_float.def:1110
1110 FORMAT_FLOAT(4,)
(gdb) list
1105 }
1106
1107 switch (kind)
1108 {
1109 case 4:
1110 FORMAT_FLOAT(4,)
1111 break;
1112
1113 case 8:
1114 FORMAT_FLOAT(8,)

So, you need to go read FORMAT_FLOAT in write_float.def.

--
steve

[jfh]

Thank you Steve. As you said elsewhere, -ffpe-trap=inexact is the modern equivalent of -ffpe_trap=precision and they are not often worth using. Gfortran does not have a bug here.

[Robin Vowels]

On Monday, October 25, 2021 at 8:13:05 AM UTC+11, feen...@gmail.com wrote:
> Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
>
> program test
> write(*,*) 'a'
> write(*,*) 50.
> stop
> end
>
> and run it under gdb with a breakpoint at the first executable statement::
>
> Breakpoint 1, test () at test.for:2
> 2 write(*,*) 'a'
> (gdb) n
> a
> 3 write(*,*) 50.
> ( gdb) n
>
> Program received signal SIGFPE, Arithmetic exception.
> 0x00007ffff7ba5a21 in ?? () from /lib64/libgfortran.so.3
>
> It appears that the attempt to write a small integer triggers an arithmetic exception.

.
There are no small integers in this program.

[gah4]

On Sunday, October 24, 2021 at 2:13:05 PM UTC-7, feen...@gmail.com wrote:
> Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
>
> program test
> write(*,*) 'a'
> write(*,*) 50.
> stop
> end

Some (many!) years ago, I was at a Kahan seminar on floating point
arithmetic, where he mentioned that no processors have an exception
for precision loss. After, I asked about IBM S/360 and successors,
which do have one. But only for all precision loss. If even one bit
is left, no interrupt.

Even more, the usual way to zero a register is to subtract it from
itself, so would cause the interrupt.

But otherwise, I/O library routines should save the flag state,
and set appropriate flags for their use. Then restore before
returning.

[jfh]

But there may be situations where the user wants to know precision has been lost even in I/O library routines. (I don't know of any but that may just be my ignorance.)

[gah4]

On Saturday, October 30, 2021 at 6:32:35 PM UTC-7, jfh wrote:

(snip, I wrote)

> > But otherwise, I/O library routines should save the flag state,
> > and set appropriate flags for their use. Then restore before
> > returning.

> But there may be situations where the user wants to know
> precision has been lost even in I/O library routines.
> (I don't know of any but that may just be my ignorance.)

In general, one knows what their own code does, but not what
library routines do. I suspect that is especially true in this case.

There isn't actually a "precision" exception, but "inexact", which
is raised by most floating point operations. But consider, more
likely, that you have underflow exception enabled. And it might
be that the I/O routines, in normal operation, generate underflow
exceptions. In most cases, you wouldn't want to know that.

If the I/O routines did something that you would want to know
about, you should be able to ask specifically for them, and the
routines should specifically test for them.

I almost wrote "library routines" without the I/O in front, and
probably also mean that. In many cases, library routines
should specifically test for things, and give the appropriate
exception, instead of just letting them happen.

For example, CABS should not generate an overflow when
the result doesn't overflow, (and is correct), even if
intermediate values overflow.

But as well as I know it, you likely never want inexact enabled.

[Robin Vowels]

On Sunday, October 31, 2021 at 6:05:45 AM UTC+11, gah4 wrote:
> On Sunday, October 24, 2021 at 2:13:05 PM UTC-7, feen...@gmail.com wrote:
> > Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
> >
> > program test
> > write(*,*) 'a'
> > write(*,*) 50.
> > stop
> > end
> Some (many!) years ago, I was at a Kahan seminar on floating point
> arithmetic, where he mentioned that no processors have an exception
> for precision loss. After, I asked about IBM S/360 and successors,
> which do have one. But only for all precision loss. If even one bit
> is left, no interrupt.
>
> Even more, the usual way to zero a register is to subtract it from
> itself, so would cause the interrupt.

.
A FP register can be cleared by either loading zero into it,
or by loading zero from another register that has zero in it.

[Robin Vowels]

On Sunday, October 31, 2021 at 6:05:45 AM UTC+11, gah4 wrote:

> On Sunday, October 24, 2021 at 2:13:05 PM UTC-7, feen...@gmail.com wrote:
> > Is there somewhere a discussion of gfortran option "ffpe-trap=precision"? I tried recompiling my code with all the debugging traps enabled, but I couldn't make "Arithmetic Exception" go away. Eventually I made this test program:
> >
> > program test
> > write(*,*) 'a'
> > write(*,*) 50.
> > stop
> > end
> Some (many!) years ago, I was at a Kahan seminar on floating point
> arithmetic, where he mentioned that no processors have an exception
> for precision loss. After, I asked about IBM S/360 and successors,
> which do have one.

.
They have three.
1. Loss of significance, when the result of a float operation yields a zero
mantissa.
2. Exponent underflow.
3. Exponent overflow.

.
> But only for all precision loss. If even one bit
> is left, no interrupt.
>
> Even more, the usual way to zero a register is to subtract it from
> itself, so would cause the interrupt.

.
There are several ways that zero can be placed in a FPR.
1. Load zero from storage;
2. Load from another register that has zero in it;
3. Subtract the register from itself, but first setting the
significance mask bit to zero with SPM and restoring it
after the operation.