Complex divide problem in HP-UX 9.0 Fortran

[uunet!apollo.hp.com!johnston]

Cleve Moler writes (about HP-UX Fortran):

> The first thing that set me off today was the Fortran compiler and
> math library... Here is a test program that illustrates just
> two of my complaints.

        real*8 y
        complex*16 z
c
c     The following should produce y = IEEE floating point infinity.
c     But this message is produced at compile time:
c
c     Warning on line 9 of xxx.f: divide by zero detected; zero result inserted
c                                                          ^^^^
        y = 1.d0/0.d0
        write(6,*) y
c
c     The following should NOT produce y = 0.
c     It appears that zabs(z) is being computed as
c        sqrt(real(z)**2 + imag(z)**2)
c     without any scaling to protect against under/overflow.
c
        z = dcmplx(3.0d-200,4.0d-200)
        y = zabs(z)
        write(6,*) y
c
c     The output is y = 0.0d0
c
        end


The problem with complex division has been around for a while.  It was fixed
in our library source code in late 1992, unfortunately too late for the 9.0 release
of HP-UX.  Then came 9.0 Fortran, which inlined complex division, using the old
algorithm!  That means you can't fix the problem with just a new library, you
need a fixed compiler as well. 28

If your Fortran compiler is on maintenance, you can request S700 Fortran compiler 
patch #7 (dated Nov. 12/93) from your HP field office.  You should also ask
for a fixed math library, which is available in patch PHSS_3219 (dated Oct 12/93).

Dividing by zero (at compile time) to get zero is hard to defend; it appears to 
be an historic anomaly.  I am filing a bug report.  If you divide by zero at run 
time you will get the right thing (namely, IEEE infinity).

Both of these problems should be fixed in HP-UX 10.0.

-Ed Johnston

HP-UX Math Library Project        Email: johnstonach.hp.com  
Hewlett-Packard                   (508) 436-5987
300 Apollo Drive            
Chelmsford, MA 01824