0.0000000567 becomes 0.000000000, why?

[Baosheng_B_]

Hi, there is a parameter sigma which should be equal to 5.67*10**(-8) in my programming.  When the code is like this

          sigma= 5.67*10**(-8)

         L=0.1 !m

         a=1.0

         e=1.0

         ua=2000.0 !K

         dx=L/21.0

         t=sigma*dx*(a*ua**4-e*298.0**4)

        print*, t

the results is t=0.00000000

However, when the code is like this

         sigma= 0.0000000567

         L=0.1 !m

         a=1.0

         e=1.0

         ua=2000.0 !K

         dx=L/21.0

         t=sigma*dx*(a*ua**4-e*298.0**4)

        print*, t

I get the exact result that t=4317.871

I don't know why I got two different results because of the sigma function. Could you tell me why? 

[Craig_Dedo · Accepted solution]

You don't tell us what compiler you are using, or compiler version, operating system with version, or hardware.  Therefore, I can only guess as to what is going on.  Right now I don't have the time to run your source code through a compiler to see what, if anything, changes.

It also appears that you have little understanding of how Fortran represents numbers or how Fortran evaluates arithmetic expressions.

Here is my guess.  It appears that the root cause of your problem is roundoff in operations on floating-point numbers.  With numbers that can include a fractional component, the computer cannot represent most real numbers exactly, so they are represented by the closest approximation.  When you perform arithmetic operations on such numbers, the computer also needs to use the closest approximation for the result.  Hence, with every operation, there is a small amount of error between the exact result and the closest approximation that the computer can generate.

All Fortran compilers must provide at least two methods of representing floating-point, also known as real, numbers.  These are known as single precision and double precision.  Many Fortran compilers provide more than two.  For example, the GFortran compiler on Windows provides four different precisions.  Most, but not all, Fortran compilers use the IEEE definitions for single and double precision.  With these definitions, single precision is 7 decimal digits and an exponent range of 38, double precision is 15 decimal digits and an exponent range of 308.

By rule, in Fortran, default real numbers are single precision.  If you want any other precision, you need to specify it using a kind type parameter when you define a variable or you use a numeric constant.  In your example, all of the numeric constants do not have kind type parameters after them, so they are single precision.  Most likely, they are accurate to only 7 decimal digits.

In the first version, the program is computing sigma by a series of arithmetic operations.  The compiler may be computing a value that is too small to be represented, so it may be assigning a value of zero to sigma.  Again, this is only a guess.

I need to point out that in Fortran, and many other languages, using notation like 10**(-8) to specify exponential notation is not correct.  So instead of writing:

[fortran]

sigma = 5.67*10**(-8)

[/fortran]

you should write:

[fortran]

sigma = 5.67E-08

[/fortran]

The "E" in the notation stands for 10 raised to the power of the exponent after the "E".

I strongly recommend that programmers use double precision for calculating results whenever possible because you get far less roundoff with each arithmetic operation.  I have a module of commonly used named constants that I use in all programs and procedures.  Among many other things I define the named constants for all of the precisions that I need for real number.  E.g., I use these definitions:

[fortran]

Integer, Parameter :: SP = Selected_Real_Kind (P=6,R=35)

Integer, Parameter :: DP = Selected_Real_Kind (P=15,R=300)

[/fortran]

Using these definitions you would define sigma as double precision by writing:

[fortran]

sigma = 5.67E-08_DP

[/fortran]

Hope this helps.

[mecej4]

Baosheng B. wrote:

          0.0000000567 becomes 0.000000000, why?

Try this very short Fortran program, and consult Fortran textbooks and manuals to understand why the result printed is correct, even if it is not what you expect.

[fortran]program real
real sigma
sigma=10**(-1)
write(*,*)sigma
end[/fortran]

[Baosheng_B_]

Thank you Craig Dedo, it works. The fortran I'm using is fortran 90, intel compiler.

[John_Campbell]

Baosheng,

