CODEWITHSUNDEEP
fortranprecisionfortran90
user3225087
user3225087

Reputation: 225

d0 when taking roots of numbers

So in general, I understand the difference between specifying 3. and 3.0d0 with the difference being the number of digits stored by the computer. When doing arithmetic operations, I generally make sure everything is in double precision. However, I am confused about the following operations:

64^(1./3.) vs. 64^(1.0d0/3.0d0)

It took me a couple of weeks to find an error where I was assigning the output of 64^(1.0d0/3.0d0) to an integer. Because 64^(1.0d0/3.0d0) returns 3.999999, the integer got the value 3 and not 4. However, 64^(1./3.) = 4.00000. Can someone explain to me why it is wise to use 1./3. vs. 1.0d0/3.0d0 here?

Upvotes: 0

Views: 255

Answers (2)

agentp
agentp

Reputation: 6989

this is a peculiar fortuitous case where the lower precision calculation gives the exact result. You can see this without the integer conversion issue:

 write(*,*)4.d0-64**(1./3.),4.d0-64**(1.d0/3.d0)

 0.000000000 4.440892E-016

In general this does not happen, here the double precision value is "better"

 write(*,*)13.d0-2197**(1./3.),13.d0-2197**(1.d0/3.d0)

 -9.5367E-7 1.77E-015

Here, since the s.p. calc comes out slightly high it gives you the correct value on integer conversion, while the d.p. result will get rounded down, hence be wrong, even though the floating point error was smaller.

So in general, no you should not consider use of single precision to be preferred.

in fact 64 and 125 seem to be the only special cases where the s.p. calc gives a perfect cube root while the d.p. calc does not.

Upvotes: 0

M. S. B.
M. S. B.

Reputation: 29391

The issue isn't so much single versus double precision. All floating point calculations are subject to imprecision compared to true real numbers. In assigning a real to an integer, Fortran truncates. You probably want to use the Fortran intrinsic nint.

Upvotes: 2

Related Questions