function g = lacz_gamma(a, cg=7)
  p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \
        771.32342877765313, -176.61502916214059, 12.507343278686905, \
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ];
  x=a;
  if ( x < 0.5 )
    g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) );
  else
    x = x - 1.0;
    t = p(1);
    for i=1:(cg+1)
      t = t + p(i+1)/(x+double(i));
    endfor
    w = x + double(cg) + 0.5;
    g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t;
  endif
endfunction


for i = 1:10
  printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));
endfor