double factrl(int n)
{
    double a[33],temp;
	int ntop,j;
    ntop = 0;
    a[1] = 1.0;
    if( n < 0 )
	{
        cout<<"negative factorial";
        exit(1);
	}
    else if (n <= ntop)
	{
        temp = a[n + 1];
	}
    else if (n <= 32)
	{
        for( j = ntop + 1;j<=n;j++)
            a[j + 1] = j * a[j];
        ntop = n;
        temp = a[n + 1];
	}
    else
        temp = exp(gammln(n + 1.0));
	return temp;
}