# numdif.s                         J F Monahan, June 2001
 # demonstration of numerical differentiation
 #  using log(gamma(x))
 x <- .5
 k <- -4*(1:15)
 h <- 2**k
 xph <- x + h
 xmh <- x - h
 fx <- lgamma(x)
 fxph <- lgamma(xph)
 fxmh <- lgamma(xmh)
 fp1 <- (fxph - fx) / h
 fp2 <- (fxph - fxmh) / (2*h)
 fpp <- (fxph -2*fx+fxmh) / (h*h)

 tabl <- paste('k',format(k),' 2**k',format(h),'fp1',format(fp1),
           'fp2',format(fp2),'fpp',format(fpp) )
 cat(tabl,sep="\n")
                 
 rm(list=ls())