Fortran COMPLEX calculates different from C++

Question

I have completed a port from Fortran to C++ but have discovered some differences in the COMPLEX type. Consider the following codes:

PROGRAM CMPLX
    COMPLEX*16 c
    REAL*8 a
    c = (1.23456789, 3.45678901)
    a = AIMAG(1.0 / c)
    WRITE (*, *) a
END

And the C++:

#include <complex>
#include <iostream>
#include <iomanip>

int main()
{
    std::complex<double> c(1.23456789, 3.45678901);
    double a = (1.0 / c).imag();

    std::cout << std::setprecision(15) << " " << a << std::endl;
}

Compiling the C++ version with clang++ or g++, I get the output: -0.256561150444368 Compiling the Fortran version however gives me: -0.25656115049876993

I mean, doesn't both languages follow the IEEE 754? If I run the following in Octave (Matlab):

octave:1> c=1.23456789+ 3.45678901i
c =  1.2346 + 3.4568i
octave:2> c
c =  1.2346 + 3.4568i
octave:3> output_precision(15)
octave:4> c
c =  1.23456789000000e+00 + 3.45678901000000e+00i
octave:5> 1 / c
ans =  9.16290109820952e-02 - 2.56561150444368e-01i

I get the same as the C++ version. What is up with the Fortran COMPLEX type? Am I missing some compiler flags? -ffast-math doesn't change anything. I want to produce the exact same 15 decimals in C++ and Fortran, so I easier can spot porting differences.

Any Fortran gurus around? Thanks!

Is there a type with the assignment statement? Should it be `double a = (1.0 / c.imag());`? — Thomas Matthews, Jan 12 '15 at 19:32

High Performance Mark · Accepted Answer · 2015-01-12T19:42:45.757

5

In the Fortran code replace

c = (1.23456789, 3.45678901)

with

c = (1.23456789d0, 3.45678901d0)

Without a kind the real literals you use on the rhs are, most likely, 32-bit reals and you probably want 64-bit reals. The suffix d0 causes the compiler to create 64-bit reals closest to the values you provide. I've glossed over some details in this, and there are other (possibly better) ways of specifying the kind of a real number literal but this approach should work OK on any current Fortran compiler.

I don't know C++ very well, I'm not sure if the C++ code has the same problem.

If I read your question correctly the two codes produce the same answer to 8sf, the limit of single precision.

As for IEEE-754 compliance, that standard does not cover, so far as I am aware, the issues of complex arithmetic. I expect the f-p arithmetic used behind the scenes produces results on complex numbers within expected error bounds in most cases, but I'm not aware that they are guaranteed as error bounds on f-p arithmetic are.

edited Jan 12 '15 at 19:42

answered Jan 12 '15 at 19:36

High Performance Mark

77,191
7
105
161

d0 suffix solved this simple example, thanks for that. But, I need to check for more similar problems tomorrow. Please hang around :) – user1284878 Jan 12 '15 at 19:52
3

Btw, In C++ all float literals are double by default. That's what got me. – user1284878 Jan 12 '15 at 19:55
I tried d0 at work with more complicated Fortan code and yes, it does the trick. Scary because they obviously was using float-precision for constants. Not any more :) Tanks for the quick answer! – user1284878 Jan 13 '15 at 09:45

score 0 · Answer 2 · answered Jan 12 '15 at 19:34

0

I would propose to change all Fortran contants to DP

1.23456789_8 (or 1.23456789D00) etc

and use DIMAG instead of AIMAG

answered Jan 12 '15 at 19:34

Severin Pappadeux

18,636
3
38
64

1

Note that `dimag` is not Fortran, and that use of `_8` is not portable. – francescalus Jan 12 '15 at 19:55

Fortran COMPLEX calculates different from C++

2 Answers2