#include <stdio.h>
#include <math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf_gamma.h>

/* gcc -o loggamma loggamma.c -lgsl -lm */

void main()
{
  double zr,zi,Zr,Zi;
  gsl_sf_result lnr,arg;
  gsl_complex logg,expg;
  int i,j;
  FILE *fh,*fh2;
  char buffer[1024];

  fh=fopen("loggamma.asc","w");
  fh2=fopen("loggamma2.asc","w");

  for(i=0;i<=100;i++)
  {
    zr=1;
    zi=i/10.0;
    gsl_sf_lngamma_complex_e(zr,zi,&lnr,&arg);
    GSL_SET_COMPLEX(&logg,lnr.val,arg.val);
    expg=gsl_complex_exp(logg);

    Zr=GSL_REAL(expg);
    Zi=GSL_IMAG(expg);
//  printf("%.8e %.8e  %.12f %.12f  %.8e %.8e\n",zr,zi,lnr.val,arg.val,Zr,Zi);
    printf("%.8e %.8e  %.12f %.12f  %.12f %.12f\n",zr,zi,lnr.val,arg.val,lnr.val*cos(arg.val),lnr.val*sin(arg.val));
    sprintf(buffer, "%.8e, %.8e\n", Zr, Zi);
    for(j=0; buffer[j]!='\0'; j++)
    {
      if(buffer[j]=='e') buffer[j]='\'';
    }
    fprintf(fh, "%s", buffer);
    sprintf(buffer, "%.9e, %.9e\n", lnr.val, arg.val);
    for(j=0; buffer[j]!='\0'; j++)
    {
      if(buffer[j]=='e') buffer[j]='\'';
    }
    fprintf(fh2, "%s", buffer);
  }
  fclose(fh);
}