square root equals zero for qud precision?

[Darrell]

Compare the output of square root of the value below
using Intel Fortran version 11.1.5.5 and 12.1.102.102 update A runtime libraries
Compare using real, real*8 (double), and real*16 quadruple precision floating point values

Our calculation is melting down with a NaN due to dividing by zero from the sqrt until we upgrade intel fortran runtime libraries on our customers computers.

sqrt(0) is zero and you should never get zero any other time, right?

x = 8512183.56904001394286751747131348

print *,sqrt(x)
print *,dsqrt(x)
print *,x**0.5
print *,x**(1./2.)
print *,x**(1.e0/2.e0)
print *,x**(1.d0/2.d0)
print *,x**(1.q0/2.q0)

this only happened in the Release version and not in the Debugger--it was quite the frustrating experience in tracking it down. While you can switch the compiler version in Visual Studio, the runtime libraries have the same name so it's which version is in the system PATH. We spent a lot of time with dependency walker tracking this down.

what is so special about this number?
why does Intel keep messing with how square root works?

we've spent a lot of time on the floating point compiler options:

precise vs fast
speed vs accuracy

We're using Windows XP (32-bit) we tried Microsoft C++ 2005, 2008, and 2010 too just to check if Fortran calls the C sqrt?

I'll go compare the results with the Gnu and Portland Group compilers then try the same thing with C++ and Java...

[Steven_L_Intel1]

We are always looking to improve our math libraries, especially in cases where they do not return "correctly rounded infinite precision" results.

If your program is linked to the DLL libraries, then you will get whatever the latest version is installed on the system. If you link to the static libraries, it will not change based on what is on the target system.

The various /fp options don't have much if any effect on quad-precision, which is all done in software.

The code you posted will not compile unless x is declared real(8), and even then the constant you assign to x is single-precision. I don't understand what you are trying to show here. Can you elaborate? When I declare x to be real(8), build and run this, I get:

[plain]   2917.56466407173
   2917.56466407173
   2917.56466407173
   2917.56466407173
   2917.56466407173
   2917.56466407173
   2917.56466407173458507466422908363[/plain]

which all looks perfectly reasonable. What is the problem?

[Darrell]

LIBIFCOREMD.DLL v11.1.5.5 (Beta1) calls LIMMMD.DLL v11.0.11.2 where sqrt lives as far as I can tell using dependency walker. Just replacing the runtime libraries with the version from 12.1 makes the problem go away. The offending line of code is marked below. My original post oversimplified the problem, sorry.

Weve spent a lot of time on the compiler floating point settings and DLL Hell. The double (not quad) version works fine with VFv6.

Here are the settings we're using looking for accurate floating point calculations vs fast:
/nologo /fpe:0 /fp:source /fpconstant /Qfp-speculation=safe /module:"Release\" /object:"Release\" /libs:dll /threads /c

[fxfortran]      program main
c
      integer debug ! 0=off; 1=debug prints on
      real*16 g4    ! quad precision input
      real*8 :: z1  ! double precision output 1
      real*8 :: z2  ! double precision output 2
      real*8 :: z3  ! double precision output 3
      real*8 :: z4  ! double precision output 4
c
      debug = 1
      g4 = 2.837394523013337981D+06
      call eightroot(debug,g4,z1,z2,z3,z4)
c
      print *,'g4 = ',g4
      print *,'z1 = ',z1
      print *,'z2 = ',z2
      print *,'z3 = ',z3
      print *,'z4 = ',z4
      print *,'breakpoint here so console stays open in debugger'
c
      end program main
c
      subroutine eightroot(debug,g4,z1,z2,z3,z4)
      implicit none
      integer debug
      real*16 g4
      real*8 :: z1,z2,z3,z4
      real*16 :: x1,y1,cc,dd,ee,x2,y2
      real*16 :: rtdd
      real*16 :: arg1,arg2,arg3,arg4
      arg1= 3.e0*g4
      arg2= 8.e0+3.e0*g4
      arg3= sqrt(arg1) ! arg3.ne.0 or NaN--move line up diff answer !
      arg4= sqrt(arg2)
      cc= arg3*arg4
      dd=-3.e0*g4+cc
      rtdd= sqrt(dd)
c
      ee= - rtdd**5 + g4*(24.e0*rtdd - 12.e0*rtdd**3 +
     . ((54.e0-36.e0*rtdd)*g4 + 6.e0*cc))
c
      if (ee .lt. 0.e0) then
         ee= 4.e0
      end if
      x1= 0.5e0*( 2.e0-rtdd )
      y1= 1.e0/(2.e0*sqrt(6.e0*g4))*sqrt(ee)
      ee= g4*(54.e0*g4+6.e0*cc-rtdd*(24.e0-36.e0*g4))
     .   +dd**1.5e0*(12.e0*g4+dd)
      x2= 0.5e0*( 2.e0+rtdd )
      y2= 1.e0/(2.e0*sqrt(6.e0*g4))*sqrt(ee)
      ee= sqrt(x1**2.e0+y1**2.e0)
      z1= y1/sqrt(2.e0)/sqrt(ee-x1)
      z2= sqrt(ee-x1)/sqrt(2.e0)
      ee=sqrt(x2**2.e0+y2**2.e0)
      z3=y2/sqrt(2.e0)/sqrt(ee-x2)
      z4=sqrt(ee-x2)/sqrt(2.e0)
c
      return
      end subroutine eightroot
[/fxfortran]

[Steven_L_Intel1]

I know we've fixed bugs in the quad-precision libraries over the years. It loioks as this is what you ran into. But even with single and double, we're always looking for places where we can do better.