function fdist
for n=1:10, for m=1:10, f(n,m)=fa(n,m,0.05); end, end
disp([' m n',sprintf(' %7.0f',1:10)]);
for m=1:10, disp(sprintf('%2d: %s',m,sprintf(' %7.3f',f(:,m)))); end

% Fa(n,m) : f-distribution critical value --------------------------------------
function f=fa(n,m,a)

% set range of search
r=[0,10]; while a<fq(n,m,r(2)), r(2)=r(2)*2; end

% binary search
while 1
    f=(r(1)+r(2))/2; p=fq(n,m,f);
    if abs(p-a)<a*1E-4, break, elseif p<a, r(2)=f; else r(1)=f; end
end

% f-distribution upper probability ---------------------------------------------
function p=fq(n,m,x)
p=1-betainc(n*x/(m+n*x),n/2,m/2);
