1

I'm translating a code from Fortran. I get a weird behavior during assignment. I know that adding the code could be very helpful but I can't add the exact code (I'm not authorized) and I was not successful to replicate it.

The lines are the following (qk was predefined and qk1 was not):

   print*,"qk",qk     
   print*,"qk1",qk1
   QK1=QK
   print*,"qk",qk
   print*,"qk1",qk1

and I get these values printed:

    qk   21909779.000000000     
    qk1   6.44842193E+32
    qk   21909779.000000000     
    qk1   21909780.0

The point is that I would expect to get qk1 equal to qk... why are they different? When I try to replicate it, obviously I get the same values printed.

Since I didn't add the code I do not expect a precise answer... does anyone have any idea about what to check?

francescalus
  • 30,576
  • 16
  • 61
  • 96
scana
  • 111
  • 1
  • 11
  • 1
    are qk and qk1 defined as exactly the same type and precision? The number of decimal 0's may be indicating a difference in precision – rvbarreto Jan 19 '20 at 14:25
  • 1
    Probably not relevant, but the first line shown causes the value of a variable called `qt` to be printed. I'm always suspicious that what may be a careless typo in posting the code to SO is a result of carelessness in coding ... – High Performance Mark Jan 19 '20 at 14:31
  • @rvbarreto qt is "double precision" and qk1 is not defined – scana Jan 19 '20 at 14:31
  • @HighPerformanceMark yes you are right, the code is right but I made a mistake while copying. I'll correct it now. – scana Jan 19 '20 at 14:34
  • @rvbarreto that was it! Once I defined qk1 as "double precision" I get the same values. While translating I took it for granted that qk1 was defined, but it was not. Ok I solved the problem, many thanks to both of you. Nevertheless, I'm still curious... why fortran "added one"? – scana Jan 19 '20 at 14:40
  • 2
    So there is a probably a difference between the data type of qk and qk1 (real versus double precision), but to be able to tell this add a [minimal working example (MWE)](https://tex.meta.stackexchange.com/q/228) that illustrates your problem. – albert Jan 19 '20 at 14:40
  • 1
    No time to check at moment but you are close to the max sig figures that "single precision" real can provide, my suspicion is that due to this you can't represent the number exactly, and so you are getting the nearest approximation, which is 1 away. see https://stackoverflow.com/questions/56968179/exponentiation-in-fortran-gfortran-to-high-precision/56969033#56969033 and https://stackoverflow.com/questions/59581069/borwein-s-algorithm-for-the-calculation-of-pi-in-fortran-is-converging-too-fast/59581284#59581284 for something similar – Ian Bush Jan 19 '20 at 14:52
  • 1
    Your compiler tried to fix the mistake for you and used a different precision the number. Take a look at wikipedia article for 'Single-precision floating-point format' and also for 'Double-precision floating-point format' and pay attention to how IEEE 754 defines number representation. This may help to understand the limitations behind the precision you choose and also behind conversions – rvbarreto Jan 19 '20 at 14:55
  • Also, there are flags that you can use in your compiler to prevent this. If you are using gfortran, take a look at https://staff.washington.edu/rjl/uwamath583s11/sphinx/notes/html/gfortran_flags.html at Wconversion and other flags – rvbarreto Jan 19 '20 at 14:58
  • If you check the value of `spacing(qk1)` you will find it is equal to `2`. In other words, single precision floating point values in this region are either all odd or all even. In this case it seems they are all even. – RussF Jan 20 '20 at 02:56
  • 1
    An IEEE-binary32 floating point number, can represent all integers between ±2^24 accurately (±16777216), the number you try to represent is bigger and hence rounding will occur. `qk` is most likely a double precision number which can hold this value, but `qk1` is of lower precision. Hence the rounding. – kvantour Jan 20 '20 at 13:47

1 Answers1

6

The reason is that as indicated in the comments qk is single precision, qk1 double, and at the values required the spacing between single precision reals is 2:

Program one_out

  Use, Intrinsic :: iso_fortran_env, Only : real32, real64

  Implicit None

  Real( real64 ) :: qk64
  Real( real32 ) :: qk32

  qk64 = 21909779.0_real64
  qk32 = qk64

  Write( *, * ) 'qk64 = ', qk64
  Write( *, * ) 'qk32 = ', qk32

  Write( *, * ) '64 spacing ', Spacing( qk64 )
  Write( *, * ) '32 spacing ', Spacing( qk32 )

End Program one_out

ian@eris:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -std=f2008 one_out.f90 
one_out.f90:11:9:

   qk32 = qk64
         1
Warning: Possible change of value in conversion from REAL(8) to REAL(4) at (1) [-Wconversion]
ian@eris:~/work/stack$ ./a.out
 qk64 =    21909779.000000000     
 qk32 =    21909780.0    
 64 spacing    3.7252902984619141E-009
 32 spacing    2.00000000

The most important lesson here is always, always use Implicit None - and the second is that compiler warning are useful, turn them on and work out what they mean!

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
  • I wonder whether an explanation/link, at the level of the questioner, as to what "spacing" is (or why it matter in this case) would help? – francescalus Jan 20 '20 at 11:06
  • OK - be-meetinged most of the rest of the day but will do, or feel free to edit the answer yourself – Ian Bush Jan 20 '20 at 11:16
  • 2
    I'd probably write something like "Loosely, the value 21909779 falls in the gap between the nearby _representable numbers_. This gap is the spacing referred to: with a spacing of 2 between 21909778 and 21909780 it isn't possible to assign exactly the value 21909779". But that's a bit rubbish without a decent link. – francescalus Jan 20 '20 at 15:07