Sure, telling gfortran to use glibc as the libm for everything would be reasonable. All I’m saying is that at the end of the day, it’s gfortran’s responsibility to make sure that gfortran doesn’t introduce bugs into user code, and currently it’s failing to do that.
yeah, I’ve just re-opened the issues you mentioned. They should be relatively easy to fix.
That is a really weird specification, since it means that an implementation of arctan that isn’t processor-dependent doesn’t follow the spec. Also, if you read “processor-dependent” to mean that the approximation can only depend on the processor (and X obviously), gfortran is still possibly breaking the standard by also returning a different result based on whether the value is evaluated at compile time or runtime.
Ah, that makes more sense. I still think gfortran is not standards compliant in this case since the same compiler gives 2 different answers.
The 06 value is 0x3ff1b6e192ebbe45.
Indeed, the ...06 value is the nearest of ...04
Curious. The .s file seems to show why, and the compiler code would show how it got there, but wondering what this variant shows in that programming environment …
program test
use,intrinsic :: iso_fortran_env, only : wp => real64
use,intrinsic :: iso_fortran_env, only : stderr=>error_unit,stdin=>input_unit,stdout=>output_unit
real(wp),parameter :: x = 2.0_wp
character(len=*),parameter :: g='(*(g0,1x))'
real(wp) :: y
character(len=20) :: mode
write(stdout,g)'1.10714871779409050301706546017853704007004764540143264667653920743 ...'
y = x
write(stdout,g)'nearest ', nearest (atan(x),+1.0_wp) ! Nearest representable number
write(stdout,g)'nearest ', nearest (atan(x),-1.0_wp) ! Nearest representable number
write(stdout,g)'spacing ', spacing (atan(x)) ! Smallest distance between two numbers of a given type
write(stdout,g)'digits ', digits (atan(x)) ! Significant digits function
write(stdout,g)'epsilon ', epsilon (atan(x)) ! Epsilon function
write(stdout,g)'precision ', precision (atan(x)) ! Decimal precision of a real kind
inquire(unit=stdout,round=mode)
write(stdout,g)'round ', mode
write(stdout,g) atan(x)
write(stdout,g) atan(2.0_wp)
write(stdout,g) atan(y)
write(stdout,g) atan(x).eq.atan(y)
write(stdout,g) atan(2.0_wp) == atan(y)
write(stdout,g) atan(x-epsilon(x))
write(stdout,g) atan(x+epsilon(x))
write(stdout,g) atan(nearest(x+epsilon(x),+1.0_wp))
write(stdout,g) atan(nearest(x-epsilon(x),-1.0_wp))
write(stdout,g,round='up') atan(nearest(x+epsilon(x),+1.0_wp)), atan(nearest(x-epsilon(x),-1.0_wp)), 'UP'
write(stdout,g,round='down') atan(nearest(x+epsilon(x),+1.0_wp)), atan(nearest(x-epsilon(x),-1.0_wp)), 'DOWN'
write(stdout,g,round='zero') atan(nearest(x+epsilon(x),+1.0_wp)), atan(nearest(x-epsilon(x),-1.0_wp)), 'ZERO'
write(stdout,g,round='nearest') atan(nearest(x+epsilon(x),+1.0_wp)), atan(nearest(x-epsilon(x),-1.0_wp)), 'NEAREST'
write(stdout,g,round='compatible') atan(nearest(x+epsilon(x),+1.0_wp)), atan(nearest(x-epsilon(x),-1.0_wp)), 'COMPATIBLE'
write(stdout,'(b0)') atan(x), atan(y)
write(stdout,'(b0)') x, y
end program test
I don’t have any problem with .58 ULP. For functions like atan, I think anything under 1 ULP is justifiable. That, said, I think it’s a fairly major problem for the result to depend on whether gfortran decides to constant prop the answer.
I’m not seeing where the 0.42 and 0.58 numbers are coming from. When I do the arithmetic, it looks like the ***6 result is closer than the ***4 result, not the other way around. The former is high by +.097 in the last decimal digit, while the latter is low by -.103. However, that is an odd way to compare errors, since arbitrary list-directed i/o was used to print the values. What if extra digits had been printed in the first place? That is why I suggested the hex format to look at the actual bits, that eliminates one extra level of noise in the output.
BTW, as an aside, it requires an extension to the language to print a real number with a Z format. Either that or tricks with transfer() or equivalence.
In this case, we have an extended precision result to compare to, the compiler does not have that information. I personally would not call this a bug in the compiler, or even an error in the compile time or the runtime libraries. This is just an artifact of floating point arithmetic. Programmers are supposed to know that these kinds of differences can happen. As I said before, they can be caused by rounding modes, algorithm choices, truncation errors, using extended registers for intermediate results, optimization levels, and on and on. Certainly anyone who has programmed fortran with a cross-compiler would not be surprised with seeing such small differences in compile-time and run-time evaluations.
I don’t have that document handy right now for comparison, but ironically, I was just reading a post of yours (kargl) in c.l.f. from a few years ago where you were explaining how to use equivalence in order to print a double precision value with a Z format. In any case, if this is fixed now in f2018, that is a good thing. I personally probably use Z formats more for printing real values than integer values, and I always liked the fact that gfortran supported that as an extension (before it was standardized).
Speaking of f2018, isn’t there now a standard format for reading/printing floating point values in a hex interchange format, more specifically with a decimal exponent and a hex mantissa? I think this takes care of the hidden bit issue too, which is sometimes confusing.
As for the tan() values, isn’t it correct that the printed value ***6 is closer to the extended precision value than the ***4 value? I guess what is happening is that the printed values are enough different from the actual binary values to flip the errors the other way around. That is why looking at the bits is important.
1 Like
julia> atan(2.0)
1.1071487177940904
julia> atan(big(2))-atan(2.0)
9.40447137356637993148121741102451826466765392e-17
julia> atan(big(2))-nextfloat(atan(2.0))
-1.2799989118936750876991415950791887985332346e-16
The 04 value is closer. The printing happens the way it does because most languages default to printing the shortest decimal expansion that uniquely identifies a floating point number.
1 Like
You might argue the processor here does not abide by the principle of least surprise. Is that a bug? It is up to the gfortran users to decide and follow-up collectively, as appropriate.
The EX descriptor serves for that. But it seems not be supported yet by gfortran, at least ver. 11.3:
program hex
real :: x=1.5, y=atan(2.0)
print '(2ex20.8)', x,y
end program hex
>gfortran hexreal.f90
hexreal.f90:3:18:
3 | print '(2ex20.8)', x,y
| 1
Error: Period required in format specifier E at (1)
>ifort hexreal.f90 && ./a.out
0XC.00000000P-3 0X8.DB70D000P-3
FYI: here is the clang results on the same Mac:
#include <stdio.h>
#include <math.h>
int main()
{
double x,y;
x = 2.0;
y = atan(x);
printf("%21.18le %21.18le\n", x, y);
return 0;
}
clang -O0 test.c -o test
Result is:
2.000000000000000000e+00 1.107148717794090631e+00
1.1071487177940906
The decimal precision of 8-bytes real is, as returned by precision(1d0), 15. We are discussing the 17th digit.
1 Like