dim c(12)

sub gamma(z)
    local accm, k, k1_factrl

    accm = c(1)
    if accm=0 then
        accm = sqrt(2*PI)
        c(1) = accm
        k1_factrl = 1
        for k=2 to 12
            c(k) = exp(13-k)*(13-k)^(k-1.5)/k1_factrl
            k1_factrl = k1_factrl * -(k-1)
        next
    end if
    for k=2 to 12
        accm = accm + c(k)/(z+k-1)
    next
    accm = accm * exp(-(z+12))*(z+12)^(z+0.5)
    return accm/z
end sub

sub si(x)
    print x using "%18.13f"
end sub


for i = 0.1 to 2.1 step .1
    print i, " = "; : si(gamma(i))
next