4

I have translated some code from Fortran to C++ and both codes give me the same result for a given input with the exception of two data points in the middle of my data set.

My code calculates the distance between points and does some interesting things with that information. Two points in the C++ code are found to be at one distance from each other and at a different distance in Fortran. The code is long, so I won't post it.

This strikes me as weird because the two "strange points" are right in the middle of my code, whereas all of the other 106 points behave the same.

I have already read the Goldberg paper, and it makes me believe that real and float ought to be the same on my 32-bit system.

nbro
  • 15,395
  • 32
  • 113
  • 196

2 Answers2

1

A real in Fortran may be float (which is kind 4) or double (kind 8) in C++.

It also may depend on your compiler options (i.e. math extensions, optimization, square root implementation, etc).

nbro
  • 15,395
  • 32
  • 113
  • 196
Pixelchemist
  • 24,090
  • 7
  • 47
  • 71
0

In most C/C++ implementations you'll encounter, float corresponds to REAL*4, and double corresponds to REAL*8.


This StackOverflow answer (and related comments) describe Fortran 90's types: Fortran 90 kind parameter.


Differences in floating point computations may arise due to different evaluation order. Floating point arithmetic is very sensitive to evaluation order, especially where addition and subtraction among values with a wide dynamic range is involved.

Also, C/C++ math and math libraries default to double precision in many surprising places, unless you explicitly ask for a float. For example, the constant 1.0 is a double precision constant. Consider the following code fragment:

float x;
x = 1.0 + 2.0;
x = x + 3.0;

The expression 1.0 + 2.0 is computed at double precision, and the result cast back to single precision. Likewise, the second statement x + 3.0 promotes x to double, does the arithmetic, and then casts back to float.

To get single precision constants and keep your arithmetic at single precision, you need to add the suffix f, as follows:

float x;
x = 1.0f + 2.0f;
x = x + 3.0f;

Now this arithmetic will all be done at single precision.

For math library calls, the single-precision variant also usually has an f suffix, such as cosf or sqrtf.

nbro
  • 15,395
  • 32
  • 113
  • 196
Joe Z
  • 17,413
  • 3
  • 28
  • 39