Mecej4 indicated the problem with your coding, perhaps the following will more clearly show the problem.

[fortran]    real*8 sigma, s
    sigma = 5.8*10**(-8)
    s     = 5.8*10.**(-8)
    write (*,*) sigma,s
    end[/fortran]

The following perhaps better shows what is required to improve the precison due to the implied type and kind of each operation.

[fortran]    real*8 sigma, s, s8,ss8
    sigma = 5.8   * 10   **(-8)
    s     = 5.8   * 10.  **(-8)
    s8    = 5.8d0 * 10.  **(-8)
    ss8   = 5.8d0 * 10.d0**(-8)
    write (*,*) sigma
    write (*,*) s
    write (*,*) s8
    write (*,*) ss8
    write (*,*) 5.8d-8
    end[/fortran]

Craig,

I can understand the use of kind qualifiers, as in:
    Integer, Parameter :: I8 = Selected_Int_Kind  (R=10)
    write (*,*) 10_i8**10, 10**10

But, Why "5.67E-08_DP" !!

John

[Baosheng_B_]

Hi Mecej4 I tried what you said, but it didn't work,  sigma is still written as 0.00000 on the screen. Thank  you all the same.

Hi John, I tried what you said. they worked except this one:

03	sigma = 5.8   * 10   **(-8)

Without a "." or "d0" after the values, it still shows 0.0000. Thank you very much. 

[DavidWhite]

Baosheng,

I wonder whether this last case is a matter of integer arithmetic.  You are raising an integer to an integer power, and so the result is an integer.  In this case, 1E-8 as an integer is zero.

David

[John_Campbell]

The problem is integer arithmetic, due to the precidence of **

5.8*10**(-8)  becomes   5.8 * ( 10**(-8) ) and 10**(-8) becomes integer = 0

The problem with ( 10.**(-8) ) is that this becomes real(4) so 1.e-8 with only 7 figures of accuracy.

To further explain the alternatives I provided, to identify both type and kind problems, when comparing 10.d0 ** (-8) to 10. ** (-8), even though 10._4 converts to the same value as 10._8 ( the extra ~32 bits are all zero), the resulting value becomes real(4) and 1.e-8 differs from 1.d-8, as the extra ~32 bits of 1.d-8 are not zero. I'd expect that the same 64-bit calculation would be performed, but for 10._4 the result returned would be converted/truncated to 4 bytes. ( I think there are compiler options to retain the register values which may retain the precision in calculating the s8 example I provided)

Similarly s = 5.8_4 and s = 5.8_8 are different, while s = 58._4 and s = 58._8 would result in the same value for s, as the value 58. is a "rational" binary number.

John

[Baosheng_B_]

Hi David,

Thank you for your reminding.  it's a matter of integer arithmetic. I tried other examples, let sigma1=30**(-3), sigma2=30.**(-3), sigma3=30**(-3.0), sigma4=30.**(-3.0), respectively. Only does sigma1 equal to zero, the rest of other three sigmas get the same right value.

[Baosheng_B_]

Hi John,

I see. If I want to get a real number, not a integer or zero, there should be at least a real number(following a dot "." or something) before or after "**".

Baosheng

[Xj_Kong]

Tips: Whenever you need to represent any "exact" number (say sigma = 5.67E-08, tau=0.99991200) in fortran, just type sigma = 5.67D-8 and tau=0.999912d0 for short and you get the accurate recognition by compilers. 

[Craig_Dedo]

David and Baosheng:

Yes, it is a matter of integer arithmetic.  I completely overlooked this very important issue in my rather lengthy answer yesterday.  My bad.  Still, everything I said in my post of yesterday is correct, as far as it goes.

In the assignment statement, sigma = 5.67*10**(-8), exponentiation has a higher precedence than multiplication or division.  Therefore, the exponentiation is performed first.  In exponentiation between two integers, the Fortran standard requires the use of integer arithmetic.  Therefore, 10**(-8) is evaluated as 1 / 100000000, using integer division, not floating-point division.  Thus the result is zero, not 1.0E-08.  This is specifically required by the rules in sections 7.1.5.2.1 and 7.1.5.2.2 of the Fortran 2008 standard.  If you want the quotation of these rules, I will be happy to provide them in a separate message.

