Different Precision in C++ and Fortran

Any platform where Fortran `double precision` is different from C `double`?

Depends on your idea of "practical compiler", I guess.

Neither the C standard nor the Fortran standard require IEEE conformance, but most[citation needed] compilers adhere to it (as do most processors, so this is as much a matter of optimization as it is of conformance). Cray and VAX processors are the notable exceptions, so you may find a compiler for those systems that conforms to the processor's implementation of float instead. (e.g., on a Cray, a float is 8 bytes).

For more: http://www.quadibloc.com/comp/cp0201.htm

Why the c and fortran versions of this same program produce different results?

The difference is due to the order of evaluation combined with the limited precision of the floating point type.

If you change the Fortran version to

t = t + (a(i,k) * a(i,k) + b(k,j) * b(k,j))

i.e. add parenthesis around the terms with a and b, you get the same result for both languages. The C version already uses this order of evaluation due to the use of the += assignment operator.

As mentioned in the comments, this is expected behavior at the limits of the available precision.

Why is the same floating-point constant printed differently in Fortran and in C?

By default, Fortran's REAL number constants are single-precision; In C, however, floating point literals have double precision.

When you translate 0.8804418 to single precision and then print it as a double in C, you get 0.8804417849 printed (demo)

float x = 0.8804418f;
printf("%.10g\n", x);

Fortran's printout appears to be the same number rounded up.

Fortran's syntax for double-precision REAL numbers uses suffix d:

print *, 0.8804418d+0

This prints 0.88044180000000005 (demo).

Discrepancy between the values computed by Fortran and C++

In Fortran, by default, floating point literals are single precision, whereas in C/C++ they are double precision.

Thus, in your Fortran code, the expression for calculating NUMERATOR is done in single precision; it is only converted to double precision when assigning the final result to the NUMERATOR variable.

And the same thing for the expression calculating the value that is assigned to the DENOMINATOR variable.

Fortran accuracy and speed vs. C

Modifying the Fortran program such that it corresponds to the C version by making all calculations in double precision:

program calcpi
  implicit none
  integer :: i
  integer, parameter :: p = selected_real_kind(15)
  real(p) :: pi

  pi=0.0_p
  do i = 0,1000000000
     pi = pi + 1.0_p/(4.0_p*i+1.0_p)
     pi = pi - 1.0_p/(4.0_p*i+3.0_p)
  end do

  pi = pi * 4.0_p

  write(*,*) pi

end program calcpi

Compiling with -O2 using GCC 4.4.3 on x86_64-linux-gnu on a Xeon X3450 (2.67 GHz) I get the following timings:

$ time ./calcpi_c 
Pi=3.141593

real    0m13.903s
user    0m13.860s
sys 0m0.010s
$ time ./calcpi_fort 
   3.1415926530880767     

real    0m13.876s
user    0m13.840s
sys 0m0.000s

IOW, they are more or less indistinguishable. Which is about what one would expect for such a simple example.

Parameter precision in fortran

Note: this is essentially the same as Francescalus' answer, but with some extra fluf.

The issue at hand here is related to the IMPLICIT statement.

The Fortran standard makes the following statements :

5.5 IMPLICIT statement

1. In a scoping unit, an IMPLICIT statement specifies a type, and
   possibly type parameters, for all implicitly typed data entities whose
   names begin with one of the letters specified in the statement.
   Alternatively, it may indicate that no implicit typing rules are to
   apply in a particular scoping unit.

2. <snip>

3. In each scoping unit, there is a mapping, which may be null, between
   each of the letters A, B, ..., Z and a type (and type parameters).
   An IMPLICIT statement specifies the mapping for the letters in its
   letter-spec-list. IMPLICIT NONE specifies the null mapping for all
   the letters. If a mapping is not specified for a letter, the
   default for a program unit or an interface body is default integer if
   the letter is I, J, ..., or N and default real otherwise, and the
   default for an internal or module procedure is the mapping in the host
   scoping unit.

So in short, by default everything is of type default REAL unless it starts with I,J,...,N, then it is of type default INTEGER.

In the example in the question, the variables n1 and nT1 are hence default INTEGER unless specified otherwise. Thus the following might be a solution :

subroutine skip(IUNIT,NBYTS)
   implicit double precision (A-H,O-Z)
   character c*1
   integer(kind=integer_kind), parameter :: n1 = 1024, nT1 = 8*n1
   ...
end subroutine skip

As a general remark on variable declaration I would like to make the following remarks:

• use implicit none as default. It makes debugging easier
• avoid star notation, it is not part of the standard for anything but characters, and for characters it is declared obsolescent.
• make use of kind parameters declared in the intrinsic module iso_fortran_env
• be aware that double precision does not necessarily mean "double precision" in the literal or IEEE sense of the word. It just means, twice the storage of a real and a higher precision than real. (a higher precision could be one-digit).

---

Related Topics