# incomplete gamma function #
MODE R = LONG LONG REAL;
REAL th = (
  PROC sm = (REAL lt)REAL: IF R(1.0) < (R(1.0) + lt) THEN sm(0.3 * lt) ELSE lt FI;
  sm(1));
# ; sum(k, t) returns summation from k to infinity where tk is (s+k)(kth term). #
PROC gamma = (R s, REAL z)R: (
  PROC sum = (INT k, R tk)R: (ABS tk < th | 0 | tk/(s+k) + sum(k+1, -tk*z/(k+1)));
  sum(0, long long exp(s*long long ln(z))));
(1=1 | FOR i TO 50 DO REAL x = i*0.5; print((x, gamma(2, x), newline)) OD)
#print((gamma (2, 2)))
5.93994150290161924318001515082546789777105362271272355595523381 e -1#
#print((gamma (1, 1)))
6.32120558828557678404476229838539132554188868968232165492163198 e -1#
#print((gamma (2, 1)))
2.64241117657115356808952459677078265108377737936464330984326397 e -1#
#print((gamma (2, 25)))
9.99999999638913459510935464538794142601741720329606230121161191e e -1#