function F = t_cdf (t, k)
% F = t_cdf (t, k)
% Compute the t distribution CDF efficiently (without calling betainc)
% for positive scalar integer degree of freedom k,  up to 1e4

switch k                % degree of freedom
  case 1                % one degree of freedom 
    F = 0.5 + atan(t)/pi;
  case 2                % two degrees of freedom 
    F = 0.5 + t ./ (2 * sqrt(2 + t .^ 2));
  otherwise             % more than two degrees of freedom  k
    ts = t ./ sqrt(k);
    ttf = 1 ./ (1 + ts .^ 2);
    u = s = 1;
    if mod (k, 2)    % odd degree of freedom k
      m = (k - 1) / 2;
      for ii = 2:m
        u .*= (1 - 1/(2*ii - 1)) .* ttf;
        s += u;
      end
      F = 0.5 + (ts .* ttf .* s + atan(ts)) / pi;
    else             % even degree of freedom k
      m = k / 2;
      for ii = 1:(m - 1)
        u .*= (1 - 1/(2*ii)) .* ttf;
        s += u;
      end
      F = 0.5 + (ts .* sqrt(ttf) .* s) / 2.0;
    end
end

% ---------------------------------------------------------
