The Euler Gamma Function: the first shot

The easiest way to implement the Euler gamma function is to call it from a library provided by the vendor. At this stage, however, the program is no longer portable, so in the long run it is better to implement your own function.

Solaris Fortran provides two versions of the Euler gamma function: a double precision version called d_lgamma(x) and a single precision version called r_lgamma(x). Both functions return the logarithm of the Euler gamma function. Under IRIX 6.3 Fortran does not have an easy access to the Euler gamma function. There is a definition traditionally provided in libm.a, but that definition works for C only, and we would have to write a special wrapper to call it from a Fortran program.

So let us do it non-portably under Solaris first. Then we'll implement our own version and test it against that provided by Solaris.

Solaris based definition looks as follows:

  real(kind=long) function gamma(a)
    real (kind=long), intent(in) :: a

    interface
       function d_lgamma(x)
         double precision :: d_lgamma
         double precision, intent(in) :: x
       end function d_lgamma
    end interface

    gamma = exp(d_lgamma(a))

  end function gamma

Because in this case we can't include a definition of d_lgamma via a module interface, since it is not provided, we have to provide that definition ourselves, so that the compiler would know what conversions it may have to perform before passing data to d_lgamma and when retrieving data from it. This is done by the interface  statement, in which d_lgamma and its argument are declared to be of type double precision. A double precision number in Fortran  is a floating point number with twice the precision of a standard ``word''. For example, on 32-bit architectures a double precision floating point number is a 64-bit number. A ``double precision'' specification is much less precise than a specification obtained by using selected_real_kind(9,99), although the latter may often default to double precision. Fortran vendors are not obliged to provide any precision that a user may request. So, for example, if you were to ask for, say,

integer, parameter :: long = selected_real_kind(99,999)

Solaris Fortran would return long=-3 (this happens during compilation, not during the runtime), and then the compiler would complain that statements such as 2.0_-3 are illegal. You can easily check that our specification

integer, parameter :: long = selected_real_kind(9,99)

returns long=8.