Next: Definition of Up: The Goodness of Fit Previous: Defining a Module and

### 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`.

Next: Definition of Up: The Goodness of Fit Previous: Defining a Module and
Zdzislaw Meglicki
2001-02-26