double gammln(double xx)
{
	int j;
	float temp;
    double cof[6],stp,half,one,fpf,x,tmp,ser;
    cof[1] = 76.18009173;
    cof[2] = -86.50532033;
    cof[3] = 24.01409822;
    cof[4] = -1.231739516;
    cof[5] = 0.00120858003;
    cof[6] = -0.00000536382;
    stp = 2.50662827465;
    half = 0.5;
    one = 1.0;
    fpf = 5.5;
    x = xx - one;
    tmp = x + fpf;
    tmp = (x + half) * log(tmp) - tmp;
    ser = one;
    for (j = 1;j<=6;j++)
	{
        x = x + one;
        ser = ser + cof[j] / x;
    }
    temp = tmp + log(stp * ser);
	return temp;
}