numerical evaluation of quotients of gamma functions

Maybe I'll avoid overflows by using bigfloats--the cancellation in log_gamma(float(x)) - log_gamma(float(y)))
might be too large for some calculations I might do:

  (%i1) x : 10^7+1;
  (%o1) 10000001

  (%i2) y: 10^7;
  (%o2) 10000000

  (%i3) rel(a,b) := abs((a-b)/a)$

set fpprec : 16 . This is about IEEE binary64 numbers with larger exponent range 

  (%i4) gamma(bfloat(x)) / gamma(bfloat(y)), fpprec : 16;
  (%o4) 9.99999999999898b6

Loose about 3 digits 

  (%i5) rel(y,%);
  (%o5) 1.019798219203949b-13

Using binary64 numbers and the log_gamma function, we loose about 8 digits:

  (%i6) exp(log_gamma(float(x)) - log_gamma(float(y)));
  (%o6) 1.0000000152118005*10^7

  (%i7) rel(y,%);
  (%o7) 1.521180048584938*10^-8

--Barton

s.edu/mailman/listinfo/maxima