#define N 30

function gamma(z)
{
    local i;			// int
    local g;			// real

    g = 1.0;
    for(i = 1; i < N; i++)
	g *= (1.0 + 1.0 / i)^z / (1.0 + z / i);
    return g / z;
}

gopen(500, 500, "blue", 1);
gwindow(-5, 4, 4, -5);

for(x = -5; x < 4; x += 0.5)
    gline(x, -5, x, 5);
for(y = -5; y < 5; y += 0.5)
    gline(-5, y, 4, y);

gpen(0);
for(x = -5; x <= 4; x += 0.004)
    gpset(x, gamma(x));

gpause();