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 gammaBecause 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 integer, parameter :: long = selected_real_kind(99,999)Solaris Fortran would return

`long=-3`

(this happens during
compilation, `2.0_-3`

are illegal.
You can easily check that our specificationinteger, parameter :: long = selected_real_kind(9,99)returns

`long=8`

.