CMPLX: Fortran 2008 complex intrinsics on branch cuts and minus zero

Download CMPLX Download CMPLX
Download CMPLX
Download CMPLX

Fortran 2003 standard added support for complex arguments in intrinsic functions LOG and SQRT. Fortran 2008 standard added support for complex arguments in intrinsic functions ACOS, ACOSH, ASIN, ASINH, ATAN, ATANH. Fortran 2003 added support for IEEE 754, the floating point standard, which defines minus zero and rules for its use.

Although seemingly an esoteric topic, the values of complex functions on branch cuts are very important, e.g. in applications in elasticity, where the most important data is always on boundaries. Another example is the error function, erf(z), used extensively e.g. in acoustics. erf(z) is calculated using SQRT of complex arguments, including on branch cuts.

The aim of this test is validate correctness of Fortran 2008 intrinsics LOG, SQRT, ACOS, ACOSH, ASIN, ASINH, ATAN, ATANH on branch cuts with the use of minus zero.

CMPLX package is distributed under BSD license

Feedback and support

With any questions, bug reports, feedback and any other comments please submit a ticket.

Latest release

1-FEB-2017: Release 2.1.

This release contain new functionality. The images can now be plotted with gnuplot. The report can now be made with groff. The pgplot and latex options have been improved. Several typos in the comments have been fixed.

How to build

(Refer to README for more info.) The top level directory has several Makefiles. Simply choose one for your platform. The top level Makefile will then descend in all subdirecories and build the code there. At present, there 5 subdirectories:

common       - Fortran code common to all modules and programs.
doc          - build a PostScript and a PDF with analytical
               derivations for the test values on the branch cuts.
terminal     - a terminal test program, which reports in plain text
               the values of the complex intrinsics for special
               values on branch cuts and branch points.
with_gnuplot - a graphical report with conformal maps of the branch
               cuts and the rest of the complex plane. The results
               are saved to file, the gnuplot instructions are also
               saved to file. Then gnuplot is called to produce the
               the plots. The Latex is used to prepare the report in
               PostScript and PDF.
with_pgplot  - a graphical report with conformal maps of the branch
               cuts and the rest of the complex plane. PGPLOT library
               is used extensively to produce the plots directly from
               the executable program. Then Latex is used to prepare
               the report in PostScript and PDF.

The precision of complex kinds is controlled by named constant fk in fun.f90.

Brief description of contents of CMPLX package

There are 4 main parts to the CMPLX package:

  1. A simple terminal program term.f90, which calculates values of all 8 above intrinsic functions at some special complex arguments and compares against the analytical solutions. The output is in a form of a text table in a terminal. Only a Fortran 2008 conforming processor is required to build this program.

    Each mterm* module contains a single public subroutine with checks for several special arguments, e.g. mtermacos.f90 has:

    module mtermacos
    use fun
    use msetfmt
    implicit none
    private
    public :: termacos
    contains
    subroutine termacos( total, pass )
    

    with check of this type:

    !*** test 1 ***
     total = total + 1
        z = cmplx( - hugefk, zerop, kind=fk )
      res = acos(z)
       re = real( res, kind=fk )
       im = aimag( res )
    ! test
      if ( im .lt. 0 .and. re .gt. 0 .and. ieee_is_finite(im) ) then
        pass = pass + 1
        word = "PASS"
      else
        word = "FAIL"
      end if
    ! write out
      write (*,fmt) "z=(-HUGE,+0)=", z
      write (*,fmt) "ACOS(z)=", res, word
    

    Module fun.f90 defines the named constants fk, zerop, zerom, hugefk etc. (It defines some other functions which are not used ATM.)

    Module msetfmt.f90 defines the output formats based on precision and range of the floating point kind - real32, real64 or real128.

  2. Folder with_gnuplot contain program main.f90 used to calculate and plot conformal maps with gnuplot, including of the branch cuts. The resulting PostScript plots are included in a Latex document - gnuplot res.tex On Cray the report is made with groff from res.1.
  3. Folder with_pgplot contain program main.f90 used to calculate and plot conformal maps with PGPLOT, including of the branch cuts. The resulting PostScript plots produced with PGPLOT are included in a Latex document res.tex.
  4. Folder doc contains detailed but simple analytical derivations to justify the analytical results used in term.f90 as reference values. The top level file in this directory is doc.1. A groff installation is required to build this document.

Refer to README for more details.

Fundamentals of inverse complex trigonometric and hyperbolic functions

The most autoritative information on the mathematics is probably DLMF, speficically Sections 4.2(i) The Logarithm for LOG, 4.23 Inverse Trigonometric Functions for ACOS, ASIN, ATAN and 4.37 Inverse Hyperbolic Functions for ACOSH, ASINH, ATANH.

For an in depth analysis of numerical altorithms for calculating these fuctions refer to W. Kahan, Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit, in The State of the Art in Numerical Analysis, A. Iserles, M. J. D. Powell Eds., Clarendon Press, Oxford, 1987. Here's a PDF of this paper.

Results

This is work in progress. The analytical derivations are in doc.pdf. Conformal maps calculated with gfortran7 and plotted with PGPLOT are in res_pgplot.pdf, and plotted with gnuplot are in res_gnuplot.pdf. Conformal maps from Cray are in res_groff.pdf.

These results are preliminary! The tests need more thorough checking!

The tests are very simple. We check that the signs is correct for finite and infinite values and that the values are finite or infinite. The accuracy of floating point results is not assessed.

real32

real64

real128

Other Fortran related projects

validate this page