function [A, detA] = ipjfact(n, k)
%IPJFACT Hankel matrix with factorial elements.
%   [A,D] = GALLERY('IPJFACT',N,K) returns A, an N-by-N 
%   Hankel matrix and D, the determinant of A, which is 
%   known explicity.
%
%   If K = 0 (the default), then the elements of A are
%      A(i,j) = (i+j)!    
%   If K = 1, then the elements of A are
%      A(i,j) = 1/(i+j)!  
%
%   Note that the inverse of A is also known explicitly.

%   Acknowledgement: Suggested by P. R. Graves-Morris.
%
%   Reference:
%   [1] M.J.C. Gover, The explicit inverse of factorial Hankel 
%         matrices, Dept. of Mathematics, University of Bradford, 
%         1993.
%
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.6 $  $Date: 1997/11/21 23:29:39 $

if nargin == 1, k = 0; end

c = cumprod(2:n+1);
d = cumprod(n+1:2*n) * c(n-1);

A = hankel(c, d);

if k == 1
   A = ones(n)./A;
end

if nargout == 2
   d = 1;
   if k == 0
      for i=1:n-1
          d = d*prod(1:i+1)*prod(1:n-i);
      end
      d = d*prod(1:n+1);
   else
      for i=0:n-1
          d = d*prod(1:i)/prod(1:n+1+i);
      end
      if rem(n*(n-1)/2,2), d = -d; end
   end
   detA = d;
end
