Procedures to Avoid Loss of Precision: ln(1+x), etc.
These procedures compute expressions involving logarithm, exponential, sine,
cosine, hyperbolic sine, hyperbolic cosine, and the gamma function, some
having arguments that are expressions, more accurately than is possible
using the Fortran intrinsic functions, or the usual gamma function.
This software is part of the MATH77/mathc90 collection here.