fpprec: 30$

gamma_coeff(n) := block([a: makelist(1, n)],
   a[2]: bfloat(%gamma),
   for k from 3 thru n do
      a[k]: bfloat((sum((-1)^j * zeta(j) * a[k - j], j, 2, k - 1) - a[2] * a[k - 1]) / (1 - k * a[1])),
   a)$

poleval(a, x) := block([y: 0],
   for k from length(a) thru 1 step -1 do
      y: y * x + a[k],
   y)$

gc: gamma_coeff(20)$

gamma_approx(x) := block([y: 1],
   while x > 2 do (x: x - 1, y: y * x),
   y / (poleval(gc, x - 1)))$

gamma_approx(12.3b0) - gamma(12.3b0);
/* -9.25224705314470500985141176997b-15 */