Thus, if you want to use the form of 10**exponent instead of the standard exponential notation, you will need to include a decimal point in the base and possibly also in the exponent.  E.g., you would write sigma = 5.67*10.0**(-8) or sigma = 5.67*10.0**(-8.0).  Of course, if you did not specify any kind type parameter, you would by default get single precision, which on most Fortran compilers is 7 decimal digits and an exponent range of 38.  It is still very much preferred to use the standard method of exponential notation using E+nn or E-nn.

By the way, in case you are wondering, in an operation of raising a integer base to a real or complex exponent, the base is first converted to real or complex.  Raising a negative-valued real base to a real power is prohibited.

Hope this helps.

[Craig_Dedo]

John Campbell wrote:

Craig,

I can understand the use of kind qualifiers, as in:
    Integer, Parameter :: I8 = Selected_Int_Kind  (R=10)
    write (*,*) 10_i8**10, 10**10

But, Why "5.67E-08_DP" !!

John

Using constant_kind, where constant is a numeric, character, or logical constant and kind is the kind type parameter, is the standard method of specifying kind type parameters for constants since Fortran 90 was published in June 1991.  Thus, 5.67E-08_DP means exactly the same thing as 5.67D-08.  Without a trailing kind parameter, if you use E for exponent notation, you get single precision and if you use D for exponent notation, you get double precision.  In order to specify kinds other than single or double precision, you have to use the kind parameter at the end.  I find that using kind parameters for everything is neater and more consistent that sticking with the older practice from Fortran 77.

[John_Campbell]

We've travelled a long way on this post, where the answer was probably just: use  "sigma = 5.67E-8"

[FortranFan]

John Campbell wrote:

We've travelled a long way on this post, where the answer was probably just: use  "sigma = 5.67E-8"

Why do you compound the confusion now by dropping the "_DP" part; Craig Dedo's post was a highly valuable lesson in its usefulness.

[John_Campbell]

FortranFan,

I suggested sigma = 5.67e-8 as there was no initial indication of the kind of sigma, although t was reported to 7 figures.

While 5.67E-08_DP may provide a valuable lesson in the definition of precision, I regard this as an "overloaded" constant, where 5.67d-8 is more typical. I don't think we should adopt 5.67E-08_8 ? Can I describe this as the more typical way of defining a real(8) constant ?

Craig Dedo's post is a valuable lesson as 5.67E-08_16 is a necessary way of defining a real*16 constant, which should probably be provided as 5.67E-08_QP.

John

[FortranFan]

Craig Dedo wrote:

...

Among many other things I define the named constants for all of the precisions that I need for real number.  E.g., I use these definitions:

[fortran]

...
   Integer, Parameter :: DP = Selected_Real_Kind (P=15,R=300)

[/fortran]

...

Note the above-mentioned SELECTED_REAL_KIND function standardized with Fortran 95 has at least two parts: the precision part with the P argument and the range part with the R argument.

Also, note most FORTRAN 77 coders are used to "DOUBLE PRECISION" attribute to their data types for real numbers requiring more than single precision.  Now, note on most systems (IA32 processors, etc.), types declared with "DOUBLE PRECISION" attribute have a precision of 15 digits and range of 307 (roughly speaking, the type can represent numbers between -1E+307 and -1E-307 and +1E-307 and +1E+307).

So if one wants to be completely consistent with code written with DOUBLE PRECISION attribute in FORTRAN 77 (e.g., procedure calls involving older libraries/DLLs) , one should use KIND and RANGE functions on a given system to determine the precision and range of such types.  One may find a fully consistent named constant is

[fortran]

   INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(P=15, R=307)

[/fortran]

However, one can also take a shortcut as

[fortran]

   INTEGER, PARAMETER :: DP = KIND(1D0)

