Two small number subtract get zero in Fortran

Question

I am writing a fortran code, which works with high precision. Here I am trying to use ISO_FORTRAN_ENV to achieve that. I find out that if I subtract two small number, the answer I get will be zero. But if I add them it would be fine.

Here is my program:

program test

  USE ISO_FORTRAN_ENV, ONLY : REAL32

  IMPLICIT NONE
  REAL(REAL32), PARAMETER :: RE_L=-1.7499576837060936950743114606, &
                            RE_R=-1.7499576837060933110499019727
  REAL(REAL32) :: A 

  A=(RE_R)-(RE_L)
  PRINT 20, A
  20 FORMAT(f50.40)


end program test

The result I got is 0.00000000000000000000000000000. And I do not know why.

We have many related questions asked already, I suggest readimg them https://stackoverflow.com/questions/10016261/precision-problems-of-real-numbers-in-fortran https://stackoverflow.com/questions/44052886/fortran-floating-point-equality https://stackoverflow.com/questions/25715548/assigning-a-lower-precision-number-to-a-higher-precision-in-fortran90 — Vladimir F Героям слава, Feb 07 '19 at 07:12

score 2 · Accepted Answer · answered Feb 07 '19 at 02:49

2

You are bumping into Single Precision significant digits in the literals, and in the REAL32 range.

Bump to REAL64, add a double specifier on the literals 1.7499576837060933110499019727d+0

prompt$ gfortran toosmallorg.f90
prompt$ ./a.out
        0.0000000000000004440892098500626161694527

Even with REAL64 d literals, you only get 16 digits after the decimal point, and the displayed result is wrong. With REAL32 only 7 digits after the decimal are recognized (under the default settings this Fedora machine has with gfortran 7.3). REAL128 might give enough wiggle room with the literals.

Fair warning; Not an expert regarding Fortran or the numeric precisions allowed

answered Feb 07 '19 at 02:49

Brian Tiffin

3,978
1
24
34

Right! Now I got the result! Thank you! – Shiqi He Feb 07 '19 at 15:23
By the way, so d is a double specifier, what about f? It cannot be taken as a specifier? – Shiqi He Feb 07 '19 at 15:40
@ShiqiHe, for floating point literals, the two types of Exponent literals are E for single, D for double. In the above `d+0` means treat as double precision but don't shift the magnitude. – Brian Tiffin Feb 09 '19 at 12:16

score 0 · Answer 2 · answered Feb 07 '19 at 02:46

0

This is a long story, you can search google for some key works like " What Every Programmer Should Know About Floating-Point Arithmetic", some very good relevant pages would show up.

In short, float point format, float or double, can only represent the number up to certain precision. The algorithm should be carefully designed in certain cases to get reasonable accuracy, for example, the minimum float point x that makes (1.0+x)-1.0!=0.0 is defined as FLT_EPSILON, it's value is surprisingly large than you may expected.

answered Feb 07 '19 at 02:46

KL-Yang

381
4
8

1

There is no FLT_EPSILON in Fortran and no `!=` in Fortran. And the book is What should every Computer scientist ... and I recommend it *only* to real computer scientists. It is *way* to advanced for an average programmer like me, especialy if they suppose I should really *know* all that stuff by heart. – Vladimir F Героям слава Feb 07 '19 at 07:00
no matter how many digits fortran print out, only about the first 7 digits in 1.234567Exx is valid for float, all other value is not accurate. you can try that there is no single float point number between 1.0 and 1.0000001192092895507812500000000000000000. – KL-Yang Feb 07 '19 at 07:20

Two small number subtract get zero in Fortran

2 Answers2