forums.silverfrost.com

[mecej4]

This short program runs when it is compiled with FTN95 8.70 as a 32-bit EXE, and prints -5.000000E+07.

When the source is then compiled with /64 and run, the run is aborted after reporting "floating point division by zero".

[JohnCampbell]

Ver 8.62 produces an "Access violation" for /64 /lgo
Ver 8.62 produces an "Floating point division by zero" for /64 /debug /lgo

Gfortran 10.2 produces :
GCC version 10.2.0
Infinity

[mecej4]

Thanks, John.

The test program is an offshoot of an attempt to understand some odd behavior (in the Hydrus6 software from USGS). Really old Fortran did not have hyperbolic functions as intrinsic functions, so the authors used an ASF to implement it, and large parts of these old codes are in single precision.

Here is a different attempt on the issue (approximate 1/tanh(x) - 1/x for small x).

[JohnCampbell]

Your results are unusual, but I can not reproduce your results.
I changed the code slightly to

[PaulLaidler]

This is just a general comment...

In the first program z is small and tanh(z) is approximately z. So wetab is the difference of two large values that are almost equal. When the calculation is successful it will be highly suseptible to round-off error.

32 bit extended FTN95 precision real(3) using the intrinsic tanh may give the best result. 32 bit double FTN95 precision real(2) might use extended precision internally for part of the calculation.

Single precision is not an option here and double precision probably gives a very poor result.

[mecej4]

Thanks again, John.

I obtained the peculiar line-11 with Gfortran 9.3-64-bit on Cygwin, 4.8.5-64-bit as well as 6.2.0-64-bit on WSL1/Ubuntu14.04.1, and 7.1.1 on a Linux cloud service.

I'll have to install a newer 10.x version and see.

[JohnCampbell]

mecej4,

I only use a pre-built gFortran for windows from equation.com, which is what clearwin_64 uses.
Don't have 9.3, 4.8.5 or 7.1.1 but do have 6.2.0 which produces on i5-2300

[mecej4]

John, there IS something strange happening here, but it seems to affect only Gfortran, and probably depends on how various versions of Gfortran react to the precision loss that Paul commented upon.

On my Dell XPS-17 laptop with an i7-2720QM, running W10-64, Gfortran 9.3 from Cygwin gives for y the strange value 1.0000E+00 on Line-11, but Gfortran 9.2 from Equation.com gives 0.0000E+00. I even asked the two compilers to generate an assembler file. The two assembler files were identical except for the lines containing

.ascii "GCC version 9.3.0"

in one and

.ascii "GCC version 9.2.0"

in the other. Therefore, the differences must be caused by the Gfortran RTLs and how those differ in the two distributions of GCC.

[JohnCampbell]

Could you include:
print '(i3,4ES12.4)',i,x,y,z,tanh(x)
or
t = tanh(x)
print '(i3,4ES12.4)',i,x,y,z,t

It is a strange value of x to produce a different result !

[mecej4]

John, I did what you asked. Starting from line-5, the values of tanh(x) were the same as that of x to the five significant digits printed. Here are lines 10-12, as produced by Gfortran 9.3.0 (Cygwin):

[JohnCampbell]

y = 1/x - 1/tanh = (1-x/tanh)/x

So x/tanh is not 1.0 for 1.024e-7 but is for 5.12e-7 and 2.048e-8 ?

tanh - x must be non-zero in your compiler, which implies an error in the TANH function interpolation.
See if the following identifies the non-zero "tanh - x" (due to rounding ?)

[mecej4]

John, even with those versions of Gfortran that gave the two 1.0000 s on Line-11, I find that replacing

x = x*0.2

by

x = x/5.0

restores expected/normal results.

Strange to find the difference between 0.2 and 1/5.0 to cause a change in a result from 0.0E0 to 1.0E0. Perhaps, a striking example for use in a course on programming?

[JohnCampbell]

I am thinking that the actual value x = 1.0240E-07, that is generated using x = x*0.2 is causing a roundoff error in the interpolation formula for tanh (x),
so that tanh(x) /= x,
so y /= 0

The value of x is very close to epsilon(x), so that a roundoff error would produce y = 1

If tanh(x) is returned as x*(1-epsilon) ; the nearest representable value less than x, ie roundoff error in the calculation of tanh

y = one/x - one/tanh
= 1/x - 1 / ( x*(1-epsilon) )
= ( x*(1-eps) - x ) / ( x * x*(1-eps) )
= x*eps / ( x * x )
= 1 when x = epsilon (x)

you have "stumbled" on a unique value of x where tanh (x) fails due to roundoff, which is only present in some versions of the gFortran tanh library.

[mecej4]

John, the following program may help us get closer to the heart of the problem. It does the same calculation as before, but with values of x close to the problematic value. With versions of Gfortran on Windows that use the implementation of tanh() in the Microsoft VC RTL, such as Mingw, the last column is always zero. With Cygwin on Windows, MS WSL-1 and on Linux, we can speculate that the underlying GNU math library causes different behavior, with the last column giving either 1.00000E+00 or 0.00000E+00.

[JohnCampbell]

I used 10 format(i2,ES16.7,3x,Z8,ES16.7,ES15.5)

and got on i5-2300