Program to calculate one- and two-way probabilities for t-tests

DEF FNt(x) = (1 + x^2/v)^(exponent)
DIM n9(16), t(136)

INPUT "enter n   (v=n-2 will be computed) = ",n
v = n-2                    
x = (v+1)/2
GOSUB gamma
g1 = gamma
x = v/2
GOSUB gamma
g2 = gamma
g2 = SQR(v*pi)*g2
constant = g1/g2
exponent = -(v+1)/2
GOSUB t:

gamma:
pi = 3.14159
y = x + 2
g = SQR(2*pi/y) * EXP(y * LOG(y) + (1 - 1/(30*y^2))/(12*y) - y)
gamma = g/(x * (x+1))
RETURN

t:
lower = 0
INPUT "t = ",tstat
upper = ABS(tstat)
tol = .001
GOSUB integration
PRINT
PRINT "probability of"
PRINT
IF tstat < 0 THEN 
  PRINT "     t  <=";
  PRINT tstat;
  PRINT USING "   =  #.######"; .5 - sum*constant
ELSE
TABLE 10-4  (CONTINUED)

  PRINT "     t  >=";
  PRINT tstat;
  PRINT USING "   =  #.######"; .5 - sum*constant
END IF
PRINT
PRINT "     ";:
PRINT -tstat;
PRINT "  <   t   <  ";
PRINT tstat;
PRINT USING "   = #.######"; (.5 - sum*constant)*2
a$ = INPUT$(1)
END

integration: 
panels = 1
n9(1) = 1
delta = (upper - lower)/panels
c = (FNt(lower) + FNt(upper))/2
t(1) = delta * c
n = 1
n8 = 2
sum = c
again:
  n = n+1
  f3 = 4
  n9(n) = n8
  panels = panels*2
  l = panels - 1
  delta = (upper - lower)/panels
  FOR ii = 1 TO (l+1)/2
    i = ii*2 - 1
    x = lower + i*delta
    sum = sum + FNt(x)
    NEXT ii
  t(n8) = delta*sum
  n7 = n9(n-1)
  k = n-1
  FOR m = 1 TO k
    j = n8 + m
    n6 = n9(n-1) + m - 1
    t(j) = (f3*t(j-1) - t(n6))/(f3 - 1)
    f3 = f3*4
    NEXT m
  IF n<=4 THEN more
  IF t(n8 + 1) = 0 THEN more
  IF ABS(t(n7 + 1) - t(n8 + 1)) <= ABS(t(n8+1)*tol) THEN out
  IF ABS(t(n8-1) - t(j)) = ABS(t(j)*tol) THEN out
  IF n > 15 THEN fail
  TABLE 10-4  (CONTINUED)

more:  n8 = j + 1
GOTO again
out: sum = t(j)
RETURN
fail: PRINT "no convergence"
RETURN