code: public static double digamma(double x)
      {
          if(x == 0)
              return Double.NaN;//complex infinity
          else if(x < 0)//digamma(1-x) == digamma(x)+pi/tan(pi*x), to make x positive
          {
              if(Math.rint(x) == x)
                  return Double.NaN;//the zeros are complex infinity
              return digamma(1-x)-PI/tan(PI*x); 
          }
          
          /*
           * shift over 2 values to the left and use truncated approximation 
           * log(x+2)-1/(2 (x+2))-1/(12 (x+2)^2) -1/x-1/(x+1), 
           * the x+2 and x and x+1 are grouped sepratly 
           */
          double xp2 = x+2;
          
          return log(xp2)-(6*x+13)/(12*xp2*xp2)-(2*x+1)/(x*x+x);
      }
docs: Computes the digamma function of the input
  @param x the input value
  @return &psi;(x)