The value small=2**-(p+3), where p is the precision, can be determine from

The value small=2**-(p+3), where p is the precision, can be determine from
lgamma(x) = -log(x) - log(1+x) + x*(1-g) + x**2*P(x) with g = 0.57...
being the Euler constant and P(x) a polynomial.  Substitution of small
into the RHS shows that the last 3 terms are negligible in comparison to
the leading term.  The choice of 3 may be conservative.

The value large=2**(p+3) is detemined from Stirling's approximation
lgamma(x) = x*(log(x)-1) - log(x)/2 + log(2*pi)/2 + P(1/x)/x
Again, substitution of large into the RHS reveals the last 3 terms
are negligible in comparison to the leading term.

Move the x=+-0 special case into the |x|<small block.

In the ld80 and ld128 implementaion, use fdlibm compatible comparisons
involving ix, lx, and llx.  This replaces several floating point
comparisons (some involving fabsl()) and also fixes the special cases
x=1 and x=2.

While here
  . Remove unnecessary parentheses.
  . Fix/improve comments due to the above changes.
  . Fix nearby whitespace.

* src/e_lgamma_r.c:
  . Sort declaration.
  . Remove unneeded explicit cast for type conversion.
  . Replace a double literal constant by an integer literal constant.

* src/e_lgammaf_r.c:
  . Sort declaration.

* ld128/e_lgammal_r.c:
  . Replace a long double literal constant by a double literal constant.

* ld80/e_lgammal_r.c:
  . Remove unused '#include float.h'
  . Replace a long double literal constant by a double literal constant.

Requested by:	bde