n=10000000; sn=0; for k=[n:-1:1] sn=sn+1.0/k; end display(sprintf("gamma(%d) = %.13f",n,sn-log(n)))