cgamma <- function(z){
  n<-length(z);
  res <- array(0, n);
  for (kk in 1:n){
    x=Re(z[kk]); y=Im(z[kk]); #gr <- 0; gi <- 0;
    a <- array(0,10);
    x1 <- 0;
    a <- c(8.333333333333333e-02, -2.777777777777778e-03, 7.936507936507937e-04, -5.952380952380952e-04, 8.417508417508418e-04, -1.917526917526918e-03, 6.410256410256410e-03, -2.955065359477124e-02, 1.796443723688307e-01, -1.39243221690590e+00);
    if (y == 0 & x == trunc(x) & x <= 0){
      gr <- 1e400; gi <- 1e400; 
    }
    else if (x < 0){
      x1 <- x;  y1 <- y;
      x <- -x;  y <- -y;
    }
    x0 <- x;
    if(x <= 7.0){
      na<-trunc(7-x);
      x0<-x+na;
    }
    z1<-sqrt(x0*x0+y*y);
    th<-atan(y/x0);
    gr <- (x0-0.5)*log(z1)-th*y-x0+0.5*log(2*pi);
    gi <- th*(x0-0.5)+y*log(z1)-y;
    for (k in 1:10){
      t <- z1^(1-2*k);
      gr <- gr+a[k]*t*cos((2*k-1)*th);
      gi <- gi-a[k]*t*sin((2*k-1)*th);
    }
    k<- 10+1;
    if(x <= 7){
      gr1=0;
      gi1=0;
      for (j in 0:(na-1)){
        gr1 <- gr1 + 0.5*log((x+j)^2+y*y);
        gi1 <- gi1 + atan(y/(x+j));
        end;
        }
      j <- na-1+1;
      gr <- gr-gr1;
      gi <-gi-gi1;
    } 
    if(x1 < 0){
      z1 <- sqrt(x*x+y*y);
      th1 <- atan(y/x);
      sr <- -sin(pi*x)*cosh(pi*y);
      si <- -cos(pi*x)*sinh(pi*y);
      z2 <- sqrt(sr*sr+si*si);
      th2 <- atan(si/sr);
      if(sr < 0){
        th2 <- pi+th2;
      }
      gr <- log(pi/(z1*z2))-gr;
      gi <- -th1-th2-gi;
      x <- x1;
      y <- y1;
    }
    g0 <- exp(gr);
    gr <- g0*cos(gi);
    gi <- g0*sin(gi);
    res[kk] <-gr + 1i*gi;
    #print(c(x+1i*y, res[k]))
  }
  res; 
}