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
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_lgammavia 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_lgammaand when retrieving data from it. This is done by the
interfacestatement, in which
d_lgammaand 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_-3are illegal. You can easily check that our specification
integer, parameter :: long = selected_real_kind(9,99)returns