function a = ijfact3_inverse ( n )

%*****************************************************************************80
%
%% ijfact3_inverse() returns the inverse of the IJFACT3 matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 December 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Karen Habermann,
%    An explicit formula for the inverse of a factorial Hankel matrix,
%    arXiv:1808.02880v2 [math.NA] 14 Aug 2018
%
%  Input:
%
%    integer N, the order of A.
%
%  Output:
%
%    real A(N,N), the matrix.
%
  a = zeros ( n, n );
  for i = 1 : n
    for j = 1 : n
      s = 0.0;
      for k = 0 : i - 1
        if ( j - 1 <= n - i + k )
          s = s + nchoosek ( n - i + k, j - 1 ) * nchoosek ( n + k - 1, k );
        end
      end
      c1 = ( -1.0 ) ^ mod ( n + i + j + 1, 2 );
      c2 = factorial ( i - 1 );
      c3 = factorial ( j );
      c4 = nchoosek ( n - 1, i - 1 );
      c5 = nchoosek ( n + j - 1, j );
      a(i,j) = c1 * c2 * c3 * c4 * c5 * s;
    end
  end

  return
end
