Floating point math accuracy, c++ vs fortran

user_12962785 0
Viewed 531

I have tried to implement a recursive function, in both C++ and Fortran, which calculates the value of the n'th Legendre polynomial, at x. In Fortran I have

recursive function legendre(n, x) result(p)
  integer, intent(in) :: n
  real(8), intent(in) :: x
  real(8) :: p

  if(n == 0) then
    p = 1.0
  else if (n == 1) then
    p = x
  else
    p = (2.0*real(n,dp)-1.0)*x*legendre(n-1,x)-(real(n,dp)-1.0)*legendre(n-2,x)
    p = p / real(n,dp)
  end if 
end function legendre

and then in C++ I have

double legendre(int n, double x) {
  double p;
  if(n == 0) return 1.0;
  else if(n == 1) return x;
  else {
    p = (2.0*(double)n - 1.0)*x*legendre(n-1,x)-((double)n - 1.0)*legendre(n-2,x);
    p /= (double)n;
    return p;
  }
}

These two functions seem to be exactly the same to me, both using double precission, but the result from the Fortran function is substantial different from the C++ result. For example, legendre(7,-0.2345) = 0.28876207107499049178814404296875 according to WolframAlpha. The two codes above, when compiled with no optimizations produce

Fortran : 0.28876206852410113

C++ : 0.28876207107499052285

I know that the answers should not be the same due to floating point arithmetic, but the difference in value here for double precision seems somewhat large to me. What is the reason that the Fortran value is so far off from the other two ?

fortran c++
2 Answers

Although the variables in your FORTRAN function are defined as double-precision (8 bytes), the constants you have specified are default (single-precision, 4-byte) values.

According to this discussion, that means the arithmetic is performed to single-precision accuracy:

Even if the variable that you are assigning the result to is defined to be DP, the Fortran standard requires that the arithmetic on the constants be performed using SP. That is because you are using default real constants, since you do not have any kind type parameter at the end of the constants. By rule, default real is SP.

And, further on in the same discussion:

...Starting with Fortran 90, published in June 1991, this practice of "promoting" SP constants to DP is prohibited.

So, in order to force double-precision maths, specify the constants as DP: instead of, for example, 1.0, specify 1.0D0 (and so forth for the others).

user_108710730

Thanks to Adrian's response, I was able to fix the problem. Nothing was wrong with the code in the Fortran function, the issue was the value of x which I was passing to it. Even though in the main program, I had defined x as real(8) and latter assigned it the value with a simple

x = -0.2345

which I thought should be double precision. It should actually be

x = -0.2345_dp

This results in the two functions having the same answer. I believe it is likely due to the reason that Adrian pointed out.

user_129627850
Related