function q = chisqq(p,v)
%CHISQQ	Quantiles of the chi-square distribution.
%	Q = CHISQQ(P,V) satisfies P(X < Q) = P, where X follows a
%	chi-squared distribution on V degrees of freedom.
%	V must be a scalar.
%
%	See also CHISQP.

%	Gordon K Smyth, University of Queensland, gks@maths.uq.edu.au
%	26 July 1998

%	Reference:  Johnson and Kotz (1970). Continuous Univariate
%	Distributions, Volume I. Wiley, New York.

if v <= 0, error('V must be positive'); end
if v == 1, q = normq((1-p)/2).^2; return; end
if v == 2, q = -2*log(1-p); return; end

q = NaN * ones(size(p));

k = find(p <= 0);
if any(k), q(k) = 0; end

k = find(p >= 1);
if any(k), q(k) = Inf; end

k = find((p > 0) & (p < 0.05));
if any(k)
%  Polynomial approximation to density
   q(k) = 2*exp( 2/v*(gammaln(v/2+1)+log(p(k))) );
end

k = find((p > 0.05) & (p < 1));
if any(k)
%  Wilson-Hilferty approximation
   c = 2/9/v;
   q(k) = v*( normq(p(k))*sqrt(c) + 1 - c ).^3;
end

k = find((p > 0) & (p < 1));
if any(k)
%  Newton iteration
   a = v/2*log(2) + gammaln(v/2);
   Q = q(k);
   P = p(k);
   for i=1:3,
      Q = Q - (chisqp(Q,v) - P) ./ exp((v/2-1)*log(Q) - Q/2 - a);
   end
   q(k) = Q;
end
