Quad precision (real*16) accuracy problems

[yozohida]

Hi,

The quadruple precision available in ifort-9.0 is supposed to
be a IEEE-style floating point format with 113 bits of accuracy, right?

It seems that sometimes it doesn't seem to have the full 113 bits
of accuracy; the following program finds the smallest i such that
fl(2^i + 1) - 2^i != 1, but if fails around i = 81 already.

Is there a but in the quad precision emulation somewhere?
 

      program test
        real*16 a, c, x
        integer i

        a = 1
        c = 1
        i = 0

        do while (c .eq. 1)
          i = i + 1
          a = 2 * a
          x = a + 1
          c = x - a
        end do
        print *, 'i = ', i
        print *, 'a = ', a
      end program

It should print out "i = 113" or there abouts, but it
prints out i = 81. It works fine for real*8 case, printing
out i = 53. Various flags such as -O0, -mp, etc doesn't seem
to have much effect.

This is on a P4 system.
Funny thing is that when compiled on a P3 system, it works
just as expected ... perhaps something to do with SSE2 instructions?
Similar problem seem to exist in 8.1 compilers as well.

This was brought to my attention while debugging LAPACK compiled
to work with quad precision (in dlamch.f).

-Yozo (UC Berkeley)

[Steven_L_Intel1]

Exactly which 9.0 version do you have? There have been bugs fixed in the quad precision support over the past year. I get 113 with a compiler from October.

[yozohida]

The problem was showing up with both 9.0.021 and 9.0.031.

However, it turned out that the dynamic executable
was picking up libimf.so from the C++ compiler (icc-9.0.023)
installation. When I set the LD_LIBRARY_PATH to have
Fortran library path first, then everything worked.

More importantly, upgrading to the newest C++ compiler (9.0.030)
also fixed the problem.

So, the summary is: the bug seems to exist in 9.0.021, but fixed in
9.0.31 on fortran side, and similarly on the C++ side.

Sorry for the confusion, I should have checked for latest
updates before asking... (but never realized the C++ installation
affected this :-)

Thanks,
-Yozo

[TimP]

Several of the same libraries appear in the C and Fortran installations. Prior to the 8.0 compilers, Intel C and Fortran shared the same library, in a default installation. This could cause problems if both were not upgraded together. The separate libraries don't solve the problem, if builds from one compiler are run with the other library first on LD_LIBRARY_PATH, or in mixed language builds.

[V_Koch1]

Hi,

I am having a similar problem. Compiling your program with ifort 9.1.036 seems fine as it prints the correct value i=113. The following code

program test2
real*16 t

t = 2.q0 ** (-63)
print *, '1 + 2**(-63) =', 1 + t
t = 2.q0 ** (-64)
print *, '1 + 2**(-64) =', 1 + t
print *, '(1 + 2**(-64)) - 1 =', (1 + t)-1
end program

should print something like
1 + 2**(-63) = 1.0000000000000000001084202172485504
1 + 2**(-64) = 1.0000000000000000000542101086242752
(1 + 2**(-64)) - 1 = 5.4210108624275221700372640043497086Q-020

but I get (a wrong second line)
1 + 2**(-63) = 1.00000000000000000010842021724855
1 + 2**(-64) = 1.00000000000000000000000000000000
(1 + 2**(-64)) - 1 = 5.421010862427522170037264004349709E-0020

This is on a Core Duo System using SuSE Linux 10.1.

-Volker

[Steven_L_Intel1]

From what I can tell, the computation is correct but the formatted output is not. I'll report this to development.

[jimdempseyatthecove]

I have seen similar problems with real(8) computations that involved constants in expressions where the constant demoted the expression.

Experiment with:

program test
real*16 a, c, x, One, Two
integer i
One = 1
Two = 2
a = One
c = One
i = 0

do while (c .eq. One )
i = i + 1
a = Two* a
x = a + One
c = x - a
end do
print *, 'i = ', i
print *, 'a = ', a
end program

If the above works properlythen this would deffinately point to the source of the problem.

Also, until you get the compiler update as a work aroundyou can change your code to use variables in lieu of literal values.

Jim Dempsey

[Steven_L_Intel1]

As I said, the actual computation is fine. You can see the correct value if you dump the variable in hex. It's the formatted output that is getting it wrong.