gamma functions should be computed in log

test.c

#include <stdio.h>
#include <math.h>
#include <complex.h>

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>

static double complex gamma_fftlog(double complex z) {
  /* Lanczos coefficients for g = 7 */
  static double p[] = {
    0.99999999999980993227684700473478,
    676.520368121885098567009190444019,
    -1259.13921672240287047156078755283,
    771.3234287776530788486528258894,
    -176.61502916214059906584551354,
    12.507343278686904814458936853,
    -0.13857109526572011689554707,
    9.984369578019570859563e-6,
    1.50563273514931155834e-7
  };

  if(creal(z) < 0.5)
    return M_PI / (csin(M_PI*z)*gamma_fftlog(1. - z));
  z -= 1;
  double complex x = p[0];
  for(int n = 1; n < 9; n++)
    x += p[n] / (z + (double)(n));
  double complex t = z + 7.5;
  return sqrt(2*M_PI) * cpow(t, z+0.5) * cexp(-t) * x;
}

static double complex lngamma_fftlog(double complex z) {
  return clog(gamma_fftlog(z));
}

static double complex lngamma_fftlog_gsl(double complex z) {
  gsl_sf_result lnr, phi;
  gsl_sf_lngamma_complex_e(creal(z), cimag(z), &lnr, &phi);
  return lnr.val + I*phi.val;
}


int main(int argc, char *argv[]) {
  double complex vals[5] = {
    -10 + I * 1.0,
    1e-300,
    10 + I * 1.0,
    100,
    1000};
  int i;

  for (i=0;i<5;++i) {
    printf("val: %60.50e +  %60.50ei\n", creal(vals[i]), cimag(vals[i]));
    printf("    ours: %60.50e + %60.50ei\n",
           creal(lngamma_fftlog(vals[i])),
           cimag(lngamma_fftlog(vals[i])));
    printf("    gsl : %60.50e + %60.50ei\n",
           creal(lngamma_fftlog_gsl(vals[i])),
           cimag(lngamma_fftlog_gsl(vals[i])));
  }
}