[/fortran]

which is a style I'm not fond of.

[FortranFan]

John Campbell wrote:

FortranFan,

I suggested sigma = 5.67e-8 as there was no initial indication of the kind of sigma, although t was reported to 7 figures.

While 5.67E-08_DP may provide a valuable lesson in the definition of precision, I regard this as an "overloaded" constant, where 5.67d-8 is more typical. I don't think we should adopt 5.67E-08_8 ? Can I describe this as the more typical way of defining a real(8) constant ?

Craig Dedo's post is a valuable lesson as 5.67E-08_16 is a necessary way of defining a real*16 constant, which should probably be provided as 5.67E-08_QP.

John

A sidebar: as you know, asterisk (*) declarations such as INTEGER*2, REAL*4, REAL*8, CHARACTER* are marked for obsolescence in the newer standards even though most compilers still support them.  I prefer not to mention them in communications at all since it's a convention I don't like but which is too sticky.

Re: "5.67d-8 is more typical. I don't think we should adopt 5.67E-08_8", well most of our code adopts the convention explained by Craig above and use either "_DP/_dp".  Use of "D" (0D0, etc.) is generally avoided.  For generic code intended to work with a variety of real kinds, we also use "_WP/_wp" convention denoting working precision.  You'll also find it used widely in a recent book, "Numerical Computing with Modern Fortran" by Hanson and Hopkins that was reviewed by Steve in his Dr Fortran blog.

[John_Campbell]

Fortranfan,

I am not a fan of SELECTED_REAL_KIND, as I find what it implies to be very misleading.

The concept that you can request a precision or exponent range that is required for the numerical calculation to be carried out is far removed from reality, as all you have is real*4, real*8 and either real*10 or real*16. Whenever I use SELECTED_REAL_KIND, I have to go and check what are the values for P= or R= that are consistent with the real*x I am going to use.

Even then; real*4 is never good enough, real*8 should work and if you require real*16, then you should rewrite your algorithm. How do you check a real*16 calculation ?

Have you ever developed a numerical computation by deciding I need 13 significant figures so set P=13 ? You don't consider the precision of the inputs, but the expected precision of the result and estimate what *4, *8 or better might provide. Given that consideration, you select the KIND required and then have to check which P will get it right. My background is numerical methods computation, so I don't have an interest in a verbose definition that does not define the accuracy of the result.

Enough of my runblings, I'll let you get back to Computing with Modern Fortran.

John

[Steven_L_Intel1]

John, you misunderstand what SELECTED_REAL_KIND is for. The problem it solves is that there are no guarantees what the precision and range of various real kinds are, across implementations. Also, syntax such as real*4, etc., is non-standard. You use SELECTED_REAL_KIND to specify a minimum acceptable precision and range and let the implementation give you the kind that best meets your needs.

[John_Campbell]

Steve,

My apologies for continuing this discussion, as the standard and it's direction is not going to change, but to me SELECTED_REAL_KIND misunderstands what Fortran is for. There are few number precisions available on the existing hardware platforms in use, and the consideration of what precision to use is more based on what number formats are available and not the P= or R=, as byte size is often the primary consideration.
Having experienced most hardware types over the last 40 years, and received and converted code from many sources, the specification of real*8 is a much preferable syntax that any of the other offerings. Certainly, when this is not used in a code, it flags that there can be many problems ahead in using this code. While use of "DP = Selected_Real_Kind (P=15,R=300)" may be elegant, often it can be difficult to find with code listings, as it can be located in modules removed from the main code. Look at all the KIND definitions associated with ifort.
That there was no guidance for the value of kind in the F90 standard was a shortcoming of the standard, as with an intrinsic function to return the byte count, although I suppose the existence of votes from hardware representatives not based on 8-bit bytes didn't help.
I start from the position that Fortran is used for numerical computation and the standard should support that. I think that the recent standards have moved away from this and what Fortran is used for has either been forgotten or has changed significantly.

John