function theResult = inv_fisher(p, dof1, dof2, tol)

% inv_fisher -- Inverse Fisher F probability.
%  inv_fisher(p, dof1, dof2, tol) returns the Fisher F value
%   for which the integral from f to Inf is probability p,
%   for dof1 and dof2 degrees-of-freedom, computed to
%   tolerance tol (default = 1e-8).
 
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 25-Jan-1999 14:06:53.

MAXIMUM_ITERATIONS = 100;
SUB_SAMPLES = 100;
DEFAULT_TOLERANCE = 1.0e-8;

if nargin < 3, help(mfilename), return, end
if nargin < 4, tol = DEFAULT_TOLERANCE; end

if ischar(p), p = eval(p); end
if ischar(dof1), dof1 = eval(dof1); end
if ischar(dof1), dof2 = eval(dof2); end
if ischar(tol), tol = eval(tol); end

if length(p) > 1 | length(dof1) > 1 | length(dof2) > 1
	if length(p) > 1
		result = zeros(size(p));
		if length(dof1) == 1, dof1 = dof1 .* ones(size(p)); end
		if length(dof2) == 1, dof2 = dof2 .* ones(size(p)); end
	elseif length(dof1) > 1
		result = zeros(size(dof1));
		if length(p) == 1, p = p .* ones(size(dof1)); end
		if length(dof2) == 1, dof2 = dof2 .* ones(size(dof1)); end
	elseif length(dof2) > 1
		result = zeros(size(dof2));
		if length(p) == 1, p = p .* ones(size(dof2)); end
		if length(dof1) == 1, dof1 = dof1 .* ones(size(dof2)); end
	end
	for i = 1:prod(size(p))
		result(i) = inv_fisher(p(i), dof1(i), dof2(i), tol);
	end
	if nargout > 0
		theResult = result;
	else
		disp(result)
	end
	return
end

a = 0;
b = 1;

while fisher(b, dof1, dof2) > p
	a = b;
	b = 2 * b;
end

result = [];

for i = 1:MAXIMUM_ITERATIONS
	delta = (b - a) ./ SUB_SAMPLES;
	a = a:delta:b;
	a(end) = b;
	c = fisher(a, dof1, dof2);
	k = 1:length(c)-1;
	f = find(c(k) >= p & c(k+1) <= p);
	f = f(1);
	b = a(f+1);
	a = a(f);
	c = (c(f) + c(f+1))/2;
	err = c - p;
	if abs(err) < tol
		result = (a + b)/2;
		break;
	end
end

if nargout > 0
	theResult = result;
else
	disp(result)
end
